Gaussian Processes (GPs) are among the most powerful tools in probabilistic machine learning. They provide not just predictions, but complete predictive distributions—means and variances that quantify epistemic uncertainty. This makes GPs invaluable for applications where knowing "how sure are we?" is as important as the prediction itself: active learning, Bayesian optimization, safety-critical systems, and more.

But GPs have a problem: they don't scale. The exact GP posterior requires $O(n^3)$ computation and $O(n^2)$ memory—precisely the limitations we've been discussing throughout this module.

Sparse Gaussian Processes solve this by introducing inducing points—a small set of $m$ pseudo-observations that summarize the entire dataset. Unlike Nyström (which also uses representative points), sparse GPs work within a principled Bayesian framework, optimizing the inducing point locations and values through variational inference.

The result is a GP that scales to millions of points while retaining full uncertainty quantification.

The core idea of sparse GPs is to introduce inducing variables—function values at a set of $m$ inducing locations $\mathbf{Z} = {\mathbf{z}_1, \ldots, \mathbf{z}_m}$ (also called inducing inputs or pseudo-inputs).

The setup:

Let $\mathbf{f} = (f(\mathbf{x}_1), \ldots, f(\mathbf{x}_n))^T$ be the function values at training inputs, and $\mathbf{u} = (f(\mathbf{z}_1), \ldots, f(\mathbf{z}_m))^T$ be the function values at inducing locations.

Under the GP prior:

$$\begin{bmatrix} \mathbf{f} \ \mathbf{u} \end{bmatrix} \sim \mathcal{N}\left( \mathbf{0}, \begin{bmatrix} \mathbf{K}{nn} & \mathbf{K}{nm} \ \mathbf{K}{mn} & \mathbf{K}{mm} \end{bmatrix} \right)$$

where:

• $\mathbf{K}_{nn}$: $n \times n$ kernel matrix among training points
• $\mathbf{K}{nm} = \mathbf{K}{mn}^T$: $n \times m$ kernel matrix between training and inducing points
• $\mathbf{K}_{mm}$: $m \times m$ kernel matrix among inducing points

The conditional distribution:

From standard Gaussian conditioning, the distribution of $\mathbf{f}$ given $\mathbf{u}$ is:

$$p(\mathbf{f} | \mathbf{u}) = \mathcal{N}(\mathbf{K}{nm}\mathbf{K}{mm}^{-1}\mathbf{u}, \mathbf{K}{nn} - \mathbf{K}{nm}\mathbf{K}{mm}^{-1}\mathbf{K}{mn})$$

Define:

• $\boldsymbol{\mu}f = \mathbf{K}{nm}\mathbf{K}_{mm}^{-1}\mathbf{u}$ — the mean (function of $\mathbf{u}$)
• $\mathbf{Q}{nn} = \mathbf{K}{nm}\mathbf{K}{mm}^{-1}\mathbf{K}{mn}$ — the Nyström approximation to $\mathbf{K}_{nn}$
• $\boldsymbol{\Sigma}f = \mathbf{K}{nn} - \mathbf{Q}_{nn}$ — the residual covariance

The graphical model:

Sparse GPs assume conditional independence given $\mathbf{u}$:

$$p(\mathbf{f}, \mathbf{u}) = p(\mathbf{f} | \mathbf{u}) p(\mathbf{u})$$

This structure means:

1. All information from $\mathbf{f}$ to predictions passes through $\mathbf{u}$
2. Given $\mathbf{u}$, different training points become conditionally independent (or nearly so)
3. The inducing points act as a "bottleneck" summarizing the data

Different sparse GP methods make different assumptions about how $\mathbf{f}$ relates to $\mathbf{u}$. The two most important are FITC (Fully Independent Training Conditional) and VFE (Variational Free Energy).

FITC (Snelson & Ghahramani, 2006):

FITC assumes the off-diagonal elements of the residual covariance $\boldsymbol{\Sigma}_f$ are zero:

