Hyperparameter Learning in Gaussian Processes

With the marginal likelihood objective and its gradients established, we now tackle the optimization itself. While the mathematical machinery looks clean, practical optimization of GP hyperparameters presents significant challenges:

1. Non-convexity: The log marginal likelihood surface has multiple local optima corresponding to qualitatively different models.

2. Scale sensitivity: Hyperparameters span many orders of magnitude—length-scales from 0.001 to 1000, variances from $10^{-6}$ to $10^6$.

3. Ill-conditioning: The Hessian can have eigenvalues differing by factors of $10^{10}$, causing gradient descent to oscillate wildly.

4. Parameter interactions: Signal variance and noise variance trade off; multiple length-scales interact in complex ways.

5. Computational cost: Each function evaluation requires $O(n^3)$ operations, limiting how many optimization iterations are practical.

This page develops strategies to navigate these challenges reliably.

Several optimization algorithms are suitable for GP hyperparameter learning. Each has trade-offs between per-iteration cost, convergence speed, and robustness.

First-Order Methods (Gradient-Based):

These methods use only gradient information $\nabla_\theta \mathcal{L}$:

Quasi-Newton Methods:

These methods approximate second-order information from gradient history:

Second-Order Methods:

These use explicit Hessian information:

Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) is a quasi-Newton method that approximates the inverse Hessian using information from the last $m$ gradient evaluations (typically $m = 5$ to $20$).

Algorithm Overview:

1. Compute gradient: $\mathbf{g}k = -\nabla\theta \mathcal{L}(\theta_k)$ (negative because we maximize)

2. Compute search direction: $\mathbf{d}_k = \mathbf{H}_k \mathbf{g}_k$ where $\mathbf{H}_k$ approximates $\nabla^2 \mathcal{L}^{-1}$

3. Line search: Find step size $\alpha$ satisfying Wolfe conditions:

   • Armijo (sufficient decrease): $\mathcal{L}(\theta_k + \alpha \mathbf{d}_k) \geq \mathcal{L}(\theta_k) + c_1 \alpha \mathbf{g}_k^\top \mathbf{d}_k$
   • Curvature: $|\nabla\mathcal{L}(\theta_k + \alpha \mathbf{d}_k)^\top \mathbf{d}_k| \leq c_2 |\mathbf{g}_k^\top \mathbf{d}_k|$

4. Update: $\theta_{k+1} = \theta_k + \alpha \mathbf{d}_k$

5. Update inverse Hessian approximation: Store $(\mathbf{s}_k, \mathbf{y}_k)$ pairs where:

   • $\mathbf{s}k = \theta{k+1} - \theta_k$
   • $\mathbf{y}k = \nabla\mathcal{L}(\theta{k+1}) - \nabla\mathcal{L}(\theta_k)$

Why L-BFGS Excels for GPs:

1. Handles scale differences: The inverse Hessian approximation automatically adapts to different parameter scales.

2. Memory efficiency: Only stores $m$ vector pairs, regardless of problem dimensionality.

3. Superlinear convergence: Achieves better-than-linear convergence near the optimum.

4. Robust line search: The Wolfe conditions ensure progress even with approximate Hessians.

Key Parameters:

The marginal likelihood surface is non-convex. Different local optima can correspond to genuinely different explanations of the data:

Example: Two Valid Models

Consider noisy observations of a moderately varying function:

• Model A: Smooth function ($\ell = 1.0$) + substantial noise ($\sigma_n^2 = 0.5$)

  • Interpretation: 'The true function is smooth; variability is noise.'

• Model B: Wiggly function ($\ell = 0.1$) + minimal noise ($\sigma_n^2 = 0.01$)

  • Interpretation: 'The true function is complex; observations are accurate.'

Both can achieve similar marginal likelihood values, yet they make very different predictions for new points (especially in uncertainty estimates).

Strategy: Multiple Random Restarts

The most reliable approach to finding the global optimum:

1. Generate diverse initial points: Sample initial hyperparameters from a broad prior.

2. Run optimization from each: Apply L-BFGS (or chosen optimizer) independently.

3. Select best: Choose the result with highest marginal likelihood.

4. Check consistency: If runs converge to very different solutions with similar LML, investigate further—the data may genuinely support multiple explanations.

How Many Restarts?

Good initialization accelerates convergence and increases the probability of finding the global optimum.

Data-Informed Initialization:

Rather than random sampling, use data statistics to guide initialization:

Initialization for Random Restarts:

When using multiple restarts, diversify initial points:

1. One data-informed point: Use the strategy above for at least one restart.

2. Latin Hypercube Sampling: Generate well-spread points in hyperparameter space.

3. Log-uniform sampling: Sample $\log \theta \sim \text{Uniform}(\log \theta_{\min}, \log \theta_{\max})$ to cover multiple orders of magnitude.

4. Perturbations of data-informed: Add Gaussian noise to the data-informed initialization.

How we represent hyperparameters matters as much as how we optimize them. Poor parameterization causes ill-conditioning, constraint violations, and slow convergence.

The Natural (Constrained) Parameterization:

Hyperparameters have natural constraints:

