Bayesian Deep Gaussian Process

• Bayesian Deep Gaussian Process is a hierarchical model where layered GP mappings capture latent representations and propagate uncertainty.
• The methodology employs variational inference and sparse inducing points to handle non-Gaussian posterior complexities in multi-layer GP compositions.
• Applications range from optimization and emulation to structured output modeling, emphasizing nonstationary behavior and hierarchical feature learning.

A Bayesian Deep Gaussian Process (DGP) is a hierarchical probabilistic model formed by composing multiple Gaussian-process-distributed latent functions. In contrast to a shallow Gaussian process, which places a prior directly on a single input–output map, a DGP introduces latent representations layer by layer, so that each transformation is itself stochastic and uncertainty is propagated through the hierarchy. This construction is used to obtain hierarchical feature learning, input warping, and effectively nonstationary behavior while retaining a Bayesian treatment of functions and predictive uncertainty (Jakkala, 2021, Damianou et al., 2012).

1. Core formulation

The standard GP regression model places a prior

$f \sim \mathcal{GP}(m(\cdot), k(\cdot,\cdot)),$

with observation model

$$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$

A Bayesian DGP replaces this single latent function by a composition of GP-distributed mappings,

$$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$

with $\mathbf{h}^{(0)}=\mathbf{x}$ and the output generated from the final layer. Because the hidden outputs of one layer become the inputs of the next, the resulting model is generally not itself a GP, and its posterior is typically non-Gaussian and analytically intractable (Jakkala, 2021).

The original variational DGP formulation of Damianou and Lawrence is a latent-variable hierarchy in which observed outputs $Y$, intermediate latent layers $X_h$, and a top latent parent $Z$ are linked by successive GP mappings. In the unsupervised case emphasized there, the model evidence takes the form

$$\log p(Y) = \log \iint p(Y\mid X)\,p(X\mid Z)\,p(Z)\,dX\,dZ,$$

which is already intractable because latent variables enter nonlinearly inside covariance matrices. That paper used this Bayesian hierarchy together with a strict variational lower bound on the marginal likelihood for model selection over the number of layers and nodes per layer (Damianou et al., 2012).

A common supervised two-layer specification is

$$\mathbf{h}^{(1)} = f^{(1)}(\mathbf{x}) + \epsilon^{(1)}, \qquad \epsilon^{(1)} \sim \mathcal{N}\bigl(\mathbf{0},\sigma_1^2\mathbf{I}\bigr),$$

$$\mathbf{y} = f^{(2)}\bigl(\mathbf{h}^{(1)}\bigr) + \epsilon^{(2)}, \qquad \epsilon^{(2)} \sim \mathcal{N}\bigl(\mathbf{0},\sigma_2^2\mathbf{I}\bigr).$$

Here the first GP layer produces a latent embedding and the second GP layer maps that embedding to outputs. In the materials-design setting of a recent application, the concrete implementation used 2 layers, 32 inducing points per layer, and 10 latent GPs in the first layer, with the first layer acting as an input-dependent warping of the feature space (Alvi et al., 17 Sep 2025).

2. Bayesian inference and posterior approximation

Exact Bayesian inference in DGPs is difficult for two coupled reasons: latent outputs of one layer are inputs to the next, and the posterior over hidden representations becomes non-Gaussian. The standard workaround is sparse variational inference with inducing variables. In the survey formulation, a doubly stochastic variational objective is

$$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$0

The “doubly stochastic” terminology refers to stochasticity from minibatch optimization and from sampling through latent layers. This inference family became a major practical step because earlier variational schemes could force excessive independence across layers (Jakkala, 2021).

In concrete two-layer implementations, the evidence lower bound is often written as

$$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$1

with Gaussian inducing posteriors

$$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$2

This is the formulation used in recent DGP-based Bayesian optimization for materials design, where the model was trained with Adam, learning rate $$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$3, and early stopping after 50 epochs (Alvi et al., 17 Sep 2025).

Alternative approximate Bayesian routes have also been developed. One line uses stochastic expectation propagation and probabilistic backpropagation. In that setting, the posterior over inducing variables is approximated by a tied-factor Gaussian

$$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$4

and layerwise Gaussian moment propagation is used to approximate the tilted normalizer needed for EP-style updates. This provides a minibatch-compatible Bayesian training scheme with FITC sparsification and $$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$5 site memory under SEP rather than $$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$6 under standard EP (Bui et al., 2015).

Another line replaces explicit variational densities by implicit posterior samplers. Neural Operator Variational Inference (NOVI) defines an implicit posterior $$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$7 over inducing variables through a neural generator and trains it by minimizing a Regularized Stein Discrepancy against the true DGP posterior. With the optimal discriminator, the objective becomes

$$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$8

so the method reduces to minimizing a scaled Fisher divergence. The same paper gives a predictive bias bound

