Scalable Gaussian Processes

The previous page established an uncomfortable truth: exact Gaussian Process inference scales as O(N³), rendering it impractical for datasets beyond tens of thousands of points. The fundamental challenge is the N×N kernel matrix K—computing it, storing it, and inverting it all scale prohibitively with dataset size.

Sparse Gaussian Process methods offer a compelling resolution: what if we could compress the information from N training points into a much smaller set of M representative points, where M << N? If successful, we would replace the N×N matrix operations with M×M operations, dramatically reducing computational burden while preserving the essential GP predictions.

This seemingly simple idea—approximating a large GP with a smaller one—conceals remarkable mathematical depth. The challenge is not merely computational; it's representational. How do we choose which M points best summarize N points? What information is lost in the compression? How do we quantify the approximation error? And critically, how do we ensure the resulting approximation remains a valid Gaussian Process with meaningful uncertainty estimates?

This page develops the foundational sparse GP framework that addresses these questions, establishing the theoretical apparatus for all modern scalable GP methods.

At the heart of all sparse GP methods lies the concept of inducing points (also called pseudo-inputs, support points, or basis points). Instead of conditioning the GP directly on all N training observations, we introduce a set of M << N synthetic locations $Z = \{\mathbf{z}_1, ..., \mathbf{z}_M\}$ in the input space, along with corresponding function values $\mathbf{u} = [f(\mathbf{z}_1), ..., f(\mathbf{z}_M)]^\top$.

The key insight is that under the GP prior, the function values at any finite set of points are jointly Gaussian. In particular, the training function values $\mathbf{f} = [f(\mathbf{x}_1), ..., f(\mathbf{x}_N)]^\top$ and inducing values $\mathbf{u}$ satisfy:

$$\begin{bmatrix} \mathbf{f} \\ \mathbf{u} \end{bmatrix} \sim \mathcal{N}\left( \mathbf{0}, \begin{bmatrix} K_{ff} & K_{fu} \\ K_{uf} & K_{uu} \end{bmatrix} \right)$$

where:

• $K_{ff}$ is the N×N kernel matrix over training points
• $K_{uu}$ is the M×M kernel matrix over inducing points
• $K_{fu}$ is the N×M cross-covariance matrix between training and inducing points
• $K_{uf} = K_{fu}^\top$ is the transpose

The conditional distribution of $\mathbf{f}$ given $\mathbf{u}$ is also Gaussian:

$$p(\mathbf{f} | \mathbf{u}) = \mathcal{N}(K_{fu}K_{uu}^{-1}\mathbf{u}, \, K_{ff} - K_{fu}K_{uu}^{-1}K_{uf})$$

This conditioning relationship is exact under the GP prior—no approximation has been made yet. The magic of sparse GPs happens when we make strategic approximations to this conditional structure.

The graphical model perspective:

Sparse GPs can be understood through the lens of probabilistic graphical models. In the exact GP, all training points $\mathbf{f}$ are correlated through the full covariance matrix $K_{ff}$. In sparse GPs, we introduce inducing variables $\mathbf{u}$ as intermediate latent variables that mediate the dependence:

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

Different sparse GP approximations correspond to different assumptions about the conditional $p(\mathbf{f}|\mathbf{u})$—specifically, which dependencies are retained and which are discarded. This graphical model viewpoint clarifies that we're not merely subsampling data, but intelligently restructuring the probabilistic model.

The simplest approach to scaling GPs is to ignore most of the data entirely. Subset of Data (SoD) selects M representative training points from the full dataset of N points and performs exact GP inference using only this subset.

Let $(X_M, \mathbf{y}_M)$ denote the selected M points. Inference proceeds exactly as for standard GPs:

$$\mu_* = \mathbf{k}{*M}^\top (K{MM} + \sigma_n^2 I)^{-1} \mathbf{y}M$$ $$\sigma^2 = k(\mathbf{x}_, \mathbf{x}*) - \mathbf{k}{*M}^\top (K_{MM} + \sigma_n^2 I)^{-1} \mathbf{k}_{*M}$$

where $K_{MM}$ is the M×M kernel matrix over the selected subset.

Subset selection strategies:

