Gaussian processes from first principles

Suppose you are walking a path through hilly country, and you know the height of the ground at five places along it, from signposts or a map. Two of the places are close together; between the others there are long stretches with nothing marked. You stop in one of those stretches. How high is the ground under your feet, and how sure can you be?

A Gaussian process is one way to answer that question. For every place you haven't measured, it gives a best guess, plus a band around it that shows how far off the guess could plausibly be. The same reasoning works when the input is not a place on a path but the setting of an experiment or of a slow computer simulation. That makes Gaussian processes useful as stand-ins for slow simulations, an idea that goes back at least to Sacks, Welch, Mitchell and Wynn (1989), and the band tells you where the next measurement would teach you the most before you take it.

This article builds a Gaussian process from two numbers up. You need to know what a bell curve is, and what a mean and a standard deviation are. Everything else is explained on the way. The figures are live, so move their sliders.

Before any math, make a guess. Picture the five known heights as points on a graph, with the position along the path across and the height up: two points close together, a wide gap in the middle. Someone asks you to draw a smooth profile of the ground through them and to shade the region where the true profile could plausibly run. Where would your shading be widest: right at a known point, in a gap between points, or beyond the last point? Hold on to your answer. Figure 3 checks it.

Two numbers that move together#

Start smaller than a whole path, with two uncertain numbers: the heights of the ground at two points, A and B, a few steps apart. You don't know either height exactly. But if the ground at A turns out higher than you expected, the ground at B very likely is higher too: the two values move together. Two points a few steps apart are nearly the same height; two points a long way apart can be completely unrelated.

To keep the arithmetic simple, measure each height in standard deviations from the average height of the area. Before you know anything, each one then has a mean of 0 and a standard deviation of 1. How strongly they move together is their correlation, written $\rho$ (rho). A correlation of 0 means knowing one tells you nothing about the other (for values that follow a joint bell curve, as here); a correlation of 1 means they are always equal.

Now measure the height at A and find the value $a$. If the two heights follow a joint bell curve (a two-dimensional Gaussian distribution), the standard rule for what you should now believe about B is short (Rasmussen & Williams, 2006, eq. A.6):

$$\begin{aligned} \text{best guess for } B &= \rho\, a, \\ \text{spread of } B &= \sqrt{1-\rho^2}. \end{aligned}$$

The best guess for B moves toward what you measured at A, by a fraction set by the correlation. The spread shrinks from 1 to $\sqrt{1-\rho^2}$, which is small only when $\rho$ is close to 1. With $\rho = 0.8$ the spread drops to 0.6; with $\rho = 0.3$ it barely moves, to 0.95.

Correlation between A and B, ρ

0.80

Measured height at A

1.5

Best guess for B

ρ × a = 0.80 × 1.5 = 1.20

Spread of B

√(1 − ρ²) = 0.60, was 1

95% range for B

0.02 to 2.38

Set $\rho$ to 0, and the measurement at A changes nothing about B: the two points are too far apart to be related. Set it to 0.99, and B is pinned almost exactly to A, as it would be one step away. This rule, applied to many numbers at once, is all a Gaussian process does.

From two numbers to a whole curve#

The profile of the ground along the path is a long list of numbers: the height at $x = 0$, at $x = 0.1$, at $x = 0.2$, and so on. An uncertain profile is a long list of uncertain numbers. To describe it the way we described A and B, we need a bell curve for each height, which is easy (mean 0, spread 1), and a correlation for every pair of heights, which is the interesting part.

For ground that rises and falls smoothly, nearby points should have similar heights, and points far apart can be unrelated. One rule that says exactly this makes the correlation fall with the distance $d$ between two points:

$$\rho(d) = \exp\!\left(-\frac{d^2}{2\ell^2}\right).$$

