Standard Gaussian Process (GP) methods offer a powerful, non-parametric Bayesian approach for tasks like regression and classification, providing not just predictions but also principled uncertainty estimates. However, their practical application is often hindered by computational demands, especially as dataset sizes grow.

The Computational Barrier in Standard Gaussian Processes

Recall from the section on GP regression formulation that the core computations involve the $N \times N$ covariance matrix $K_{NN}$ , constructed from evaluating the kernel function between all pairs of the $N$ training data points. Key operations include:

Training (Hyperparameter Optimization): Finding optimal kernel hyperparameters typically involves maximizing the marginal likelihood, which requires computing the inverse $K_{NN}^{-1}$ and the determinant $|K_{NN}|$ . Matrix inversion generally scales as $O(N^3)$ .
Prediction: Computing the predictive mean and variance for a new point $x_*$ requires solving linear systems involving $K_{NN}$ or multiplying by $K_{NN}^{-1}$ . These operations typically scale between $O(N)$ and $O(N^2)$ per test point, depending on whether $K_{NN}^{-1}$ is precomputed and reused.

The $O(N^3)$ complexity for training makes standard GPs computationally infeasible for datasets where $N$ grows beyond a few thousand data points. Modern datasets frequently contain hundreds of thousands or millions of samples, necessitating approximation techniques to make GPs practical at scale.

Strategies for Scaling Gaussian Processes

The main computational bottleneck is the manipulation of the dense $N \times N$ covariance matrix $K_{NN}$ . Approximation methods aim to circumvent this bottleneck, primarily by reducing the effective number of points the model depends on or by introducing structural assumptions into the covariance matrix. The most successful and widely used family of techniques relies on sparse approximations.

Sparse methods introduce a smaller set of $M$ inducing points ( $M \ll N$ ), denoted by $Z = \{z_1, ..., z_M\}$ , which act as representative points or a summary of the full dataset. The core idea is that the function's behavior across the entire input space can be adequately captured by its values $u = f(Z)$ at these $M$ locations. The inference and prediction process is then reformulated to depend computationally on $M$ rather than $N$ .

Sparse Approximations via Inducing Variables

Instead of conditioning directly on the full dataset $\mathbf{y}$ , sparse GP methods condition the GP $f$ on the function values $u$ at the inducing locations $Z$ . The approximation hinges on how the relationship between $u$ and $\mathbf{y}$ is handled. The inducing locations $Z$ themselves can be chosen in various ways: as a subset of the training inputs $X$ , placed on a grid, selected randomly, or, most effectively, optimized as part of the model training process.

Let's examine some common sparse approximation techniques.

Subset of Datapoints (SoD)

This is the most straightforward baseline: simply select a random subset of $M$ data points from the original $N$ points and perform exact GP inference using only this subset.

Pros: Easy to implement. Reduces complexity to $O(M^3)$ for training.
Cons: Discards potentially useful information from the $N-M$ ignored data points. The choice of subset significantly impacts performance, and a random subset is often suboptimal.

Subset of Regressors (SoR)

The SoR method uses $M$ inducing points $Z$ (not necessarily a subset of $X$ ) and bases predictions on these points. The predictive mean is approximated as $\mathbb{E}[f_* | \mathbf{y}] \approx K_{*M} K_{MM}^{-1} u_M$ , where $u_M$ represents the GP values at $Z$ . Various ways exist to estimate $u_M$ from the training data $\mathbf{y}$ , often involving approximations themselves.

Pros: Reduces prediction cost to $O(M^2)$ (after computing $K_{MM}^{-1}$ ). Can perform better than SoD if $Z$ is chosen well.
Cons: Less theoretically grounded than variational methods. Optimal selection or optimization of $Z$ is not inherent to the original formulation.

Fully Independent Training Conditional (FITC)

FITC introduces a specific structural assumption: given the inducing variables $u$ , the training outputs $y_i$ are conditionally independent. This approximates the true covariance matrix $K_{NN}$ with a sum of a low-rank term and a diagonal matrix: $K_{NN} \approx Q_{NN} + \text{diag}(\Lambda)$ , where $Q_{NN} = K_{NM} K_{MM}^{-1} K_{MN}$ .

Pros: Allows for faster computation, with training complexity often around $O(NM^2)$ . Offers better performance than SoR in many cases.
Cons: The conditional independence assumption can be inaccurate, potentially leading to poor uncertainty estimates.

Variational Free Energy (VFE) / Sparse Variational GP (SVGP)

Modern sparse GP methods are predominantly based on variational inference (VI). Instead of directly approximating the covariance matrix, these methods approximate the posterior distribution $p(f | \mathbf{y})$ . They introduce a variational distribution $q(u)$ over the inducing variables $u = f(Z)$ , where $q(u)$ is typically chosen to be a multivariate Gaussian $\mathcal{N}(m, S)$ . The goal is to minimize the KL divergence between the approximate posterior $q(f) = \int p(f|u) q(u) du$ and the true posterior $p(f | \mathbf{y})$ .