$$y = f(x) + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\beta^{-1}).$$9

This suggests a more expressive posterior family than mean-field Gaussian approximations, at the cost of a minimax optimization problem (Xu et al., 2023).

3. Architectural variants and structured-output extensions

Bayesian DGPs have been extended by altering either the layerwise kernel structure or the latent architecture. In computer vision, Convolutional DGPs replace global kernels by convolutional kernels. If an image $$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$0 is decomposed into patches $$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$1, the image-level function is

$$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$2

which induces covariance

$$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$3

A weighted variant introduces patch weights $$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$4,

$$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$5

This imports local spatial structure and some translation robustness into the deep GP hierarchy while retaining sparse variational inference as the training mechanism (Kumar et al., 2018).

For sequential data, utterance-level SRU-DGP models replace the affine transformations in a Simple Recurrent Unit by GP-distributed latent functions,

$$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$6

with analogous GP replacements for the reset gate and hidden-state update. The recurrent dependence remains in the lightweight SRU state update, while GP regression stays parallelizable across frames. That work also reformulated the ELBO at the utterance level and used low-rank random-feature approximations for whole-utterance covariance sampling (Koriyama et al., 2020).

Several extensions target structured outputs rather than only structured inputs. For transient multivariate outputs, a deep matrix-variate Gaussian process emulator replaces independent-output DGP layers by matrix-normal layers,

$$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$7

$$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$8

The output covariance $$\mathbf{h}^{(l)} = f^{(l)}(\mathbf{h}^{(l-1)}), \qquad f^{(l)} \sim \mathcal{GP}(m^{(l)}, k^{(l)}),$$9 captures dependence across output coordinates, which in that contaminant-localization application were time points of a transient concentration curve (Park et al., 2018).

For multi-fidelity regression, conditional DGPs use low-fidelity data to define the latent lower layer directly: $\mathbf{h}^{(0)}=\mathbf{x}$0 Under an outer squared-exponential kernel, moment matching yields an effective high-fidelity kernel

$\mathbf{h}^{(0)}=\mathbf{x}$1

making uncertainty propagation from the lower-fidelity posterior explicit in the induced kernel (Lu et al., 2020).

Other structured approximations include GP-DRF, which combines an exact GP input layer with deep random-feature layers so that trees, graphs, and sequences can be mapped first into a fixed-dimensional latent vector and then processed by a scalable DGP-like stack (Laradji et al., 2019). A complementary weight-space view shows that squared-exponential-kernel DGPs can be translated into deep bottlenecked trigonometric networks; in the wide limit, these reproduce the same effective covariance functions as the corresponding function-space DGPs and expose a neural tangent kernel for the model class (Lu et al., 2022).

4. Applications in optimization, emulation, and scientific computing

A major use of Bayesian DGPs is as uncertainty-aware surrogates inside sequential design and emulator pipelines. In cost-aware batch Bayesian optimization for complex materials design, a DGP surrogate is used inside a loop that alternates between cheap heterotopic single-task queries and expensive isotopic multi-objective batches. The single-objective acquisition uses

$\mathbf{h}^{(0)}=\mathbf{x}$2

and the batch multi-objective step uses

$\mathbf{h}^{(0)}=\mathbf{x}$3

In that refractory high-entropy alloy campaign, the DGP surrogate was described as hierarchical and multi-output, and the overall framework converged to optimal formulations in fewer iterations with cost-aware queries than conventional GP-based BO (Alvi et al., 17 Sep 2025).

In likelihood-free inference, DGPs have been used as discrepancy surrogates in BOLFI when the target surface is multimodal, non-stationary, or heteroscedastic. The model there is a latent-variable DGP with one latent-variable layer and two GP layers, and Bayesian optimization is adapted to its non-Gaussian predictive uncertainty by quantile-conditioned moments

$\mathbf{h}^{(0)}=\mathbf{x}$4

leading to the acquisition

$\mathbf{h}^{(0)}=\mathbf{x}$5

On multimodal targets, these DGP surrogates outperformed ordinary GP surrogates, while remaining comparable on unimodal cases (Aushev et al., 2020).

Multi-fidelity engineering optimization is another recurrent setting. For simulated chemical reactors, a five-layer DGP was used as an end-to-end surrogate across five discrete mesh fidelities $\mathbf{h}^{(0)}=\mathbf{x}$6. Design selection used a highest-fidelity UCB rule,

$\mathbf{h}^{(0)}=\mathbf{x}$7

and fidelity selection used cost-weighted uncertainty,

$\mathbf{h}^{(0)}=\mathbf{x}$8

The paper emphasized uncertainty propagation across fidelities and the budgetary advantage of mixing cheap and expensive simulations (Savage et al., 2022).

For large computer experiments, fully Bayesian DGP inference has been combined with Vecchia approximations at every GP layer. A two-layer Vecchia-DGP takes the form