The number $\ell$ is the length scale: the distance over which heights stay strongly related. At a distance of one length scale the correlation is 0.61. At three length scales it is 0.01, and the two heights have practically nothing to say about each other. For the ground, the length scale sets how wide the hills are. This rule is called the squared exponential kernel, and Rasmussen and Williams (2006, eq. 2.31) use it as their running example. A kernel, also called a covariance function, is the rule that gives the correlation (more generally, the covariance) for any pair of inputs.

That is enough to draw random profiles. Pick 121 points between 0 and 10, fill in the correlation of every pair with the rule above, and draw one sample from the resulting 121-dimensional bell curve. Each sample is a whole profile.

Length scale, ℓ

1.00

New draws

The rule behind the draws: correlation against distance

Move the length scale. Short length scales give bumpy ground, because heights only a little apart are already nearly unrelated. Long ones give broad, gentle hills. Nothing else changes: the figure reuses the same random numbers for every setting, and the rule alone decides how they are tied together.

Those draws come from a Gaussian process. The textbook defines one as "a collection of random variables, any finite number of which have a joint Gaussian distribution" (Rasmussen & Williams, 2006, p. 13). Pick any points you like, two or 121 or a million, and the heights there follow a bell curve in that many dimensions, with correlations given by the kernel. Two things specify a Gaussian process completely: a mean function (here 0 everywhere, our guess before any data) and the kernel.

In standard notation

With a mean function $m(x)$ and a kernel $k(x, x')$, one writes $f \sim \mathcal{GP}(m, k)$. For any inputs $x_1, \dots, x_n$, the vector of outputs is Gaussian, with mean vector $m_i = m(x_i)$ and covariance matrix $K_{ij} = k(x_i, x_j)$. The kernel used here, with an output scale $\sigma_f$ that this article fixes at 1, is

$$k(x, x') = \sigma_f^2 \exp\!\left(-\frac{(x - x')^2}{2\ell^2}\right).$$

Not every rule will do: a kernel must give a valid covariance matrix, one that is positive semidefinite, for any choice of inputs (Rasmussen & Williams, 2006, sec. 4.1). The squared exponential does.

Measuring a few points#

Now use the five known heights. In the language of the last two sections, five of the uncertain heights become known, and we want to update all the others. It is the same move as in Figure 1, with more numbers. Each unknown height is correlated with each known one, by an amount set by their distance, and the rule for jointly Gaussian values combines all those correlations into a new best guess and a new spread.

Recall your guess from the start, then reveal the result.

Show the Gaussian process

In the starting view, the band closes at each known point. If a height is known exactly, nothing is left to be unsure about there. Every drawn profile passes through every dot.

Between known points, the band opens by an amount that depends on the gap. Between the two close points at $x = 2.4$ and $x = 3.1$, the band is about ±0.09 at its widest. Both neighbors are at most 0.35 away, much less than the length scale of 1.2, so they are strongly correlated with everything in between. In the middle of the wide gap, at $x = 4.65$, the nearest known points are 1.55 away, the correlation with each is only 0.43, and the band grows to about ±1.29.

Beyond the last known point, the band returns to where it started. At $x = 10$ the nearest known point is 2.7 away, the correlation is 0.08, and the band is ±1.95: practically the ±1.96 it had before any data. Far from everything measured, the model knows only what it knew at the start, and its best guess drifts back to the average height, 0 on this scale.

So the uncertainty is widest far beyond the data, and within the data it is widest in the largest gaps. The length scale decides what counts as large: with a length scale of 5, the same gap would barely open.

The formulas behind Figure 3

Collect the measured inputs in $X$ and the measured outputs in $\mathbf{y}$. Let $K$ be the matrix of kernel values between the measured inputs, and $\mathbf{k}_*$ the vector of kernel values between a new input $x_*$ and each measured input. Conditioning the joint Gaussian on the measurements gives (Rasmussen & Williams, 2006, eq. 2.19)