• $\ell > 0$ (length-scale positive)
• $\sigma_f^2 > 0$ (signal variance positive)
• $\sigma_n^2 > 0$ (noise variance positive)

Optimizing directly with constrained optimization is possible but adds complexity.

The Log Parameterization (Unconstrained):

Map constrained parameters to unconstrained via: $$\bar{\theta} = \log \theta \quad \Rightarrow \quad \theta = \exp(\bar{\theta})$$

Advantages:

1. Unconstrained optimization: Use any optimizer without bounds.
2. Scale invariance: Parameters spanning orders of magnitude become comparable.
3. Numerical stability: Small positive values don't underflow.

Gradient Transformation:

With $\theta = \exp(\bar{\theta})$, the chain rule gives:

$$\frac{\partial \mathcal{L}}{\partial \bar{\theta}} = \frac{\partial \mathcal{L}}{\partial \theta} \cdot \frac{\partial \theta}{\partial \bar{\theta}} = \frac{\partial \mathcal{L}}{\partial \theta} \cdot \theta$$

This multiplication by $\theta$ naturally scales gradients—large parameters have proportionally larger gradients, appropriate for their scale.

Alternative Parameterizations:

Knowing when to stop optimization is crucial—stopping too early leaves performance on the table; stopping too late wastes computation.

Standard Convergence Criteria:

1. Gradient norm: $|\nabla\mathcal{L}| < \epsilon_g$

   • Typical: $\epsilon_g = 10^{-5}$ to $10^{-8}$
   • Indicates: First-order optimality condition nearly satisfied

2. Function change: $|\mathcal{L}^{(k)} - \mathcal{L}^{(k-1)}| < \epsilon_f$

   • Typical: $\epsilon_f = 10^{-9}$ to $10^{-12}$
   • Indicates: Marginal likelihood barely improving

3. Parameter change: $|\theta^{(k)} - \theta^{(k-1)}| < \epsilon_x$

   • Typical: $\epsilon_x = 10^{-6}$ to $10^{-9}$
   • Indicates: Hyperparameters stabilized

4. Maximum iterations: Fallback to prevent infinite loops

Practical Considerations:

Recommended Composite Criterion:

Stop when ANY of:

• Gradient norm $< 10^{-5}$ AND function change $< 10^{-9}$ for 3 consecutive iterations
• Function change $< 10^{-12}$ for 5 consecutive iterations
• Maximum iterations reached

The 'consecutive iterations' requirement prevents premature stopping in flat regions.

The marginal likelihood alone can sometimes drive hyperparameters to pathological values. Adding bounds or regularization prevents such behavior.

Hyperparameter Bounds:

Jitter for Numerical Stability:

Always add a small 'jitter' term to the noise variance: $$\mathbf{K}_y = \mathbf{K} + (\sigma_n^2 + \epsilon) \mathbf{I}$$

with $\epsilon \approx 10^{-6}$ to $10^{-8}$. This ensures $\mathbf{K}_y$ is positive definite even if $\sigma_n^2$ becomes very small.

Hyperprior Regularization:

An alternative to hard bounds is adding a log-prior to the objective:

$$\tilde{\mathcal{L}}(\theta) = \mathcal{L}(\theta) + \log p(\theta)$$

This corresponds to MAP estimation of hyperparameters rather than maximum likelihood.

Common Hyperpriors:

The gradient of the log-prior simply adds to the marginal likelihood gradient.

Each marginal likelihood evaluation costs $O(n^3)$ for $n$ training points. For large datasets, this dominates total computation.

Budget Allocation:

Strategies for Large Datasets:

1. Subsampling for initialization: Optimize on a random subset (e.g., 1000 points) to find good starting hyperparameters, then fine-tune on full data.

2. Warm starting: If training multiple related models, use previous optimum as starting point.

3. Early stopping in restarts: If a restart's LML is below the current best after a few iterations, abort early.

4. Parallel restarts: Run multiple restarts simultaneously on different cores.

5. Sparse approximations: Use inducing point methods (covered in Module 5) where the cost is $O(nm^2)$ with $m \ll n$ inducing points.

Approximating Gradients:

For very large $n$, consider:

• Stochastic gradients using minibatches (requires careful variance reduction)
• Matrix-free methods that avoid forming full $\mathbf{K}$
• Kronecker or Toeplitz structure exploitation for gridded data

When optimization doesn't converge or produces unreasonable hyperparameters, systematic debugging is essential.

Common Failure Modes and Solutions:

Debugging Checklist:

1. Verify gradients numerically (as covered in Page 1)
2. Check data for issues: Duplicates, huge outliers, constant features
3. Visualize LML surface: Plot 1D or 2D slices to understand landscape
4. Examine intermediate values: Log $\mathbf{K}_y$ condition number, gradient norms
5. Test on synthetic data: Generate data from known GP, verify recovery
6. Compare with known implementation: GPyTorch, GPflow, scikit-learn

We've developed a comprehensive toolkit for optimizing GP hyperparameters. Let's consolidate the key strategies:

Looking Ahead:

The next page examines cross-validation as an alternative to marginal likelihood for hyperparameter selection. Understanding both approaches allows you to choose the right tool for each situation and use cross-validation as a sanity check when marginal likelihood behavior is unexpected.