The performance of SoD depends critically on how the M points are selected:

1. Random sampling: Simple but can miss important regions; high variance in performance
2. Uniform/Grid sampling: Works for low-dimensional, evenly distributed data
3. K-means clustering: Select cluster centers to cover the input space
4. Farthest point sampling: Greedily select points that maximize minimum distance to existing selections
5. Leverage score sampling: Weight points by their influence on the GP posterior

Despite its limitations, SoD serves as an important baseline. If a sophisticated sparse GP method cannot outperform well-selected SoD, its additional complexity is not justified. Moreover, SoD reveals the information-theoretic limits of what M points can capture about N points—setting expectations for all sparse methods.

Subset of Regressors (SoR), also known as the Nyström approximation in the kernel methods literature, takes a fundamentally different approach. Rather than ignoring data, SoR uses all N observations but approximates the kernel matrix $K_{ff}$ with a low-rank structure.

The key approximation is:

$$K_{ff} \approx Q_{ff} = K_{fu}K_{uu}^{-1}K_{uf}$$

Here, $Q_{ff}$ is the Nyström approximation of the full kernel matrix—a rank-M approximation that can be computed and stored efficiently. The intuition is that the kernel function can be approximated by its projection onto the span of kernel functions centered at the M inducing points.

Mathematical derivation:

Consider the kernel as defining a feature map $\phi(\mathbf{x})$ in a (possibly infinite-dimensional) Hilbert space:

$$k(\mathbf{x}, \mathbf{x}') = \langle \phi(\mathbf{x}), \phi(\mathbf{x}') \rangle$$

The Nyström approximation projects this feature map onto the subspace spanned by $\{\phi(\mathbf{z}_1), ..., \phi(\mathbf{z}_M)\}$. The projected feature map is:

$$\tilde{\phi}(\mathbf{x}) = K_{uu}^{-1/2} \mathbf{k}_u(\mathbf{x})$$

where $\mathbf{k}_u(\mathbf{x}) = [k(\mathbf{z}_1, \mathbf{x}), ..., k(\mathbf{z}_M, \mathbf{x})]^\top$. The approximate kernel becomes:

$$\tilde{k}(\mathbf{x}, \mathbf{x}') = \langle \tilde{\phi}(\mathbf{x}), \tilde{\phi}(\mathbf{x}') \rangle = \mathbf{k}u(\mathbf{x})^\top K{uu}^{-1} \mathbf{k}_u(\mathbf{x}')$$

Predictive equations under SoR:

Substituting the Nyström approximation into the GP predictive equations:

$$\mu_* = \mathbf{k}{*u}^\top (K{uu} + K_{uf}\Sigma^{-1}K_{fu})^{-1} K_{uf}\Sigma^{-1}\mathbf{y}$$

$$\sigma_^2 = \mathbf{k}{*u}^\top (K{uu} + K_{uf}\Sigma^{-1}K_{fu})^{-1} \mathbf{k}_{u}$$

where $\Sigma = \sigma_n^2 I$ is the noise covariance.

Using the Woodbury identity, these can be rewritten as:

$$\mu_* = \mathbf{k}{*u}^\top K{uu}^{-1} K_{uf} (Q_{ff} + \Sigma)^{-1} \mathbf{y}$$

The term $(Q_{ff} + \Sigma)^{-1}$ still involves an N×N matrix, but the low-rank structure of $Q_{ff}$ enables efficient computation via the matrix inversion lemma:

$$(Q_{ff} + \Sigma)^{-1} = \Sigma^{-1} - \Sigma^{-1}K_{fu}(K_{uu} + K_{uf}\Sigma^{-1}K_{fu})^{-1}K_{uf}\Sigma^{-1}$$

The Deterministic Training Conditional (DTC) approximation, also known as the Projected Latent Variables approximation, provides a probabilistically coherent interpretation of the Nyström/SoR method. Rather than viewing the approximation as a matrix factorization trick, DTC interprets it through the lens of conditional Gaussian distributions.

Recall the exact joint distribution:

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

The exact conditional $p(\mathbf{f}|\mathbf{u})$ has covariance $K_{ff} - K_{fu}K_{uu}^{-1}K_{uf}$, which is non-diagonal and captures residual correlations between training points not explained by the inducing points.

DTC's key approximation is to replace this conditional with a deterministic mapping:

$$p_{\text{DTC}}(\mathbf{f}|\mathbf{u}) = \mathcal{N}(K_{fu}K_{uu}^{-1}\mathbf{u}, \, \mathbf{0})$$

That is, given the inducing values $\mathbf{u}$, the training function values are exactly determined with no residual uncertainty. This is a dramatic simplification—we assume the inducing points capture all variability in the function.

Consequences of the DTC assumption:

With observation noise $\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$, the likelihood becomes:

$$p(\mathbf{y}|\mathbf{u}) = \mathcal{N}(K_{fu}K_{uu}^{-1}\mathbf{u}, \, \sigma_n^2 I)$$

The marginal likelihood, integrating out $\mathbf{u}$:

$$p(\mathbf{y}) = \int p(\mathbf{y}|\mathbf{u})p(\mathbf{u})d\mathbf{u} = \mathcal{N}(\mathbf{0}, \, Q_{ff} + \sigma_n^2 I)$$

where $Q_{ff} = K_{fu}K_{uu}^{-1}K_{uf}$ is exactly the Nyström approximation.

The posterior over inducing points:

$$p(\mathbf{u}|\mathbf{y}) = \mathcal{N}(\mu_u, \Sigma_u)$$

where: $$\Sigma_u = (K_{uu}^{-1} + \sigma_n^{-2}K_{uf}K_{fu})^{-1} = K_{uu}(K_{uu} + \sigma_n^{-2}K_{uf}K_{fu})^{-1}K_{uu}$$ $$\mu_u = \sigma_n^{-2}\Sigma_u K_{uf}\mathbf{y}$$

The Fully Independent Training Conditional (FITC) approximation addresses DTC's systematic variance underestimation by retaining diagonal residual variance—the part of each training point's variance not explained by the inducing points.

Recall the exact conditional:

$$p(f_i | \mathbf{u}) = \mathcal{N}(\mathbf{k}{iu}K{uu}^{-1}\mathbf{u}, \, k_{ii} - \mathbf{k}{iu}K{uu}^{-1}\mathbf{k}_{ui})$$

The conditional variance $k_{ii} - \mathbf{k}{iu}K{uu}^{-1}\mathbf{k}_{ui}$ represents the function uncertainty at point $i$ that is not captured by the inducing points. This is non-negative and can be substantial in regions far from inducing points.

FITC makes the approximation that, conditioned on $\mathbf{u}$, the training function values are independent (but not deterministic):

$$p_{\text{FITC}}(\mathbf{f}|\mathbf{u}) = \prod_{i=1}^N \mathcal{N}(f_i; \mathbf{k}{iu}K{uu}^{-1}\mathbf{u}, \, \lambda_i)$$

where $\lambda_i = k_{ii} - \mathbf{k}{iu}K{uu}^{-1}\mathbf{k}_{ui}$ is the residual variance at point $i$.

The FITC covariance structure:

Under FITC, the prior covariance over training points becomes:

$$K_{\text{FITC}} = Q_{ff} + \text{diag}(K_{ff} - Q_{ff}) = Q_{ff} + \Lambda$$

where $\Lambda = \text{diag}(\lambda_1, ..., \lambda_N)$ is a diagonal matrix of residual variances.

This structure is crucial: FITC says "training points are correlated through the inducing points (via $Q_{ff}$), but have independent residual noise beyond that (via $\Lambda$)." The diagonal structure of $\Lambda$ enables efficient computation.

The FITC marginal likelihood:

$$\log p(\mathbf{y}) = \log \mathcal{N}(\mathbf{0}, Q_{ff} + \Lambda + \sigma_n^2 I)$$

Defining $D = \Lambda + \sigma_n^2 I$ (a diagonal matrix), this becomes:

$$\log p(\mathbf{y}) = \log \mathcal{N}(\mathbf{0}, Q_{ff} + D)$$

The Woodbury formula gives:

$$(Q_{ff} + D)^{-1} = D^{-1} - D^{-1}K_{fu}(K_{uu} + K_{uf}D^{-1}K_{fu})^{-1}K_{uf}D^{-1}$$

Since $D$ is diagonal, $D^{-1}$ is trivially O(N). The computational bottleneck remains the M×M matrix $(K_{uu} + K_{uf}D^{-1}K_{fu})$, preserving O(NM²) complexity.

PITC: Partially Independent Training Conditional:

A natural extension is PITC, which partitions training points into blocks and models block-wise dependencies while ignoring cross-block dependencies. PITC interpolates between FITC (every point is its own block) and the exact GP (all points in one block). It offers improved approximation at the cost of O(NB²) additional computation where B is the block size. In practice, FITC often provides sufficient improvement over DTC, and modern variational methods (covered in subsequent pages) supersede PITC.

Sparse GP methods introduce approximation error. Understanding when this error is acceptable—and when it corrupts results—is essential for practitioners. Several perspectives illuminate approximation quality:

The projection interpretation:

The Nyström approximation $Q_{ff} = K_{fu}K_{uu}^{-1}K_{uf}$ is a rank-M approximation to the full kernel matrix. It is optimal (in Frobenius norm) among rank-M approximations that pass through the inducing columns of $K$. The approximation error is:

$$\|K_{ff} - Q_{ff}\|F^2 = \sum{i=M+1}^N \lambda_i^2$$

where $\lambda_i$ are eigenvalues of $K_{ff}$ not captured by the inducing point projection. If the kernel matrix has rapidly decaying spectrum (common for smooth kernels), a small M captures most variance.

The information-theoretic view:

From an information theory perspective, the inducing points $\mathbf{u}$ should maximize the mutual information $I(\mathbf{f}; \mathbf{u})$—the information they provide about the full function values. This leads to the objective:

$$\max_Z I(\mathbf{f}; \mathbf{u}) = \frac{1}{2}\log|K_{ff}| - \frac{1}{2}\log|K_{ff} - Q_{ff}|$$

Equivalently, we minimize the residual entropy $H(\mathbf{f}|\mathbf{u})$. Greedy algorithms for inducing point selection often optimize this criterion.

When sparse approximations fail:

Sparse GPs struggle in specific scenarios:

1. Rapidly varying functions: If the true function has length-scale much shorter than inducing point spacing, the approximation misses important structure

2. Non-stationary functions: Standard kernels assume uniform length-scale; local structure may require many inducing points

3. High dimensions: Covering input space requires exponentially many inducing points (curse of dimensionality)

4. Outlier regions: Data points far from inducing points have less predictive influence than they should

Efficient implementation of sparse GPs requires careful attention to numerical linear algebra. Several techniques are essential for production-quality code:

Cholesky-based computation:

Never explicitly compute matrix inverses. Instead, use Cholesky factorization:

$$K_{uu} = L_{uu}L_{uu}^\top$$

where $L_{uu}$ is lower-triangular. Then $K_{uu}^{-1}\mathbf{v}$ is computed by solving $L_{uu}L_{uu}^\top\mathbf{x} = \mathbf{v}$ via forward and back substitution. This is more numerically stable and equally fast.

Memory-efficient cross-covariance:

The N×M cross-covariance matrix $K_{fu}$ dominates memory for large N. For prediction, we often don't need the full matrix—only its product with certain vectors. When possible, compute $K_{fu}\mathbf{v}$ incrementally without materializing the full matrix.

Incremental updates:

When adding new inducing points (e.g., during greedy selection), use rank-1 updates to existing Cholesky factors rather than recomputing from scratch.

GPU acceleration:

Matrix-matrix products dominate sparse GP computation. These are highly parallelizable, making GPUs effective even at moderate N. The key is to batch operations and minimize CPU-GPU transfers.

We have developed the fundamental framework for sparse Gaussian Processes—methods that reduce complexity from O(N³) to O(NM²) by compressing the GP through M inducing points. The key insights:

What's next:

The methods presented here—DTC, FITC, and their variants—are foundational but somewhat ad-hoc. In the next page, we examine inducing point methods more systematically, exploring strategies for selecting and even learning the optimal inducing point locations. Following that, we develop the elegant variational sparse GP framework that unifies and extends these classical approximations with principled Bayesian inference.