This minimization is achieved by maximizing the Evidence Lower Bound (ELBO) on the marginal likelihood $\log p(\mathbf{y})$ :

\mathcal{L}( \theta, Z, m, S) = \mathbb{E}_{q(f)}[\log p(\mathbf{y} | f)] - \text{KL}[q(u) || p(u)]

Here, $p(u) = \mathcal{N}(0, K_{MM})$ is the GP prior evaluated at the inducing points $Z$ . $\theta$ represents the kernel hyperparameters. A significant advantage of this variational framework is that the inducing point locations $Z$ can be treated as variational parameters and optimized alongside the kernel hyperparameters $\theta$ and the parameters of $q(u)$ (i.e., $m$ and $S$ ) by maximizing the ELBO using gradient-based methods. This allows the model to learn the optimal placement of inducing points to best summarize the data.

Stochastic Variational Gaussian Processes (SVGP)

The VFE framework naturally leads to stochastic optimization, making it highly scalable. The first term in the ELBO, the expected log-likelihood, can be written as a sum over the $N$ data points: $\sum_{i=1}^N \mathbb{E}_{q(f_i)}[\log p(y_i | f_i)]$ . This structure is amenable to the techniques of Stochastic Variational Inference (SVI) covered in Chapter 3.

Instead of computing the expectation over the full dataset at each step, we can use mini-batches $B \subset \{1, ..., N\}$ to compute noisy but unbiased estimates of the ELBO's gradient:

\nabla \mathcal{L} \approx \frac{N}{|B|} \sum_{i \in B} \nabla \mathbb{E}_{q(f_i)}[\log p(y_i | f_i)] - \nabla \text{KL}[q(u) || p(u)]

These stochastic gradients are then used to update the parameters $(\theta, Z, m, S)$ using optimizers like Adam or RMSprop.

Pros: Scales GP training to very large datasets (millions of points). The per-iteration cost depends on the mini-batch size $|B|$ and the number of inducing points $M$ , typically $O(|B|M^2 + M^3)$ . Provides a robust theoretical foundation via variational inference. Allows joint optimization of inducing locations, variational parameters, and hyperparameters.
Cons: Requires careful tuning of optimization parameters (learning rate, batch size). Convergence can sometimes be sensitive to initialization. The quality of approximation still depends on $M$ .

Choosing the Number of Inducing Points M

The number of inducing points, $M$ , is a critical hyperparameter in sparse GP methods. It directly controls the balance between computational efficiency and the fidelity of the approximation to the full GP.

Computational Cost: Increases polynomially with $M$ . For SVGP, the per-iteration cost is roughly $O(M^3 + |B|M^2)$ , and prediction cost is $O(M^2)$ per point. Memory usage also scales with $M^2$ .
Approximation Accuracy: Generally improves as $M$ increases. A larger $M$ allows the variational distribution $q(u)$ to better capture the complexities of the true posterior $p(u|\mathbf{y})$ . As $M$ approaches $N$ , the VFE approximation theoretically converges to the exact GP.

In practice, $M$ is often determined by the available computational budget (time and memory) and monitored via performance on a validation set. Common choices range from $M=100$ to $M=1000$ or slightly more, depending on the dataset size and complexity. Visualizing the learned inducing point locations $Z$ can sometimes provide insights.

Comparison between an exact GP fit and a sparse approximation (SVGP) using M=5 optimized inducing points (purple crosses) on a 1D regression task. The sparse GP effectively approximates the mean and uncertainty of the exact GP at a fraction of the computational cost by learning the optimal positions for its inducing points.

Summary and Trade-offs

Approximation methods, especially sparse variational techniques like SVGP, are indispensable for applying Gaussian Processes to large-scale problems.

Scalability: They break the $O(N^3)$ barrier of exact GPs, often achieving complexity that scales linearly with $N$ (using stochastic optimization) and polynomially with the number of inducing points $M$ .
Approximation Quality: The fidelity depends on the method and the number of inducing points $M$ . Variational methods provide a principled optimization framework (maximizing the ELBO) that balances model fit with complexity control, often yielding high-quality approximations. Optimizing inducing point locations $Z$ is a significant advantage of methods like VFE/SVGP.
Practical Use: Libraries such as GPflow (TensorFlow), GPyTorch (PyTorch), and the GP module in scikit-learn (offering some basic approximations) provide implementations of these scalable methods, making them readily available.

While exact GPs are preferred for smaller datasets where computation is not limiting, sparse approximations empower practitioners to use the flexible, non-parametric, and uncertainty-aware nature of Gaussian Processes on the large datasets common in contemporary machine learning. The choice among sparse methods typically involves considering the dataset size, desired accuracy, and available computational resources, with SVGP being a strong default for very large datasets.