$$p_{\text{FITC}}(\mathbf{f} | \mathbf{u}) = \prod_{i=1}^{n} p(f_i | \mathbf{u}) = \mathcal{N}(\boldsymbol{\mu}_f, \text{diag}(\boldsymbol{\Sigma}_f))$$

This makes training points conditionally independent given $\mathbf{u}$. The resulting approximate likelihood is:

$$p_{\text{FITC}}(\mathbf{y} | \mathbf{u}) = \mathcal{N}(\boldsymbol{\mu}_f, \text{diag}(\boldsymbol{\Sigma}_f) + \sigma^2 \mathbf{I})$$

Pros: Fast, captures heteroscedastic residual variance Cons: Approximation is heuristic, not principled; can underestimate uncertainty

VFE (Titsias, 2009):

VFE takes a variational inference approach. Instead of modifying the model, it finds the best approximation to the true posterior within a constrained family.

The variational objective (ELBO) is:

$$\mathcal{L} = \log \mathcal{N}(\mathbf{y} | \mathbf{0}, \mathbf{Q}{nn} + \sigma^2 \mathbf{I}) - \frac{1}{2\sigma^2} \text{tr}(\mathbf{K}{nn} - \mathbf{Q}_{nn})$$

The first term is the marginal likelihood under the Nyström approximation. The second term is a trace penalty that measures how much variance is lost by the approximation.

Key insight: The trace term $\text{tr}(\mathbf{K}{nn} - \mathbf{Q}{nn})$ penalizes poor inducing point placement. If the inducing points span the training data well, $\mathbf{Q}{nn} \approx \mathbf{K}{nn}$ and the penalty is small. If they miss important regions, the penalty is large.

The VFE objective in detail:

Maximizing the ELBO:

$$\mathcal{L} = -\frac{n}{2}\log(2\pi) - \frac{1}{2}\log|\mathbf{Q}{nn} + \sigma^2\mathbf{I}| - \frac{1}{2}\mathbf{y}^T(\mathbf{Q}{nn} + \sigma^2\mathbf{I})^{-1}\mathbf{y} - \frac{1}{2\sigma^2}\text{tr}(\mathbf{K}{nn} - \mathbf{Q}{nn})$$

Using the Woodbury identity, we can evaluate this in $O(nm^2 + m^3)$ instead of $O(n^3)$:

$$|\mathbf{Q}{nn} + \sigma^2\mathbf{I}| = \sigma^{2n} |\mathbf{K}{mm}|^{-1} |\mathbf{K}{mm} + \sigma^{-2}\mathbf{K}{mn}\mathbf{K}_{nm}|$$

Similarly, $(\mathbf{Q}_{nn} + \sigma^2\mathbf{I})^{-1}$ can be computed efficiently using matrix inversion lemmas.

While VFE reduces complexity from $O(n^3)$ to $O(nm^2)$, we still need to process all $n$ data points. For truly massive datasets (millions of points), even $O(nm^2)$ is prohibitive.

Stochastic Variational Gaussian Processes (SVGP) (Hensman et al., 2013, 2015) enable minibatch training by reformulating the variational objective to decompose over data points.

The key innovation:

SVGP introduces an explicit variational distribution over the inducing variables:

$$q(\mathbf{u}) = \mathcal{N}(\mathbf{m}, \mathbf{S})$$

where $\mathbf{m} \in \mathbb{R}^m$ and $\mathbf{S} \in \mathbb{R}^{m \times m}$ (positive definite) are variational parameters to be optimized.

The ELBO becomes:

$$\mathcal{L} = \sum_{i=1}^{n} \mathbb{E}_{q(f_i)}[\log p(y_i | f_i)] - \text{KL}(q(\mathbf{u}) | p(\mathbf{u}))$$

where $q(f_i)$ is the marginal of $f_i$ under $q(\mathbf{u})$.

Computing the ELBO:

1. KL divergence between Gaussians has a closed form:

