Gaussian Process Regression

Chapter 7 ended with a prior: a Gaussian process that assigns a probability to every function before any data arrives. This chapter adds the data. We evaluate the unknown function at a few inputs, and ask what the prior, together with those values, says about the function everywhere else.

The answer needs no new machinery. A Gaussian process says that any finite set of function values is jointly Gaussian, and Section 4.5 showed how to condition a joint Gaussian on some of its coordinates. Gaussian process regression is that one formula, applied with the observed inputs on one side and the inputs we want to predict on the other. Everything else in the chapter is about reading the result: what the posterior mean does between and beyond the data, why the uncertainty depends on where we looked but not on what we saw, what changes when observations are noisy, and how to compute all of it without numerical trouble.

The posterior uncertainty is the quantity the rest of the book spends. Every acquisition function in Chapter 12 is a rule for turning it into the next query, so it is worth knowing its shape well.

8.1 Conditioning on observations #

Write $f$ for the unknown function on an input domain $\X$, and suppose it has a Gaussian process prior with mean zero and kernel $k$:

We observe $f$ at $n$ inputs $\vx_1, \dots, \vx_n$, collected in a set $X$, and for now the observations are exact: $y_i = f(\vx_i)$. Stack the values into a vector $\vy = (y_1, \dots, y_n)^\T$. We want the distribution of $f$ at $m$ new inputs $X_*$, whose unknown values we stack into $\vf_*$.

The smallest case was worked in Example 4.3: two function values with correlation 0.8, one of them observed at 1.2, and a belief about the other that moved to mean 0.96 and standard deviation 0.6. What follows is that computation with $n$ observed values, any number of unobserved ones, and a kernel supplying the correlations.

By the definition of a Gaussian process, $\vy$ and $\vf_*$ are jointly Gaussian. Their covariance is built entry by entry from the kernel, so it has four blocks: covariances among the observed inputs, between observed and new inputs, and among the new inputs.

Here $[\mK]_{ij} = k(\vx_i, \vx_j)$ holds the covariances among observed inputs, $[\mK_*]_{ij} = k(\vx_i, \vx_{*j})$ those between observed and new inputs, and $[\mK_{**}]_{ij} = k(\vx_{*i}, \vx_{*j})$ those among new inputs.

Conditioning on $\vy$ is now the operation from Section 4.5.

So the posterior over the new values is

Nothing in this derivation depended on how many new inputs we chose, or which. The posterior is again a Gaussian process: for any finite set of new inputs, the predictions are jointly Gaussian with the mean and covariance above. That closure is what makes the method practical. The data are absorbed once, and the result can be queried anywhere.

For a single new input $\vx$, the blocks become a vector and two numbers. Write $\vk(\vx) = (k(\vx, \vx_1), \dots, k(\vx, \vx_n))^\T$ for the covariances between $\vx$ and the observed inputs. Then the posterior mean and variance are

These two lines are the working form of the whole chapter. The rest of the book writes $\mu_n(\vx)$ and $\sigma_n(\vx)$, always with their argument, when the number $n$ of observations matters. They are not to be confused with the noise standard deviation $\sigma_n$ of Section 8.3, which has no argument and whose subscript stands for noise.

8.2 Reading the posterior #

The figure below computes Equation (8.3) on a grid of 160 inputs. The blue line is the posterior mean; the shaded band covers the mean plus and minus 1.96 posterior standard deviations, which holds 95% of the posterior probability at each input. The strip underneath plots the standard deviation by itself.

A few experiments make the equations concrete.

Add a point far from the others. The band pinches to nearly nothing at the new input and opens again on either side. How quickly it opens is set by the lengthscale $\ell$ of the kernel: the RBF kernel $k(x, x') = \exp\!\left(-(x - x')^2 / 2\ell^2\right)$ has fallen to about $0.14$ two lengthscales away, so an observation says little about inputs more than two lengthscales from it.

Watch the mean between and beyond the data. Between nearby observations the mean interpolates smoothly. Far from all observations it returns to zero, the prior mean, and the band returns to the prior width. The model does not extrapolate trends; it reverts to what it believed before seeing data.

Shrink the lengthscale. The mean starts to wiggle back to zero between points, and the band balloons in every gap. Lengthen it and the mean becomes a stiff curve that may miss the points entirely if they disagree. Choosing the lengthscale is the subject of Chapter 9.

The case of a single observation shows the structure without any matrices.

8.2.1 The mean is a sum of bumps #

The single-observation case generalizes. Define the weights $\bm{\alpha} = \mK^{-1}\vy$. Then the mean in Equation (8.3) is

a weighted sum of $n$ kernels, one centered on each observation. Turning on the kernel bumps in the figure draws each term. When two observations are close, their kernels overlap and the weights must compensate for each other, which is why a weight can be much larger than the value it helps to fit, or of the opposite sign.

Equation (8.4) is also the prediction of kernel ridge regression, a method with no probabilistic reading, when its regularization strength equals the noise variance introduced in the next section. The Gaussian process adds the variance, which ridge regression does not have, and which Bayesian optimization needs (Kanagawa et al., 2018).

8.3 Noisy observations #

Real evaluations are rarely exact. A training run with a different random seed gives a different accuracy; a person walking in an exoskeleton has good and bad strides. The standard model adds independent Gaussian noise to each observation:

Because the noise is independent of $f$ and across observations, it adds $\sigma_n^2$ to the variance of each observation and nothing to any covariance. The joint distribution Equation (8.1) keeps its shape with $\mK$ replaced by $\mK + \sigma_n^2 \mI$, and so does the derivation:

Two things change in the picture. Move the noise slider in Figure 8.1 and the mean stops passing through the points: it now trades fit against smoothness, the way ridge regression does. The band also stops pinching to zero, because a noisy observation cannot pin down $f$ exactly.

The noise variance is usually not known. It is a hyperparameter, fitted along with the lengthscale in Chapter 9. Too little noise makes the model chase every fluctuation; too much makes it ignore real structure. The two can also trade off against each other: a short lengthscale with little noise and a long lengthscale with much noise can explain the same wiggly data, which is one reason fitted hyperparameters deserve a skeptical look.

8.4 Computing it #

The formulas contain a matrix inverse, but a careful implementation never forms one. The matrix $\mK + \sigma_n^2 \mI$ is symmetric and positive definite, so it has a Cholesky factorization $\mL\mL^\T$ with $\mL$ lower triangular (Section 3.5). Solving a triangular system costs $O(n^2)$ and is numerically stable; inverting a nearly singular matrix is neither.

Step 4 is the variance formula in disguise: $\mathbf{v}^\T\mathbf{v} = \vk^\T \mL^{-\T}\mL^{-1}\vk = \vk^\T(\mK + \sigma_n^2\mI)^{-1}\vk$.

The cost splits into a part paid once per data set and a part paid per prediction. The factorization in step 1 takes $O(n^3)$ time and $O(n^2)$ memory. After that, each mean costs $O(n)$ and each variance $O(n^2)$. A laptop factorizes a matrix with a few thousand rows in well under a second, and a Bayesian optimization run rarely has more than a few hundred observations, so the cubic cost is seldom the bottleneck in this book. Larger data sets need approximations that summarize the data with a smaller set of inducing points (Quiñonero-Candela and Rasmussen, 2005; Titsias, 2009).

import numpy as np

def rbf(a, b, ell=0.12):
    return np.exp(-0.5 * (a[:, None] - b[None, :]) ** 2 / ell**2)

def gp_posterior(x, y, xs, noise=1e-4, ell=0.12):
    L = np.linalg.cholesky(rbf(x, x, ell) + noise * np.eye(len(x)))
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    Ks = rbf(x, xs, ell)                 # n x m
    mean = Ks.T @ alpha
    v = np.linalg.solve(L, Ks)           # n x m
    var = 1.0 - np.sum(v**2, axis=0)     # k(x, x) = 1 for this kernel
    return mean, var

Even a well-posed kernel matrix can fail to factorize in floating point when two inputs are nearly identical, because two rows become nearly equal and the smallest eigenvalue rounds to zero or below. Implementations add a small "jitter", a constant such as $10^{-6}\sigma_f^2$, to the diagonal and retry with a larger one if the factorization still fails. With observation noise the problem rarely arises, since $\sigma_n^2$ already plays that role.

8.5 Posterior samples #

The mean and the band summarize the posterior one input at a time. They do not show what a single plausible function looks like, because neighboring values are strongly correlated. To see whole functions, draw samples from the joint posterior Equation (8.2) on a grid of inputs. With $\mSigma_*$ the posterior covariance and $\mL_*$ its Cholesky factor, each sample is

the sampling recipe of Section 4.3. Turn on the samples in Figure 8.1. Each violet curve passes through (or near) every observation, stays inside the band most of the time, and is as smooth as the kernel allows. Each is a function the model considers possible.

Samples are not only a visualization. Thompson sampling, one of the acquisition rules in Section 12.5, draws one posterior sample and evaluates the objective where that sample is largest. On a grid of $m$ points, sampling costs $O(m^3)$ for the factorization, which limits grids to a few thousand points; for larger or continuous domains, samples can be drawn as functions using random features or pathwise updates (Rahimi and Recht, 2007; Wilson et al., 2020).

8.6 Pitfalls #

Three habits prevent most surprises in practice.

Standardize the outputs. The zero prior mean and unit amplitude assume the function's values are centered near zero with spread near one. A function whose values sit around 1000 would be pulled toward zero away from the data. Subtract the mean of the observations and divide by their standard deviation before fitting, and undo the transformation on the predictions. Libraries such as BoTorch do this with an outcome transform (Balandat et al., 2020).

Scale the inputs. A single lengthscale assumes all input directions vary on comparable scales. Map each input to $[0, 1]$ first, and give each dimension its own lengthscale when they matter differently (Section 9.2).

Distrust extrapolation. Beyond the data, the posterior returns to the prior by construction. If the objective has a trend that continues past the observed range, a stationary kernel will not predict it. That is usually acceptable in optimization over a bounded domain, but it is a reason to make the domain no larger than it needs to be.

8.7 Exercises #

Further reading #

• Rasmussen and Williams (2006), chapter 2, is the standard derivation, in both the weight-space and function-space views, with the algorithm used here.
• Garnett (2023), chapters 2 to 4, develops Gaussian processes with Bayesian optimization in mind, including the inference choices this chapter treats as fixed.
• Görtler et al. (2019) is an interactive visual introduction that complements the figures in this chapter.
• Williams and Rasmussen (1996) introduced Gaussian process regression to machine learning. The same predictor had long been used in geostatistics as kriging (Krige, 1951; Matheron, 1963).
• Kanagawa et al. (2018) sets out the exact correspondences between Gaussian process regression and kernel methods such as kernel ridge regression.