Gaussian Processes and Deep Kernel Learning

Gaussian processes occupy a distinctive position in probabilistic deep learning. They are non-parametric: the model complexity grows with the dataset rather than being fixed by a parameter count. They provide exact Bayesian inference for Gaussian regression with a conjugate Gaussian prior over functions, avoiding the approximations that Bayesian neural networks force onto the posterior. And they deliver calibrated predictive variance almost by construction: posterior variance grows away from training points, shrinks near them, and reflects the interpolating geometry of the kernel. The canonical reference (Rasmussen and Williams 2006) is the standard GP textbook and the reference this part builds on.

The main obstacle to GP deployment at deep-learning scales is the $\mathcal{O}(N^3)$ cost of inverting the $N \times N$ kernel matrix and the $\mathcal{O}(N^2)$ memory cost of storing it¹. Four scaling strategies have dominated the literature. Inducing-point variational methods, introduced by (Titsias 2009) and scaled to stochastic gradients by (Hensman et al. 2013), replace the full dataset with a small set of $M \ll N$ pseudo-inputs whose values serve as sufficient statistics for the posterior, reducing cost to $\mathcal{O}(NM^2)$ per evaluation and $\mathcal{O}(M^3)$ for the inducing-point covariance; choosing $M$ and the inducing-point locations is itself a tuning problem². Structured kernel interpolation (the KISS-GP line of (Wilson and Nickisch 2015)) exploits Kronecker and Toeplitz structure in the kernel matrix when inputs lie on grids (or are interpolated to grids) to reduce the inversion to $\mathcal{O}(N)$ or $\mathcal{O}(N \log N)$. Conjugate-gradient Lanczos methods (Gardner et al. 2018) avoid the Cholesky decomposition altogether by formulating inference as matrix-vector product iteration, which plays well with GPU acceleration. Random Fourier features (Rahimi and Recht 2007) take the orthogonal kernel-side route: by Bochner's theorem any shift-invariant kernel admits the integral representation $k(x-x') = \mathbb{E}_{\omega}[\cos(\omega^\top(x-x'))]$ for some spectral density, so averaging $D$ Monte-Carlo samples $z_i(x) = \sqrt{2/D}\cos(\omega_i^\top x + b_i)$ yields an explicit $D$-dimensional feature map that reduces kernel ridge regression to a $D \times D$ primal system at $\mathcal{O}(N D^2)$ cost; the construction is still the standard scaffolding behind modern primal-form GP approximations.

Expressive kernel design has co-evolved with scalable inference. Deep GPs (Damianou and Lawrence 2013) stack GPs as compositional priors: the output of one GP feeds the input of the next, yielding a non-Gaussian prior over functions whose inductive bias is richer than any single stationary kernel; the practical depth ceiling is shallow because intermediate layers tend to collapse to near-deterministic mappings without identity initialization or skip-connections³, and beyond two or three layers the marginal-likelihood improvements saturate with mild overfit⁴. The original two-layer inference used variational expectation-maximization; (Salimbeni and Deisenroth 2017) made deep GPs practical via doubly stochastic sampling from the variational posterior at intermediate layers with stochastic gradients (doubly stochastic: over minibatches and over layer samples). Deep kernel learning (DKL) (Wilson et al. 2016a, 2016b) takes a complementary path: a neural feature extractor $h_\theta$ feeds a fixed-form GP kernel $k_{\text{base}}$, yielding $k_{\text{DKL}}(x, x') = k_{\text{base}}(h_\theta(x), h_\theta(x'))$ trained jointly by maximizing the marginal likelihood. (Ober et al. 2021) documents the failure modes: the feature extractor maps training inputs to coincident latent points and predictive variance collapses⁵, and the marginal-likelihood signal that supposedly regularizes the joint optimization can itself overfit on hyperparameters with the deep parameterization⁶.

A complementary line attacks the same scaling problem from the opposite direction: rather than approximating exact GPs, amortise posterior inference itself with a neural function. Neural Processes (Garnelo et al. 2018) take a meta-learning view. At training time the model sees many regression tasks. For each task it receives a context set $\mathcal{C} = \{(x_i, y_i)\}_{i=1}^{n_C}$ and a target set $\mathcal{T} = \{(x^*, y^*)\}$. A shared encoder $h$ maps each context pair to a representation $r_i = h(x_i, y_i)$, which are aggregated into a permutation-invariant summary $r = \frac{1}{n_C}\sum_i r_i$. A latent variable $z$ is sampled from a distribution $q(z \mid r)$, and a decoder $g(x^*, z)$ outputs a predictive distribution over $y^*$. The objective is an ELBO: the decoder's log-likelihood term rewards accurate predictions at target points; the KL term regularises $q(z \mid \mathcal{C})$ toward $q(z \mid \mathcal{C} \cup \mathcal{T})$, which is the GP-analogue of the posterior contracting around the full data. At test time, inference is a single forward pass through the encoder and decoder: $\mathcal{O}(n_C)$ in the context size, not $\mathcal{O}(n_C^3)$.

The Conditional Neural Process (CNP) variant drops $z$ entirely and conditions the decoder deterministically on the aggregated $r$. This makes training simpler (no ELBO, direct log-likelihood) but removes the ability to sample diverse function realisations from the same context; predictive variance comes only from the decoder's output distribution, not from uncertainty over $z$.

The Neural Process family inherits two GP properties (calibrated predictive variance and a function-space prior) and replaces the $\mathcal{O}(n_C^3)$ kernel inversion with $\mathcal{O}(n_C)$ amortised inference at the price of an offline meta-training stage. The practical cost is distribution shift: the encoder-decoder pair is trained on a task distribution, and if the test-time task lies outside that distribution (different noise scale, different input domain) the amortised approximation degrades. Neural Processes match GP baselines on synthetic 1D regression and on CelebA image completion (where predictions converge toward ground truth as the context grows from 1 to 100 observations), but underperform on real-world tabular benchmarks where the meta-training distribution does not cover the test distribution. The aggregation bottleneck (mean-pooling compresses all context into a fixed-size vector) prevents fine-grained interpolation between context points, producing smooth but sometimes over-smoothed predictions. Attentive variants of the family replace mean-pooling with cross-attention from target inputs to context pairs, recovering most of the GP-accuracy gap, and are the contemporary baseline for the family.

The final strand in this part is the NNGP-NTK bridge. (Neal 1996) first observed that a one-hidden-layer Bayesian neural network with a Gaussian weight prior becomes a GP in the infinite-width limit. (Lee et al. 2018) and (Matthews et al. 2018) extended this to arbitrarily deep fully connected networks with a recursive kernel formula whose two hyperparameters $(\sigma_w^2, \sigma_b^2)$ must be tuned to the edge of chaos to avoid degenerate fixed points⁷. (Novak et al. 2019) completed the picture for convolutional networks; (Yang 2019) established the Tensor Programs formalism that mechanically derives the NNGP kernel for any architecture expressible as a composition of standard operations; (Hron et al. 2020) applied the same machinery to multi-head attention. On the parallel optimization-theoretic track, (Jacot et al. 2018) showed that infinite-width networks trained by gradient descent with squared loss behave like kernel regression with the Neural Tangent Kernel (NTK), a cousin of the NNGP kernel that captures the linearized training dynamics rather than the Bayesian posterior. The NNGP/NTK limits formalize what "Bayesian neural network" and "trained neural network" mean in the infinite-width regime, but they pay a price: infinite-width Bayesian inference does not learn features, so the CIFAR-10 gap to SGD-trained finite-width CNNs persists⁸. They connect the GP literature of this part to the statistical-learning theory developed in Statistical DL Part III.

The table groups the GP literature into five families, walked through in chronological order. Exact GPs (Full-GP, GPyTorch's Cholesky path) deliver the gold-standard posterior at $\mathcal{O}(N^3)$ cost⁹, leaving everything else as an approximation budget; the kernel-design choice is itself a non-trivial inductive-bias commitment that does not transfer between datasets¹⁰. Sparse GPs (VFE / Titsias, Hensman SVGP, KISS-GP) reduce that cost to $\mathcal{O}(NM^2)$ or $\mathcal{O}(N)$ by representing the posterior through $M$ pseudo-inputs or a structured-grid interpolation; natural-gradient updates (Salimbeni, Eleftheriadis, et al. 2018) and orthogonally-decoupled parameterizations (Salimbeni, Cheng, et al. 2018) further accelerate convergence. The $M$-budget remains hand-tuned¹¹, and Hensman's SVGP is what made GPs Adam-trainable on GPU minibatches.

Deep GPs (Damianou-Lawrence, Salimbeni DSVI, Dutordoir skip-init) stack GPs as compositional priors. The 2-layer / 3-layer / 5-layer rows in the table tell a single story: the first composition step lifts the held-out NLL noticeably; the third composition step lifts it slightly; further depth saturates and risks mild overfit¹². Without identity-init or residual connections the intermediate latent collapses to a near-deterministic mapping¹³, which is why Dutordoir's skip-init recipe matters disproportionately to its incremental NLL gain. Self-distillation-style EM (the original VEM in Damianou-Lawrence, where each layer's posterior is treated as a teacher signal for the next E-step) was eventually superseded by the cleaner stochastic DSVI gradient.

Deep kernel learning (Wilson DKL, SV-DKL, Bayesian DKL of Ober et al.) puts a neural feature extractor in front of a GP kernel. On CIFAR-10 both DKL papers report their headline accuracy alongside the same backbone trained as a plain DCNN, and the reported gap is a fraction of a point, so the DKL gain over the deterministic baseline is real but small. The trouble is calibration: the feature extractor compresses training inputs onto coincident latents and predictive variance collapses¹⁴, and the ML-II marginal likelihood used for joint hyperparameter selection can itself overfit with a deep parameterization¹⁵. Ober et al.'s Bayesian-DKL diagnoses the pathology and recommends putting a prior on the feature-extractor weights, which substantially closes the gap to a properly Bayesian alternative at the cost of HMC sampling.

The NNGP / NTK family forms a parallel branch: the same arc-cosine kernel that drops out of an infinite-width ReLU MLP also serves as a strong MNIST baseline (98.8% top-1 with depth 2-20), and Novak's CNN-GAP variant lifts CIFAR-10 to 77.4% with no SGD at all. Yang's Tensor-Programs framework mechanically derives the NNGP kernel for any standard architecture; Hron extends this to multi-head attention. The persistent ceiling is the feature-learning gap¹⁶: the same paper's SGD-trained finite CNN with pooling sits 7.7 points above the NNGP-CNN-GAP number, which is the usual reason given for preferring practical neural networks over their infinite-width Bayesian limits despite the calibration story being on the GP side. Tuning $(\sigma_w^2, \sigma_b^2)$ onto the edge of chaos is required to avoid degenerate fixed points¹⁷, and the resulting kernel inherits the kernel-design problem in a different costume¹⁸.

A diagonal cut across the table separates the closed-form methods (Full-GP, NNGP, NTK; Gaussian likelihood or kernel-only) from the approximate-inference methods (SVGP, Deep-GP, DKL; non-Gaussian likelihoods, deep parameterizations, or both)¹⁹. Three readings. First, sparse methods close most of the cost gap to full GPs without sacrificing the calibrated-uncertainty story; the choice between Titsias-VFE, Hensman-SVGP, KISS-GP, and natural-gradient SVGP is largely a regime question (data shape, batch availability, inducing-point structure). Second, deep GPs help modestly and saturate quickly²⁰, so depth in the GP world is qualitatively different from depth in the deterministic-network world. Third, deep kernel learning is competitive on CIFAR-10 but the original papers underreported the variance-collapse pathologies that Ober et al.'s Bayesian-DKL diagnoses. The fix, full Bayes over feature-extractor weights, restores calibration at the cost of HMC sampling and gives up the original DKL argument for using a GP read-out.

References