Hyperparameter Learning in Gaussian Processes

Gaussian Processes define probability distributions over functions. But a GP is not fully specified until we choose its hyperparameters—the parameters of the mean function $m(\mathbf{x})$ and covariance function $k(\mathbf{x}, \mathbf{x}')$ that determine the prior beliefs about the functions we're modeling.

Consider the ubiquitous Squared Exponential (RBF) kernel:

$$k(\mathbf{x}, \mathbf{x}') = \sigma_f^2 \exp\left(-\frac{|\mathbf{x} - \mathbf{x}'|^2}{2\ell^2}\right)$$

This kernel has two hyperparameters: the signal variance $\sigma_f^2$ (controlling the amplitude of functions) and the length-scale $\ell$ (controlling smoothness). Add the observation noise variance $\sigma_n^2$, and we have three hyperparameters demanding careful selection.

How should we choose these? Unlike model parameters which can be learned through direct optimization of a training objective, hyperparameters govern the shape of our prior—and poor choices lead to catastrophically wrong uncertainty estimates.

In a fully Bayesian treatment, hyperparameters $\boldsymbol{\theta}$ (comprising kernel parameters, noise variance, and mean function parameters) would be integrated out rather than point-estimated. Given training data $\mathcal{D} = {(\mathbf{x}i, y_i)}{i=1}^n$, the predictive distribution would be:

$$p(y_* | \mathbf{x}*, \mathcal{D}) = \int p(y* | \mathbf{x}_*, \mathcal{D}, \boldsymbol{\theta}) , p(\boldsymbol{\theta} | \mathcal{D}) , d\boldsymbol{\theta}$$

This integral marginalizes over all possible hyperparameter values, weighted by their posterior probability given the data. However, this integral is typically intractable for GP hyperparameters, as the posterior $p(\boldsymbol{\theta} | \mathcal{D})$ doesn't have a closed form.

Practical Resolution: The standard approach is to find a point estimate $\hat{\boldsymbol{\theta}}$ by maximizing the marginal likelihood (also called the evidence or type-II maximum likelihood):

$$\hat{\boldsymbol{\theta}} = \arg\max_{\boldsymbol{\theta}} , p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta})$$

This effectively selects hyperparameters that make the observed data most probable under the GP prior.

Why not just maximize fit to training data?

A naive approach would be to choose hyperparameters that minimize training error or maximize the standard likelihood. This leads to severe overfitting:

• Setting the length-scale $\ell$ very small allows the GP to interpolate training points exactly
• Setting the noise variance $\sigma_n^2$ to near-zero forces the GP through every data point
• The resulting model has excellent training fit but terrible generalization

The marginal likelihood naturally penalizes such overly complex solutions, as we'll see in the derivation.

Let's derive the marginal likelihood from first principles. In GP regression with Gaussian observation noise, we model:

$$y_i = f(\mathbf{x}_i) + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$$

where $f$ is drawn from a GP prior: $f \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))$.

Step 1: Joint Distribution of Observations

Given the GP prior, the function values at training points $\mathbf{f} = [f(\mathbf{x}_1), \ldots, f(\mathbf{x}_n)]^\top$ follow a multivariate Gaussian:

$$\mathbf{f} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{K})$$

where $\boldsymbol{\mu} = [m(\mathbf{x}_1), \ldots, m(\mathbf{x}n)]^\top$ and $K{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$.

The observations $\mathbf{y}$ are obtained by adding independent Gaussian noise:

$$\mathbf{y} | \mathbf{f} \sim \mathcal{N}(\mathbf{f}, \sigma_n^2 \mathbf{I})$$

Step 2: Marginalizing Over Function Values

The marginal likelihood is obtained by integrating out the latent function values:

$$p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta}) = \int p(\mathbf{y} | \mathbf{f}) , p(\mathbf{f} | \mathbf{X}, \boldsymbol{\theta}) , d\mathbf{f}$$

Since both $p(\mathbf{y} | \mathbf{f})$ and $p(\mathbf{f} | \mathbf{X}, \boldsymbol{\theta})$ are Gaussian, this integral has a closed-form solution. The marginal distribution of $\mathbf{y}$ is:

