---
title: GradInf for Deep Gaussian Processes
url: https://www.emergentmind.com/topics/gradinf
type: topic
---

# GradInf for Deep Gaussian Processes

Searching arXiv for the primary paper and closely related uses of “GradInf” to ground the article in current literature.
GradInf, in the context of Gaussian-process-based surrogate modeling, denotes gradient information in Gaussian-process-based surrogate modeling. In “Deep Gaussian Processes with Gradients” [2512.18066], the term is associated with a full Bayesian framework for deep Gaussian processes (DGPs) that supports both gradient-enhancement and gradient posterior predictive distributions. The framework is designed for expensive nonstationary computer experiments, where one or more latent Gaussian processes warp the input space into a plausibly stationary regime before outer Gaussian-process regression is applied, and it is implemented in the `deepgp` package on CRAN with optional Vecchia approximation to circumvent cubic computational bottlenecks.

## 1. Motivation and problem setting

Many computer simulators are expensive and nonstationary: correlation structure and local smoothness change dramatically across the input space, including ignition cliffs, plateaus, and phase transitions [2512.18066]. A stationary kernel in a single GP imposes a single correlation and smoothness scale everywhere, which is often mismatched. DGPs address this by composing Gaussian processes so that latent warpings transform the original inputs into a domain where a stationary outer regression is more plausible.

Gradient information is increasingly available via adjoint solvers and automatic differentiation. In the single-layer GP setting, gradients can sharpen learning of local curvature, improve sample efficiency, and enhance uncertainty quantification, especially when function evaluations are scarce. In the DGP setting, gradients additionally constrain both the latent warping and the outer regression, reducing ambiguity in the composition and helping capture steep regions without inflating variance elsewhere [2512.18066].

The reported advantages are threefold. First, gradient data improve sample efficiency because derivatives constrain the latent warping and the outer regression simultaneously. Second, derivative cross-covariances tie local slope information to function predictions, lowering posterior variance where gradients are informative and yielding better calibrated intervals. Third, under deep composition the gradient of the map responds to both outer and inner layer behavior via the chain rule, which is particularly relevant for complex nonstationarity [2512.18066].

## 2. Single-layer Gaussian processes with gradients

