REVIEW 3 major objections 4 minor 1 cited by
Gradient-enhancement and Gradient Predictions for Deep Gaussian Process Modeling of Expensive Computer Experiments
T0 review · 3 major / 4 minor · reviewed 2026-08-03 · deepseek-v4-flash
Pith's one-line read This paper shows a two-layer deep Gaussian process can be upgraded to use and predict gradients via the chain rule, beating gradient-enhanced GPs and plain DGPs on nonstationary surfaces.
desk verdict Useful, well-engineered gradient-enhancement for deep GPs with a real Jacobian omission in the MCMC sampler; conditional accept after fixing the likelihood and the abstract mismatch. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing object is the multivariate chain rule written as the linear system ∇_x y = (∇_x W)(∇_w y) (Eq. 13), in which ∇_x W is the D×D matrix of partial derivatives of the latent warp W with respect to the original inputs and ∇_w y collects the outer layer's derivatives with respect to the warped coordinates. The identity does double duty: during training, observed input-gradients are pushed backwards through the warp to form pseudo-observations of the latent-space gradients that enter the outer-layer likelihood; at prediction, inferred latent-space gradients are pushed forwards to produce posterior predictive distributions of the input-gradients. Around this identity the paper build
What would settle it
On a one-dimensional nonstationary example such as the step-like function in the paper, record det(∇_x W) at each training location across the MCMC iterations: if any determinant is zero or the linear solve in Eq. (13) is numerically unstable, gradient-enhancement is not computable there. Then run the sampler twice — once with the Eq. (16) acceptance ratio and once with a |det(∇_x W)| factor accounting for the change of variables — and compare the posterior predictive distributions of the response and its gradients; any visible divergence shows the omitted Jacobian term is not negligible.
Extended reading notes
Core claim
The central claim is that a deep Gaussian process can be made fully gradient-aware while staying within exact MCMC-based Bayesian inference, and that doing so pays off precisely on the nonstationary surfaces where DGPs are chosen in the first place. Because the derivatives of a Gaussian process are themselves GPs, the paper upgrades both layers of the two-layer DGP to carry derivative information: the inner layer samples each latent node jointly with its input-gradients, and the outer layer carries the response together with its gradients with respect to the warped inputs. The two are linked by the multivariate chain rule ∇_x y = (∇_x W)(∇_w y), which converts observed input-gradients into l
Load-bearing premise
Everything rests on the chain-rule conversion being sound: observed input-gradients are turned into latent-space gradients by solving ∇_x y = (∇_x W)(∇_w y), which requires the latent warp's derivative matrix ∇_x W to be invertible at every training location, while the MCMC acceptance probability in Eq. (16) treats the change of variables as carrying no Jacobian term.
Editorial extensions
If this is right
- On nonstationary response surfaces, gradient-enhanced DGPs give lower RMSE and better-calibrated uncertainty than gradient-enhanced GPs and non-enhanced DGPs, with the largest gains when training data are scarce.
- DGP posterior predictions of gradients are more accurate than GP gradient predictions on nonstationary functions, which should strengthen downstream tasks that rely on gradients, such as derivative-based optimization and active learning.
- Gradient information can be folded into a DGP without abandoning full Bayesian inference: the latent warp and its gradients are sampled jointly by MCMC, avoiding variational or moment-matching shortcuts.
- The Vecchia approximation built into the open-source software keeps gradient-enriched models computationally feasible, extending the method to larger training sets than the cubic-cost baseline would allow.
Reading between the lines
- The chain-rule bridge is compositional, so the same construction should extend to DGPs with more than one hidden layer by linking the gradient identities layer by layer; the paper restricts itself to a single latent layer, but nothing in the mechanism appears to require that.
- Because converting observed input-gradients into latent-space gradients requires solving Eq. (13) at every training location, training designs that avoid near-singular latent warps — or a regularized fallback for the solve — could prove important in applications with flat or degenerate regions.
- A direct testable extension: run the sampler with and without a |det(∇_x W)| Jacobian factor in the Eq. (16) acceptance ratio on a strongly warping function; any divergence in the resulting posterior predictive distributions would indicate the change of variables is not Jacobian-free.
- The gradient predictions the paper develops could feed gradient-based acquisition functions in Bayesian optimization, where the DGP's learned warp might locate sharp features that stationary gradient surrogates smooth over — an application the paper points to but does not implement.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper develops a Bayesian deep Gaussian process (DGP) framework that incorporates gradient information. The model places gradient-enhanced Gaussian priors on both the outer layer (response and response gradients with respect to latent warped inputs) and the inner layer (latent warping nodes and their gradients with respect to original inputs). The multivariate chain rule is used to convert between input-space gradients and warped-space gradients. Gradient-enhanced DGP training is performed by augmenting elliptical slice sampling to sample the latent warping and its gradient; DGP gradient predictions are obtained by applying GP gradient prediction at both layers and multiplying by the inferred warping Jacobian. The methodology is implemented in the CRAN package deepgp and benchmarked on three nonstationary functions and a quantum mechanics-inspired example, comparing GP, DGP, gradient-enhanced GP, and gradient-enhanced DGP in terms of RMSE and CRPS for both responses and gradients. The central empirical claim is that on nonstationary surfaces the gradient-enhanced DGP outperforms gradient-enhanced GPs and non-enhanced DGPs.
Significance. If the proposed framework is statistically valid, it fills a real gap: gradient-enhanced DGPs and DGP gradient predictions are useful for expensive nonstationary simulators, especially for downstream optimization and active learning. The paper is also accompanied by open-source software and reproducible benchmark code, which is a strength. However, the validity of the Bayesian inference and the uncertainty quantification claims rests on the correctness of the likelihood used in the MCMC scheme and on the gradient predictive moment formulas. As detailed below, both contain technical gaps that need to be addressed before the results can be taken as Bayesian posterior inference.
major comments (3)
- [Section 3.3, Eqs. (13) and (16)] The ESS acceptance probability uses the wrong observed-data likelihood. Observed input-space gradients are related to the warped-space gradients by A_i ∇_w y_i = ∇_x y_i, where A_i = ∇_x W_i is D×D. The density of the observations given W_all is therefore p(y, ∇_x y | W_all) = p(y, ∇_w y | W_all) ∏_{i=1}^n |det A_i|^{-1}. Eq. (16) evaluates the Gaussian density at the solved ∇_w y but omits the ∏ |det A_i|^{-1} factor. Because A_i changes at every ESS proposal (17), the stationary distribution is ∝ prior(W_all) p(y, ∇_w y | W_all) ∏ |det A_i|, not the model posterior. This biases the posterior over warpings toward regions with large |det ∇_x W| and can seriously distort uncertainty quantification. The required invertibility of A_i at every training location is also not discussed; near-singular A_i will make the solve in Eq. (13) unstable.
- [Section 3.2, moment formulas] The stated variance of the DGP gradient prediction is not the variance of the linear combination. For a given posterior draw, ∂y/∂x_d = Σ_i (∂w_i/∂x_d) ∂y/∂w_i, so the conditional variance is (∇_x w_d)^T Cov(∂y/∂w_1,...,∂y/∂w_D) (∇_x w_d). The displayed formula V[∂y/∂x_d] = Σ_i (∂w_i/∂x_d)^2 V[∂y/∂w_i] drops all off-diagonal terms of the predictive covariance from Eq. (5). These cross-covariances are generally nonzero for the squared-exponential kernel, so the reported gradient CRPS may not reflect the actual posterior predictive distribution. If the software uses this formula for UQ, the gradient uncertainty statements are not reliable.
- [Section 3.3, Eq. (13)] The method requires the D×D matrix ∇_x W_i to be invertible at every training location i. Under a continuous GP prior, exact singularity has probability zero, but near-singularity is possible and can produce extreme solved values of ∇_w y, especially because the Gaussian kernel gradients can be near zero or near-duplicated. The manuscript does not acknowledge this assumption, discuss numerical safeguards, or report any diagnostics (e.g., condition numbers) in the benchmark experiments.
minor comments (4)
- [Section 3.3] Eq. (16) is described as 'the gradient-enhanced version of the likelihood of the outer layer.' It is actually the density of the transformed pseudo-observations (y, ∇_w y), not the likelihood of the observed data (y, ∇_x y). The wording should be revised together with the Jacobian correction.
- [Section 3.1] In the sentence introducing Wall, 'column-binded' should be 'column-bound' or 'column-bind.' Similarly, 'intractible' in the introduction and 'surrgotes' in Section 6 are typos.
- [Figure 2 caption] The caption says 'For observed training locations, we sample these (Section 3.3)' referring to dw/dx. Section 3.3 actually solves for ∇_w y using Eq. (13); what is sampled is W_all. Please clarify the wording.
- [Section 4.1] The Vecchia ordering description 'then append each ∂y_i/∂x_d using the same ordering, in turn' is ambiguous. Does this mean all derivatives with respect to x_1 for all i, then all with respect to x_2, or derivative-by-derivative within each observation? A short pseudocode block or explicit index sequence would remove ambiguity.
Circularity Check
No circularity found: the derivation is self-contained and validated on out-of-sample benchmarks.
full rationale
The paper's central claim is that gradient-enhanced DGPs and DGP gradient predictions improve over gradient-enhanced GPs and standard DGPs. The derivation chain is self-contained: Eq. (13) is the multivariate chain rule connecting observed input-space gradients ∇_x y to latent-space gradients ∇_w y via the sampled warping Jacobian ∇_x W. Gradient-enhanced training solves this linear system for ∇_w y given observed ∇_x y; gradient prediction at test locations applies the same chain rule to GP posterior draws of ∇_x W and ∇_w y. This is a deterministic change of variables, not a fitting of the answer being predicted. Benchmarks use re-randomized LHS training designs and out-of-sample LHS test sets with known closed-form gradients, so the reported improvements are not forced by construction. The only notable technical issue is the omission of the change-of-variables Jacobian det(∇_x W) in the Eq. (16) acceptance probability, which affects Bayesian correctness of the sampler but is not a circularity: it does not make the prediction equal to the input. Self-citations (the deepgp package and prior DGP work by the same group) are used for implementation and comparison, not as load-bearing proof of the new method's validity.
Assumptions & free parameters
free parameters (4)
- kernel lengthscales θ (outer and D inner)
- output variance τ²
- nugget/jitter ε =
fixed small value
- Vecchia conditioning size m and ordering rules
assumptions (6)
- standard math The Gaussian/squared exponential kernel is twice differentiable, so derivative observations and gradient predictions are well-defined.
- domain assumption Two-layer DGP with single latent layer of dimension D, independent latent nodes, noise-free unit-variance zero-mean latent GPs.
- domain assumption The simulator is deterministic, so the nugget can be fixed at a small numerical value rather than estimated.
- ad hoc to paper The latent Jacobian ∇_x W is square and invertible at every training location so Eq. (13) can be solved for ∇_w y.
- domain assumption MCMC with ESS and Metropolis-Hastings, combined with stochastic imputation, yields samples from the correct posterior.
- domain assumption Vecchia approximation with the proposed ordering and tie-breaking provides a faithful approximation of the full GP posterior.
Cite this review
Pith. "Pith review of Gradient-enhancement and Gradient Predictions for Deep Gaussian Process Modeling of Expensive Computer Experiments." pith.science (2026). https://pith.science/paper/NUDLH63P
@misc{pith2026251218066,
author = {Pith},
title = {Pith review of: Gradient-enhancement and Gradient Predictions for Deep Gaussian Process Modeling of Expensive Computer Experiments},
year = {2026},
howpublished = {\url{https://pith.science/paper/NUDLH63P}},
note = {Machine review of arXiv:2512.18066}
}
read the original abstract
Deep Gaussian processes (DGPs) are popular surrogate models for complex nonstationary computer experiments. DGPs use one or more latent Gaussian processes (GPs) to warp the input space into a plausibly stationary regime, then use typical GP regression on the warped domain. While this composition of GPs is conceptually straightforward, the functional nature of the multi-dimensional latent warping makes Bayesian posterior inference challenging. Traditional GPs with smooth kernels are naturally suited for the integration of gradient information, but the integration of gradients within a DGP presents new challenges and has yet to be explored. We propose a novel and comprehensive Bayesian framework for DGPs with gradients that facilitates both gradient-enhancement and gradient posterior predictive distributions. Our focus is on surrogate modeling of expensive, deterministic, and nonstationary computer experiments. Gradient-enhancement is most impactful when data is limited, and gradient predictions are most useful for downstream surrogate modeling tasks like optimization and active learning. We benchmark both contributions (gradient-enhanced DGPs and DGP gradient predictions), separately and together, on a variety of nonstationary test functions as well as real quantum mechanics computer experiments that simulate molecular energy and forces as a function of atomic position. On nonstationary surfaces, our gradient-enhanced DGPs outperform gradient-enhanced GPs and non-enhanced DGPs, and our DGP gradient predictions are more effective than GP gradient predictions. We provide open-source software in the "deepgp" package on CRAN, with optional Vecchia approximation to circumvent cubic computational bottlenecks.
Figures
Figures from the paper (5 more)
Forward citations
Cited by 1 Pith paper
-
Profile Bayesian Optimization for Expensive Computer Experiments
A two-stage acquisition scheme estimates profile optima of expensive black-box simulations over a control parameter, beating standard BO and profile-EI on tested benchmarks.
Reference graph
Works this paper leans on
-
[1]
Unexpected improvements to expected improvement for bayesian optimization
Ament, S., Daulton, S., Eriksson, D., Balandat, M., and Bakshy, E. (2023). “Unexpected improvements to expected improvement for bayesian optimization.”Advances in Neural Information Processing Systems, 36, 20577–20612. Baker, E., Barbillon, P., Fadikar, A., Gramacy, R. B., Herbei, R., Higdon, D., Huang, J., Johnson, L. R., Ma, P., Mondal, A., et al. (2022...
arXiv 2023
-
[2]
MIT press Cambridge, MA. Rumsey, K. (2025).duqling: Library of UQ Test Functions. R package version 2.0.0, commit b5df65ce6691434919e00b28455bea08f64abbe7. Rumsey, K. N., Gibson, G. C., Francom, D., and Morris, R. (2025). “All Emulators are Wrong, Many are Useful, and Some are More Useful Than Others: A Reproducible Comparison of Computer Model Surrogates...
arXiv 2025
-
[3]
Multifidelity data fusion via gradient-enhanced Gaussian process regression
Deng, Y., Lin, G., and Yang, X. (2020). “Multifidelity data fusion via gradient-enhanced Gaussian process regression.”arXiv preprint arXiv:2008.01066. Dunlop, M. M., Girolami, M. A., Stuart, A. M., and Teckentrup, A. L. (2018). “How deep are deep Gaussian processes?”Journal of Machine Learning Research, 19, 54, 1–46. 17 Dwight, R. and Han, Z.-H. (2009). “...
arXiv 2020
-
[5]
On the instability issue of gradient-enhanced Gaussian process emulators for computer experiments
He, X. and Chien, P. (2018). “On the instability issue of gradient-enhanced Gaussian process emulators for computer experiments.”SIAM/ASA Journal on Uncertainty Quantification, 6, 2, 627–644. Hung, T.-H. and Chien, P. (2021). “A random Fourier feature method for emulating computer models with gradient information.”Technometrics, 63, 4, 500–509. Jacobson, ...
2018
-
[6]
Empirical assessment of deep gaussian process surrogate models for engineering problems
Rajaram, D., Puranik, T. G., Ashwin Renganathan, S., Sung, W., Fischer, O. P., Mavris, D. N., and Ramamurthy, A. (2021). “Empirical assessment of deep gaussian process surrogate models for engineering problems.”Journal of Aircraft, 58, 1, 182–196. Rasmussen, C. E. and Williams, C. K. (2006).Gaussian processes for machine learning, vol
2021
-
[11]
Vecchia-approximated deep Gaussian processes for computer experiments
Springer. Sauer, A., Cooper, A., and Gramacy, R. B. (2023a). “Vecchia-approximated deep Gaussian processes for computer experiments.”Journal of Computational and Graphical Statistics, 32, 3, 824–837. Sauer, A., Gramacy, R. B., and Higdon, D. (2023b). “Active learning for deep Gaussian process surrogates.” Technometrics, 65, 1, 4–18. Sauer, A. E. (2023). “...
2023
-
[15]
Gradient-enhanced reliability analysis of transonic aeroelastic flutter
Stanford, B., Sauer, A., Jacobson, K., and Warner, J. (2022). “Gradient-enhanced reliability analysis of transonic aeroelastic flutter.” InAIAA SciTech 2022 Forum,
2022
-
[16]
Bayesian optimization with gradients
MIT press Cambridge, MA. Wu, J., Poloczek, M., Wilson, A. G., and Frazier, P. (2017). “Bayesian optimization with gradients.” Advances in neural information processing systems,
2017
Show all 16 references
-
[17]
Sequential learning of active subspaces
Wycoff, N., Binois, M., and Wild, S. M. (2021). “Sequential learning of active subspaces.”Journal of Computational and Graphical Statistics, 30, 4, 1224–1237. Wycoff, N. B. (2021). “Gradient-Based Sensitivity Analysis with Kernels.” Ph.D. thesis, Virginia Tech. Yang, Y., Ming,...
2021 arXiv
-
[30]
J., Williams, B
Santner, T. J., Williams, B. J., Notz, W. I., and Williams, B. J. (2003).The design and analysis of computer experiments, vol
2003
-
[31]
Strictly proper scoring rules, prediction, and estimation
Gneiting, T. and Raftery, A. E. (2007). “Strictly proper scoring rules, prediction, and estimation.”Journal of the American statistical Association, 102, 477, 359–378. Gramacy, R. B. (2020).Surrogates: Gaussian process modeling, design, and optimization for the applied science...
2007
-
[65]
Deep Gaussian processes for calibration of computer models (with discussion)
Marmin, S. and Filippone, M. (2022). “Deep Gaussian processes for calibration of computer models (with discussion).”Bayesian Analysis, 17, 4, 1301–1350. McKay, M. D., Beckman, R. J., and Conover, W. J. (2000). “A comparison of three methods for selecting values of input variab...
2022
-
[282]
Global optimization of atomic structures with gradient-enhanced Gaussian process regression
Kaappa, S., Del R ´ ıo, E. G., and Jacobsen, K. W. (2021). “Global optimization of atomic structures with gradient-enhanced Gaussian process regression.”Physical Review B, 103, 17, 174114. Kang, M., Sch¨ afer, F., Guinness, J., and Katzfuss, M. (2024). “Asymptotic properties o...
2021 arXiv
-
[632]
Validated reduced computational methods for realistic RDC injection modeling
Stein, M. L. (1999).Interpolation of spatial data. Springer-Verlag. Stumbar, W., Stigliano, W., Grunenwald, J., Salek, P., Athmanathan, V., Meyer, T., Webb, A., Fugger, C., Miki, K., Perkins, D., et al. (2025). “Validated reduced computational methods for realistic RDC injecti...
1999
-
[1033]
Estimation and model identification for continuous spatial processes
Vecchia, A. V. (1988). “Estimation and model identification for continuous spatial processes.”Journal of the Royal Statistical Society Series B: Statistical Methodology, 50, 2, 297–312. Wang, C.-Z., Nagisetty, K. R., Montanari, F., and Hill, D. C. (2015). “Application of adjoi...
1988
-
[2276]
Scaling Gaussian process regression with derivatives
Eriksson, D., Dong, K., Lee, E., Bindel, D., and Wilson, A. G. (2018). “Scaling Gaussian process regression with derivatives.”Advances in neural information processing systems,
2018
Reviewed August 3, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.