$$\mathbf{y} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{K} + \sigma_n^2 \mathbf{I}) = \mathcal{N}(\boldsymbol{\mu}, \mathbf{K}_y)$$

where we define $\mathbf{K}_y = \mathbf{K} + \sigma_n^2 \mathbf{I}$ as the covariance matrix of the observations.

Step 3: The Log Marginal Likelihood

The log of the multivariate Gaussian density gives us the log marginal likelihood (LML):

$$\log p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta}) = -\frac{1}{2}(\mathbf{y} - \boldsymbol{\mu})^\top \mathbf{K}_y^{-1} (\mathbf{y} - \boldsymbol{\mu}) - \frac{1}{2}\log |\mathbf{K}_y| - \frac{n}{2}\log(2\pi)$$

For convenience (and without loss of generality), assume a zero mean function $\boldsymbol{\mu} = \mathbf{0}$:

$$\boxed{\mathcal{L}(\boldsymbol{\theta}) = \log p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta}) = -\frac{1}{2}\mathbf{y}^\top \mathbf{K}_y^{-1} \mathbf{y} - \frac{1}{2}\log |\mathbf{K}_y| - \frac{n}{2}\log(2\pi)}$$

This is the central equation for GP hyperparameter learning. Let's examine each term:

The marginal likelihood automatically implements Occam's Razor—the principle that simpler explanations should be preferred over complex ones, all else being equal. This emerges naturally from the mathematical structure, not as an imposed regularization.

Geometric Intuition:

Consider the space of all possible datasets. A simple model (e.g., a GP with large length-scale, modeling only smooth functions) concentrates its probability mass on a restricted set of datasets—those exhibiting smooth patterns. A complex model (small length-scale, capable of fitting wiggly functions) spreads its probability mass across a much larger space of possible datasets.

For any particular dataset:

• If the data is genuinely simple (smooth), the simple model assigns it high probability because it's in the model's 'sweet spot'
• The complex model assigns it lower probability (despite being able to fit it) because probability mass is diluted across all the complex patterns it could fit

Conversely:

• If the data is genuinely complex, the simple model can't fit it well and assigns low probability
• The complex model assigns reasonable probability because such data is within its modeling capability

Mathematical Manifestation:

Consider how the two non-constant terms trade off:

1. Overly simple models (large $\ell$, small $\sigma_f^2$):

   • Cannot fit the data well → data fit term is very negative
   • Low complexity → complexity penalty is mild
   • Net result: low marginal likelihood due to poor fit

2. Overly complex models (small $\ell$, large $\sigma_f^2$):

   • Can fit the data perfectly → data fit term is only mildly negative
   • High complexity → complexity penalty is severe (large $\log|\mathbf{K}_y|$)
   • Net result: low marginal likelihood due to complexity penalty

3. Appropriately complex models:

   • Fit the data reasonably well → moderate data fit term
   • Complexity matched to data → moderate complexity penalty
   • Net result: optimal marginal likelihood

The marginal likelihood is maximized by models whose complexity matches the true complexity of the data-generating process.

Let's develop deeper intuition for each term in the log marginal likelihood.

The Data Fit Term: $-\frac{1}{2}\mathbf{y}^\top \mathbf{K}_y^{-1} \mathbf{y}$

This term is the negative of a Mahalanobis distance—it measures how 'typical' the observations $\mathbf{y}$ are under the GP prior with covariance $\mathbf{K}_y$.

The inverse covariance $\mathbf{K}_y^{-1}$ (precision matrix) acts as a weighting that accounts for correlations between observations. Points that are modeled as highly correlated should have similar values; deviations are more severely penalized than for weakly correlated points.

Key insight: This is not simply a sum of squared residuals. The matrix structure means:

• If $\mathbf{K}_y$ predicts strong correlations and the data exhibits them → low Mahalanobis distance
• If $\mathbf{K}_y$ predicts correlations that aren't present (or vice versa) → high Mahalanobis distance

The Complexity Penalty: $-\frac{1}{2}\log |\mathbf{K}_y|$

The determinant $|\mathbf{K}_y|$ has a beautiful geometric interpretation: it represents the volume of the probability ellipsoid defined by the Gaussian prior.

