Deep Gaussian Processes

• DeepGP is a hierarchical composition of Gaussian process mappings where each layer warps inputs to capture increasingly abstract, non-Gaussian features.
• It employs approximate Bayesian inference and scalable sparse approximations to overcome intractable exact posterior computations.
• DeepGP models are applied in areas like geophysical retrieval, image processing, and dynamic system identification, with variants addressing multitask and convolutional challenges.

DeepGP most commonly denotes a deep Gaussian process (DGP): a hierarchical composition of Gaussian-process mappings in which the output of one GP layer becomes the input to the next, yielding a deep probabilistic nonparametric model that combines multilayer expressiveness with GP-style uncertainty quantification (Damianou et al., 2012, Havasi et al., 2018). In the original formulation, a single layer is equivalent to a standard GP or the GP latent variable model, depending on whether inputs are observed or latent, while multiple layers create latent warpings and increasingly abstract representations (Damianou et al., 2012). Across later literature, DeepGPs are used to model nonstationary, hierarchical, modular, and multi-output structure, but exact inference remains intractable, so the field is largely organized around approximate Bayesian inference and scalable sparse approximations (Bui et al., 2015, Sauer et al., 2022).

1. Origins and conceptual basis

The term emerged in the 2012 paper "Deep Gaussian Processes" (Damianou et al., 2012), which introduced deep GPs as a deep belief network based on Gaussian-process mappings. The model places observed data $Y$ at the leaves, intermediate latent layers $X_1,\dots,X_{H-1}$ in the hierarchy, and a top latent layer $Z$, with a chain of GP mappings

$$Y \leftarrow X_1 \leftarrow X_2 \leftarrow \cdots \leftarrow X_{H-1} \leftarrow Z.$$

This construction generalizes both ordinary GP regression and the GP-LVM: a standard GP corresponds to a single mapping from observed inputs to outputs, whereas a GP-LVM uses one latent input layer. DeepGP stacks such mappings recursively (Damianou et al., 2012).

The original paper framed the model as fully Bayesian and emphasized approximate variational marginalization rather than direct MAP optimization of latent variables. A strict lower bound on the marginal likelihood was used for model selection over depth and latent dimensionality. In the digit experiment discussed there, this criterion selected a five-layer hierarchy on a dataset containing only 150 examples, illustrating the paper’s central claim that deep probabilistic structure can be justified even in low-data regimes (Damianou et al., 2012).

A recurrent conceptual theme in later work is the analogy to deep neural networks. One line of work states that a GP is equivalent to a single hidden-layer neural network with infinitely many hidden units, and a DGP is the corresponding multi-layer analogue (Bui et al., 2015). Another line emphasizes that random-feature reformulations make a DGP look like a Bayesian deep neural network with structured weights, while retaining GP-derived priors and uncertainty propagation (Cutajar et al., 2016). This suggests that DeepGP occupies a hybrid position between classical kernel methods and hierarchical representation learning.

2. Canonical model and induced properties

A canonical DGP is a composition of latent GP layers. One common notation writes

$$\mathbf F^0 = \mathbf X,\qquad \mathbf F^l \sim \mathrm{GP}(\mathbf F^{l-1}),\qquad l=1,\dots,L,$$

while another expresses the output model as

$$\mathbf{y}_n \mid f^L; \mathbf{x}_n \sim p\!\left(\mathbf{y}_n \mid f^L(f^{L-1}(\dots f^1(\mathbf{x}_n)))\right),$$

with each layer prior

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

These formulations describe the same object: a hierarchy of nonlinear stochastic transformations in which uncertainty is propagated through depth (Havasi et al., 2018, Boustati et al., 2019).

Relative to a shallow GP, the induced prior over functions is non-Gaussian and substantially more expressive. The literature repeatedly attributes three advantages to this composition: automatic representation learning or kernel design through latent input warping, richer non-Gaussian function priors than a single GP, and the ability to represent complicated modular or hierarchical processes by stacking simpler GP modules (Bui et al., 2015, Svendsen et al., 2020). In geophysical retrieval, for example, DGPs are motivated precisely because radiative-transfer-type relations are described as nonlinear, hierarchical, and modular rather than well captured by a single shallow kernel (Svendsen et al., 2020).