$$\text{KL}(q(\mathbf{u}) | p(\mathbf{u})) = \frac{1}{2}\left[\text{tr}(\mathbf{K}{mm}^{-1}\mathbf{S}) + \mathbf{m}^T\mathbf{K}{mm}^{-1}\mathbf{m} - m + \log\frac{|\mathbf{K}_{mm}|}{|\mathbf{S}|}\right]$$

Cost: $O(m^3)$ for Cholesky of $\mathbf{K}_{mm}$ and $\mathbf{S}$.

2. Marginal distribution $q(f_i)$ for a training point:

$$q(f_i) = \mathcal{N}(\mu_i, \sigma_i^2)$$

where:

• $\mu_i = \mathbf{k}i^T \mathbf{K}{mm}^{-1} \mathbf{m}$
• $\sigma_i^2 = k_{ii} - \mathbf{k}i^T \mathbf{K}{mm}^{-1} \mathbf{k}i + \mathbf{k}i^T \mathbf{K}{mm}^{-1} \mathbf{S} \mathbf{K}{mm}^{-1} \mathbf{k}_i$

with $\mathbf{k}_i = [k(\mathbf{z}_1, \mathbf{x}_i), \ldots, k(\mathbf{z}_m, \mathbf{x}_i)]^T$.

Cost per point: $O(m^2)$.

3. Expected log-likelihood for Gaussian observation model:

$$\mathbb{E}_{q(f_i)}[\log p(y_i|f_i)] = -\frac{1}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}[(y_i - \mu_i)^2 + \sigma_i^2]$$

Total per-iteration cost: $O(|B| \cdot m^2 + m^3)$, independent of $n$!

The quality of sparse GP approximations depends heavily on the inducing point locations $\mathbf{Z}$. Unlike Nyström where landmarks are selected from training data, sparse GPs optimize inducing locations as continuous parameters.

Initialization strategies:

1. Random subset: Sample $m$ training points as initial locations
2. K-means centers: Use cluster centroids to cover the data space
3. Greedy maximization: Select points maximizing information gain
4. Latin hypercube: Space-filling design for low dimensions

Optimization dynamics:

Inducing points are optimized jointly with other parameters via gradient descent. The gradients push inducing points toward regions where:

• Data density is high (more information to capture)
• Residual uncertainty is large (poorly explained regions)
• Function varies rapidly (need more resolution)

In practice, inducing points often migrate during training, starting from initial positions and settling into data-adaptive configurations.

Natural gradients for variational parameters:

The variational parameters $(\mathbf{m}, \mathbf{S})$ live on a curved space (Gaussian distributions form a manifold). Standard gradient descent ignores this geometry, leading to slow convergence.

Natural gradients account for the information geometry, often converging much faster:

$$\tilde{\nabla}{\mathbf{m}} \mathcal{L} = \mathbf{S}^{-1} \nabla{\mathbf{m}} \mathcal{L}$$ $$\tilde{\nabla}{\mathbf{S}} \mathcal{L} = \mathbf{S}^{-1} \nabla{\mathbf{S}} \mathcal{L} \mathbf{S}^{-1}$$

GPyTorch and other GP libraries implement efficient natural gradient updates.

Choosing the number of inducing points:

The number $m$ controls the quality-cost tradeoff:

• Too few: Poor approximation, underfit, high uncertainty
• Too many: Slow training, diminishing returns

Practical heuristics:

• Start with $m \approx \sqrt{n}$ or $m \approx 100-1000$
• Increase until ELBO plateaus or cross-validation stops improving
• For high-frequency functions or scattered data, more inducing points are needed

Let's carefully analyze the computational and memory costs of sparse GP methods.

Exact GP baseline:

• Training: $O(n^3)$ for Cholesky decomposition
• Memory: $O(n^2)$ for kernel matrix
• Prediction: $O(n)$ per test point for mean, $O(n^2)$ for variance

VFE (full-batch):

SVGP (minibatch):

Prediction complexity:

Given a trained SVGP with variational parameters $(\mathbf{m}, \mathbf{S})$:

• Mean prediction: $\mu_* = \mathbf{k}*^T \mathbf{K}{mm}^{-1} \mathbf{m}$

  • Cost: $O(m^2)$ per test point (can precompute $\mathbf{K}_{mm}^{-1}\mathbf{m}$, then $O(m)$)

