Gaussian Processes from scratch: put a prior over functions, and a Cholesky solve gives you a band that pinches at your data

Ordinary regression fits weights. A Gaussian Process throws the weights away and keeps a distribution over whole functions — no fixed shape at all. You never pick basis functions; you pick a kernel (a notion of which inputs should give similar outputs) and the maths hands you back a mean curve, a shaded uncertainty band, and actual sampled functions. I built an interactive demo where you click to drop observations and watch the band breathe — and every curve is conditioned live in the browser with a real Cholesky solve. Here's how it works.

The kernel is the entire model

A GP is defined completely by its kernel: the covariance between the outputs at two inputs. The squared-exponential (RBF) kernel says nearby inputs are strongly correlated and that correlation decays smoothly with distance. Two knobs: σ_f², the signal variance (how far the function swings), and ℓ, the length-scale (how quickly two inputs stop mattering to each other). Slide ℓ and the fit goes from wiggly to smooth — that single choice is what makes the whole thing tick.

import numpy as np
def rbf(a, b, sf2=1.0, ell=1.0):
    d = a[:, None] - b[None, :]          # all pairwise input differences
    return sf2 * np.exp(-d**2 / (2 * ell**2))

Three kernel matrices, no design matrix

Conditioning needs three blocks of one big joint covariance: K between the observed inputs (plus observation noise on its diagonal), K* between the test inputs and the observed ones, and K** between the test inputs themselves. That's it — no weights, no basis functions written down anywhere.

K   = rbf(X, X)  + noise*np.eye(len(X))   # train-train  (N x N)
Ks  = rbf(Xs, X)                          # test-train   (M x N)
Kss = rbf(Xs, Xs)                         # test-test    (M x M)

Condition with a Cholesky — the posterior mean

The posterior mean is μ* = K*(K + σ_n²I)⁻¹y. You never form that inverse directly. Factor K = L Lᵀ once with a Cholesky — numerically stable and half the cost — and solve two triangular systems. The vector alpha is then reused for every test point. This is the exact solve the demo runs on each click and each slider move.

L     = np.linalg.cholesky(K)                          # K = L Lᵀ
alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))    # solve (K) alpha = y
mean  = Ks @ alpha                                     # posterior mean

The band pinches at data, reverts to the prior between

The posterior covariance is Σ* = K** − K*(K + σ_n²I)⁻¹K*ᵀ: start from the prior and subtract whatever the data explains. Its diagonal is the per-point variance; the square root, times two, is the shaded band. At an observation the subtracted term nearly cancels the prior, so the band pinches to a noise-limited waist; far from any point K*→0, nothing is subtracted, and it relaxes back to the prior height σ_f. That behaviour is the whole payoff — the model knows where it does not know.

v    = np.linalg.solve(L, Ks.T)              # reuse L
cov  = Kss - v.T @ v                         # K** - K* (K)^-1 K*ᵀ
std  = np.sqrt(np.clip(np.diag(cov), 0, None))
band_lo, band_hi = mean - 2*std, mean + 2*std

This is the sharpest contrast with Bayesian linear regression, which was the previous day's build. There a straight-line model's uncertainty grows without bound as you extrapolate. A GP's band instead dips to a noise floor at every observation and rises back to a flat ceiling between and beyond points — it never claims to know more than the prior where it has no data.

Sampling whole functions makes it visceral

The band summarises the spread; drawing actual functions makes it click. A sample is a draw from the joint Gaussian N(mean, Σ*) over the whole test grid — Cholesky-factor Σ* (with a whisper of jitter for stability) and push standard-normal noise through it. These are the grey curves in the demo, and every single one threads exactly through your observations. Where they fan apart is the uncertainty.

Ls = np.linalg.cholesky(cov + 1e-6*np.eye(len(Xs)))   # jitter keeps it PD
for _ in range(10):
    f = mean + Ls @ np.random.randn(len(Xs))          # one plausible function
    plt.plot(Xs, f, color="grey", alpha=0.3)

The one-liner, and where a GP earns its keep

scikit-learn ships all of this as GaussianProcessRegressor, and it goes further: hand it a kernel and it learns the length-scale, signal variance and noise by maximising the log marginal likelihood — no hand-tuning of sliders. Ask for return_std=True and you get calibrated error bars for free.

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel as C
kernel = C(1.0) * RBF(length_scale=1.0) + WhiteKernel(noise_level=0.02)
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True).fit(X[:, None], y)
mean, std = gp.predict(Xs[:, None], return_std=True)   # predictions WITH error bars

A GP is really Bayesian linear regression with an infinite basis — the kernel trick lets you pick a kernel and never write the basis down at all. The cost is the O(N³) Cholesky, so for big N you reach for sparse or inducing-point variants. Use one when data is small-to-medium, smoothness is a fair assumption, and an honest uncertainty band actually matters: forecasting, active learning, Bayesian optimisation. And that band isn't just decoration — the natural next step is to use it to decide where to sample next.

Click on the chart to drop points and watch the band pinch:
https://dev48v.infy.uk/ml/day37-gaussian-processes.html