Learn AI Create the Future — with Pranay
Main site

Classic Algorithms in Depth

Gaussian process regression

A Gaussian process predicts a value and an honest error bar together, by treating "smooth curves through my data" as the model itself.

Read these first

On this page 5
  1. Why it exists
  2. How it works
  3. A real example you have seen
  4. Remember this
  5. What to learn next

One lesson, three depths. Pick the one that fits you today — you can switch any time.

Beginner — No maths. Plain English.

Gaussian process regression predicts a value and, in the same breath, tells you how sure it is.

Imagine asking a sensible friend to guess the afternoon temperature on your roof. She has readings from the morning and the evening, but nothing from midday. Near her readings, she answers firmly: "about 20 degrees, give or take one."

For 2 pm she still answers. But she widens her hands: "somewhere between 20 and 27, honestly."

Most regression models give you the guess without the hand-widening. A Gaussian process (GP) gives both, and the widening is computed, not vibes.

Why it exists

An ordinary regression line commits to one formula and reports one number per question — equally confidently everywhere, including regions it has never seen. For a doctor's dosage model or an engineer's safety margin, confident guesses in blind spots are dangerous.

A Gaussian process starts differently. Instead of one formula, imagine every smooth curve that passes near your data points. Where the curves bunch together, the answer is settled — narrow band. Where they fan out, in the gaps between readings, the model is honestly unsure — wide band.

The only ingredient you owe it is a definition of "smooth": how strongly nearby inputs should agree. That definition is a kernel, the same similarity-scoring idea from the kernel trick.

How it works

temp │      ╭─────╮ ← the band of curves that fit the readings
     │  ●──╮│     │╭───●
     │      ╰──╮  ││        band is NARROW near readings ●
     │         ╰──╯│        band is WIDE in the midday gap
     └────────────────────── hour
        readings ●    gap    readings ●

A real example you have seen

Weather apps that say "27°, likely between 25 and 30" are doing this kind of reasoning. GPs also sit inside hyperparameter tuning tools — the software that picks a model's dials for you. They model "score as a function of settings". The uncertainty band then decides which setting to try next. That is the trick behind Bayesian optimisation in tools like Optuna.

Remember this

  • A GP predicts a band, not a point: value plus honest uncertainty.
  • The band is narrow near data, wide in gaps — exactly where it should be.
  • The kernel encodes your one assumption: how smooth the world is.

What to learn next

Developer — Code and libraries.

Setup

bash
pip install scikit-learn

Outputs verified with scikit-learn 1.7.2 on CPU.

Seven readings and a midday blind spot

gp_roof.py
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

# hour of day -> roof temperature in deg C; nothing between 8:00 and 20:00
X = np.array([[0.0], [2.0], [4.0], [6.0], [8.0], [20.0], [22.0]])
y = np.array([21.0, 20.0, 19.0, 19.5, 22.0, 26.0, 23.5])

# RBF: "hours close together have similar temperatures, fading over ~3 hours"
# WhiteKernel: "each reading carries sensor noise of about half a degree"
kernel = RBF(length_scale=3.0) + WhiteKernel(noise_level=0.5)
gp = GaussianProcessRegressor(kernel=kernel, optimizer=None, normalize_y=True).fit(X, y)

for hour in (5.0, 9.0, 14.0, 21.0):
    mean, std = gp.predict([[hour]], return_std=True)
    print(f"{int(hour):2d}:00 -> {mean[0]:5.1f} deg C, 95% range +/- {1.96 * std[0]:.1f}")
print("log marginal likelihood:", round(gp.log_marginal_likelihood_value_, 2))
Output
 5:00 ->  19.8 deg C, 95% range +/- 3.8
 9:00 ->  21.7 deg C, 95% range +/- 4.3
14:00 ->  22.1 deg C, 95% range +/- 5.5
21:00 ->  24.2 deg C, 95% range +/- 3.8
log marginal likelihood: -8.92

The walkthrough

Watch the error bars, not the means. At 5:00 and 21:00, snug between readings, the band is ±3.8. At 14:00 — six hours from the nearest reading — it swells to ±5.5. Nobody programmed "be less sure at midday". It fell out of the maths, which is the entire selling point.

The kernel is the model. RBF(length_scale=3.0) declares that temperatures 3 hours apart are strongly related and 10 hours apart barely at all. WhiteKernel(noise_level=0.5) declares the sensor itself wobbles. Change these declarations and predictions change — a GP is opinionated about smoothness because you told it what to believe.

