Learning Kernel Hyperparameters in Gaussian Processes

A Gaussian process kernel like the squared-exponential has knobs — a lengthscale $\ell$, an amplitude $\sigma^2$, and a noise level $\sigma_n^2$ — and the fitted model can look wildly different depending on where those knobs are set. Set $\ell$ too small and the GP snakes through every noisy wiggle in the training data; set it too large and it flattens into a nearly straight line. Neither extreme is chosen by hand-tuning in practice — it is chosen by directly maximizing a number the data itself hands you: the marginal likelihood.

The analogy: a scorecard that punishes both underfitting and showing off

Imagine judging two witnesses describing the same event. One gives an almost too-perfect account that fits every tiny detail — plausible, but suspiciously overfit to what actually happened, unlikely to generalize to a slightly different retelling. The other gives a vague, generic account that fits nothing specifically. A good judge wants the account that fits the specific evidence well without being needlessly elaborate to do it. The marginal likelihood is exactly that judge for a GP: it rewards hyperparameters that fit the observed data closely, but it also automatically penalizes settings that only fit by being needlessly flexible — no separate validation set required, because the penalty is built into the formula itself.

The formula, term by term

The log marginal likelihood for a GP with data $y$, covariance matrix $K$ (built from the kernel), and noise variance $\sigma_n^2$ is:

$$\log p(y) = -\frac{1}{2} y^\top (K + \sigma_n^2 I)^{-1} y \; - \; \frac{1}{2}\log\left|K + \sigma_n^2 I\right| \; - \; \frac{n}{2}\log(2\pi)$$

The first term is the data-fit term — how well the current covariance explains the observed $y$ values; it gets better as the model fits the data more tightly. The second term is the complexity penalty — the log-determinant grows when the kernel is very flexible (small $\ell$ produces a covariance matrix closer to a scaled identity, which has a larger determinant), automatically punishing over-flexible settings. The third term is a constant that does not depend on the hyperparameters. Maximizing this sum, by gradient ascent over $\ell$, $\sigma^2$, and $\sigma_n^2$, finds the setting that best balances "fits the data" against "isn't needlessly wiggly to do it."

Worked example 1: two lengthscales, two scores

Two observations at $x=0, y=1$ and $x=1, y=1.1$, comparing $\ell=0.5$ against $\ell=5$, both with $\sigma^2=1$, $\sigma_n^2 = 0.01$. At $\ell=0.5$: $k(0,1)=\exp(-1/(2\times0.25)) = 0.135$, giving a covariance matrix with small off-diagonal terms — the model treats these two points as nearly independent, so the data-fit term is unremarkable and the complexity penalty (log-determinant) is large, since the matrix is close to $\sigma_n^2 I$-scaled and flexible. At $\ell = 5$: $k(0,1) = \exp(-1/50) = 0.980$ — near-perfect correlation, matching how close $y=1$ and $y=1.1$ actually are, giving a much better data-fit term with only a small complexity cost. Here the smoother, larger-lengthscale setting wins the marginal likelihood comparison because it explains this particular pair of nearby, similar values more efficiently.

Worked example 2: overfitting made concrete

With noisy training data $y = \{1.0, 1.3, 0.9, 1.4, 1.0\}$ at evenly spaced $x = 0,1,2,3,4$ and a held-out point at $x=2.5$ with true value near $1.15$: a too-small lengthscale ($\ell=0.2$) fits the training points almost exactly (data-fit term looks great) but predicts wildly at $x=2.5$ — say $0.6$, an error of $0.55$ — because it has learned the noise, not the trend. A moderate lengthscale ($\ell=1.5$) fits training points slightly less perfectly but predicts around $1.1$ at $x=2.5$, an error of only $0.05$. The marginal likelihood favors the moderate setting even without ever seeing the held-out point, because its complexity penalty already flags the tiny-lengthscale fit as suspiciously flexible.

This is literally the optimizer used here — the hyperparameters $\ell$, $\sigma^2$, $\sigma_n^2$ are climbed via gradient ascent on the marginal likelihood surface exactly like a loss surface, complete with the risk of getting stuck at a poor local optimum.

What this means in practice

This is why GPs are prized for small-to-medium datasets in finance — volatility surface interpolation, sparse alternative-data series — where a held-out validation set would be too small to trust, but the marginal likelihood can still be computed exactly from all available points. The optimization is non-convex, though, so different starting values for $\ell$ can converge to different local optima with meaningfully different fitted behavior.