The single-layer construction begins with a zero-mean Gaussian process with a smooth kernel:
\[
f \sim GP(0, k), \qquad x, x' \in \mathbb{R}^D.
\]
For a twice differentiable kernel, the key covariance identities are [2512.18066]
\[
\mathrm{Cov}(f(x), f(x')) = k(x,x'),
\]
\[
\mathrm{Cov}\!\left(f(x), \frac{\partial f(x')}{\partial x'_j}\right) = \frac{\partial k(x,x')}{\partial x'_j},
\]
\[
\mathrm{Cov}\!\left(\frac{\partial f(x)}{\partial x_i}, \frac{\partial f(x')}{\partial x'_j}\right) = \frac{\partial^2 k(x,x')}{\partial x_i \partial x'_j}.
\]

At training inputs \(X=[x_1,\dots,x_n]\), the stacked prior over function values and all partial derivatives is
\[
\begin{bmatrix}
f(X)\\
\nabla f(X)
\end{bmatrix}
\sim
N\!\left(
0,\;
\tau^2 [K_{\bullet\bullet}(X)+\epsilon I]
\right),
\]
where \(K_{\bullet\bullet}(X)\) is the block covariance obtained by differentiating \(k\) appropriately, with blocks \(K_{00}\), \(K_{0d}\), \(K_{d0}\), and \(K_{df}\), and \(\epsilon\) is a small diagonal jitter used to stabilize numerics even for deterministic simulators [2512.18066].

The gradient-enhanced likelihood augments function observations with derivative observations:
\[
y_i = f(x_i) + \epsilon_{yi}, \qquad \epsilon_{yi}\sim N(0,\sigma_y^2),
\]
\[
g_i = \nabla f(x_i) + \epsilon_{gi}, \qquad \epsilon_{gi}\sim N(0,\Sigma_g),
\]
with \(\Sigma_g\) typically diagonal, using per-component gradient noise \(\sigma_{g,d}^2\). For deterministic experiments, the framework fixes a small jitter \(\epsilon\) across function and gradient blocks; for noisy gradients arising from numerical differentiation or adjoint tolerances, setting \(\sigma_y^2\) and \(\sigma_{g,d}^2>0\) is described as prudent [2512.18066].

Kernel smoothness determines which derivative quantities are admissible. The squared exponential kernel is infinitely differentiable and yields clean gradient covariances. Matérn-\(5/2\) is twice differentiable and valid for gradient modeling, whereas Matérn-\(3/2\) is once differentiable, so gradients are admissible but Hessians are not. The stated rule is that at least \(C^1\) smoothness is required to include gradients, and at least \(C^2\) smoothness to include gradient–gradient covariances [2512.18066].

For test inputs \(\mathcal{X}\), the posterior predictive distribution over functions and gradients is Gaussian:
\[
\begin{bmatrix}
f(\mathcal{X})\\
\nabla f(\mathcal{X})
\end{bmatrix}
\mid (X,y)
\sim N(\mu^\*,\Sigma^\*),
\]
with
\[
\mu^\* = K_{\bullet 0}(\mathcal{X},X)\,[K_{00}(X)+\epsilon I]^{-1}y,
\]
\[
\Sigma^\* =
\tau^2\!\left[
K_{\bullet\bullet}(\mathcal{X})
-
K_{\bullet 0}(\mathcal{X},X)\,[K_{00}(X)+\epsilon I]^{-1}K_{0\bullet}(X,\mathcal{X})
\right].
\]
When gradient observations are available, \(K_{00}(X)\) and \(y\) are replaced by \(K_{\bullet\bullet}(X)\) and \([y;\nabla y]\), respectively [2512.18066].

## 3. Deep composition and chain-rule propagation

The DGP construction considered for surrogate modeling is a two-layer composition
\[
z_1(x) = f_1(x), \qquad y(x) = f_2(z_1(x)).
\]
The reported focus is a one-latent-layer DGP with \(D\) latent nodes, one per input dimension, independent given \(X\):
\[
w_d = f_{1,d}(X), \qquad f_{1,d}\sim GP(0,k_{1,d}), \qquad d=1,\dots,D,
\]
\[
y\mid W \sim GP(0,\tau^2 k_2(W)),
\]
with \(W=[w_1,\dots,w_D]\). Each layer has its own lengthscales, and the outer layer has scale \(\tau^2\) [2512.18066].

The key mechanism is gradient propagation through the latent warping. For
\[
y(x)=f_2(W(x)),
\]
the gradient obeys
\[
\nabla y(x)=J_{f_2}(W(x))\cdot J_W(x),
\]
where \(J_{f_2}\) is the \(1\times D\) Jacobian \([\partial y/\partial w^1,\dots,\partial y/\partial w^D]\) and \(J_W(x)\) is the \(D\times D\) Jacobian of \(W(x)\), with entries \(\partial w_i/\partial x^d\). For \(L\) layers, the full chain rule is
\[
\nabla y(x)=J_{f_L}(z_{L-1}(x))\cdot J_{f_{L-1}}(z_{L-2}(x))\cdots J_{f_1}(x).
\]
This is the formal basis for using gradient information to constrain both the inner warping and the outer response [2512.18066].

The joint model separates inner-layer gradients with respect to \(X\) from outer-layer gradients with respect to \(W\). The inner layer jointly models \([w_d;\partial w_d/\partial x^1;\dots;\partial w_d/\partial x^D]\) via the GP-with-gradients prior. The outer layer jointly models \([y;\partial y/\partial w^1;\dots;\partial y/\partial w^D]\) via a GP-with-gradients over the warped inputs \(W\). Observed output gradients with respect to \(X\) are linked to the outer-layer gradient variables through the linear system implied by the chain rule:
\[
\begin{bmatrix}
\partial w_1/\partial x^1 & \cdots & \partial w_D/\partial x^1\\
\vdots & & \vdots\\
\partial w_1/\partial x^D & \cdots & \partial w_D/\partial x^D
\end{bmatrix}
\cdot
\begin{bmatrix}
\partial y/\partial w^1\\
\vdots\\
\partial y/\partial w^D
\end{bmatrix}
=
\begin{bmatrix}
\partial y/\partial x^1\\
\vdots\\
\partial y/\partial x^D
\end{bmatrix}.
\]
Thus, given \(\nabla_x W\) and observed \(\nabla_x y\), one solves for \(\nabla_w y\); given \(\nabla_w y\) and \(\nabla_x W\), one predicts \(\nabla_x y\) [2512.18066].

## 4. Bayesian inference framework

The prior specification uses zero means for both layers, with identity-mean on the inner layer as an optional alternative. Smooth kernels are required, with Gaussian and Matérn-\(5/2\) given as the main choices. Distinct lengthscales are assigned per inner node and one lengthscale to the outer layer. The outer scale \(\tau^2\) uses the reference prior
\[
\pi(\tau^2)\propto \frac{1}{\tau^2},
\]
which is analytically integrated in the Gaussian likelihood. Noise terms \(\sigma_y^2\) and \(\Sigma_g\) are optional; deterministic simulators default to jitter \(\epsilon\) [2512.18066].

Posterior inference proceeds by sampling the latent warpings \(W\) and their gradients \(\nabla_x W\) via Elliptical Slice Sampling (ESS) from the inner-layer GP-with-gradients prior, one latent node at a time in a Gibbs schedule. For the \(d\)-th node, the proposal over the node and its gradient stack is
\[
w_{d,\mathrm{all}}^{(t)}
=
w_{d,\mathrm{all}}^{(t-1)}\cos(\gamma)
+
w_{d,\mathrm{all}}^\star\sin(\gamma),
\qquad
\gamma\in[0,2\pi],
\]
with shrinkage until acceptance. Acceptance uses the outer-layer gradient-enhanced likelihood, where \(\nabla_w y\) is constructed from the current \(\nabla_x W\) sample and the observed \(\nabla_x y\) through the chain rule. Lengthscales are updated by Metropolis–Hastings with Gamma priors and sliding-window proposals, interleaved with ESS in the Gibbs loop. Variational inference is explicitly not used; the paper favors full Bayesian MCMC because of performance issues observed with some variational and doubly-stochastic approximations for surrogates [2512.18066].

Gradient posterior prediction in the DGP is described as a four-step procedure at each MCMC iteration. First, infer \([W(x^\*),\nabla_x W(x^\*)]\) by GP-with-gradients conditioning on \(X\) and the sampled inner-layer draws. Second, infer \([y(x^\*),\nabla_w y(x^\*)]\) by GP-with-gradients conditioning on \(W\) and the observed \(y\). Third, obtain \(\nabla_x y(x^\*)\) via the chain rule. Fourth, aggregate across MCMC iterations to obtain posterior means and variances. The paper also gives a stochastic-imputation moment propagation formula:
\[
E\!\left[\frac{\partial y}{\partial x^d}\right]
=
\sum_{i=1}^D
\left(\frac{\partial w_i}{\partial x^d}\right)
E\!\left[\frac{\partial y}{\partial w^i}\right],
\]
\[
\mathrm{Var}\!\left[\frac{\partial y}{\partial x^d}\right]
=
\sum_{i=1}^D
\left(\frac{\partial w_i}{\partial x^d}\right)^2
\mathrm{Var}\!\left[\frac{\partial y}{\partial w^i}\right],
\]
followed by averaging across iterations and adding between-iteration variance through the law of total variance [2512.18066].

Identifiability is treated as a benign but unavoidable issue. Latent warpings are not unique; reflections or monotone transforms that preserve pairwise latent distances can yield similar outer-layer likelihoods. ESS samples may therefore differ by signs or shifts. The stated interpretation is that this is benign for the outer GP because it depends on pairwise distances. Gradient data regularize the latent map by informing local stretching and compression, but overfitting can occur if gradient noise is underestimated [2512.18066].

## 5. Computation, Vecchia approximation, and software

A central computational difficulty is that naive GP likelihoods and conditioning scale cubically, and gradient enhancement increases the effective observation count from \(n\) to
\[
N = n(1 + D).
\]
For DGPs, this burden is duplicated across inner and outer layers and across many MCMC iterations [2512.18066].

The reported remedy is the Vecchia approximation, which factorizes the Gaussian density into low-order conditionals:
\[
\log p(y)\approx \sum_{i=1}^{N}\log p(y_i\mid y_{c_i}),
\qquad
c_i\subset \{1,\dots,i-1\},\quad |c_i|\le m.
\]
With an ordering and nearest-neighbor conditioning sets, the resulting sparse precision representation yields \(O(Nm^2)\) time and \(O(Nm)\) memory per layer, with \(m\) controlling the accuracy–speed trade-off [2512.18066].

For function-plus-gradient observations, the details of ordering and conditioning are part of the method. Responses are ordered first, followed by derivatives per point, for example \(y_1,\dots,y_n,\partial y_1/\partial x^1,\dots,\partial y_1/\partial x^D,\partial y_2/\partial x^1,\dots\). Conditioning uses nearest neighbors in input space, and ties are broken by preferring response values over derivatives. The same Vecchia machinery is applied to the inner-layer nodes with gradients, the outer-layer GP over \(W\) with gradients, and predictive conditionals, thereby accelerating ESS, likelihood evaluation, and posterior prediction, including joint gradient prediction [2512.18066].

The `deepgp` package on CRAN operationalizes these ideas. The reported capabilities are Bayesian DGPs with one latent layer, per-node GP priors, gradient-enhanced training for both inner and outer layers, gradient posterior prediction at test inputs via chain-rule propagation, optional Vecchia approximation for all layers and operations, full MCMC for latent warping by ESS and lengthscales by Metropolis–Hastings, jitter \(\epsilon\) for stability, and identity or zero means on the inner layer [2512.18066].

## 6. Empirical behavior

The empirical study uses nonstationary simulations with 30 Monte Carlo repetitions each and Latin hypercube designs. The three benchmarks are Squiggle, Plateau, and Ignition, chosen to represent strong nonstationarity under increasing dimension and sample size [2512.18066].

| Benchmark | Setting | Reported outcome |
|---|---|---|
| Squiggle | \(D=2,\; n=25\) | geDGP yields substantially lower RMSE and CRPS than geGP |
| Plateau | \(D=3,\; n=30\) | geDGP consistently best; DGP without gradients can beat GP with gradients |
| Ignition | \(D=6,\; n=100\) | geDGP best overall in RMSE and CRPS; Vecchia used throughout |

On Squiggle, described as strong nonstationarity with an “S”-shaped ridge, DGPs outperform GPs for \(y\) and uncertainty quantification, and gradient enhancement improves both. The paper states that geDGP yields substantially lower RMSE and CRPS than geGP, that DGP gradient predictions without gradient training data improve over GP gradient predictions, and that GEK, an RFF-based stationary GP with gradients, underperforms geDGP [2512.18066].

On Plateau, a \(D=3\) surface with flat regions separated by a steep drop, DGPs strongly outperform GPs. The paper highlights that DGP without gradients can beat GP with gradients, emphasizing the importance of nonstationarity modeling. geDGP is reported as consistently best, and occasional outliers in gradient prediction are said not to compromise function prediction quality [2512.18066].

On Ignition, a \(D=6\) problem with a steep ignition cliff, Vecchia is used throughout because \(N\) reaches up to 700 for gradient-enhanced runs. Here again DGP predicts gradients better than stationary geGP, and geDGP is best overall in RMSE and CRPS [2512.18066].

The stated empirical takeaway is that the largest gains from gradient data in DGPs occur under strong nonstationarity, limited function samples, and moderate-to-high dimensions where warping plus slope information disentangle local behavior. The same section reports that gradient-enhanced DGPs deliver better calibrated uncertainty quantification, as reflected in lower CRPS, than the alternatives considered [2512.18066].

## 7. Assumptions, limitations, and broader uses of the term

The framework assumes smooth kernels, specifically \(C^2\) smoothness for gradient–gradient covariances, and deterministic simulators with small jitter \(\epsilon\), although gradient noise can be incorporated when calibrated carefully. The principal limitations identified are deep-warping identifiability, overfitting risks when gradient noise is underestimated, and residual scalability constraints even after reduction to \(O(Nm^2)\), especially when \(nD\) is large and Jacobian products become unstable [2512.18066].

The paper also identifies open challenges. These include active learning with gradients, in which function and gradient points are acquired jointly to maximize information about latent warpings and sharp features; hybrid inducing–Vecchia schemes for massive datasets; PDE-constrained surrogates that incorporate physics-informed gradient structure or adjoint outputs directly into multi-layer GP priors; and uncertainty-aware gradient exploitation in downstream optimization [2512.18066].

The label “GradInf” also appears in other research domains. “GradInf: Gradient Estimation as Probabilistic Inference” uses the name for a formal reduction from gradient estimation to probabilistic inference, implemented as a probabilistic programming system based on coupling, factorization, inference, and automatic differentiation [2607.07840]. In causal discovery, “Trust Your \(\nabla\): Gradient-based Intervention Targeting for Causal Discovery” states that GIT embodies the GradInf principle because intervention targets are selected from gradient signals of the causal structure loss [2211.13715]. In graph explainability, “Graph-based Integrated Gradients for Explaining Graph Neural Networks” presents GB-IG as a gradient-based influence method tailored to GNNs [2509.07648]. In “Inference in Graded Bayesian Networks”, graded inference is presented as a generalized Viterbi algorithm for graded Bayesian networks [1901.01837]. This broader usage suggests that GradInf functions less as a single canonical doctrine than as a recurring label for methods that use gradients as primary carriers of information, but in the DGP literature its most specific meaning is gradient information in Gaussian-process-based surrogate modeling [2512.18066].

Source: https://www.emergentmind.com/topics/gradinf