Hyperparameter Learning in Gaussian Processes

Having established the marginal likelihood as our objective for hyperparameter learning, we now face a computational challenge: how do we efficiently optimize it?

The marginal likelihood landscape is typically non-convex with multiple local optima. Gradient-free methods (grid search, random search, evolutionary algorithms) can work for low-dimensional hyperparameter spaces but become prohibitively expensive as dimensionality increases. With multiple length-scales (one per input dimension), periodic kernels, or composite kernel structures, the hyperparameter space easily reaches 10-50 dimensions.

Gradient-based optimization is the solution. Methods like L-BFGS, conjugate gradient, or Adam can efficiently navigate high-dimensional spaces when gradients are available. For Gaussian Processes, we're fortunate: the log marginal likelihood has closed-form gradients.

Let's establish the notation we'll use throughout this derivation.

The Log Marginal Likelihood (LML):

$$\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)$$

where:

• $\mathbf{y} \in \mathbb{R}^n$ — Training observations
• $\mathbf{K}_y = \mathbf{K} + \sigma_n^2 \mathbf{I}$ — Noisy covariance matrix
• $\mathbf{K} = [k(\mathbf{x}_i, \mathbf{x}j)]{i,j=1}^n$ — Kernel matrix
• $\boldsymbol{\theta}$ — Vector of all hyperparameters (kernel parameters + noise variance)

Our Goal:

Compute $\frac{\partial \mathcal{L}}{\partial \theta_j}$ for each hyperparameter $\theta_j$. This enables gradient-based optimization:

$$\boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + \eta \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}^{(t)})$$

Key Quantities:

We'll repeatedly use:

• $\boldsymbol{\alpha} = \mathbf{K}_y^{-1} \mathbf{y}$ — The 'dual representation' of predictions
• $\frac{\partial \mathbf{K}_y}{\partial \theta_j}$ — Derivative of covariance matrix w.r.t. hyperparameter $\theta_j$

The structure of the gradient will depend on how each hyperparameter enters the covariance matrix. Fortunately, all common kernels have simple, differentiable forms.

Before diving into the gradient derivation, we need several matrix calculus identities. These are essential tools for differentiating expressions involving matrices.

Identity 1: Derivative of Matrix Inverse

For a matrix $\mathbf{A}$ that depends on a scalar $\theta$:

$$\frac{\partial \mathbf{A}^{-1}}{\partial \theta} = -\mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial \theta} \mathbf{A}^{-1}$$

Proof sketch: Differentiate $\mathbf{A}^{-1} \mathbf{A} = \mathbf{I}$ using the product rule.

Identity 2: Derivative of Log Determinant

For a positive definite matrix $\mathbf{A}$:

$$\frac{\partial \log |\mathbf{A}|}{\partial \theta} = \text{tr}\left( \mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial \theta} \right)$$

Proof sketch: Use $\log|\mathbf{A}| = \sum_i \log \lambda_i$ and the sensitivity of eigenvalues, or the identity $\frac{\partial |\mathbf{A}|}{\partial A_{ij}} = |\mathbf{A}| (\mathbf{A}^{-1})_{ji}$.

Identity 3: Derivative of Quadratic Form

For vectors $\mathbf{a}, \mathbf{b}$ independent of $\theta$:

$$\frac{\partial}{\partial \theta} \left( \mathbf{a}^\top \mathbf{A}^{-1} \mathbf{b} \right) = -\mathbf{a}^\top \mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial \theta} \mathbf{A}^{-1} \mathbf{b}$$

This follows directly from Identity 1.

Identity 4: Trace Cyclic Property

For matrices of compatible dimensions:

$$\text{tr}(\mathbf{ABC}) = \text{tr}(\mathbf{CAB}) = \text{tr}(\mathbf{BCA})$$

This allows rearranging terms inside traces, which is crucial for efficient computation.

Identity 5: Trace and Outer Products

For vectors $\mathbf{a}, \mathbf{b}$ and matrix $\mathbf{M}$:

$$\mathbf{a}^\top \mathbf{M} \mathbf{b} = \text{tr}(\mathbf{M} \mathbf{b} \mathbf{a}^\top) = \text{tr}(\mathbf{a}^\top \mathbf{M} \mathbf{b})$$

This allows converting quadratic forms to traces when convenient.

We now derive the gradient $\frac{\partial \mathcal{L}}{\partial \theta_j}$ for an arbitrary hyperparameter $\theta_j$.