Several works also stress that a DGP can act as a nonstationary model even when each layer uses a stationary kernel. The mechanism is the latent warping: the covariance of the outer GP is evaluated on learned latent coordinates rather than on the original inputs. In active-learning work for simulator surrogates, this is described as automatic warping of the input space, allowing the model to distribute modeling capacity differently across regions and to represent abrupt regime changes more effectively than a stationary GP (Sauer et al., 2020).

3. Inference, variational objectives, and scalability

The main technical obstacle is that exact posterior inference in a DGP is generally intractable because latent variables from one layer enter the kernel evaluations of later layers. The earliest DeepGP paper handled this with inducing variables, Gaussian variational distributions over latent layers, and a strict lower bound on the marginal likelihood (Damianou et al., 2012). Subsequent work developed several inference regimes rather than a single consensus method.

One influential route combined sparse GP layers with stochastic Expectation Propagation and probabilistic backpropagation. In "Training Deep Gaussian Processes using Stochastic Expectation Propagation and Probabilistic Backpropagation" (Bui et al., 2015), the posterior approximation is written as $q(u)\propto p(u)g(u)^N$, replacing per-datapoint EP factors by a single shared factor and reducing memory cost to $\mathcal{O}(M^2)$. The method uses FITC approximations, Gaussian moment propagation through layers, and minibatch Adam optimization.

A second route made DGPs practical at much larger scale by replacing each GP layer with explicit random-feature expansions and training with stochastic variational inference. "Random Feature Expansions for Deep Gaussian Processes" (Cutajar et al., 2016) reports experiments on datasets with up to 8 million observations and architectures with up to 30 hidden layers. On the large MNIST setup it reports $99.14\%$ accuracy with the RBF version, and on the airline dataset it reports $X_1,\dots,X_{H-1}$0 accuracy and $X_1,\dots,X_{H-1}$1 MNLL for the RBF model (Cutajar et al., 2016). This work shifts DGP computation from kernel-matrix decompositions toward matrix products and reparameterized stochastic optimization.

A third line improves sparse variational efficiency rather than changing the representation. "Deep Gaussian Processes with Decoupled Inducing Inputs" (Havasi et al., 2018) uses separate inducing sets for the mean and variance of each layer, with $X_1,\dots,X_{H-1}$2 for the mean and $X_1,\dots,X_{H-1}$3 for the variance. The paper compares a standard DGP with $X_1,\dots,X_{H-1}$4 to a decoupled DGP with $X_1,\dots,X_{H-1}$5, $X_1,\dots,X_{H-1}$6, and reports that the decoupled model is faster for models with at least one hidden layer while achieving better test log-likelihood and lower RMSE on almost all settings (Havasi et al., 2018).

A fourth route targets posterior expressiveness. "Implicit Posterior Variational Inference for Deep Gaussian Processes" (Yu et al., 2019) argues that mean-field Gaussian posteriors across layers are restrictive and biased, and replaces the explicit variational posterior with an implicit generator trained adversarially. The paper casts DGP inference as a two-player game and states that, at Nash equilibrium, the learned implicit posterior coincides with an unbiased posterior belief. Empirically it reports posterior sample generation about $X_1,\dots,X_{H-1}$7 faster and training about $X_1,\dots,X_{H-1}$8 faster than SGHMC in the Airline experiment (Yu et al., 2019).

For fully Bayesian surrogate modeling, "Vecchia-approximated Deep Gaussian Processes for Computer Experiments" (Sauer et al., 2022) applies the Vecchia approximation at every GP layer, using sparse local conditionals to reduce dense cubic costs to roughly $X_1,\dots,X_{H-1}$9. This paper explicitly contrasts a "big data" regime of optimization-based approximate inference with a "small data" regime of posterior integration, and uses the Vecchia construction to bridge them while preserving uncertainty quantification. Its satellite-drag case study scales to $Z$0 training runs (Sauer et al., 2022).

4. Structured and specialized DeepGP architectures

Once the canonical composition was established, the literature diversified into architectures that impose additional structure on sharing, locality, frequency content, or output dependence. These variants retain the basic GP-composition principle but alter either the layer semantics or the inducing representation.

