Distributions over Functions

The last sections of Chapter 6 borrowed a model this book had not yet built, the Gaussian process. This part builds it, starting from the last model that was built in full. Section 5.4 fitted a straight line the Bayesian way. It put a Gaussian prior on the line's two weights, observed a few noisy points, and found that the posterior over the weights was again Gaussian, available in closed form. Every quantity it produced, a prediction, its uncertainty, a plausible line drawn at random, came out of a few lines of algebra on Gaussians.

The trouble is the line. The book's running example has a wide hill on the left and a narrow, taller peak on the right (Section 1.2), and nothing guarantees that a real objective, a validation accuracy or a person's comfort, is straight either. We want to keep the algebra and drop the line.

This chapter does that in three moves. First, it gives the linear model many features, so that it can bend, and asks what the model's prior says about the function's values. Second, it notices that this prior depends on the features only through one function of two inputs, the kernel, and lets the number of features grow without bound until only the kernel is left. Third, it takes what remains as a definition: a Gaussian process, a probability distribution whose samples are whole functions. The chapter ends by drawing functions from that distribution and reading what each choice of kernel says about the function we are looking for. Chapter 8 then conditions this prior on data.

7.1 From weights to functions #

Write the straight line as $f(x) = w_1 + w_2 x$: two weights, each multiplying a fixed function of the input, the constant $1$ and $x$ itself. Bayesian linear regression places a Gaussian prior on the weight vector $\vw$ and needs nothing else from the model, because the function is linear in $\vw$. That linearity is what makes the algebra work, and nothing in it requires the two fixed functions to be $1$ and $x$.

So replace them with any $M$ fixed functions $\phi_1, \dots, \phi_M$ of the input, called features or basis functions, and keep the weights:

Here $\boldsymbol{\phi}(\vx) = (\phi_1(\vx), \dots, \phi_M(\vx))^\T$ stacks the features of one input, and $\mSigma_p$ is the prior covariance of the weights (the subscript is for prior). The model is still linear in $\vw$, so every Gaussian computation of Section 5.4 still applies. Only the features changed.

Which features? Polynomials and sines are common choices. This chapter uses bumps. For a one-dimensional input $x$, the $j$th bump is a Gaussian-shaped function centered at its own point $c_j$, with a width $s$ shared by all bumps:

A weighted sum of bumps rises where the weights are positive and dips where they are negative. With enough bumps, close enough together, it can take almost any smooth shape.

7.1.1 A prior over weights is a prior over functions #

Draw a weight vector from its prior and the model becomes one particular function. Draw again and it becomes another. The prior over weights is therefore, with no further assumption, a prior over functions, and the most direct way to see what it believes is to draw from it.

Look at the band first. With six bumps it bulges over each bump and narrows between them: the model is confident that $f$ is close to zero midway between two centers, for no reason except that no feature lives there. The functions it draws inherit the pattern. They swing over the centers and sag toward zero in the gaps. Nobody believes that about an objective function. It is an artifact of where the features happen to sit, and Section 7.2 removes it.

7.1.2 What the prior says about function values #

To make "what the prior believes" precise, ask about the function's values at a finite set of inputs $\vx_1, \dots, \vx_n$. Stack the values into a vector $\vf = (f(\vx_1), \dots, f(\vx_n))^\T$, and stack the features of those inputs into the $n \times M$ matrix $\boldsymbol{\Phi}$ whose $i$th row is $\boldsymbol{\phi}(\vx_i)^\T$. Then Equation (7.1), applied to each input in turn, says $\vf = \boldsymbol{\Phi}\vw$: the function values are a linear map of the weights.

So the values at any $n$ inputs are jointly Gaussian,

and every entry of the covariance has the same form:

Two things about this result carry the rest of the chapter. First, the weights have disappeared. The distribution of the function values is stated entirely through $k$, a function that takes two inputs and returns the covariance of the function's values there. Second, nothing restricted the number of inputs or which ones we chose. Any finite list of inputs gets a joint Gaussian, and every one of those Gaussians is built from the same $k$.