Starting Point:

$$\mathcal{L} = -\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)$$

The third term is constant with respect to hyperparameters. Let's differentiate the first two.

Gradient of the Data Fit Term:

$$\frac{\partial}{\partial \theta_j}\left( -\frac{1}{2}\mathbf{y}^\top \mathbf{K}_y^{-1} \mathbf{y} \right)$$

Using Identity 3:

$$= -\frac{1}{2} \cdot \left( -\mathbf{y}^\top \mathbf{K}_y^{-1} \frac{\partial \mathbf{K}_y}{\partial \theta_j} \mathbf{K}_y^{-1} \mathbf{y} \right)$$

$$= \frac{1}{2} \mathbf{y}^\top \mathbf{K}_y^{-1} \frac{\partial \mathbf{K}_y}{\partial \theta_j} \mathbf{K}_y^{-1} \mathbf{y}$$

Defining $\boldsymbol{\alpha} = \mathbf{K}_y^{-1} \mathbf{y}$:

$$= \frac{1}{2} \boldsymbol{\alpha}^\top \frac{\partial \mathbf{K}_y}{\partial \theta_j} \boldsymbol{\alpha}$$

Gradient of the Complexity Term:

$$\frac{\partial}{\partial \theta_j}\left( -\frac{1}{2}\log |\mathbf{K}_y| \right)$$

Using Identity 2:

$$= -\frac{1}{2} \text{tr}\left( \mathbf{K}_y^{-1} \frac{\partial \mathbf{K}_y}{\partial \theta_j} \right)$$

Combining the Terms:

$$\boxed{\frac{\partial \mathcal{L}}{\partial \theta_j} = \frac{1}{2} \boldsymbol{\alpha}^\top \frac{\partial \mathbf{K}_y}{\partial \theta_j} \boldsymbol{\alpha} - \frac{1}{2} \text{tr}\left( \mathbf{K}_y^{-1} \frac{\partial \mathbf{K}_y}{\partial \theta_j} \right)}$$

Or equivalently, using the trace-quadratic identity ($\boldsymbol{\alpha}^\top \mathbf{A} \boldsymbol{\alpha} = \text{tr}(\boldsymbol{\alpha} \boldsymbol{\alpha}^\top \mathbf{A})$):

$$\boxed{\frac{\partial \mathcal{L}}{\partial \theta_j} = \frac{1}{2} \text{tr}\left( \left( \boldsymbol{\alpha} \boldsymbol{\alpha}^\top - \mathbf{K}_y^{-1} \right) \frac{\partial \mathbf{K}_y}{\partial \theta_j} \right)}$$

Let's develop intuition for what each gradient term represents.

Term 1: Data Fit Contribution

$$\frac{1}{2} \boldsymbol{\alpha}^\top \frac{\partial \mathbf{K}_y}{\partial \theta_j} \boldsymbol{\alpha}$$

Recall $\boldsymbol{\alpha} = \mathbf{K}_y^{-1} \mathbf{y}$ represents the 'residuals in the dual space.' This term measures how the data fit changes when we perturb $\theta_j$:

• If increasing $\theta_j$ makes $\mathbf{K}_y$ 'align better' with the outer product $\mathbf{y}\mathbf{y}^\top$, this term is positive
• It encourages hyperparameters that make the prior covariance match the empirical covariance of the data

Term 2: Complexity Contribution

$$-\frac{1}{2} \text{tr}\left( \mathbf{K}_y^{-1} \frac{\partial \mathbf{K}_y}{\partial \theta_j} \right)$$

This term measures how the model complexity changes with $\theta_j$:

• If increasing $\theta_j$ increases $|\mathbf{K}_y|$ (the prior volume), this term is negative
• It penalizes hyperparameters that expand the space of functions unnecessarily

The Balance:

At the optimum, these two forces are in equilibrium:

$$\boldsymbol{\alpha}^\top \frac{\partial \mathbf{K}_y}{\partial \theta_j} \boldsymbol{\alpha} = \text{tr}\left( \mathbf{K}_y^{-1} \frac{\partial \mathbf{K}_y}{\partial \theta_j} \right)$$

For each hyperparameter, the marginal gain in data fit from increasing it equals the marginal cost in complexity. This is the differential form of Occam's Razor.

Connection to Fisher Information:

The complexity term $\text{tr}(\mathbf{K}_y^{-1} \frac{\partial \mathbf{K}_y}{\partial \theta_j})$ is related to the Fisher information of the Gaussian distribution with covariance $\mathbf{K}_y$. Hyperparameters with large Fisher information (the model is sensitive to their values) are more heavily penalized—the model 'pays more' for parameters that strongly affect predictions.

The gradient formula requires $\frac{\partial \mathbf{K}_y}{\partial \theta_j}$—the derivative of the covariance matrix with respect to each hyperparameter. This decomposes as:

$$\frac{\partial \mathbf{K}_y}{\partial \theta_j} = \frac{\partial \mathbf{K}}{\partial \theta_j} + \frac{\partial (\sigma_n^2 \mathbf{I})}{\partial \theta_j}$$

The noise derivative is simple: $\frac{\partial (\sigma_n^2 \mathbf{I})}{\partial \sigma_n^2} = \mathbf{I}$ if $\theta_j = \sigma_n^2$, and zero otherwise.

The kernel derivative $\frac{\partial \mathbf{K}}{\partial \theta_j}$ depends on the specific kernel form. Let's work through the most common cases.

Squared Exponential (RBF) Kernel:

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

Derivative w.r.t. signal variance $\sigma_f^2$: $$\frac{\partial k}{\partial \sigma_f^2} = \frac{k(\mathbf{x}, \mathbf{x}')}{\sigma_f^2}$$

Derivative w.r.t. length-scale $\ell$: $$\frac{\partial k}{\partial \ell} = k(\mathbf{x}, \mathbf{x}') \cdot \frac{r^2}{\ell^3}$$

For log-parameterization (common in practice), with $\bar{\ell} = \log \ell$: $$\frac{\partial k}{\partial \bar{\ell}} = \frac{\partial k}{\partial \ell} \cdot \ell = k(\mathbf{x}, \mathbf{x}') \cdot \frac{r^2}{\ell^2}$$

Matérn Kernels:

$$k_\nu(r) = \sigma_f^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{\sqrt{2\nu} r}{\ell} \right)^\nu K_\nu\left( \frac{\sqrt{2\nu} r}{\ell} \right)$$

The derivatives are more complex due to the Bessel function. For the commonly used special cases:

Matérn-1/2 (Exponential): $k(r) = \sigma_f^2 \exp(-r/\ell)$ $$\frac{\partial k}{\partial \ell} = k(r) \cdot \frac{r}{\ell^2}$$

Matérn-3/2: $k(r) = \sigma_f^2 \left(1 + \frac{\sqrt{3}r}{\ell}\right) \exp\left(-\frac{\sqrt{3}r}{\ell}\right)$ $$\frac{\partial k}{\partial \ell} = \sigma_f^2 \frac{3r^2}{\ell^3} \exp\left(-\frac{\sqrt{3}r}{\ell}\right)$$

Matérn-5/2: $k(r) = \sigma_f^2 \left(1 + \frac{\sqrt{5}r}{\ell} + \frac{5r^2}{3\ell^2}\right) \exp\left(-\frac{\sqrt{5}r}{\ell}\right)$ $$\frac{\partial k}{\partial \ell} = \sigma_f^2 \frac{5r^2}{3\ell^3} \left(1 + \frac{\sqrt{5}r}{\ell}\right) \exp\left(-\frac{\sqrt{5}r}{\ell}\right)$$

Automatic Relevance Determination (ARD) uses a separate length-scale for each input dimension, enabling automatic feature selection:

$$k(\mathbf{x}, \mathbf{x}') = \sigma_f^2 \exp\left( -\frac{1}{2} \sum_{d=1}^D \frac{(x_d - x'_d)^2}{\ell_d^2} \right)$$

This introduces $D$ length-scale parameters ${\ell_1, \ldots, \ell_D}$, and we need gradients for each.

Derivation:

Define the weighted squared distance: $$r_d = (x_d - x'_d)^2$$

Then: $$\frac{\partial k}{\partial \ell_d} = k(\mathbf{x}, \mathbf{x}') \cdot \frac{r_d}{\ell_d^3}$$

With log-parameterization ($\bar{\ell}_d = \log \ell_d$): $$\frac{\partial k}{\partial \bar{\ell}_d} = k(\mathbf{x}, \mathbf{x}') \cdot \frac{r_d}{\ell_d^2}$$

Interpreting ARD Gradients:

After optimization, the learned length-scales reveal feature importance:

• Large $\ell_d$: Input dimension $d$ has little effect on predictions; the function varies slowly along this axis. This dimension is effectively 'turned off.'

• Small $\ell_d$: Input dimension $d$ strongly influences predictions; the function varies rapidly along this axis. This dimension is highly relevant.

The gradient structure reveals why ARD works: for an irrelevant feature (constant along dimension $d$), the squared distance term $r_d$ is small for all pairs, making the gradient negligible. The length-scale then grows large with little marginal likelihood change. For relevant features, $r_d$ varies substantially, and the gradient drives $\ell_d$ to an appropriate value.

Computing gradients naively is expensive. Each hyperparameter's gradient requires:

$$\frac{\partial \mathcal{L}}{\partial \theta_j} = \frac{1}{2} \text{tr}\left( \left( \boldsymbol{\alpha} \boldsymbol{\alpha}^\top - \mathbf{K}_y^{-1} \right) \frac{\partial \mathbf{K}_y}{\partial \theta_j} \right)$$

Computing $\mathbf{K}_y^{-1}$ explicitly costs $O(n^3)$ and produces an $n \times n$ matrix. With $p$ hyperparameters, naive computation costs $O(pn^3)$, which is prohibitive for large datasets.

Efficient Strategy:

1. Cholesky decomposition once: Factor $\mathbf{K}_y = \mathbf{L}\mathbf{L}^\top$ in $O(n^3)$ time.

2. Solve for $\boldsymbol{\alpha}$ once: Compute $\boldsymbol{\alpha} = \mathbf{K}_y^{-1}\mathbf{y}$ via back-substitution in $O(n^2)$ time.

3. Solve for $\mathbf{K}_y^{-1}$ column-by-column: If needed, compute each column of $\mathbf{K}_y^{-1}$ by solving $\mathbf{L}\mathbf{L}^\top \mathbf{b}_i = \mathbf{e}_i$. Total: $O(n^3)$.

4. Compute each gradient via trace: Use the relation: $$\text{tr}(\mathbf{A}\mathbf{B}) = \sum_{i,j} A_{ij} B_{ji}$$

   This is element-wise multiply and sum: $O(n^2)$ per hyperparameter.

Complexity Analysis:

The $O(n^3)$ Cholesky and inverse computations dominate, but they're done once per evaluation. Additional hyperparameters add only $O(n^2)$ cost each.

Modern numerical libraries achieve peak performance through vectorization. Here we present an optimized implementation that minimizes memory allocations and maximizes cache efficiency.

Key Optimizations:

1. Avoid explicit $\mathbf{K}_y^{-1}$: In many cases, we can restructure computations to avoid forming the full inverse.

2. Hadamard product for traces: Use $\text{tr}(\mathbf{A}\mathbf{B}) = \mathbf{1}^\top (\mathbf{A} \odot \mathbf{B}^\top) \mathbf{1}$ where $\odot$ is element-wise multiplication.

3. Pre-compute shared structures: Matrices like $\boldsymbol{\alpha}\boldsymbol{\alpha}^\top$ are computed once.

4. Log-parameterization: Working with $\log \theta$ instead of $\theta$ is numerically stable and eliminates positivity constraints.

When implementing GP gradients, numerical verification is essential. Analytic gradient errors can cause silent failures—the optimization may converge to a suboptimal point without any obvious error message.

Finite Difference Approximation:

The simplest approach uses centered finite differences:

$$\frac{\partial \mathcal{L}}{\partial \theta_j} \approx \frac{\mathcal{L}(\theta_j + \epsilon) - \mathcal{L}(\theta_j - \epsilon)}{2\epsilon}$$

with $\epsilon \approx 10^{-5}$ to $10^{-7}$.

Relative Error Check:

$$\text{relative error} = \frac{|\text{analytic} - \text{numerical}|}{\max(|\text{analytic}|, |\text{numerical}|, 10^{-8})}$$

A relative error below $10^{-4}$ is typically acceptable; below $10^{-6}$ is excellent.

We've developed the complete machinery for computing gradients of the log marginal likelihood. Let's consolidate the key results:

Looking Ahead:

With gradients in hand, we can apply gradient-based optimization algorithms to find hyperparameters that maximize the marginal likelihood. The next page explores optimization strategies—which algorithms to use, how to handle multiple optima, and practical considerations for robust hyperparameter learning.