Scalable Gaussian Processes

Gaussian Processes are often described as one of the most elegant frameworks in machine learning—a method that provides not just predictions but principled uncertainty estimates, automatically adapts model complexity through marginal likelihood optimization, and offers a coherent Bayesian treatment of the learning problem. However, beneath this mathematical beauty lies a computational reality that has historically limited GP adoption: standard Gaussian Processes do not scale to large datasets.

This limitation is not a minor inconvenience or an edge case encountered only by specialists. It is a fundamental characteristic of the exact GP formulation that manifests almost immediately as datasets grow beyond thousands of examples. Understanding why GPs struggle with scale—the mathematical origins of their computational complexity—is essential before we can appreciate the ingenious approximations that modern GP practitioners have developed to overcome these barriers.

In this page, we dissect the exact GP computational pipeline, reveal the cubic bottleneck that dominates runtime, analyze memory requirements that can exhaust modern hardware, and establish the quantitative scaling relationships that define the boundary between tractable and intractable GP inference.

To understand GP scaling limitations, we must first examine what computations exact GP inference requires. Recall that for a dataset with N training points $\{(\mathbf{x}i, y_i)\}{i=1}^N$ and a Gaussian Process with mean function $m(\mathbf{x})$ and kernel $k(\mathbf{x}, \mathbf{x}')$, the predictive distribution at a test point $\mathbf{x}_*$ is:

$$p(f_* | \mathbf{y}, X, \mathbf{x}*) = \mathcal{N}(\mu, \sigma_^2)$$

where:

$$\mu_* = \mathbf{k}_*^\top (K + \sigma_n^2 I)^{-1} \mathbf{y}$$

$$\sigma_^2 = k(\mathbf{x}_, \mathbf{x}*) - \mathbf{k}^\top (K + \sigma_n^2 I)^{-1} \mathbf{k}_$$

Here $K$ is the $N \times N$ kernel matrix with entries $K_{ij} = k(\mathbf{x}_i, \mathbf{x}j)$, $\mathbf{k}* = [k(\mathbf{x}1, \mathbf{x}), ..., k(\mathbf{x}N, \mathbf{x})]^\top$ is the $N \times 1$ cross-covariance vector, and $\sigma_n^2$ is the observation noise variance.

These innocent-looking equations conceal substantial computational work. Let's decompose them:

Why Cholesky decomposition?

The term $(K + \sigma_n^2 I)^{-1}\mathbf{y}$ could theoretically be computed via direct matrix inversion, but this approach has poor numerical stability and even worse complexity. Instead, we exploit the symmetry and positive-definiteness of the kernel matrix to compute its Cholesky decomposition:

$$K + \sigma_n^2 I = L L^\top$$

where $L$ is lower-triangular. This factorization allows us to:

1. Solve linear systems via forward and back substitution in O(N²)
2. Compute log-determinants as $\log|K + \sigma_n^2 I| = 2 \sum_i \log L_{ii}$
3. Generate samples from the GP prior or posterior efficiently

The decomposition itself, however, requires O(N³/3) operations—approximately N³/3 multiplications and N³/3 additions. For N = 10,000, this is roughly 333 billion floating-point operations. On a modern CPU achieving 100 GFLOP/s, this takes approximately 50 minutes of pure computation—and this is before any hyperparameter optimization, which may require hundreds of such decompositions.

Before the O(N³) computation becomes limiting, memory requirements often prove prohibitive. The kernel matrix $K$ is an $N \times N$ dense matrix that must be stored in memory. Using double-precision floating-point numbers (64 bits = 8 bytes each), the memory requirements scale quadratically:

$$\text{Memory} = N^2 \times 8 \text{ bytes}$$

This quadratic growth rapidly exhausts available RAM:

Additional memory considerations:

The kernel matrix itself is not the only memory consumer:

• Cholesky factor L: Another N² doubles, though this can overwrite K
• Working arrays: Cholesky algorithms require additional N-length workspace
• Gradient computation: Storing derivatives for hyperparameter optimization
• Multiple kernel components: Composite kernels require intermediate matrices
• GPU memory: If using GPU acceleration, VRAM is typically 12-80GB—more limiting than CPU RAM

In practice, the rule of thumb is that exact GP inference requires approximately 3× the memory of the kernel matrix when including all auxiliary storage. This pushes the practical limit down to roughly N ≈ 30,000 on a 32GB workstation, and N ≈ 10,000-15,000 on a typical laptop.

Let's develop concrete intuition for how the O(N³) complexity manifests in wall-clock time. The Cholesky decomposition of an N×N matrix requires:

$$\text{FLOPs} = \frac{N^3}{3} + O(N^2) \approx \frac{N^3}{3}$$

Modern optimized BLAS implementations (Intel MKL, OpenBLAS, cuBLAS) achieve near-peak theoretical performance:

• CPU (single core): ~20-50 GFLOP/s per core
• CPU (multi-core, 8 cores): ~100-300 GFLOP/s
• GPU (Tesla V100): ~7,000 GFLOP/s (double precision)
• GPU (A100): ~10,000 GFLOP/s (double precision)

Using these rates, we can estimate wall-clock time:

The hyperparameter optimization multiplier:

These timings represent a single model evaluation. In practice, hyperparameter optimization via marginal likelihood maximization requires many evaluations—typically 50-500 iterations of gradient-based optimization, each requiring:

1. One Cholesky decomposition
2. One triangular solve
3. Computation of the log marginal likelihood
4. Gradient computation (requiring additional triangular solves per hyperparameter)

For a model with 3 hyperparameters (length-scale, variance, noise) optimized over 100 iterations, the total time is roughly:

$$T_{\text{total}} \approx 100 \times (T_{\text{Cholesky}} + 4 \times T_{\text{solve}}) \approx 100 \times 1.2 \times T_{\text{Cholesky}}$$

This means hyperparameter optimization for N = 50,000 on a high-end GPU could take approximately 10 minutes—still feasible, but already demanding. For N = 100,000, we're looking at an hour or more even with top-tier hardware.

The log marginal likelihood, which serves as the objective function for GP hyperparameter optimization, requires the full Cholesky decomposition:

$$\log p(\mathbf{y} | X, \theta) = -\frac{1}{2}\mathbf{y}^\top (K_{\theta} + \sigma_n^2 I)^{-1} \mathbf{y} - \frac{1}{2}\log|K_{\theta} + \sigma_n^2 I| - \frac{N}{2}\log(2\pi)$$

Each term presents computational challenges:

1. The quadratic form $\mathbf{y}^\top K^{-1} \mathbf{y}$: Requires solving $K\mathbf{\alpha} = \mathbf{y}$, which needs the Cholesky factor

2. The log-determinant $\log|K|$: Equals $2\sum_i \log L_{ii}$, requiring the Cholesky factor

3. Gradients $\frac{\partial}{\partial \theta_j} \log p(\mathbf{y}|X,\theta)$: Require computing traces of the form $\text{tr}(K^{-1} \frac{\partial K}{\partial \theta_j})$

The gradient computation is particularly expensive because it involves:

The iterative penalty:

Gradient-based optimization typically requires 50-500 iterations to converge, depending on initialization, problem conditioning, and the number of hyperparameters. Each iteration requires a full likelihood and gradient evaluation. This means:

• A single inference pass at N = 10,000 might take 2 seconds on a GPU
• Hyperparameter optimization (100 iterations) takes 3-4 minutes
• With 5 random restarts to escape local optima: 15-20 minutes
• Cross-validation for model comparison: hours

This iterative cost structure makes interactive GP workflows impractical for datasets beyond ~10,000 points, even when memory is sufficient. Users cannot rapidly experiment with different kernels, test hypotheses, or perform exploratory analysis when each model variant requires significant wait times.

Beyond pure computational cost, large-scale GP inference faces numerical precision challenges that become increasingly problematic as N grows. These issues arise from the accumulated round-off error inherent in floating-point arithmetic.

Condition number deterioration:

The condition number of a matrix, $\kappa(K) = \frac{\lambda_{\max}}{\lambda_{\min}}$, measures its sensitivity to numerical perturbations. For kernel matrices:

• Eigenvalues typically span many orders of magnitude
• Squared exponential kernels can produce condition numbers > 10^{12}
• Large N increases the probability of near-singular submatrices

When $\kappa(K) \approx 10^{16}$ (machine precision for doubles), the Cholesky decomposition may fail or produce meaningless results.

The jitter dilemma:

Practitioners routinely add a small "jitter" term to the diagonal of the kernel matrix to ensure positive definiteness:

$$K_{\text{stable}} = K + \epsilon I$$

Typical values of $\epsilon$ range from 10⁻⁶ to 10⁻³, chosen heuristically. This introduces a fundamental trade-off:

• Too small: Numerical instability persists
• Too large: Model accuracy degrades (adds artificial noise)

As N increases, larger jitter is often needed, but this increasingly corrupts the model. This is a symptom of the ill-conditioning inherent to large GP covariance matrices—another manifestation of scaling difficulties that approximation methods must address.

Combining memory, computation, and numerical considerations, we can establish practical limits for exact GP inference across different hardware profiles. These limits represent the boundary beyond which approximation methods become mandatory, not optional.

The modern data reality:

These limits must be contextualized against the scale of modern datasets:

• Computer vision: ImageNet has 1.2M images
• Natural language: Large text corpora contain billions of tokens
• Scientific computing: Simulation outputs often exceed millions of points
• Time series: IoT sensors generate millions of observations per day
• Spatial data: Geographic datasets cover millions of locations

Even the most generous exact GP limits (N ≈ 100,000 with multi-GPU) represent a tiny fraction of available data in most contemporary applications. This gap between data availability and GP tractability is the central motivation for scalable GP methods.

Domain-specific implications:

The O(N³) barrier has historically relegated GPs to roles where small sample sizes are inherent:

• Bayesian optimization: Typically N < 1,000 function evaluations
• Emulator/surrogate modeling: 100s to low 1000s of simulation runs
• Small-sample active learning: Carefully selected examples
• Hyperparameter tuning: Limited budget for costly evaluations

Expanding GPs beyond these niches—to large-scale regression, classification, and time series—requires the approximation methods we will explore in subsequent pages.

Having established the precise nature of the exact GP scalability barrier, we can now preview the approaches that form the remainder of this module. Each method attacks the bottleneck from a different angle:

Each approach trades exactness for tractability—introducing approximation error in exchange for manageable computation. The art of scalable GPs lies in understanding these trade-offs and selecting methods appropriate for specific applications.

Critically, modern scalable GP methods are not merely "good enough" approximations. Variational sparse GPs, for instance, provide formal evidence lower bounds that enable principled model comparison. Random feature methods offer provable approximation guarantees that bound worst-case error. These are not ad-hoc heuristics but mathematically grounded techniques with understood properties.

We have thoroughly examined the computational barriers that prevent standard Gaussian Processes from scaling to large datasets. The key insights from this analysis establish the foundation for all scalable GP methods:

What's next:

In the following pages, we will systematically develop the mathematical machinery that overcomes these barriers. We begin with Sparse Gaussian Processes, which achieve O(NM²) complexity by cleverly compressing the GP representation through a small set of inducing points. This approach—and its variational extensions—forms the backbone of modern scalable GP implementations.