Gaussian Processes for Machine Learning

In classical supervised learning, you choose a model family first — say, a polynomial of degree 3 — and then fit its parameters. But why degree 3? A degree-2 polynomial might miss a real curve; degree 10 might memorize noise. The modeler must guess the right complexity, and guessing wrong wastes data (via cross-validation) or produces overconfident predictions.

Neural networks have a similar problem in reverse: they have too many parameters, so regularization techniques are bolted on to prevent overfitting. The result works, but provides no principled measure of how uncertain the prediction is.

What if, instead of choosing a fixed model and fitting its finite parameters, we placed a probability distribution directly over functions and let the data tell us which functions are plausible?

A Gaussian Process is the infinite-dimensional extension of the multivariate Gaussian distribution. Instead of describing the probability of a finite vector, it describes the probability of an entire function.

The formal definition is elegant: a GP is a collection of random variables, any finite number of which have a joint Gaussian distribution. It is fully specified by two things:

• A mean function $m(x) = \mathbb{E}[f(x)]$ — your default guess before seeing data (often set to zero).

• A covariance function $k(x, x') = \text{Cov}[f(x), f(x')]$ — also called the kernel — which encodes how correlated function values at two inputs should be.

Think of the kernel as a similarity ruler: it answers "if I know $f(x)$, how much does that tell me about $f(x')$?" Nearby inputs get high covariance (so the function is smooth); distant inputs get low covariance (so the function is free to do something different).

The kernel function is where domain knowledge enters. Different kernels express different beliefs about the function you expect to see:

• The Squared Exponential (RBF) kernel $k(x,x') = \sigma^2 \exp\!\bigl(-\frac{\|x-x'\|^2}{2\ell^2}\bigr)$ encodes infinitely smooth functions. The length-scale $\ell$ controls how quickly correlation decays with distance — small $\ell$ allows rapid wiggles, large $\ell$ enforces broad trends. The signal variance $\sigma^2$ controls the vertical scale.

• The Matérn family allows rougher functions, parameterized by a smoothness parameter $\nu$. At $\nu = 1/2$ you get an Ornstein-Uhlenbeck process (very rough); at $\nu \to \infty$ it converges to the RBF.

• Periodic kernels model repeating patterns — perfect for seasonal data.

• Kernels can be composed: adding two kernels mixes their patterns; multiplying creates interactions (e.g., a locally-periodic function that changes amplitude over time).

The magic happens when data arrives. Suppose we observe $n$ training points $\{(x_i, y_i)\}$ where $y_i = f(x_i) + \varepsilon$ and $\varepsilon \sim \mathcal{N}(0, \sigma_n^2)$ is observation noise. We want the posterior distribution $p(f_* | X, y, x_*)$ at a new test input $x_*$.

Because everything is Gaussian, the posterior is also Gaussian — and we can write its mean and variance in closed form. The predictive mean is a weighted combination of training outputs, where the weights come from the kernel. The predictive variance is large far from training data and shrinks near observed points.

This is the key insight: the GP gives you error bars for free. No bootstrapping, no ensemble of models — the uncertainty is a natural output of the Bayesian update.

How do you choose the right kernel and its hyperparameters (like the length-scale $\ell$)? GPs have a principled answer: maximize the marginal likelihood — the probability of the observed data after integrating out the function.

The log marginal likelihood decomposes into three terms. Think of it as a scoring system:

• Data fit: how well the model explains the training data (wants complex models).

• Complexity penalty: how much "volume" the model assigns to functions that don't match the data (penalizes complex models).

• Normalization constant: a fixed term.

The fit term rewards models that pass near the data. The complexity term penalizes models that spread probability over too many functions. Together they implement an automatic Occam's razor: the simplest model that explains the data wins.

For classification, the output is a class label, not a continuous value. The GP models a latent function $f(x)$ that is then squashed through a sigmoid (for binary) or softmax (for multi-class) to produce class probabilities.

The challenge is that the sigmoid breaks the Gaussian conjugacy — the posterior over $f$ is no longer Gaussian. Two main approximations are used:

• Laplace approximation: fit a Gaussian to the posterior by finding the mode and using the local curvature (Hessian) to define the spread. Fast but can be inaccurate for skewed posteriors.

• Expectation Propagation (EP): iteratively refines local Gaussian approximations at each data point. More accurate than Laplace, especially for small datasets.

Despite the approximation, GP classifiers provide calibrated uncertainty — they say "I'm 60% sure this is class A" rather than forcing a hard decision.

One of the book's deepest contributions is showing that GPs are not an isolated technique but a unifying perspective:

• SVMs are GPs without uncertainty. The SVM solution is the MAP (maximum a posteriori) estimate of a GP — the single most probable function. The GP adds the full posterior: not just the best guess, but the distribution around it.

• Infinite neural networks are GPs. Neal (1996) showed that a single-layer neural network with infinitely many hidden units converges to a GP. The kernel is determined by the weight prior and the activation function. This connection has been revived in the Neural Tangent Kernel theory.

• Splines are GP MAP estimates. Cubic smoothing splines are the MAP solution of a GP with a specific kernel. Regularization networks and RBF networks also emerge as special cases.

This unification means insights from one field transfer to others — kernel design in SVMs informs GP priors, and GP uncertainty quantification enriches neural network predictions.

The elephant in the room: GP regression requires inverting the $n \times n$ covariance matrix, which costs $O(n^3)$ time and $O(n^2)$ memory. For 10,000 points this is manageable; for millions, it's prohibitive.

Several approximation strategies have been developed:

• Sparse / inducing-point methods: pick $m \ll n$ representative "inducing points" and approximate the full GP through them. Reduces cost to $O(nm^2)$. FITC and variational methods (Titsias 2009) are popular choices.

• Structured kernels: if the data lies on a grid, Kronecker and Toeplitz structure in the covariance matrix enables $O(n \log n)$ computations.

• Random features: approximate the kernel with a finite-dimensional feature map (Rahimi & Recht 2007), turning the GP into a Bayesian linear regression in feature space.

These approximations have made GPs practical for datasets of hundreds of thousands of points, though neural networks remain the tool of choice for truly massive datasets.

The GP framework gave machine learning a language for uncertainty that parametric models lacked. Bayesian optimization — the most direct descendant — uses a GP surrogate to answer the question every engineer asks: "where should I experiment next?" That question, and the principled Bayesian answer GPs provide, is why this book remains essential reading two decades after publication.