$$\begin{aligned} \bar{f}(x_*) &= \mathbf{k}_*^\top K^{-1} \mathbf{y}, \\ \operatorname{var}\, f(x_*) &= k(x_*, x_*) - \mathbf{k}_*^\top K^{-1} \mathbf{k}_*. \end{aligned}$$

With a single measurement and $\sigma_f = 1$, we have $K = 1$, $\mathbf{k}_* = \rho$ and $\mathbf{y} = a$, and these become $\rho a$ and $1 - \rho^2$: Figure 1's rule. The best guess is also a sum of bumps, one centered on each measurement and shaped like the kernel: $\bar f(x_*) = \sum_i \alpha_i\, k(x_i, x_*)$ with $\boldsymbol{\alpha} = K^{-1}\mathbf{y}$ (eq. 2.27). In code, $K^{-1}$ is never formed; one Cholesky factorization of $K$ serves both formulas (Algorithm 2.1). The code behind these figures follows the same steps.

Does the band depend on what you measured?#

Now a second guess: suppose the five heights had come out different, at the same five places. Would the band change shape? Decide, then press New values, same places in Figure 3.

The best guess moves; the band does not. Look at the variance formula above: it contains the measured inputs, but not the measured outputs. Rasmussen and Williams (2006, p. 18) put it plainly: the variance "does not depend on the observed targets, but only on the inputs; this is a property of the Gaussian distribution."

That has a practical consequence: you can see where the next measurement would teach you the most before you take it. Press Measure where the band is widest a few times and watch the band collapse, as each new measurement goes where the model is least sure.

One caveat: this holds for a fixed kernel. In practice the length scale is usually estimated from the measurements, as the next sections explain, and then the values do matter, through the length scale they suggest.

The assumptions live in the kernel#

Everything Figure 3 shows follows from one choice: the correlation rule and its length scale. Drag the length scale down to 0.4. The best guess now falls back toward 0 between known points, and the band balloons in every gap: the model assumes the ground could do almost anything a short distance from a known point. Drag it up to 3. The profile becomes stiff and the band thin: the model assumes the ground cannot bend quickly, so it trusts its interpolation.

The choice of kernel matters too. The squared exponential produces extremely smooth curves, which can be differentiated any number of times. Stein (1999) argued that this is unrealistic for many physical processes and recommended the Matérn family of kernels, which allows rougher curves (as summarized by Rasmussen & Williams, 2006, sec. 4.2.1).

When the band is wrong#

So far the band has been the model's own account of how unsure it is. Can the model be wrong about that? Figure 3 can show it, because its five heights were not chosen by hand: they were read off a pretend ground profile, itself drawn from a Gaussian process with a length scale of 1.2. Real ground would be rougher; this profile is smooth enough that the model's assumptions can be exactly right.

Press Start again and turn on Hidden profile. At a length scale of 1.2 the band is calibrated: averaged over many hidden profiles, it covers the hidden profile at about 95% of positions.

Now make a guess. Turn Hidden profile off and set the length scale to 3. The band narrows: the model now assumes broad, gentle hills. Will the hidden profile stay inside this narrow band? Decide, then turn Hidden profile on again.

It does not: toward the right end, it runs far outside the band. That is a model being confidently wrong, and nothing in the band warns you. The band is the uncertainty that follows from the kernel, not a measurement of the model's error: if the kernel is wrong, the band can be narrow and wrong.

That is why the length scale is rarely set by hand. It is usually chosen by asking which value makes the measured data most probable under the model, which is called maximizing the marginal likelihood (Rasmussen & Williams, 2006, sec. 5.2). It works well in practice, though the same book cautions that it can overfit when there are many settings to tune.

Exact and noisy measurements#

A surveyed height is the same every time you look it up. A phone's altitude reading is not: take it twice at the same spot and you get two numbers. For noisy measurements, the model treats each observation as the true value plus a random error with standard deviation $\sigma_n$, which adds $\sigma_n^2$ to the diagonal of $K$ (Rasmussen & Williams, 2006, eqs. 2.22–2.24).

