Gaussian process regression

A Gaussian process (GP) is a probability distribution over functions. Before we see data, it describes which functions we find plausible: smooth or rough, repeating or not, flat or with a trend. This belief is set by one function, the kernel $k(x, x')$. After we see data, Bayes' rule gives the posterior: the functions that agree with both the belief and the data. A GP gives more than one curve. It gives a mean prediction and an honest error bar at every $x$.

In the playground you add data points with a tap and see four things live:

• the posterior mean $\mu(x)$, the best guess for the function,
• the 95% band $\mu(x) \pm 1.96\,\sigma(x)$, where the true function lies with probability 0.95 (if the model is right),
• three posterior samples: random functions that are plausible given the data,
• the log marginal likelihood (the evidence), a single number that says how well the kernel explains the data. Two buttons tune the kernel by making this number as large as possible.

1. The playground #

Interactive figure: loading…

2. Try this #

1. From prior to posterior. Press Clear. The band is the same everywhere: $\pm 1.96\,\sigma_f$. Now add one point near $(0, 1)$. Question: how wide is the band at $x = 0$, at $x = 1$ and at $x = 2$? What to see: at the point the band shrinks to about $\pm 0.2$ (the noise $\sigma_n = 0.1$ is all that is left). At $x = 1$ the standard deviation is still $0.8$. At $x = 2$ it is $0.99$, so the band is almost the prior band again. One point teaches the GP only about its neighbourhood, and the size of the neighbourhood is the length scale $\ell = 1$.
2. The length scale. Pick 5 points and move $\ell$ from $0.3$ to $3$. Question: what happens between the points, for example at $x = 2$? Which value does the evidence prefer? What to see: with $\ell = 0.3$ the mean returns to $0$ between the points and the band is as wide as the prior: the points are too far apart to say anything about each other. With $\ell = 3$ the curve is stiff and the band is thin everywhere, but the evidence drops from about $-5.9$ to about $-20$: such a confident, smooth curve explains these points badly. Now press Gradient ascent. It ends near $\ell = 0.35$ with $\log p \approx -5.4$ (mostly by making $\sigma_f$ smaller), and it may say that the evidence is almost flat there. It is: every $\ell$ below about $0.4$ gives almost the same evidence. Five points far apart cannot tell these length scales apart. Here the evidence rules out long length scales, but not short ones.
3. Two explanations of noisy data. Pick noisy (30 points, true noise $0.25$) and press Gradient ascent. Set Band for to $f(x)$ and then to $y$. Then move the $\sigma_n$ slider all the way to the left ($\sigma_n = 0.01$) and the $\ell$ slider to about $0.15$, and press Gradient ascent again. Question: do both runs end at the same place? What to see: the first run climbs from $\log p \approx -20.4$ to $-10.2$ in about 35 steps and ends near $\ell = 1.0$, $\sigma_n = 0.19$: a smooth function plus noise (the true noise is $0.25$; 30 points give only a rough estimate). There, the band for $f$ misses 7 of the 30 points, and the band for $y = f(x) + \varepsilon$ contains all of them. The second run ends at $\ell \approx 0.36$, $\sigma_n = 0.01$ (the lower limit), $\log p \approx -14.6$: a wiggly function that passes through every point. Both are local maxima of the evidence. The data prefer the first by about 4.4 nats, a likelihood ratio of about $e^{4.4} \approx 80$.
4. Kernels are beliefs: extrapolation. Pick periodic. The data stop near $x = 0.4$. Choose the RBF kernel and press Gradient ascent. Then pick periodic under Try again (this resets the hyperparameters and selects the periodic kernel) and press Gradient ascent. Question: what does each model predict for $x$ between $2$ and $5$? What to see: the RBF model ends near $\ell = 0.18$ with $\log p \approx -18.2$. To the right of the data its mean falls back to $0$ and its band becomes the prior band: it knows nothing there. The periodic model finds the period $p \approx 1.60$ (the true one) with $\log p \approx 0.35$, and it continues the pattern to $x = 5$ with a thin band. The evidence prefers it by about 18 nats. Both runs say "not fully converged". More presses raise $\log p$ only a little: plain gradient ascent is slow on the flat ridge that is left.
5. A local maximum and grid search. Pick trend + season: a straight trend plus a pattern with period $1.5$, with the Sum kernel and a wrong start $p = 2$. Press Gradient ascent, then Grid search, then Gradient ascent again. Question: where does each step end? What to see: the first run ends at $p \approx 3.0$ with $\log p \approx 0.2$ (press again if it says "not fully converged"; then it stops at a local maximum). This is a real local maximum: a pattern that repeats every $1.5$ also repeats every $3$. Grid search tries all $36 \times 28$ values of $(p, \sigma_n)$ on the map and jumps to $p \approx 1.45$; the second gradient ascent climbs to $p \approx 1.48$ with $\log p \approx 6.8$. Now switch to Linear and to RBF. The linear model continues the trend but not the season; the RBF model continues neither.
6. A jump breaks the belief. Pick step: a function that jumps at $x = 0.5$, with noise $0.06$. Press Gradient ascent with RBF and then with Matérn 3/2. Question: what noise level do the models choose? What to see: RBF ends at $\sigma_n \approx 0.43$ and Matérn 3/2 at $\sigma_n \approx 0.37$, six to seven times the true noise. Both kernels believe that the function is equally smooth everywhere (they are stationary), so the best they can do is a long length scale plus a lot of "noise" that absorbs the jump. The rougher Matérn kernel explains the data a little better ($\log p \approx -17.6$ against $-18.1$).

3. Kernels are prior beliefs #

A GP with mean $0$ and kernel $k$ is written $f \sim \mathcal{GP}(0, k)$. It means: for any finite set of inputs $x_1, \dots, x_n$, the vector of function values is Gaussian,

So the kernel is a covariance function: $k(x, x')$ is the covariance of $f(x)$ and $f(x')$. If $k(x, x')$ is large, the values at $x$ and $x'$ move together; if it is near $0$, they are independent. Every choice of kernel is a belief about the function. The matrix $K$ must be positive semi-definite for every set of inputs, so not every function of two variables is a kernel.

Interactive figure: loading…

The kernels of the playground ($r = |x - x'|$):

Kernel	$k(x, x')$	Belief about $f$
RBF (squared exponential)	$\sigma_f^2 \exp\!\big(-r^2 / 2\ell^2\big)$	very smooth: all derivatives exist; changes on a scale $\ell$
Matérn 3/2	$\sigma_f^2 \big(1 + \sqrt3\,r/\ell\big) \exp\!\big(-\sqrt3\,r/\ell\big)$	rougher: one derivative, not two
Periodic	$\sigma_f^2 \exp\!\big(-2\sin^2(\pi r / p) / \ell^2\big)$	repeats exactly with period $p$
Linear	$\sigma_f^2 (1 + x x')$	straight lines $a + bx$ with $a, b \sim \mathcal N(0, \sigma_f^2)$
Sum	$k_{\text{per}} + \sigma_{\text{lin}}^2 (1 + x x')$	a periodic pattern plus an independent trend

• The length scale $\ell$ is the distance over which the function changes noticeably. For the RBF kernel, the correlation between $f(x)$ and $f(x + \ell)$ is $e^{-1/2} \approx 0.61$.
• The signal standard deviation $\sigma_f$ sets the typical size of $f$: for the first three kernels $k(x, x) = \sigma_f^2$ at every $x$.
• Stationary kernels depend only on $x - x'$, so the belief is the same everywhere. The linear kernel is not stationary: its variance $\sigma_f^2(1 + x^2)$ grows with $|x|$.
• Kernels can be combined. If $f_1 \sim \mathcal{GP}(0, k_1)$ and $f_2 \sim \mathcal{GP}(0, k_2)$ are independent, then $f_1 + f_2 \sim \mathcal{GP}(0, k_1 + k_2)$. Products of kernels are kernels too.

4. The math in short #

The posterior #

We observe $y_i = f(x_i) + \varepsilon_i$ with independent noise $\varepsilon_i \sim \mathcal N(0, \sigma_n^2)$. Write $K = k(X, X)$ for the $n \times n$ matrix of the data inputs, $K_y = K + \sigma_n^2 I$, and $k_* = k(X, x_*)$ for the vector of covariances between the data and a new input $x_*$. Under the prior, $\mathbf y$ and $f(x_*)$ are jointly Gaussian:

Conditioning a Gaussian on part of its vector (a Schur complement) gives the posterior, again Gaussian:

• The mean is a weighted sum of kernel bumps: $\mu(x_*) = \sum_i \alpha_i\, k(x_i, x_*)$ with $\alpha = K_y^{-1}\mathbf y$.
• The variance is the prior variance minus what the data explain. It does not depend on $\mathbf y$, only on where the data are. Far from all data $k_* \approx 0$, so $\sigma^2 \approx k(x_*, x_*)$: back to the prior.
• This is the band for $f(x_*)$. A new observation $y_* = f(x_*) + \varepsilon_*$ has variance $\sigma^2(x_*) + \sigma_n^2$ (the setting Band for: y).
• Posterior samples on a grid come from the full posterior covariance $\Sigma = K_{**} - K_{*X} K_y^{-1} K_{X*}$: draw $z \sim \mathcal N(0, I)$ and compute $\mu + Sz$ with $SS^\top = \Sigma$.

Computing it. Never form $K_y^{-1}$. Factor once with the Cholesky decomposition $K_y = LL^\top$, then use triangular solves:

The evidence and how to tune the kernel #

The hyperparameters $\theta = (\ell, \sigma_f, \sigma_n, \dots)$ are chosen by maximizing the log marginal likelihood: the log probability of the observed $\mathbf y$ when $f$ is integrated out. Since $\mathbf y \sim \mathcal N(0, K_y)$,

The two terms pull in opposite directions. A very short length scale lets $f$ pass through every point, so the data-fit term is good. But then $K_y$ is close to $(\sigma_f^2 + \sigma_n^2) I$: the model spreads its probability over a huge set of wiggly functions, $\log\det K_y$ is large, and the penalty is large. A very long length scale is the opposite. The evidence picks the simplest model that explains the data: an automatic form of Occam's razor. This is why the evidence in experiment 2 rejects $\ell = 3$, and why it can have several local maxima (experiment 3).

The gradient has a closed form. With $\alpha = K_y^{-1}\mathbf y$ and using $d\log\det K = \operatorname{tr}(K^{-1} dK)$:

The playground optimizes $u = \log\theta$ (so every hyperparameter stays positive and all of them move on the same scale) with plain gradient ascent and a simple line search: $u \leftarrow u + \eta\, \nabla_u \log p$. A step is kept only if the evidence increases; otherwise $\eta$ is made smaller and the step is tried again. After a good step $\eta$ grows by a factor $1.5$, and no log value moves more than $0.25$ in one step. Each hyperparameter stays inside its slider range. The run stops when the evidence rises by less than $10^{-4}$ in 50 steps; a test with the second derivatives (the Hessian) then says whether this is a maximum or an almost flat ridge. Plain gradient ascent is slow when the evidence is much more curved in one direction than in another. For the periodic kernel, the evidence is very sharp in $p$ and flat in $\ell$, so it can need several presses. Gradient ascent finds only a local maximum near the start. Grid search evaluates the evidence at all $36 \times 28$ points of the map (log-spaced values of two hyperparameters, the others fixed) and takes the best one: slow, but it cannot be trapped. Libraries such as scikit-learn and GPyTorch use L-BFGS or Adam with several random restarts.

5. The cost: $O(n^3)$ #

The Cholesky factorization of the $n \times n$ matrix $K_y$ needs about $n^3/3$ floating-point operations, and the matrix needs $n^2$ numbers of memory. After that, the mean at a new point costs $O(n)$ and the variance $O(n^2)$. Every new value of the hyperparameters needs a new factorization, so each step of gradient ascent costs $O(n^3)$, and a grid search costs (number of grid points) $\times\, O(n^3)$.

Interactive figure: loading…

This is why exact GPs are used for up to about $10^4$ points. For larger data sets there are approximations:

• Sparse GPs (inducing points): summarize the data with $m \ll n$ pseudo-points. The cost drops to $O(nm^2)$.
• Iterative solvers: conjugate gradients with fast matrix–vector products (GPyTorch). Each iteration costs $O(n^2)$ or less.
• Structure: in one dimension, a Matérn GP can be written as a state-space model, and the Kalman filter solves it in $O(n)$.
• Random features: replace the kernel by an explicit finite feature map, which gives Bayesian linear regression.

6. Related lessons #

• Positive definite matrices and Cholesky: why $K$ must be positive semi-definite, the Cholesky factor, and the Schur complement behind Gaussian conditioning.
• Factorizations and LU: solving with a factorization instead of an inverse, and low-rank (Nyström, inducing-point) approximations.
• Matrix calculus: $d\log\det K = \operatorname{tr}(K^{-1}dK)$ and $d(K^{-1}) = -K^{-1}\,dK\,K^{-1}$, the two rules behind the gradient of the evidence.
• Floating point and conditioning: why $K$ with a long length scale is badly conditioned, and why a small "jitter" on the diagonal helps.
• Eigenvalue algorithms and Krylov methods: conjugate gradients for large GPs.
• Gaussian integrals and Laplace's method: the marginal likelihood as a Gaussian integral.
• Line search and quasi-Newton methods: the optimizers used to tune GP hyperparameters.
• State-space models and the Kalman filter and the Kalman filter playground: a Matérn GP in time is a state-space model.
• Generalization in deep learning: the neural tangent kernel and the link between wide networks and kernel methods.
• ML for Science: a full lesson on Gaussian processes, scoring rules and calibration. ML Fundamentals: kernels and support vector machines. Probability and Statistics: the multivariate Gaussian and Bayesian inference.