Why optimizer=None? By default, scikit-learn tunes kernel settings by maximising the log marginal likelihood — a built-in score of how well the kernel explains the data. On seven points that automation can chase extremes (for instance, deciding the sensor is noiseless) and emit a ConvergenceWarning. Fixing the kernel keeps this demo honest and repeatable. With ~30+ points, drop optimizer=None and let it tune; compare fits via gp.log_marginal_likelihood_value_ — higher is better.

normalize_y=True centres and scales the targets internally. Without it, a GP quietly assumes your values hover around zero, and temperatures near 22 confuse the prior.

Common mistakes

Reading the band as a guarantee. ±5.5 is the model's uncertainty given its assumptions. A midday sun-spike more violent than the kernel's idea of smoothness will fall outside the band. Garbage assumptions, confident garbage bands.

No noise term on noisy data. Omit the WhiteKernel (or alpha) and the GP treats every reading as gospel, threading the curve through each wobble — overfitting with error bars attached. Real sensors always deserve a noise term.

Using exact GPs on big data. Training cost grows with the cube of the number of points. A thousand points is fine; fifty thousand freezes your machine. For large data, use sparse approximations (GPyTorch, GPflow) or a different model family.

Forgetting to scale inputs in multi-feature GPs. One length_scale shared across features with wildly different units makes "nearby" meaningless. Scale features, or give each its own length scale via RBF(length_scale=[1.0, 1.0]) — which also enables automatic relevance detection.

Try it yourself

Add a midday reading — 12.0 hours, 31.0 degrees — and rerun. Predict beforehand: what happens to the 14:00 band? Then change length_scale to 0.5 and to 10.0, and describe in one sentence what each setting believes about roof temperature.

What to learn next

Researcher — Mathematics and papers.

Definition

A Gaussian process is a collection of random variables, any finite subset of which is jointly Gaussian. It is specified by a mean function $m(x)$ (usually 0 after centring) and covariance kernel $k(x, x')$:

$$ f \sim \mathcal{GP}(m, k) $$

For training inputs $X$ (with targets $y$, noise variance $\sigma_n^2$) and test point $x_*$, the posterior predictive is Gaussian with:

$$ \mu_* = k_^\top (K + \sigma_n^2 I)^{-1} y, \qquad \sigma_^2 = k(x_, x_) - k_^\top (K + \sigma_n^2 I)^{-1} k_ $$

Where:

  • $K$ — the $n \times n$ Gram matrix, $K_{ij} = k(x_i, x_j)$.
  • $k_$ — the vector of covariances between $x_$ and each training point.
  • $\sigma_n^2$ — observation noise (the WhiteKernel term).
  • $\mu_, \sigma_^2$ — predictive mean and variance.

Note $\sigma_*^2$ does not depend on $y$ at all — uncertainty is a function of where you have data, not of the values measured there.

Hyperparameters and the marginal likelihood

Kernel parameters $\theta$ (length scale $\ell$, signal variance, noise) are set by maximising the log marginal likelihood:

$$ \log p(y \mid X, \theta) = -\tfrac{1}{2} y^\top (K_\theta + \sigma_n^2 I)^{-1} y - \tfrac{1}{2} \log |K_\theta + \sigma_n^2 I| - \tfrac{n}{2} \log 2\pi $$

The first term rewards fitting the data; the second penalises complex (large-determinant) covariances — an automatic Occam's razor, no validation set required. The surface is non-convex; restarts matter (n_restarts_optimizer).

Cost

Exact inference: $O(n^3)$ time for the Cholesky factorisation, $O(n^2)$ memory, then $O(n)$ per predictive mean and $O(n^2)$ per variance. Scaling routes: sparse inducing-point methods — FITC (Snelson and Ghahramani, 2006), variational SVGP (Titsias, 2009; Hensman et al., 2013) at $O(nm^2)$ for $m \ll n$ inducing points — and structure-exploiting solvers in GPyTorch (Gardner et al., 2018).

Connections and standing

  • Textbook: Rasmussen and Williams (2006), Gaussian Processes for Machine Learning — free online, still the reference.
  • Kernel ridge regression computes the same mean without the variance; the GP is the Bayesian reading of the same algebra.
  • A single-hidden-layer network with infinite width converges to a GP (Neal, 1996); deep versions birthed the NNGP and neural-tangent-kernel literature (Lee et al., 2018; Jacot et al., 2018).
  • Workhorse of Bayesian optimisation (Snoek et al., 2012): the GP posterior drives acquisition functions like expected improvement — the machinery behind modern hyperparameter tuners.

What to learn next