The strip under Figure 7.1 plots $k(x_0, x)$ against $x$ for one fixed $x_0$: the covariance of $f(x_0)$ with the value at every other input. Drag $x_0$ onto the center of a bump and the curve is tall and narrow. Drag it between two centers and the curve shrinks. In this model, how strongly two values are related depends on where they are, not only on how far apart they are.

7.2 The kernel #

The bump model has two defects, and the example has just named the second. Its beliefs depend on where the bumps were put, as the pinched band in Figure 7.1 shows. And with $M$ features it has only $M$ degrees of freedom, so with more than $M$ inputs its covariance matrix is singular: the prior is certain about combinations of values it has no reason to be certain about. Both defects have the same cure. Use more bumps, packed more densely, all the way to a continuum.

That cure sounds expensive. Computing $\boldsymbol{\phi}(x)$ for a million features costs a million evaluations per input. But Equation (7.2) never needs the features themselves, only the inner products in Equation (7.3). If those can be computed directly, the number of features stops mattering. For bumps they can, in closed form.

7.2.1 Infinitely many bumps #

Space the centers evenly, $\Delta$ apart, along the whole real line, and give the weights the same variance $\sigma_w^2$, independently, so that $\mSigma_p = \sigma_w^2 \mI$. Each new bump adds variance, so as the spacing shrinks the weight variance must shrink in proportion, or the functions would grow without bound. Choose $\sigma_w^2 = \Delta / (s\sqrt{\pi})$; the constant is the one that gives the result below unit variance.

Step 3 used the fact that the product of two Gaussian bumps is again a Gaussian bump; Section B.4 states the general rule. Write $\ell = \sqrt{2}\,s$, and multiply every weight variance by a constant $\sigma_f^2$, which multiplies $k$ by $\sigma_f^2$. The result is