For diagonal covariance (independent observations): $$|\mathbf{K}y| = \prod{i=1}^n (K_{y,ii}) = \prod_{i=1}^n (k(\mathbf{x}_i, \mathbf{x}_i) + \sigma_n^2)$$

Each factor represents the prior variance at one observation point. The product gives the volume of the n-dimensional hyperrectangle within which functions are likely.

For full covariance (correlated observations): Correlations reduce the volume—if knowing $y_1$ tells us something about $y_2$, the effective 'space' of possibilities shrinks. The determinant captures this reduction.

Connection to eigenvalues: $$\log|\mathbf{K}y| = \sum{i=1}^n \log \lambda_i$$

where $\lambda_i$ are eigenvalues of $\mathbf{K}_y$. More complex models have larger eigenvalues (each eigenvalue represents variance along a principal direction), leading to larger determinants and more severe penalties.

Understanding how each hyperparameter influences the marginal likelihood provides crucial intuition for debugging GP models and interpreting fitted values.

Length-scale $\ell$:

The length-scale controls the 'wiggliness' of functions sampled from the GP.

Signal Variance $\sigma_f^2$:

The signal variance controls the overall amplitude of functions.

Noise Variance $\sigma_n^2$:

The noise variance is particularly interesting because it affects both terms in opposite directions.

The noise variance essentially controls the signal-to-noise ratio the model assumes. Too low, and the GP overfits noise as signal. Too high, and it ignores genuine patterns.

An alternative to marginal likelihood optimization is selecting hyperparameters via cross-validation (CV). Both approaches have merit, and understanding their differences illuminates what the marginal likelihood is fundamentally doing.

Leave-One-Out Cross-Validation (LOO-CV):

For GP regression, LOO-CV has a remarkable closed-form expression:

$$\text{LOO-CV} = \sum_{i=1}^n \log p(y_i | \mathbf{X}, \mathbf{y}_{-i}, \boldsymbol{\theta})$$

where $\mathbf{y}_{-i}$ denotes all observations except the $i$-th. This measures how well the model predicts held-out points when trained on the rest.

Using the GP predictive equations and the matrix inversion lemma, this can be computed efficiently:

$$\text{LOO-CV} = -\frac{1}{2}\sum_{i=1}^n \left( \log\sigma_i^2 + \frac{(y_i - \mu_i)^2}{\sigma_i^2} \right)$$

where $\mu_i$ and $\sigma_i^2$ are the LOO predictive mean and variance for point $i$.

When do they differ?

In practice, marginal likelihood and LOO-CV often select similar hyperparameters. They diverge when:

1. Model misspecification is severe: If the true data-generating process is far from any GP prior, LOO-CV may find hyperparameters that minimize prediction error even if they don't correspond to a 'good' probabilistic model.

2. Small datasets: The marginal likelihood's complexity penalty is calibrated to the Bayesian framework; LOO-CV can underpenalize complexity when data is scarce.

3. Heteroscedastic or non-Gaussian noise: LOO-CV focuses on prediction accuracy regardless of how the noise model is specified; marginal likelihood requires correct noise modeling for valid results.

As a principled default, maximizing the marginal likelihood is preferred for GP hyperparameter learning. Cross-validation serves as a valuable sanity check and alternative when model assumptions are suspect.

Evaluating the log marginal likelihood requires computing:

1. $\mathbf{K}_y^{-1} \mathbf{y}$ — A linear system solve
2. $\log |\mathbf{K}_y|$ — The log determinant

Both operations are dominated by inverting or factorizing the $n \times n$ covariance matrix.

Cholesky Decomposition Approach:

The standard approach is to compute the Cholesky decomposition:

$$\mathbf{K}_y = \mathbf{L} \mathbf{L}^\top$$

where $\mathbf{L}$ is lower triangular. This decomposition costs $O(n^3)$ time and $O(n^2)$ memory.

Once $\mathbf{L}$ is available:

• Solve $\mathbf{K}_y^{-1} \mathbf{y}$: Forward substitution for $\mathbf{z} = \mathbf{L}^{-1}\mathbf{y}$, then backward substitution for $\boldsymbol{\alpha} = \mathbf{L}^{-\top}\mathbf{z}$. Cost: $O(n^2)$.

