Gaussian Processes for Machine Learning

Authors: Carl Edward Rasmussen, Christopher K. I. Williams
Published: 2006
Publisher: MIT Press
Link: Book Website

Summary

This comprehensive textbook establishes Gaussian Processes (GPs) as a principled Bayesian approach to nonparametric regression and classification. It provides the connection between GPs and kernel methods through the covariance function, unifying probabilistic and RKHS perspectives.

Key Contributions

1. GP-Kernel Connection

A Gaussian Process is a distribution over functions: $f \sim \mathcal{GP}(m(\cdot), k(\cdot, \cdot))$

where:

• $m(x) = \mathbb{E}[f(x)]$ is the mean function
• $k(x, x') = \text{Cov}(f(x), f(x'))$ is the covariance (kernel) function

2. Posterior Distribution

Given observations $(X, y)$, the posterior at test point $x_*$ is: $f_* | X, y, x_* \sim \mathcal{N}(\bar{f}_*, \text{var}(f_*))$

where: $\bar{f}_* = k(x_*, X)(K + \sigma^2 I)^{-1}y$ $\text{var}(f_*) = k(x_*, x_*) - k(x_*, X)(K + \sigma^2 I)^{-1}k(X, x_*)$

3. Marginal Likelihood

The log marginal likelihood for hyperparameter selection: $\log p(y|X, \theta) = -\frac{1}{2}y^T(K + \sigma^2 I)^{-1}y - \frac{1}{2}\log|K + \sigma^2 I| - \frac{n}{2}\log 2\pi$

GP-RKHS Connection

Equivalence

The GP posterior mean equals the RKHS solution: $\bar{f} = \arg\min_{f \in \mathcal{H}_k} \sum_{i=1}^n (f(x_i) - y_i)^2 + \lambda \|f\|_{\mathcal{H}}^2$

with $\lambda = \sigma^2$.

Different Perspectives

Aspect	GP View	RKHS View
Function	Random variable	Fixed unknown
Prior	$\mathcal{GP}(0, k)$	$\|f\|_{\mathcal{H}} < \infty$
Inference	Bayesian posterior	Regularized ERM
Uncertainty	Full posterior	Point estimate

Common Kernels

Squared Exponential (RBF)

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

Properties: Infinitely differentiable, very smooth functions.

Matérn Class

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

Special cases:

• $\nu = 1/2$: Exponential (Ornstein-Uhlenbeck)
• $\nu = 3/2$: Once differentiable
• $\nu = 5/2$: Twice differentiable
• $\nu \to \infty$: Squared exponential

Periodic Kernel

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

For periodic functions with period $p$.

Implementation

import numpy as np
from scipy.linalg import cholesky, cho_solve

class GaussianProcess:
    def __init__(self, kernel, noise_var=1e-6):
        self.kernel = kernel
        self.noise_var = noise_var
        
    def fit(self, X, y):
        self.X_train = X
        self.y_train = y
        
        # Compute kernel matrix
        K = self.kernel(X, X)
        K += self.noise_var * np.eye(len(X))
        
        # Cholesky decomposition for numerical stability
        self.L = cholesky(K, lower=True)
        self.alpha = cho_solve((self.L, True), y)
        
    def predict(self, X_test, return_var=False):
        # Cross-covariance
        K_star = self.kernel(X_test, self.X_train)
        
        # Posterior mean
        mean = K_star @ self.alpha
        
        if return_var:
            # Posterior variance
            v = cho_solve((self.L, True), K_star.T)  # shape: (n_train, n_test)
            K_ss = self.kernel(X_test, X_test)
            var = np.diag(K_ss) - np.sum(K_star.T * v, axis=0)
            return mean, var
        
        return mean
    
    def log_marginal_likelihood(self):
        """Compute log marginal likelihood for hyperparameter optimization."""
        n = len(self.y_train)
        
        # Data fit term
        data_fit = -0.5 * self.y_train.T @ self.alpha
        
        # Complexity penalty (log determinant via Cholesky)
        complexity = -np.sum(np.log(np.diag(self.L)))
        
        # Normalization
        norm = -0.5 * n * np.log(2 * np.pi)
        
        return data_fit + complexity + norm

def rbf_kernel(X1, X2, lengthscale=1.0, variance=1.0):
    """Squared exponential (RBF) kernel."""
    sq_dist = np.sum(X1**2, axis=1, keepdims=True) + \
              np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
    return variance * np.exp(-0.5 * sq_dist / lengthscale**2)

Hyperparameter Optimization

Gradient-Based

Maximize log marginal likelihood: $\frac{\partial \log p(y|X,\theta)}{\partial \theta_j} = \frac{1}{2}\text{tr}\left((\alpha\alpha^T - K^{-1})\frac{\partial K}{\partial \theta_j}\right)$

Automatic Relevance Determination (ARD)

Use separate lengthscale per dimension: $k(x, x') = \sigma_f^2 \exp\left(-\sum_d \frac{(x_d - x'_d)^2}{2\ell_d^2}\right)$

Small $\ell_d$ indicates relevant dimension.

Scalability

Computational Challenges

• Training: $O(n^3)$ for matrix inversion
• Prediction: $O(n^2)$ for new point
• Memory: $O(n^2)$ for kernel matrix

Sparse GP Methods

1. Subset of data: Use $m \ll n$ inducing points
2. FITC/PITC: Structured approximations
3. Variational: Inducing point variational inference

Applications

1. Bayesian optimization: Surrogate modeling with uncertainty
2. Spatial statistics: Kriging for geostatistics
3. Time series: Forecasting with uncertainty
4. Active learning: Acquisition functions using variance
5. Emulation: Surrogate for expensive simulations

Citation

@book{rasmussen2006gaussian,
  title={Gaussian Processes for Machine Learning},
  author={Rasmussen, Carl Edward and Williams, Christopher KI},
  year={2006},
  publisher={MIT Press}
}

Further Reading

• Williams, C. K. I. (1998). Computation with infinite neural networks
• Titsias, M. (2009). Variational learning of inducing variables in sparse GPs
• Hensman, J., et al. (2013). Gaussian processes for big data