For multitask learning, "Non-linear Multitask Learning with Deep Gaussian Processes" (Boustati et al., 2019) introduces a two-layer DGP in which the latent layer is split into shared processes $Z$1 and task-specific processes $Z$2, concatenated into a task representation $Z$3, and mapped to outputs by a task-specific GP. The paper proposes two segmentation strategies: hard-coded sharing in a multiprocess MTL-DGP, and soft sharing via ARD kernels whose learned relevance weights determine how strongly each task uses shared versus private latent dimensions. The stated motivation is to generalize linear multitask models such as ICM, LMC, and SLFM into nonlinear latent sharing (Boustati et al., 2019).

For multi-output surrogates, "Deep Intrinsic Coregionalization Multi-Output Gaussian Process Surrogate with Active Learning" (Chang et al., 22 Aug 2025) embeds ICM covariance structure at every layer, giving a two-layer deepICMGP in which both the latent layer and the output layer are coregionalized. The paper argues that this permits nonlinear and structured dependencies between outputs that ordinary multi-output GPs or scalar-output DGPs do not capture directly (Chang et al., 22 Aug 2025).

For image data, "^^^^10^^^^" (Blomqvist et al., 2018) replaces ordinary dense GP layers with convolutional GP layers that operate on local patches and build hierarchical feature maps. On MNIST it reports $Z$4 accuracy for a 2-layer DeepCGP and $Z$5 for a 3-layer DeepCGP. On CIFAR-10 it reports $Z$6 and $Z$7 for the 2-layer and 3-layer versions, respectively, and highlights an improvement of over 10 percentage points relative to prior GP approaches (Blomqvist et al., 2018).

For globally structured but nonstationary data, "Inter-domain Deep Gaussian Processes" (Rudner et al., 2020) replaces inducing-point variables with inter-domain inducing variables given by RKHS Fourier features. The paper’s claim is that local inducing points are inefficient for broad global structure, while shallow inter-domain GPs remain limited on nonstationary data; combining inter-domain features with deep composition addresses both issues (Rudner et al., 2020).

For an empirical-Bayes alternative to standard latent-zero-mean DGPs, "Conditional Deep Gaussian Processes: empirical Bayes hyperdata learning" (Lu et al., 2021) conditions intermediate GPs on learned hyperdata $Z$8, uses moment matching to obtain an effective kernel, and states equivalence with deep kernel learning in the limit of dense hyperdata in latent space. The paper presents this as a more Bayesian counterpart to deep kernel learning while also arguing that standard variational DGPs can exhibit posterior collapse or diminished latent variance (Lu et al., 2021).

5. Applications, theoretical analysis, and generalized descendants

DeepGPs have been deployed in domains where hierarchical or nonstationary structure is explicit. In geophysical parameter retrieval, "Deep Gaussian Processes for geophysical parameter retrieval" (Svendsen et al., 2020) uses DGPs to estimate surface dew point temperature from IASI infrared sounding data. The model classes DGP2–DGP4 are reported to outperform the full GP in predictive accuracy, while FITC and DGP1 trail the full GP, and the paper interprets the deeper gains as evidence that the retrieval problem has meaningful hierarchical structure (Svendsen et al., 2020).

For sequential design and simulator surrogates, "Active Learning for Deep Gaussian Process Surrogates" (Sauer et al., 2020) combines DGPs with an ESS-based Bayesian posterior sampler and ALC-style acquisition. The key claim is that latent warping lets the active learner place runs non-uniformly in input space, concentrating them in wiggly, high-information regions that a stationary GP would not isolate in the same way. The work is implemented in the deepgp package on CRAN (Sauer et al., 2020). The later Vecchia-based extension strengthens the same surrogate-modeling agenda at much larger scale (Sauer et al., 2022).

The theory of DGP priors has also advanced. "Deep Gaussian Process Priors for Bayesian Inference in Nonlinear Inverse Problems" (Abraham et al., 2023) proves posterior contraction for inverse problems when the truth has compositional structure and the prior includes a hierarchical structure prior over depth, intrinsic dimensions, active subsets, and smoothness. In Darcy’s problem, the paper additionally proves that Whittle–Matérn-type Gaussian priors that ignore compositional structure can contract at a polynomially slower rate than the DGP prior for certain generalized additive diffusivities (Abraham et al., 2023). This suggests that deep compositional priors are not merely heuristic expressivity devices but can have asymptotic advantages when the target itself is compositional.