• Compute $\log|\mathbf{K}_y|$: Since $|\mathbf{K}y| = |\mathbf{L}|^2$ and $|\mathbf{L}| = \prod_i L{ii}$, we have $\log|\mathbf{K}y| = 2\sum_i \log L{ii}$. Cost: $O(n)$.

The final expression using the Cholesky factor:

$$\mathcal{L}(\boldsymbol{\theta}) = -\frac{1}{2}|\mathbf{L}^{-1}\mathbf{y}|^2 - \sum_i \log L_{ii} - \frac{n}{2}\log(2\pi)$$

Understanding the shape of the marginal likelihood as a function of hyperparameters is crucial for successful optimization.

General Characteristics:

1. Non-convexity: The marginal likelihood is in general non-convex in hyperparameters. Multiple local maxima can exist, corresponding to qualitatively different models (e.g., 'smooth function + high noise' vs. 'wiggly function + low noise').

2. Log-scale parameterization: Hyperparameters like $\ell$, $\sigma_f^2$, $\sigma_n^2$ are strictly positive and often span multiple orders of magnitude. Working with $\log \ell$, $\log \sigma_f^2$, $\log \sigma_n^2$ instead is standard practice.

3. Flat regions: In sparsely sampled regions of input space, the marginal likelihood can be insensitive to length-scale (any scale larger than the data spacing has similar effect).

4. Sharp ridges: The noise-signal variance trade-off often creates sharp ridges in the marginal likelihood surface, making optimization challenging.

Common Pathologies:

Local vs. Global Optima:

Different local optima can correspond to genuinely different explanations of the data:

• Explanation 1: The true function is smooth with little noise; observations are nearly exact.
• Explanation 2: The true function is complex with much noise; observations are noisy samples.

Both might achieve reasonable marginal likelihood values. Domain knowledge or external validation can help distinguish between such alternatives.

The marginal likelihood connects to several other model selection criteria, providing a unifying perspective.

Connection to AIC and BIC:

The Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) approximate the marginal likelihood under certain assumptions:

$$\text{AIC} = -2\log p(\mathbf{y} | \hat{\boldsymbol{\theta}}) + 2k$$ $$\text{BIC} = -2\log p(\mathbf{y} | \hat{\boldsymbol{\theta}}) + k \log n$$

where $k$ is the number of parameters. For GP kernels, $k$ typically refers to hyperparameters. The marginal likelihood can be viewed as a more principled version that doesn't require ad hoc counting of 'effective parameters.'

Connection to Minimum Description Length (MDL):

MDL selects models that maximize compression of the data. The marginal likelihood has an MDL interpretation: $-\log p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta})$ represents the number of bits needed to encode the data using the model. Better models provide shorter codes; the marginal likelihood seeks the shortest description.

Connection to Regularization:

Maximizing the log marginal likelihood can be rewritten as:

$$\max_{\boldsymbol{\theta}} \left[ \text{Data Fit}(\boldsymbol{\theta}) - \frac{1}{2}\log|\mathbf{K}_y(\boldsymbol{\theta})| \right]$$

The second term functions as a data-dependent regularizer. Unlike fixed regularizers (L1, L2), this regularizer adapts to the data: more complex datasets justify more complex models without excessive penalty.

Connection to PAC-Bayes Bounds:

In learning theory, PAC-Bayes bounds provide generalization guarantees based on the KL divergence between posterior and prior. The marginal likelihood appears naturally in these bounds:

$$\text{Generalization Error} \leq \text{Empirical Error} + \sqrt{\frac{\text{KL}(\text{posterior} | \text{prior}) + \log n}{n}}$$

Maximizing marginal likelihood implicitly minimizes this bound under certain conditions, providing theoretical justification for its use in model selection.

We've developed a comprehensive understanding of the marginal likelihood as the cornerstone of GP hyperparameter learning. Let's consolidate the key insights:

Looking Ahead:

With the marginal likelihood established as our objective, the next page addresses how to actually optimize it. We'll derive the gradients of the marginal likelihood with respect to hyperparameters—the mathematical machinery that enables efficient gradient-based optimization for GP hyperparameter learning.