This is the RBF kernel (for radial basis function), also called the squared exponential kernel, and it is the kernel Chapter 8 works with. The same formula is the similarity measure of the RBF support vector machine, a classifier whose kernel width Chapter 22 tunes with Bayesian optimization. The construction is a standard one (Rasmussen and Williams, 2006, sec. 4.2.1). For an input with several coordinates, bumps in every direction give the same formula with $(x - x')^2$ replaced by the squared distance $\lVert \vx - \vx' \rVert^2$.

Of the whole feature model, two numbers remain. The lengthscale $\ell$ says how far apart two inputs must be before their values become nearly unrelated: the kernel falls to $0.61$ of its peak at distance $\ell$, to $0.14$ at $2\ell$, and to $0.011$ at $3\ell$. The amplitude $\sigma_f$ says how large the function is: $k(x, x) = \sigma_f^2$ is the prior variance of $f$ at any single input. A kernel with $\sigma_f = 1$, the setting of most of Chapter 8, is said to have unit amplitude.

The figure below is the model of Figure 7.1 with the number of features as a control. The readout gives the largest gap, over all inputs, between the covariance the bumps induce and the limit Equation (7.4).

A few experiments show how the limit is approached.

Start at 4 features. The band swells and pinches, and the solid curve in the strip changes shape as you drag $x_0$ across a bump. The gap to the RBF kernel is as large as the kernel itself, or larger.

Step up to 16. By 12 features the band is nearly flat and the gap is about 0.1; at 16 it is below 0.01, and the functions look like those drawn at ∞, where no features exist at all.

Shorten the lengthscale to 0.05. The bumps narrow, 16 of them are no longer enough, and the band pinches again. The number of bumps needed is roughly the width of the domain divided by the width of one bump.

That last experiment is the practical case for working with the kernel. On $[0, 1]$ with $\ell = 0.1$ a few dozen bumps would do. In $d$ input dimensions the bumps must tile a $d$-dimensional grid, so the count is raised to the power $d$: 16 bumps per axis in 10 dimensions is $16^{10} \approx 10^{12}$ features. The kernel Equation (7.4) costs $d$ subtractions and one exponential per pair of inputs, whatever the dimension.

Exchanging an explicit list of features for a function that computes their inner product directly is called the kernel trick (Rasmussen and Williams, 2006, sec. 2.1.2). It runs through the rest of the book. Section 8.1 will show that predictions from data, too, need the kernel and nothing else.

7.2.2 Which functions are kernels #

Once we stop writing features down, we need to know which functions $k$ may serve as covariances. Two conditions are necessary. The function must be symmetric, $k(\vx, \vx') = k(\vx', \vx)$, because covariance is. And for every finite set of inputs and every vector of coefficients $\mathbf{a}$,

because the left side is the variance of the weighted sum $\sum_i a_i f(\vx_i)$, and no variance is negative. A matrix with this property is positive semidefinite (Section 3.3), and a function whose matrices all have it is a positive semidefinite kernel.

Every inner product of features passes the test, since $\mathbf{a}^\T\boldsymbol{\Phi}\mSigma_p\boldsymbol{\Phi}^\T\mathbf{a} = (\boldsymbol{\Phi}^\T\mathbf{a})^\T\mSigma_p(\boldsymbol{\Phi}^\T\mathbf{a}) \ge 0$ for a covariance $\mSigma_p$. The converse is what licenses designing kernels directly: under mild technical conditions, every symmetric positive semidefinite kernel is the inner product of some list of features, possibly infinitely many, a result known as Mercer's theorem (Rasmussen and Williams, 2006, secs. 2.2 and 4.3), which Section 10.3 states with its conditions and uses to build the features from the kernel's eigenfunctions. We may therefore propose any function that passes the test and never write its features down.

The test is not a formality. Similarities that look reasonable can fail it. The "sigmoid" kernel $\tanh(a + b\,\vx^\T\vx')$, sometimes proposed by analogy with neural networks, is never positive definite, so it is not a valid covariance (Rasmussen and Williams, 2006, sec. 4.2.3). Exercise 7.2 shows a simpler failure and what it would mean: a combination of function values with negative variance.

7.3 Definition of a Gaussian process #

The previous sections arrived at a function $k$, and at the claim that it alone fixes a joint Gaussian for the function's values at any finite set of inputs. Taken as the starting point instead of a consequence, that claim is the definition.

A software engineer can read the definition as an interface. A Gaussian process is an object that accepts any finite list of inputs and answers with a mean vector and a covariance matrix, computed entry by entry from $m$ and $k$. It never has to produce the function at all inputs at once, and no computation in this book asks it to. In that sense it is a lazily evaluated random function: defined everywhere, materialized only where we look.

The answers to different queries must agree with each other. Query $\vx_1$ and $\vx_2$ together, ignore the second value, and the distribution of $f(\vx_1)$ must be what querying $\vx_1$ alone would have returned. For a Gaussian, ignoring coordinates means reading off a sub-block of the mean vector and the covariance matrix (Section 4.4). Each entry of $\mK$ depends only on its own two inputs, so that sub-block is exactly what the smaller query builds. Consistency is automatic for any covariance specified entry by entry through a kernel (Rasmussen and Williams, 2006, sec. 2.2).

The mean function is the guess before any data: where $f$ should be, input by input. In Bayesian optimization it is almost always zero or a constant, after the observed values have been standardized to mean zero and unit spread (Section 8.6). That is a weaker assumption than it sounds. Together with the kernel's amplitude, it says only that before any evaluation $f(\vx)$ lies within about $\pm 2\sigma_f$ of zero at each input, and the data move the posterior wherever they disagree. A mean function that is not constant is useful when there is real knowledge of a trend, such as a physical model whose errors the Gaussian process should learn. Unless stated otherwise, the rest of the book takes $m = 0$.

The kernel carries everything else: how strongly values at different inputs move together, and through that, how smooth, how wiggly, and how large the functions are. It must be symmetric and positive semidefinite (Section 7.2.2), and nothing more is required. The line kernel of Example 7.1 and the RBF kernel of Equation (7.4) both qualify.

Two words are easy to confuse. A Gaussian distribution describes a vector of fixed length. A Gaussian process describes a function, a collection of random variables indexed by the inputs, and any finite slice of it is a Gaussian distribution. The general name for a random collection indexed by a set is a stochastic process, and the set is its index set; in this book the index set is always the domain being optimized over. One draw of the whole function is called a sample path, or simply a sample.

7.4 Drawing functions #

A computer cannot draw a whole function. It can draw the function's values on a fine grid of inputs and connect them, which on a screen amounts to the same thing. By Equation (7.5) those values form a Gaussian vector with mean $\mathbf{m}$ and covariance $\mK$, so drawing them uses the recipe of Section 4.3: factor $\mK = \mL\mL^\T$ by Cholesky (Section 3.5), draw a vector $\vz$ of independent standard normal numbers, and return $\mathbf{m} + \mL\vz$. The result has covariance $\mL\,\E[\vz\vz^\T]\,\mL^\T = \mL\mL^\T = \mK$, as required.

The jitter in step 2 is the one of Section 3.4.4: a kernel matrix on a fine grid is nearly singular. The perturbation it adds to each sample is of order $\sqrt{\varepsilon}$, far below what a plot can show. The factorization costs $O(N^3)$ time and is paid once; each sample then costs $O(N^2)$.

import numpy as np

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

xs = np.linspace(0, 1, 200)
K = rbf(xs, xs)
L = np.linalg.cholesky(K + 1e-8 * np.eye(len(xs)))
rng = np.random.default_rng(0)
F = L @ rng.standard_normal((len(xs), 3))   # three samples, one per column

The figure runs Algorithm 7.1 on 201 grid points. It keeps the three standard normal vectors $\vz$ fixed while you change the kernel, so each change deforms the same three functions instead of drawing new ones.

Each control is worth a minute.

Drag the lengthscale from 0.3 down to 0.03. The three functions compress horizontally, like an accordion. The expected number of zero crossings in the readout grows in proportion to $1/\ell$, and the counts in the draws scatter around it.

Drag the amplitude. The functions stretch vertically and nothing else changes: not their shapes, not their crossings, not the correlation strip. Multiplying the kernel by $\sigma_f^2$ multiplies $\mL$ by $\sigma_f$, and nothing more.

Switch from RBF to Matérn 5/2 to Matérn 1/2 at a fixed lengthscale. These are kernels whose functions are rougher than the RBF's; Section 7.5.3 defines them. The Matérn 5/2 draws are hard to tell from the RBF ones except for a slightly rougher texture. The Matérn 1/2 draws are jagged at every scale, and the readout no longer gives an expected number of crossings.

Choose the periodic kernel. Every function repeats exactly, once per period, and the correlation strip returns to 1 at $x_0 \pm p$. If the lengthscale was below 0.4, it jumps to 0.8 when you switch, because for this kernel it is measured against the period (Section 7.5.4).

Move $x_0$. Where the strip is close to 1, each function's value near $x_0$ tracks its value at $x_0$, marked with a dot. Where the strip falls to zero, the two values are unrelated.

7.5 What the kernel encodes #

Every control in Figure 7.3 is a statement about the unknown function, made before any evaluation. The kernel's own parameters, the lengthscale, the amplitude, and the period, are its hyperparameters in the sense of Section 5.6: settings of the prior, kept apart from the function values the model is about. Chapter 9 fits them to data. This section reads what each one says.

7.5.1 Stationarity #

The RBF kernel depends on its two inputs only through their difference $x - x'$. A kernel with this property is stationary: shifting every input by the same amount leaves the prior unchanged, so the prior looks the same in every part of the domain. A stationary kernel can be written as a function of one argument, the difference $r = x - x'$, and we will write $k(r)$ when that is convenient. All four kernels in Figure 7.3 are stationary, which is why their bands have constant width. The six-bump model of Figure 7.1 was not stationary, because its centers made some inputs special. Nor is the line kernel $\sigma_0^2 + \sigma_1^2 x x'$ of Example 7.1, whose variance grows away from zero.

Stationary kernels are the default in Bayesian optimization. One consequence appears in Section 8.2: far from all data, the posterior returns to the same prior everywhere.

7.5.2 Lengthscale and amplitude #

For a stationary kernel the lengthscale is a frequency in disguise. A classical result in the theory of random processes gives the expected number of times a zero-mean stationary Gaussian process crosses zero upward on an interval of unit length (Rasmussen and Williams, 2006, sec. 4.1):

where $k''(0)$ is the second derivative of $k(r)$ at $r = 0$. The formula needs only the curvature of the kernel at zero distance, that is, how fast the correlation of two nearby values drops as they move apart. For the RBF kernel, $k(r) = \sigma_f^2 \exp(-r^2/2\ell^2)$ has $k''(0) = -\sigma_f^2/\ell^2$, so

With $\ell = 0.1$ a draw crosses zero upward 1.6 times on $[0, 1]$ on average; with $\ell = 0.03$, 5.3 times. The amplitude cancels, as the amplitude experiment in Figure 7.3 showed.

The lengthscale is measured in the units of the input. Measure time in milliseconds instead of seconds and the right lengthscale becomes a thousand times larger. That is why Section 8.6 recommends scaling every input to $[0, 1]$ before fitting: it makes a lengthscale of 0.1 mean the same thing in every problem. When the input has several coordinates that matter differently, each can have its own lengthscale (Section 9.2).

Even with scaled inputs, the same lengthscale does not mean the same thing in every dimension. Section 3.1.2 showed that two inputs drawn at random from the unit cube $[0, 1]^d$ are typically $\sqrt{d/6}$ apart: about $0.41$ in one dimension, $1.0$ in six, and $2.9$ in fifty. With $\ell = 0.2$, the RBF kernel at that distance, $\exp(-(d/6)/2\ell^2)$, is $0.12$ in one dimension, $4 \times 10^{-6}$ in six, and below $10^{-45}$ in fifty. In six dimensions, then, this prior already treats the values at two random inputs as unrelated, and each observation informs only a small neighborhood around itself, which is what Figure 6.3 showed from the side of information. Scaling the lengthscale with $\sqrt{d}$, here $\ell = 0.2\sqrt{d}$, keeps that correlation at $0.12$ in every dimension. In 2024, standard Bayesian optimization was shown to perform well in high dimensions once its lengthscale prior was scaled with the dimension (Hvarfner et al., 2024), and this arithmetic is one way to see why such scaling helps (inference). Section 9.5 and Section 30.1 tell that story.

The amplitude is simpler. $k(x, x) = \sigma_f^2$ is the prior variance of $f$ at each input, so before any data the model puts $f(x)$ within $\pm 1.96\,\sigma_f$ of the mean with 95% probability. Scaling a Gaussian process by a constant $c$ scales its kernel by $c^2$ (Exercise 7.3), which is why the amplitude only stretched the draws. With standardized outputs, an amplitude near 1 is the natural starting point. Libraries often store it as a separate "output scale" hyperparameter that multiplies a kernel of unit amplitude.

7.5.3 Smoothness #

The lengthscale says how fast a function varies. Smoothness says how rough it is at the smallest scales, and the two are independent: a function can wander slowly and still be jagged up close, like a coastline seen from a plane. Smoothness is set by the kernel's behavior very near zero distance. If $k(r)$ is smooth at $r = 0$, values at nearby inputs are almost perfectly correlated and the draws are smooth. If $k(r)$ has a corner at $r = 0$, nearby values decorrelate quickly and the draws are rough.

The RBF kernel is infinitely differentiable at zero, and its draws have derivatives of every order (in the mean-square sense, the one the theory of random processes uses) (Rasmussen and Williams, 2006, sec. 4.2.1). That is a strong assumption. The Matérn family, named by Stein (1999) after earlier work of Matérn, relaxes it with a smoothness parameter $\nu$: its draws are $q$ times differentiable exactly when $\nu > q$ (Rasmussen and Williams, 2006, sec. 4.2.1). The values used in practice are half-integers, for which the kernel is an exponential times a polynomial. The roughest, $\nu = 1/2$, is

Its corner at $r = 0$ produces draws that are continuous but nowhere differentiable. They are the paths of the Ornstein-Uhlenbeck process, first introduced as a model of the velocity of a particle buffeted by random collisions (Rasmussen and Williams, 2006, sec. 4.2.1). Then $\nu = 3/2$ gives once-differentiable draws, $\nu = 5/2$ twice-differentiable ones, and as $\nu \to \infty$ the Matérn kernel becomes the RBF. Section 9.1 writes out the family.

Equation (7.6) needs $k''(0)$, which the Matérn 1/2 kernel does not have. Its draws do not cross zero a finite expected number of times: a draw that crosses zero once crosses it infinitely often nearby (Rasmussen and Williams, 2006, sec. 4.1). That is why the readout in Figure 7.3 gives no expectation for this kernel, and why the count it reports depends on how fine the grid is.

Which smoothness to assume is a real modeling choice. Stein (1999) argued that the smoothness the RBF kernel assumes is unrealistic for many physical processes and recommended the Matérn class instead (Rasmussen and Williams, 2006, sec. 4.2.1). For Bayesian optimization, Snoek et al. (2012) made the same argument: draws from the RBF kernel "are unrealistically smooth for practical optimization problems", and they proposed Matérn 5/2, whose twice-differentiable draws match the assumption made by quasi-Newton methods, optimizers that estimate second derivatives from how the gradient changes between steps. They also tested the choice. Tuning three hyperparameters of a structured support vector machine that finds motifs in protein DNA sequences, about 40,000 of them, they repeated the optimization 100 times with each of several kernels and found that the choice of kernel "significantly affects performance": the RBF kernel's assumption of infinite differentiability was "too restrictive for this problem" (Snoek et al., 2012). In the other direction, values of $\nu$ above 5/2 are hard to tell apart from each other and from the RBF using finite, noisy data, which is one reason 3/2 and 5/2 are the values seen most often (Rasmussen and Williams, 2006, sec. 4.2.1).

7.5.4 Periodicity #

Some objectives repeat: an angle, a time of day, a hue on the color wheel. The periodic kernel builds the repetition into the prior,

with period $p$. A construction explains the formula. Map each input onto a circle, $\mathbf{u}(x) = \left(\cos(2\pi x/p),\, \sin(2\pi x/p)\right)$, and apply the RBF kernel to the mapped points. The squared distance between two mapped points is $4\sin^2(\pi(x - x')/p)$, which turns the RBF formula into Equation (7.7) (Rasmussen and Williams, 2006, sec. 4.2.3). Inputs exactly one period apart land on the same point of the circle, so their values have correlation 1, and every draw repeats exactly.

Because the RBF kernel is applied on the circle, the lengthscale $\ell$ of the periodic kernel is measured relative to the circle, not in the units of $x$. Near zero distance the kernel behaves like an RBF kernel with lengthscale $p\ell/(2\pi)$, so Equation (7.6) gives $1/(p\ell)$ upward crossings per unit length. That is why Figure 7.3 moves the lengthscale when you switch to this kernel. The color-preference figure of Chapter 19 uses a periodic kernel because hue wraps around.

7.5.5 What a stationary kernel cannot say #

A kernel is a strong statement, and some beliefs cannot be expressed with the kernels above. A stationary kernel cannot express a trend that continues beyond the data, because far from the data its posterior returns to the prior mean (Section 8.6). It cannot express an abrupt change, such as a setting beyond which a training run diverges, because its draws are continuous; even the Matérn 1/2 draws, rough as they are, have no jumps. And a single lengthscale cannot express a function that is flat in one region and wiggly in another. When these matter, practitioners add or multiply kernels (Section 9.1), transform the inputs before applying a kernel as the periodic kernel does, or turn to models that are not stationary.

The prior also matters more in Bayesian optimization than in most regression. A model fitted to thousands of points is shaped mostly by its data. An optimizer that can afford twenty evaluations of a ten-dimensional function has data almost nowhere, so the prior decides what it believes almost everywhere. That is why the next two chapters spend so long on the kernel: Chapter 8 shows how the prior and the data combine, and Chapter 9 shows how to let the data choose the hyperparameters.

7.6 Exercises #

Further reading #

• Rasmussen and Williams (2006), chapter 2, sets the weight-space and function-space views side by side, with the kernel trick and the definition used here; chapter 4 catalogs kernels, including the bump construction of the RBF kernel, the zero-crossing formula, and the Matérn class.
• Garnett (2023), chapters 2 and 3, develops Gaussian processes and the choice of kernel with Bayesian optimization in mind.
• Neal (1996), chapter 2, shows that neural networks with random weights become Gaussian processes as they grow infinitely wide.
• Stein (1999) makes the case, from the theory of spatial prediction, that the smoothness a kernel assumes matters, and argues for the Matérn class.
• Görtler et al. (2019) is an interactive essay in which the reader changes kernels and watches prior samples change, a companion to Figure 7.3.