Gaussian Processes

Distribution over functions, kernel selection, GP regression, sparse GPs, Bayesian optimisation

You have 5 noisy measurements of a physical process at time points {1, 2, 5, 7, 10}. You want to predict the value at time 3. A standard regression model gives you a point prediction — but how confident should you be? The measurements are noisy, and time 3 sits between two observations. A model that says "value = 4.2" with no indication of uncertainty is only half useful.

A Gaussian Process solves this by treating the unknown function itself as a random variable. A GP defines a distribution over functions: any finite collection of points {f(x₁), ..., f(xₙ)} follows a jointly Gaussian distribution. The kernel function k(x, x') defines the covariance between any two points — points close in time covary strongly, points far apart covary weakly. You specify the prior over functions by choosing the kernel. The posterior at time 3 conditions the joint Gaussian on your 5 observations, yielding a Gaussian predictive distribution: a mean (best estimate) and a variance (uncertainty). Near the data, the variance collapses — the GP knows what it knows. Far from the data, the variance expands — the GP correctly reports ignorance.

The mechanism that produces this behavior is the posterior update: μ* = K(X*, X)[K(X, X) + σ²I]⁻¹ y, and Σ* = K(X*, X*) - K(X*, X)[K(X, X) + σ²I]⁻¹ K(X, X*). The variance term Σ* shrinks at observed points and expands away from them — this is what no other regression method gives you automatically. At time 20, far beyond your training data, the posterior variance will be large: the GP correctly says it doesn't know.

The cost of this principled uncertainty is computational: the matrix inversion [K(X, X) + σ²I]⁻¹ costs O(n³) and storing K costs O(n²). At n = 10,000, this hits a hard wall on standard hardware. Everything in the sparse GP literature exists to circumvent this wall.

NOT this. "GPs are too slow for real data." Exact GPs are O(n³) and impractical for n > 10,000. But inducing point approximations (sparse GPs) reduce this to O(nm²) where m << n. For hyperparameter tuning via Bayesian optimization — where n is the number of function evaluations, typically under 200 — exact GPs are standard, fast, and the default choice in tools like BoTorch. The O(n³) wall is real, but it is the wall for exact GPs at large n. Sparse GPs, SKI, and deep kernel learning each offer a different strategy for breaking through it.

Key points

• A GP is fully specified by a mean function m(x) and a kernel k(x, x'). Any finite collection f(x₁), ..., f(xₙ) follows a multivariate Gaussian: f ~ N(m, K) where Kᵢⱼ = k(xᵢ, xⱼ). The kernel encodes structural beliefs about the function — its smoothness, length-scale, and periodicity. This is the prior over functions. A bad kernel is a bad prior: it doesn't just affect fit, it determines what shapes of function the GP can even consider.
• GP regression posterior: given noisy observations y = f(X) + ε, ε ~ N(0, σ²I), the posterior at new points X* is Gaussian with mean μ* = m(X*) + K(X*,X)[K(X,X)+σ²I]⁻¹(y-m(X)) and variance Σ* = K(X*,X*) - K(X*,X)[K(X,X)+σ²I]⁻¹K(X,X*). The variance term is what other regression methods don't give you: it collapses near observed data (the model knows what it knows) and balloons in unexplored regions (the model knows what it doesn't know).
• Defaulting to RBF is the most common GP mistake. RBF (squared exponential) is infinitely differentiable — it assumes the function is smoother than almost any real-world process. Sensor readings, financial returns, and experimental measurements are not infinitely differentiable. Matérn kernels control differentiability via ν: Matérn 3/2 (once differentiable) and 5/2 (twice differentiable) are almost always better defaults. The kernel is the most consequential modelling decision you make with a GP — it encodes what kinds of functions can explain the data.
• Hyperparameters (length-scale ℓ, signal variance σ², noise σ²_n) are learned by maximising the log marginal likelihood: log p(y|X,θ) = -½yᵀ(K+σ²I)⁻¹y - ½ log|K+σ²I| - n/2 log(2π). The first term rewards fit; the log-determinant penalises complexity — an overly flexible kernel can explain any data but pays a large penalty in the log-determinant term. This is automatic Occam's razor: the simplest kernel consistent with the data wins.
• O(n³) is the central engineering constraint. Inverting the n×n kernel matrix costs O(n³) compute and O(n²) memory. For n > 10,000, this is simply infeasible on standard hardware. This single fact drives the entire sparse GP literature — every method there is a different strategy for approximating or avoiding that matrix inversion.
• Sparse GPs introduce m << n inducing points Z that summarise the training data. SVGP (Stochastic Variational GP) places a variational distribution q(u) over inducing outputs u=f(Z) and optimises the ELBO with minibatch SGD — O(m³) per step regardless of n. With m=500 and minibatches, you can train on millions of points. The tradeoff: the inducing point approximation introduces bias — the posterior is not the true GP posterior, and predictive uncertainty may be underestimated between inducing points.
• Bayesian optimisation is where GPs earn their production keep. Use a GP as a surrogate for an expensive black-box function (hyperparameter tuning, drug discovery, materials science). The GP posterior gives you a mean prediction and uncertainty at any candidate point. Acquisition functions use both to decide where to query next: Expected Improvement (EI) = E[max(0, f(x) - f(x+))]; UCB = μ(x) + κσ(x) where κ tunes the explore-exploit tradeoff. BO with a GP is orders of magnitude more sample-efficient than grid or random search for expensive evaluations.
• GPs belong in your toolkit for: small datasets (n < 10,000) where calibrated uncertainty is the core deliverable; Bayesian optimisation (the standard surrogate in tools like BoTorch and Ax); scientific regression where domain knowledge belongs in the kernel. They don't belong in: high-dimensional unstructured inputs (images, raw text) where no meaningful kernel exists, or n > 100,000 without significant approximation infrastructure and the engineering budget to support it.
• Always normalise outputs to zero mean and unit variance before fitting a GP. GPs are sensitive to output scale in ways that quietly break the model if you skip this: the noise variance and signal variance hyperparameters are optimised assuming a certain output scale, and wrong assumptions propagate through to miscalibrated uncertainty estimates.

Recap

• GP = distribution over functions: any finite set $f(x_1)...f(x_n) sim N(m, K)$; kernel $k(x,x')$ IS the prior.
• Posterior variance collapses at data, expands away from it — the one thing other regression gives you for free.
• $O(n^3)$ compute, $O(n^2)$ memory from inverting $[K+σ^2I]$ — hard wall at n > 10,000. The central engineering fact.
• Kernel = prior over functions: RBF assumes infinite smoothness (usually wrong); Matérn 5/2 (twice-diff) is the better default.
• Hyperparameters learned by maximising log marginal likelihood — fit term vs log-det complexity penalty = automatic Occam's razor.
• Sparse GPs (SVGP): m ≪ n inducing points, $O(m^3)$/step; tradeoff = underestimated uncertainty between inducing points.
• Bayesian optimisation is where GPs earn their keep: GP surrogate + EI/UCB acquisition, far more sample-efficient than grid/random.

Check your understanding