Raise the measurement noise in Figure 3. The band no longer closes at the dots, and the best guess no longer passes through them: a single noisy measurement is not the true value, so the profile is allowed to miss it. The band in the figure is for the underlying profile; a band for a new noisy measurement would be wider still. For exact data, zero noise is the honest setting. In code, a tiny value is often kept on the diagonal anyway, for numerical stability.

The same idea, for a bridge#

Nothing in the mathematics needed the input to be a position on the ground. It needed only a distance between inputs, and the assumption that nearby inputs give similar outputs. Suppose a finite-element model of a footbridge takes an hour to run, and you have run it five times, for five different spans. You need the peak deflection for a sixth span. Read span for position and peak deflection for height, and everything above carries over: the five runs play the part of the five known heights, and the band says how far off the guess for the sixth span could plausibly be, under the kernel's assumptions. A simulation run twice with the same input gives the same output, so the noise is 0, which is how Sacks et al. (1989) treated computer experiments.

The band's practical consequence matters most here, where each new point costs an hour. Refined versions of measuring where the band is widest sit under methods engineers use for expensive models. Efficient global optimization chooses each new run by weighing the best guess against the uncertainty, to find a minimum in few runs (Jones, Schonlau & Welch, 1998). In structural reliability, AK-MCS adds runs where the Gaussian process is least sure whether the structure fails or survives, to keep the number of calls to the expensive model small (Echard, Gayton & Lemaire, 2011).

With several inputs, such as span, depth and load, each input usually gets its own length scale, and a length scale much longer than an input's range suggests that the output barely changes along that input (Rasmussen & Williams, 2006, sec. 5.1).

Where it stops working#

A sharp change, such as a cliff on the path or a switch from one failure mode to another in a structure, breaks the smoothness the kernel assumes. The Gaussian process will round the change off, and its band may not cover it.

The kernel here also uses one length scale everywhere, treating every part of the input range alike. Ground that is flat in one region and steep and broken in another, or a response that changes slowly in one region and fast in another, forces a compromise length scale that suits neither.

And the cost grows quickly with the data. Factoring the $n \times n$ matrix $K$ takes about $n^3/6$ operations (Rasmussen & Williams, 2006, Algorithm 2.1), so doubling the number of measurements makes this step eight times slower. Large data sets need approximations (Rasmussen & Williams, 2006, ch. 8). For a few dozen known heights or expensive simulations, the cost is negligible next to a single simulation run.

Why it is also called kriging#

In engineering papers, a Gaussian process used as a surrogate model is often called kriging. The name comes from geostatistics. Danie Krige (1951) studied how to value gold ore on the Witwatersrand from a limited number of samples: the problem of the path, underground, of estimating a value between the places where it was measured. Georges Matheron, who developed the theory, named the method after him (Rasmussen & Williams, 2006, sec. 2.8). The prediction theory itself goes back at least to the 1940s, to Wiener and Kolmogorov. Sacks et al. (1989) brought it to computer experiments, where each observation is the output of a slow simulation.

Back to the path#

You stopped in a long stretch between known heights. The Gaussian process gives a best guess for the height under your feet: a weighted combination of the known heights, in which the nearby ones usually count most. It also gives a band. The band is narrow next to a known point and widest in the largest gaps and beyond the last known point, and its width depends on where the known points are, not on their heights. Where it is widest is where one more measurement would teach you the most.

The band is not a guarantee. It is the uncertainty that follows from the kernel: from the assumption that the ground changes smoothly, over about one length scale. If the ground is rougher than that, or has a cliff, the band can be narrow and wrong. The same holds for the sixth span of the footbridge.

In one sentence#

A Gaussian process fixes, before any data, how strongly the outputs at any two inputs are related; each measurement then moves the best guess and shrinks the uncertainty in proportion to that relation, which is why the band is narrow near what you measured and wide far from it.

References#