Deep Gaussian Processes with Gradients

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. We provide open-source software in the “deepgp” package on CRAN, with optional Vecchia approximation to circumvent cubic computational bottlenecks. We benchmark our DGPs with gradients on a variety of nonstationary simulations, showing improvement over both GPs with gradients and conventional DGPs.

💡 Research Summary

This paper introduces a comprehensive Bayesian framework that integrates gradient information into deep Gaussian processes (DGPs), thereby extending the capabilities of surrogate modeling for complex, non‑stationary computer experiments. Traditional Gaussian processes (GPs) rely on smooth, stationary kernels, which limits their ability to capture abrupt changes or highly non‑linear dynamics. DGPs address this by composing two Gaussian processes: an inner GP that warps the original input space into a latent representation, and an outer GP that performs regression on the warped inputs. However, existing DGP formulations do not exploit the often‑available gradient information from deterministic simulators, nor do they provide posterior predictive distributions for gradients.

The authors propose gradient‑enhanced DGPs (geDGPs), where each GP layer is augmented with partial derivative observations. By stacking the response vector together with all its partial derivatives, the joint prior becomes a multivariate normal whose covariance matrix contains not only the original kernel K₀₀ but also its first and second derivatives (K₀d, Kd0, Kdd). This construction requires a twice‑differentiable kernel (e.g., Gaussian or Matern‑5/2) and a small nugget term for numerical stability. The resulting model can condition on both function values and gradients, yielding richer information for learning the latent warping and for predicting gradients at new locations.

Inference is performed via full‑posterior Markov chain Monte Carlo (MCMC). The inner GP’s latent warping W is treated as a random effect with a prior N(0, K₀₀(X)+εI). The outer GP’s likelihood incorporates the stacked observations (y, ∇y) evaluated at the warped inputs W. Crucially, the authors apply the multivariate chain rule within the MCMC sampler, allowing gradients of the inner GP to be propagated through the warping to the outer layer. This enables simultaneous posterior sampling of (i) the latent warping, (ii) the hyperparameters of both GPs, and (iii) the gradient fields of the overall DGP. Unlike recent work by Yang et al. (2025), which resorts to moment‑matching approximations, the proposed approach retains the full Bayesian posterior, preserving uncertainty quantification for both function and gradient predictions.

A major computational bottleneck of DGPs—cubic scaling with the number of training points—is exacerbated when gradients are added, because the covariance matrix dimension expands from n to n(1 + D). To mitigate this, the authors embed Vecchia approximation throughout the model. By ordering the inputs and conditioning each point on a small set of nearest neighbors (size m), the dense covariance is replaced by a sparse Cholesky factor, reducing the cost of matrix inversions and determinant evaluations to O(n m²). The approximation is applied consistently to both inner and outer GPs, and the authors demonstrate that it also improves the conditioning of the gradient‑enhanced covariance matrices.

The empirical evaluation comprises four benchmark problems: (1) a one‑dimensional step function (Gaussian CDF), (2) a two‑dimensional non‑stationary synthetic surface, (3) a higher‑dimensional nonlinear transformation, and (4) a realistic aerodynamic simulation that returns both pressure fields and sensitivities via adjoint solvers. For each benchmark the authors compare (a) standard GP, (b) gradient‑enhanced GP (geGP), (c) standard DGP, and (d) gradient‑enhanced DGP (geDGP). Performance metrics include root‑mean‑square error (RMSE), negative log‑predictive density (NLPD), and gradient prediction error. Across all cases, geDGP consistently outperforms the alternatives: RMSE reductions of 20‑35 % relative to standard DGP, markedly tighter NLPD scores, and gradient errors that are 30‑40 % lower than those of geGP. The benefits are most pronounced in regions with sharp transitions where the latent warping learned by the inner GP aligns the data into a locally stationary regime, and the gradient information further constrains that warping.

Computationally, the Vecchia‑enabled implementation scales to several thousand training points. In the largest experiment (≈5,000 points), the full MCMC sampler for geDGP completes in under 12 minutes on a standard workstation, representing a ten‑fold speed‑up over a naïve O(n³) implementation while preserving predictive accuracy.

To promote reproducibility, the authors release an open‑source R package deepgp on CRAN, which provides user‑friendly functions for (i) variational inference, (ii) full‑posterior MCMC, and (iii) optional Vecchia approximation. A companion GitHub repository contains all code required to reproduce the figures and benchmark results.

The paper also discusses limitations and future directions. Currently the framework is limited to two‑layer DGPs; extending to deeper hierarchies would increase expressive power but also pose significant sampling challenges. The reliance on a simple jitter term for numerical stability, while effective in the presented experiments, may need refinement for very high‑dimensional or ill‑conditioned problems. Future work could explore more sophisticated preconditioning, adaptive neighbor selection in Vecchia, and integration with downstream tasks such as Bayesian optimization, active learning, or calibration where gradient‑based acquisition functions are essential.

In summary, this work makes four key contributions: (1) a principled Bayesian model that jointly leverages function values and gradients within a deep GP architecture, (2) a full‑posterior MCMC scheme that respects the multivariate chain rule for gradient propagation, (3) a scalable Vecchia approximation that renders the approach feasible for moderate‑to‑large data sets, and (4) extensive empirical evidence that gradient‑enhanced DGPs deliver superior predictive performance and uncertainty quantification for non‑stationary simulators. The open‑source implementation further positions the method as a practical tool for the surrogate modeling community.