More recent descendants relax either the prior family or the likelihood family. "Deep Q-Exponential Processes" (Chang et al., 2024) replaces GP layers with Q-exponential-process layers, where $Z$9 recovers the GP case and smaller $$Y \leftarrow X_1 \leftarrow X_2 \leftarrow \cdots \leftarrow X_{H-1} \leftarrow Z.$$0 gives stronger $$Y \leftarrow X_1 \leftarrow X_2 \leftarrow \cdots \leftarrow X_{H-1} \leftarrow Z.$$1-style regularization. "Generalized and Scalable Deep Gaussian Process Emulation" (Ming et al., 25 Mar 2026) keeps latent DGP structure but adds a likelihood layer supporting heteroskedastic Gaussian, Poisson, negative binomial, and categorical responses, with implementation in dgpsi. These developments preserve the deep compositional template while broadening the class of observations and regularization regimes (Chang et al., 2024, Ming et al., 25 Mar 2026).

Two additional strands target dynamical structure and derivative information. "dynoGP: Deep Gaussian Processes for dynamic system identification" (Benavoli et al., 8 Feb 2025) constructs DeepGPs by cascading dynamic GPs, equivalent to stochastic linear time-invariant systems, with static GPs for instantaneous nonlinearities. "Distribution of Deep Gaussian process Gradients and Sequential Design for Simulators with Sharp Variations" (Yang et al., 20 Mar 2025) and "Deep Gaussian Processes with Gradients" (Booth, 19 Dec 2025) derive gradient means, covariance approximations, and gradient-aware sequential design or gradient-enhanced training, extending a classical strength of smooth GPs into the deep setting (Benavoli et al., 8 Feb 2025, Yang et al., 20 Mar 2025, Booth, 19 Dec 2025).

6. Limitations, controversies, and nomenclatural ambiguity

Several recurring limitations are explicit in the literature. First, exact inference is intractable, and different approximations make different compromises. IPVI argues that deterministic or mean-field Gaussian variational approximations can yield a biased posterior belief, while SGHMC-style stochastic alternatives reduce bias at substantially higher computational cost (Yu et al., 2019). Conditional DGPs similarly identify posterior collapse and diminished latent variance as failure modes of standard variational DGP inference (Lu et al., 2021). These papers do not reject DGPs; they identify inference quality as a central open issue.

Second, greater depth is not uniformly beneficial. The random-feature paper demonstrates architectures with up to 30 hidden layers but also states that deeper models train more slowly, even if 10- or 20-layer models can outperform shallower ones after convergence (Cutajar et al., 2016). IPVI reports that a 2-layer DGP often performs best on average for regression (Yu et al., 2019). In active-learning surrogates, three-layer DGPs can be stronger on average than shallower models but also show larger variability and overfitting tendencies, motivating deliberately parsimonious architectures (Sauer et al., 2020). A plausible implication is that depth in DeepGP is a statistical and computational design choice rather than a monotone path to better performance.

Third, the term itself is overloaded. In most of the cited literature, DeepGP denotes a deep Gaussian process. However, "Reading a Ruler in the Wild" (Pan et al., 9 Jul 2025) uses DeepGP for a completely different object: a lightweight feed-forward module inside RulerNet that regresses geometric-progression parameters $$Y \leftarrow X_1 \leftarrow X_2 \leftarrow \cdots \leftarrow X_{H-1} \leftarrow Z.$$2 from noisy ruler marks. That DeepGP is a custom 1D U-Net with only 1.6M parameters, trained on 1.23 billion randomly generated samples and benchmarked against optimization-based scale recovery; it is not a Gaussian process model (Pan et al., 9 Jul 2025). For that reason, the safest interpretation of the term depends on context: in probabilistic modeling and kernel-method literature it usually denotes deep Gaussian processes, but some recent application papers reuse the name for unrelated learned modules.

Taken together, the DeepGP literature describes not a single algorithm but a family of hierarchical stochastic-process models organized around GP composition. Its enduring themes are latent warping, uncertainty propagation, and the tension between posterior fidelity and scalability. Subsequent variants—multitask, convolutional, inter-domain, conditional, generalized, dynamic, and gradient-aware—retain that core while adapting it to specific structural assumptions or computational constraints.