ML8.5 · LESSON 5 OF 6

Gaussian processes

Update a prior, compare posterior means with MAP estimates and interpret predictive uncertainty.

THE IDEA

A prior over functions

So far the unknown was one number. Now it is a whole function f(x). A Gaussian process says: pick any set of inputs, and the function values there are jointly Gaussian with a mean (here 0) and a covariance given by a kernel k(x, x'). The RBF kernel is k = exp(-(x - x')^2 / (2 l^2)). Inputs closer than the length-scale l have values that move together; far apart, they are nearly independent.

The kernel is an assumption about smoothness: it is where you decide that nearby inputs should have similar outputs. A short length-scale allows wiggly functions, a long one forces gentle curves.

Conditioning on data

Given observations y at inputs X, with noise variance sigma^2, the posterior at test inputs has

GP regression

K   = k(X, X) + (sigma^2 + jitter) I        # training covariance
k*  = k(X, x*)                               # test point against training
mean = k*^T K^-1 y
var  = k(x*, x*) - k*^T K^-1 k*              # latent function variance

In code you never form K^-1: you solve a linear system, usually with a Cholesky factorisation. The small jitter (1e-8 here) keeps K numerically invertible, for example when two training inputs coincide and K without noise would be singular. For a single noise-free observation y = 1 at x = 5 with length-scale 1, the exact posterior is: at x = 5, mean 1 and sd 0.0001; one step away (x = 6), mean 0.6065 and sd 0.7951; at x = 9 the sd is 1, back to the prior's 1.

posterior mean and band of 2 standard deviationsobservationsprior mean 0

Click the plot to add an observation (x snaps to 0.5, y to 0.25), or use the sliders below.

length-scale1new x5new y1

noise sd

2 observations, length-scale 1, noise 0.1 solve (K + (0.1^2 + jitter) I) with jitter = 1e-8 at x = 5: mean -0.352, sd 0.797 at x = 9: mean -0.007, sd 1

Click to add observations. The band shows plus or minus 2 standard deviations of the latent function. Change the length-scale and the noise.

What to notice

The band pinches at observations and returns to the prior (width set by the kernel) far from data.
With noise 0 the mean passes through every point. With noise above 0 it smooths, and the band stays open at the points.
A longer length-scale carries information further: one point at x = 3 can move the mean at x = 6.

Predict first. You add 5 more observations at the same x where you already have one. What happens to the band at that x (noise sd 0.3)?

Chapter context

Probabilistic modeling represents uncertainty explicitly instead of presenting every fitted number as certain. Priors, likelihoods and approximation choices are assumptions that should be exposed rather than hidden behind a library call.

TRY IT YOURSELF

Part 1. Do it hands-on

Add observations by clicking or with the sliders. The latent-function sd and mean in the tasks are computed with the same Cholesky solve that you see plotted.

posterior mean and band of 2 standard deviationsobservationsprior mean 0

Click the plot to add an observation (x snaps to 0.5, y to 0.25), or use the sliders below.

length-scale1new x5new y1

noise sd

0 observations, length-scale 1, noise 0 solve (K + (0^2 + jitter) I) with jitter = 0 at x = 5: mean 0, sd 1 at x = 9: mean 0, sd 1

To do: With length-scale 1 or less, make the uncertainty small at x = 5 (sd below 0.25) but still large at x = 9 (sd above 0.9)
To do: With exactly one observation at x = 3, make the posterior mean at x = 6 exceed 0.5 by choosing the length-scale and the height
To do: Set noise to 0.5 and add two conflicting observations at x = 5, heights 1.5 and -1.5: the mean there stays within 0.1 of 0

Part 2. Do it in code

rbf_kernel(x1, x2, length_scale=1.0, variance=1.0): matrix of variance * exp(-(x1_i - x2_j)^2 / (2 length_scale^2)), shape (len(x1), len(x2)). Raise ValueError unless length_scale > 0 and variance > 0.
gp_posterior(x_train, y_train, x_test, length_scale=1.0, noise=0.0, jitter=1e-8): return (mean, var) of the latent function at x_test. Solve with np.linalg.solve on K + (noise^2 + jitter) I, never an explicit inverse, and clip the variance at 0. With no training points return the prior (zeros and the signal variance 1).
predictive_interval(mean, var, noise, z=1.96): the interval for a new noisy observation, mean -/+ z * sqrt(var + noise^2). Return (lower, upper) arrays.

import numpy as np

def rbf_kernel(x1, x2, length_scale=1.0, variance=1.0):
    pass

def gp_posterior(x_train, y_train, x_test, length_scale=1.0, noise=0.0, jitter=1e-8):
    pass

def predictive_interval(mean, var, noise, z=1.96):
    pass

K = rbf_kernel(np.array([0.0, 1.0]), np.array([0.0, 1.0]))
print(K)

Ctrl/⌘ + Enter runs it

Stuck, or something not working? You can come back to it any time.

CHECK YOUR UNDERSTANDING

Does a GP kernel encode assumptions about how function values relate?

AYesBNo

Chapter reference: math, worked example & code

This section belongs to the whole chapter, so it is the same on every lesson in it and may cover a different lesson's topic. The code for the experiment you just ran is inside the experiment itself.

Math and notation

For Bernoulli success probability theta and a Beta prior:

p(\theta\mid D)=\operatorname{Beta}(a+s,b+f),\quad E[\theta\mid D]=\frac{a+s}{a+b+s+f}.

For Beta(alpha,beta) with both parameters greater than one, the MAP is $(\alpha-1)/(\alpha+\beta-2)$ .

For a zero-mean Gaussian-process regression model with independent noise variance sigma-n squared, a test prediction has

\mu_*=k_*^T(K+\sigma_n^2I)^{-1}y,

v_*=k_{**}-k_*^T(K+\sigma_n^2I)^{-1}k_*.

K is the training covariance matrix, k-star connects the test input with training inputs, and k-star-star is its prior variance. This is latent-function variance; observation noise adds another term for a noisy predicted observation. Practical code solves systems instead of explicitly computing an inverse.

Worked example

Start with Beta(2,2), then observe three successes and one failure. The posterior is Beta(5,3). Its mean is 5/8=0.625. Its MAP is 4/6≈0.6667. Both summarize the same distribution but answer different questions.

Reference code

The following example uses the Python standard library and does not call a remote service.

from math import isclose
a, b = 2+3, 2+1
assert (a,b) == (5,3)
assert a/(a+b) == .625
assert isclose((a-1)/(a+b-2), 2/3)

Python examples are provided for study and local execution; this page does not execute Python.

Common mistakes

A posterior estimate is conditional on the chosen model. MAP is not the same as a posterior mean. A small uncertainty under a misspecified model need not mean a reliable real-world answer.