Gaussian Processes

Almost every model puts a prior on parameters and derives a distribution over functions from it. A Gaussian process reverses the order: put the prior directly on functions, and never mention parameters at all. The payoff is a posterior with exact, closed-form uncertainty that widens automatically where there are no observations, which is why GPs remain the default surrogate model in Bayesian optimisation despite scaling badly.

The definition is a collection of random variables, any finite subset of which is jointly Gaussian. A GP is specified by a mean function \(m(x)\), usually taken as zero after centring, and a covariance function \(k(x, x')\):

The posterior is one block-matrix identity

Given observations \(y = f(X) + \varepsilon\) with \(\varepsilon \sim \mathcal{N}(0, \sigma_n^2 I)\), the joint distribution of training outputs and a test output is Gaussian, so conditioning is the standard Gaussian conditional:

Two features of these formulas are worth dwelling on. The mean is a linear combination of the training targets with weights set by the kernel, which makes a GP a kernel smoother with a principled weighting. And the variance does not depend on \(y\) at all: predictive uncertainty is a function of where you observed, never of what you observed. That is a genuine property of the Gaussian model and a real limitation, since it means a GP cannot become more uncertain because the data were surprising.

Away from the training inputs, \(k_*\) goes to zero and the posterior variance returns to the prior variance. Extrapolation reverts to the prior with honest, wide error bars, which is the behaviour that makes GPs safe inside an acquisition function.

The kernel is the model

Choosing \(k\) is choosing the hypothesis space. The squared exponential kernel, \(k(x,x') = \sigma_f^2\exp(-\|x-x'\|^2/2\ell^2)\), produces infinitely differentiable sample paths, which is usually far smoother than reality. The Matérn family with \(\nu = 3/2\) or \(5/2\) gives finitely differentiable paths and is the more defensible default for physical processes (Rasmussen and Williams, 2006, Gaussian Processes for Machine Learning). Periodic kernels encode a repeating structure; a sum of kernels expresses additive structure and a product expresses interaction, so kernel construction is genuinely a modelling language.

Hyperparameters, the lengthscale \(\ell\), signal variance \(\sigma_f^2\) and noise \(\sigma_n^2\), are typically fitted by maximising the log marginal likelihood

whose first term rewards fitting the data and whose log-determinant term penalises model complexity. Occam's razor is built in rather than added, which is one of the more elegant properties in the area. Automatic relevance determination gives each input dimension its own lengthscale, and the fitted values then rank feature relevance directly.

When it breaks

Cubic cost is the defining constraint. Inverting \(K + \sigma_n^2 I\) is \(O(n^3)\) in time and \(O(n^2)\) in memory. At \(n = 10{,}000\) that is a 100-million-entry matrix and a heavy factorisation; at \(n = 100{,}000\) it is out of reach exactly. Sparse approximations using \(m \ll n\) inducing points bring this to \(O(nm^2)\), and stochastic variational GPs extend it to minibatches, at the cost of an approximation whose quality depends on inducing-point placement.

Uncertainty is not data-dependent. As noted above, \(\sigma_*^2\) ignores \(y\). A GP with homoscedastic noise cannot express "the noise is larger in this region", which is common in real measurements. Heteroscedastic GPs model \(\sigma_n^2(x)\) with a second process and lose the closed form.

High dimensions dissolve the kernel. Distance-based kernels suffer the same concentration problem as any distance method; beyond roughly 10 to 20 dimensions all pairwise distances converge and the kernel stops discriminating. This is why GP-based Bayesian optimisation is standard for tens of hyperparameters and not for thousands.

The marginal likelihood is multimodal. Optimising hyperparameters can land on a long-lengthscale, high-noise explanation, where everything is noise around a flat function, or a short-lengthscale, low-noise one, where the function wiggles through every point. Both are local optima of the same objective and they imply completely different models. Multiple restarts are necessary, and inspecting the fitted lengthscale against domain knowledge is not optional.

Gaussian noise is assumed. Counts, binary outcomes and heavy-tailed residuals break the conjugacy that gives the closed form, and the posterior then needs approximation by Laplace, expectation propagation or MCMC, which removes much of the reason to choose a GP in the first place.