$\mathbf{h}^{(0)}=\mathbf{x}$9

and reduces cubic GP costs to $Y$0 for fixed neighbor size $Y$1. That work reported practical fits up to $Y$2 and argued that posterior integration and uncertainty quantification could be preserved at that scale (Sauer et al., 2022).

Bayesian DGPs have also been adapted to correlated functional outputs in cosmology. In that setting, the latent true matter power spectrum $Y$3 is modeled through a warping layer $Y$4,

$Y$5

so that the DGP prior accounts for nonstationarity through stochastic monotone input deformation and the observation layer accounts for correlated multi-fidelity functional realizations (Walsh et al., 24 Jul 2025).

5. Theory, model selection, and posterior contraction

A defining feature of the original Bayesian DGP program is that approximate marginal likelihoods are used not only for training but also for model selection. The 2012 formulation derived a strict lower bound on the marginal likelihood and used it to compare architectures by the number of layers and latent dimensions. In that paper, the variational evidence was reported to justify a five-layer hierarchy even on a handwritten-digit dataset with only 150 examples, illustrating the model’s intended “automatic Occam’s razor” over depth and width (Damianou et al., 2012).

PAC-Bayesian theory has supplied a complementary interpretation of variational DGP training. For a class of variational DGPs, the paper on PAC-Bayesian bounds shows that choosing the PAC-Bayes prior and posterior to match the variational prior and variational posterior yields

$Y$6

so maximizing the variational lower bound is equivalent to minimizing a PAC-Bayesian empirical objective at $Y$7. The same work proves consistency of order $Y$8 and gives an oracle-type inequality comparing the variational predictor with the raw DGP predictor (Föll et al., 2019).

In nonlinear inverse problems, DGP priors have been used to formalize adaptation to compositional structure. The prior there is hierarchical: $Y$9 with $X_h$0 indexing graph structure and layerwise smoothness, and latent parameter draws of the form

$X_h$1

Under regularity, Lipschitz, and stability assumptions on the forward operator, the posterior induced by this DGP prior contracts around a compositional truth at rates determined by the intrinsic dimensions $X_h$2 and smoothnesses $X_h$3 of the latent layers. In Darcy’s problem, the same paper proves a lower bound showing that ordinary Whittle–Matérn GP priors contract at a polynomially slower rate than the DGP prior for certain generalized additive truths, making the structural difference between shallow and deep priors explicit (Abraham et al., 2023).

6. Strengths, limitations, and recurring misconceptions

The main strengths attributed to Bayesian DGPs are hierarchical feature learning, input warping, varying length scales, multi-scale structure, non-Gaussian predictive behavior, and uncertainty propagation through successive latent layers. These properties motivate their use when shallow kernels are too restrictive, especially in the presence of nonstationarity, correlated tasks, or latent compositional structure (Jakkala, 2021, Alvi et al., 17 Sep 2025).

Several limitations are equally well established. Inference is hard because hidden variables are nested through multiple GP layers and exact marginalization is unavailable. Approximate inference can be restrictive: Gaussian variational families may miss multimodal posterior structure, and sampling-based alternatives such as SGHMC are computationally demanding and difficult to tune. Early variational DGPs could also suffer degeneracies in which some layers effectively “get turned off” (Jakkala, 2021).

A common misconception is that DGPs uniformly dominate shallow GPs. The empirical record in the materials-design study is more nuanced: DGP variants helped particularly on high-dimensional, highly nonlinear, or nonconvex problems with complex latent structure, but shallow GPs remained competitive on smoother, moderate-dimensional cases, and final hypervolume could be similar on relatively smooth CALPHAD-derived landscapes (Alvi et al., 17 Sep 2025). A related misconception is that deeper uncertainty propagation automatically implies better calibrated uncertainty. The same study reports that “DGP variance estimation through doubly stochastic inference can suffer from poor calibration and underestimation of uncertainty,” motivating hybrid variants that combine DGP means with GP or multitask-GP variance estimates (Alvi et al., 17 Sep 2025).

Open technical issues recur across the literature. Output-dimension factorization is still common even when outputs are correlated. Initialization and hyperparameter tuning remain delicate. Richer posterior approximations beyond Gaussian variational families are still being explored. Scaling to highly structured datasets and to much deeper architectures remains unresolved. The survey also notes that DGP studies usually consider up to around ten layers, far shallower than the depths used in state-of-the-art vision networks (Jakkala, 2021).

Taken together, these results suggest a balanced characterization. A Bayesian DGP is neither merely a deeper GP nor simply a kernelized analogue of a neural network. It is a hierarchical prior over compositions of latent functions whose appeal lies in uncertainty-aware representation learning, but whose practical value depends strongly on inference quality, architectural choice, and the structural complexity of the target problem.