• Variance prediction: $\sigma_^2 = k_{**} - \mathbf{k}_^T \mathbf{K}{mm}^{-1} \mathbf{k}* + \mathbf{k}*^T \mathbf{K}{mm}^{-1} \mathbf{S} \mathbf{K}{mm}^{-1} \mathbf{k}*$

  • Cost: $O(m^2)$ per test point

Compare to exact GP: $O(n)$ for mean, $O(n^2)$ for variance. With $m << n$, predictions are dramatically faster.

Sparse GPs connect deeply with other kernel approximation methods we've studied.

Sparse GP vs. Nyström:

Nyström approximates the kernel matrix: $\tilde{\mathbf{K}} = \mathbf{K}{nm}\mathbf{K}{mm}^{-1}\mathbf{K}_{mn}$.

VFE uses the same low-rank structure but adds:

1. A trace penalty for the approximation error
2. Principled uncertainty quantification
3. Optimization of inducing locations

In fact, if we use fixed inducing points at training locations and ignore the trace term, VFE reduces to Nyström-based regression!

Sparse GP vs. RFF:

Both provide $m$-dimensional features, but with key differences:

The unified view:

All three methods can be understood through the lens of feature-space approximation:

1. RFF: $\phi(\mathbf{x}) = \sqrt{\frac{2}{D}}[\cos(\boldsymbol{\omega}j^T\mathbf{x} + b_j)]{j=1}^D$ — random trigonometric features

2. Nyström: $\phi(\mathbf{x}) = \mathbf{K}_{mm}^{-1/2} [k(\mathbf{z}j, \mathbf{x})]{j=1}^m$ — kernel-based features centered at landmarks

3. Sparse GP: Same features as Nyström, but with learned $\mathbf{Z}$ and full posterior inference over weights

Sparse GPs add the Bayesian layer: instead of point-estimating the weights, they maintain distributions, enabling uncertainty quantification.

Successfully deploying sparse GPs requires attention to several practical details.

1. Initialization matters:

Poor initialization can lead to:

• Local minima (inducing points stuck in suboptimal configurations)
• Collapsed inducing points (multiple points at the same location)
• Slow convergence

Best practices:

• Initialize $\mathbf{Z}$ with k-means centers or random data subset
• Initialize $\mathbf{m}$ to zeros or prior mean
• Initialize $\mathbf{S}$ to prior covariance $\mathbf{K}_{mm}$
• Warm-start by freezing $\mathbf{Z}$ initially

2. Hyperparameter optimization:

Sparse GPs still have kernel hyperparameters (length scales, variance, noise). Options:

• Joint optimization: Optimize with ELBO alongside variational parameters
• Two-stage: First fit hyperparameters with subset, then full SVGP
• Marginal likelihood: For small datasets, use exact GP for hyperparameters, then transfer to SVGP

3. Beyond regression:

SVGP extends to classification and multi-output settings:

• Classification: Use Bernoulli likelihood; expected log-likelihood no longer closed-form, but tractable via Gauss-Hermite quadrature or Monte Carlo
• Multi-output GP: Share inducing points across outputs with correlation modeling
• Deep GPs: Stack multiple GP layers with inducing points at each level

4. Software:

Mature SVGP implementations exist in:

• GPyTorch (PyTorch, highly optimized)
• GPflow (TensorFlow)
• GPy (NumPy/SciPy, slower but pure Python)
• GPML (MATLAB)

Sparse Gaussian Processes represent the most principled approach to scaling GPs while retaining their full Bayesian character—uncertainty quantification, principled model comparison, and coherent probabilistic predictions.

What's next:

We've now covered the three pillars of scalable kernel methods: Random Fourier Features, Nyström approximation, and Sparse GPs. In the final page of this module, we synthesize these techniques into a comprehensive scalability strategies guide—when to use which method, how to combine them, and how to approach real-world large-scale kernel learning problems.