Pith. sign in

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 →

arxiv 2512.18066 v2 pith:NUDLH63P submitted 2025-12-19 stat.ME

classification stat.ME MSC 62F1562M3060G15
keywords deepGaussianprocessesgradientenhancementsurrogatemodelingcomputerexperimentsnonstationaryemulationpredictionsVecchiaapproximationBayesianMCMC
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Deep Gaussian processes (DGPs) — two-layer compositions of Gaussian processes that warp a nonstationary input space into a smoother one — are popular surrogates for expensive computer experiments, but they had not been integrated with derivative information. This paper tries to establish that both gradient-related capabilities are achievable inside a fully Bayesian DGP: absorbing observed gradient data into training ('gradient enhancement') and returning posterior predictive distributions of gradients at unobserved inputs. The bridge is the multivariate chain rule, which factorizes the response's input-gradient into the latent warp's derivative matrix times the outer layer's gradient in the warped space, so the model can convert between observed input-gradients and latent-space gradients. The paper argues, with benchmarks on nonstationary test functions and, per the abstract, a real quantum-mechanics computer experiment simulating molecular energy and forces as a function of atomic position, that gradient-enhanced DGPs outperform gradient-enhanced GPs and non-enhanced DGPs, and that DGP gradient predictions are more effective than GP gradient predictions. A practitioner cares because many physics and engineering simulators already return gradients at little extra cost, so a surrogate that spends that information wisely delivers more accuracy per expensive simulation run.

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.

Watch

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

Editorial extensions of the paper, not claims the author makes directly.

  • 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.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

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)
  1. [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.
  2. [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.
  3. [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)
  1. [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.
  2. [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.
  3. [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.
  4. [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

0 steps flagged · score 0.0 of 10

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 4 free parameters · 6 assumptions · 0 invented entities

The central method rests on standard GP/DGP building blocks and a set of tuning parameters (lengthscales, nugget, Vecchia settings) plus structural assumptions about the latent warping. No new physical or theoretical entities are introduced. The most fragile element is the unstated invertibility of the latent Jacobian and the implied but omitted change-of-variables factor in the likelihood.

free parameters (4)
  • kernel lengthscales θ (outer and D inner)
    One isotropic lengthscale per GP; sampled via Metropolis-Hastings with Gamma priors (Section 4.2). The warping and outer correlation depend on these values.
  • output variance τ²
    Outer-layer variance; given reference prior π(τ²)∝1/τ² and integrated out (Section 4.2).
  • nugget/jitter ε = fixed small value
    Fixed at a small value for numerical stability on a deterministic simulator; added to entire diagonal in gradient-enhanced covariances (Sections 2 and 3.3).
  • Vecchia conditioning size m and ordering rules
    Hand-specified ordering (responses first) and nearest-neighbor conditioning sets with tie-breaking toward responses; paper says 'found to work well' and defers thorough exploration to future work (Section 4.1).
assumptions (6)
  • standard math The Gaussian/squared exponential kernel is twice differentiable, so derivative observations and gradient predictions are well-defined.
    Used throughout Sections 2-3; derivative kernels Eq. (4) require differentiability.
  • domain assumption Two-layer DGP with single latent layer of dimension D, independent latent nodes, noise-free unit-variance zero-mean latent GPs.
    Section 3; these are the DGP structural choices that make the MCMC sampler tractable; if the true simulator's nonstationarity needs more layers or correlated nodes, the method may not apply.
  • domain assumption The simulator is deterministic, so the nugget can be fixed at a small numerical value rather than estimated.
    Section 2.3 and 4; the stated focus is deterministic, expensive computer experiments.
  • 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.
    Section 3.3 uses this without statement or checking; no tie-breaking or regularization for singular Jacobians is given.
  • domain assumption MCMC with ESS and Metropolis-Hastings, combined with stochastic imputation, yields samples from the correct posterior.
    Sections 3.1-3.3 and 4.2; unless the sampler converges and is correctly derived, posterior summaries are not valid.
  • domain assumption Vecchia approximation with the proposed ordering and tie-breaking provides a faithful approximation of the full GP posterior.
    Section 4.1; paper explicitly defers thorough investigation of ordering choices to future work.

how reviews work

0 comments
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 reproduced from arXiv: 2512.18066 by the authors.

Figure 1
Figure 1. Standard and gradient-enhanced GP and DGP predictions of a simple step function (top) and its [PITH_FULL_IMAGE:figures/full_fig_p004_1.png] view at source ↗
Figure 2
Figure 2. Left: 100 ESS samples of latent w for a DGP fit to the step function from [PITH_FULL_IMAGE:figures/full_fig_p008_2.png] view at source ↗
Figure 3
Figure 3. A random sample from the gradient-DGP of Eq. ( [PITH_FULL_IMAGE:figures/full_fig_p009_3.png] view at source ↗
Figures from the paper (5 more)
Figure 4
Figure 4. Figure 4: Heatmap of the 2d squiggle function (left) from [PITH_FULL_IMAGE:figures/full_fig_p012_4.png]
Figure 5
Figure 5. Figure 5: Vecchia conditioning sets with n = 5, d = 1, and m = 3 [PITH_FULL_IMAGE:figures/full_fig_p013_5.png]
Figure 6
Figure 6. Figure 6: Simulation results for the squiggle function with [PITH_FULL_IMAGE:figures/full_fig_p015_6.png]
Figure 7
Figure 7. Figure 7: Simulation results for the plateau function with [PITH_FULL_IMAGE:figures/full_fig_p016_7.png]
Figure 8
Figure 8. Figure 8: Simulation results for the ignition function with [PITH_FULL_IMAGE:figures/full_fig_p016_8.png]

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

  1. Profile Bayesian Optimization for Expensive Computer Experiments

    stat.ME 2025-12 conditional novelty 6.0 of 10

    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

16 extracted references · 5 linked inside Pith · cited by 1 Pith paper

  1. [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...

  2. [2]

    All Emulators are Wrong, Many are Useful, and Some are More Useful Than Others: A Reproducible Comparison of Computer Model Surrogates

    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...

  3. [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). “...

  4. [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, ...

  5. [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

  6. [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). “...

  7. [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,

  8. [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,

Show all 16 references
  1. [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,...

  2. [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

  3. [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...

  4. [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...

  5. [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...

  6. [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...

  7. [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...

  8. [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,

Pith tools

Reviewed August 3, 2026 · model on record in the stance chip above.