Gaussian Processes for Machine Learning
Gaussian processes (GPs) provide a principled, nonparametric Bayesian approach to regression and classification. They model distributions over functions, enabling uncertainty quantification alongside predictions.
1. Gaussian Process Definition
1.1 Formal Definition
Definition: A Gaussian process is a collection of random variables such that any finite subset has a multivariate Gaussian distribution:
A GP is completely specified by its mean function and covariance kernel .
1.2 Mean and Kernel Functions
The mean function captures prior beliefs about the average function value:
The kernel function encodes assumptions about function smoothness, periodicity, and other structural properties:
1.3 Finite-Dimensional Consistency
A GP is well-defined if every finite-dimensional marginal is consistent:
This is guaranteed by the Kolmogorov extension theorem when the kernel is positive definite.
2. Kernel Functions
2.1 RBF (Squared Exponential) Kernel
- : signal variance (controls amplitude)
- : length scale (controls smoothness)
- Infinitely differentiable functions
- Spectral density:
2.2 Matérn Kernels
Special cases:
- : exponential kernel (rough, non-differentiable)
- : (once differentiable)
- : (twice differentiable)
- : converges to RBF kernel
where and is the modified Bessel function.
2.3 Periodic Kernel
- : period
- Captures exactly periodic patterns
- Can be combined with RBF for locally periodic behavior
2.4 Matérn-ARMA Connection
The Matérn class arises as the solution to the stochastic differential equation:
where is white noise. This connects GPs to Gaussian Markov random fields and ARMA processes.
3. GP Regression
3.1 Model
Observation model: where .
Joint distribution of training outputs and test outputs :
where , , .
3.2 Predictive Distribution
Conditioning on training data:
where:
3.3 Numerical Stability
Direct computation of is and numerically unstable. Use Cholesky decomposition:
where is lower triangular. Then:
- Solve for
- Solve for
Complexity: for decomposition, per prediction.
4. Marginal Likelihood
4.1 Definition
The marginal likelihood (evidence) is:
where represents kernel hyperparameters.
4.2 Analytical Form
For Gaussian likelihood:
Three terms:
- Data fit: measures how well the model explains the data
- Complexity penalty: penalizes complex models
- Normalization: constant
4.3 Automatic Relevance Determination (ARD)
For input-dependent length scales:
Each input dimension has its own length scale . Large indicates dimension is irrelevant.
4.4 Hyperparameter Optimization
Maximize log marginal likelihood w.r.t. :
Gradient computation (using Cholesky):
5. GP Classification
5.1 Latent GP Model
For binary classification with labels :
where is the logistic function (or probit, etc.).
The posterior is non-Gaussian, requiring approximation.
5.2 Laplace Approximation
Find the MAP estimate :
Linearize around :
where and .
5.3 Expectation Propagation
EP iteratively refines Gaussian approximations by matching moments:
- Cavity: remove factor from approximation
- Tilt: combine cavity with true factor
- Projection: project back to Gaussian
Converges to better approximation than Laplace but more expensive.
6. Sparse GP Approximations
6.1 Inducing Point Methods
Introduce inducing points . The variational free energy:
where and .
6.2 FITC (Fully Independent Training Conditional)
Approximate the joint GP by assuming independence between test points:
Complexity: instead of .
6.3 Variational Free Energy
The variational lower bound:
Optimize w.r.t. inducing points and variational parameters.
7. Multi-Output GPs
7.1 Linear Model of Coregionalization (LMC)
where are independent GPs. The kernel:
7.2 Intrinsic Coregionalization Model (ICM)
Special case where for all :
where is the coregionalization matrix.
8. Deep Gaussian Processes
8.1 Composition
A deep GP is a composition of GPs:
where the input to is the output of .
8.2 Variational Inference
Use inducing points at each layer and variational inference:
8.3 Expressiveness
Deep GPs can represent:
- Non-stationary functions
- Functions with varying length scales
- Hierarchical feature representations
9. Code: Gaussian Process Implementation
import numpy as np
from scipy.optimize import minimize
from scipy.linalg import cho_solve, cho_factor
from scipy.special import kv, gamma as gamma_fn
class GaussianProcessRegressor:
"""Full Gaussian Process Regression with hyperparameter optimization."""
def __init__(self, kernel='rbf', noise=1e-2, length_scale=1.0,
signal_variance=1.0, optimize_hyperparameters=True):
self.kernel_type = kernel
self.noise = noise
self.length_scale = length_scale
self.signal_variance = signal_variance
self.optimize_hyperparameters = optimize_hyperparameters
self.X_train_ = None
self.y_train_ = None
self.L_ = None
self.alpha_ = None
self.log_marg_lik_ = None
def _rbf_kernel(self, X1, X2, length_scale, signal_var):
"""RBF kernel."""
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return signal_var * np.exp(-0.5 * sq_dists / length_scale**2)
def _matern32_kernel(self, X1, X2, length_scale, signal_var):
"""Matérn 3/2 kernel."""
r = np.sqrt(np.sum((X1[:, None] - X2[None, :])**2, axis=2))
r = np.maximum(r, 1e-10) # Numerical stability
return signal_var * (1 + np.sqrt(3) * r / length_scale) * \
np.exp(-np.sqrt(3) * r / length_scale)
def _matern52_kernel(self, X1, X2, length_scale, signal_var):
"""Matérn 5/2 kernel."""
r = np.sqrt(np.sum((X1[:, None] - X2[None, :])**2, axis=2))
r = np.maximum(r, 1e-10)
return signal_var * (1 + np.sqrt(5) * r / length_scale +
5 * r**2 / (3 * length_scale**2)) * \
np.exp(-np.sqrt(5) * r / length_scale)
def _periodic_kernel(self, X1, X2, length_scale, signal_var, period=1.0):
"""Periodic kernel."""
r = np.abs(X1[:, None] - X2[None, :])
return signal_var * np.exp(-2 * np.sin(np.pi * r / period)**2 / length_scale**2)
def _compute_kernel(self, X1, X2, params=None):
"""Compute kernel matrix with given parameters."""
if params is None:
length_scale = self.length_scale
signal_var = self.signal_variance
else:
length_scale = np.exp(params[0]) # Log-transformed
signal_var = np.exp(params[1])
if self.kernel_type == 'rbf':
return self._rbf_kernel(X1, X2, length_scale, signal_var)
elif self.kernel_type == 'matern32':
return self._matern32_kernel(X1, X2, length_scale, signal_var)
elif self.kernel_type == 'matern52':
return self._matern52_kernel(X1, X2, length_scale, signal_var)
elif self.kernel_type == 'periodic':
return self._periodic_kernel(X1, X2, length_scale, signal_var)
else:
raise ValueError(f"Unknown kernel: {self.kernel_type}")
def _compute_cholesky(self, K):
"""Compute Cholesky decomposition with jitter."""
K_jitter = K + 1e-6 * np.eye(len(K))
try:
L = np.linalg.cholesky(K_jitter)
return L, True
except np.linalg.LinAlgError:
# Add more jitter
for jitter in [1e-5, 1e-4, 1e-3]:
try:
L = np.linalg.cholesky(K + jitter * np.eye(len(K)))
return L, True
except np.linalg.LinAlgError:
continue
return None, False
def _negative_log_marginal_likelihood(self, params):
"""Compute negative log marginal likelihood for optimization."""
try:
K = self._compute_kernel(self.X_train_, self.X_train_, params)
noise = np.exp(params[2])
K += noise * np.eye(len(self.X_train_))
L, success = self._compute_cholesky(K)
if not success:
return 1e10
# Solve for alpha
alpha = cho_solve((L, True), self.y_train_)
# Log marginal likelihood
lml = -0.5 * self.y_train_ @ alpha - \
np.sum(np.log(np.diag(L))) - \
len(self.y_train_) / 2 * np.log(2 * np.pi)
return -lml # Negative for minimization
except Exception as e:
return 1e10
def fit(self, X, y):
"""Fit GP regression model."""
self.X_train_ = X
self.y_train_ = y
if self.optimize_hyperparameters:
# Initial parameters: [log_length_scale, log_signal_var, log_noise]
x0 = np.log([self.length_scale, self.signal_variance, self.noise])
result = minimize(self._negative_log_marginal_likelihood, x0,
method='L-BFGS-B',
options={'maxiter': 100})
# Update parameters
self.length_scale = np.exp(result.x[0])
self.signal_variance = np.exp(result.x[1])
self.noise = np.exp(result.x[2])
self.log_marg_lik_ = -result.fun
# Compute final Cholesky
K = self._compute_kernel(self.X_train_, self.X_train_)
K += self.noise * np.eye(len(self.X_train_))
self.L_, success = self._compute_cholesky(K)
if success:
self.alpha_ = cho_solve((self.L_, True), self.y_train_)
else:
raise RuntimeError("Cholesky decomposition failed")
return self
def predict(self, X, return_std=True, return_cov=False):
"""Predict at new points."""
K_s = self._compute_kernel(X, self.X_train_)
K_ss = self._compute_kernel(X, X)
# Mean
mu = K_s @ self.alpha_
if return_cov:
# Full covariance
v = np.linalg.solve(self.L_, K_s.T)
cov = K_ss - v.T @ v
return mu, cov
if return_std:
# Variance
v = np.linalg.solve(self.L_, K_s.T)
var = np.diag(K_ss) - np.sum(v**2, axis=0)
var = np.maximum(var, 0)
return mu, np.sqrt(var)
return mu
def sample_prior(self, X, n_samples=5):
"""Sample from the prior GP."""
K = self._compute_kernel(X, X)
L = np.linalg.cholesky(K + 1e-8 * np.eye(len(X)))
samples = []
for _ in range(n_samples):
f = L @ np.random.randn(len(X))
samples.append(f)
return np.array(samples)
def sample_posterior(self, X, n_samples=5):
"""Sample from the posterior GP."""
K_s = self._compute_kernel(X, self.X_train_)
K_ss = self._compute_kernel(X, X)
# Posterior covariance
v = np.linalg.solve(self.L_, K_s.T)
cov = K_ss - v.T @ v
cov = (cov + cov.T) / 2 # Ensure symmetry
# Add small jitter for numerical stability
L_cov = np.linalg.cholesky(cov + 1e-8 * np.eye(len(X)))
# Mean
mu = K_s @ self.alpha_
samples = []
for _ in range(n_samples):
f = mu + L_cov @ np.random.randn(len(X))
samples.append(f)
return np.array(samples)
class GaussianProcessClassifier:
"""Gaussian Process Classification using Laplace approximation."""
def __init__(self, kernel='rbf', length_scale=1.0, signal_variance=1.0):
self.kernel_type = kernel
self.length_scale = length_scale
self.signal_variance = signal_variance
self.X_train_ = None
self.y_train_ = None
self.f_star_ = None
def _compute_kernel(self, X1, X2):
"""RBF kernel."""
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return self.signal_variance * np.exp(-0.5 * sq_dists / self.length_scale**2)
def _logistic(self, f):
"""Logistic sigmoid function."""
return 1 / (1 + np.exp(-f))
def _newton_updates(self, K, y, max_iter=20):
"""Find MAP estimate using Newton's method."""
n = len(y)
f = np.zeros(n)
for _ in range(max_iter):
# Probabilities
pi = self._logistic(f)
# Gradient
grad = y - pi
# Hessian (negative definite)
W = np.diag(pi * (1 - pi))
B = np.eye(n) + W @ K
# Solve for Newton direction
try:
L = np.linalg.cholesky(B)
b = W @ f + grad
v = np.linalg.solve(L, b)
delta = np.linalg.solve(L.T, v)
delta = f - K @ delta
f = delta
except np.linalg.LinAlgError:
break
return f
def fit(self, X, y):
"""Fit GP classifier."""
self.X_train_ = X
self.y_train_ = y
K = self._compute_kernel(X, X)
# Find MAP estimate
self.f_star_ = self._newton_updates(K, y)
return self
def predict_proba(self, X):
"""Predict probabilities."""
K_s = self._compute_kernel(X, self.X_train_)
K_ss = self._compute_kernel(X, X)
# Linear approximation around f_star_
pi = self._logistic(self.f_star_)
W = np.diag(pi * (1 - pi))
B = np.eye(len(self.X_train_)) + W @ self._compute_kernel(self.X_train_, self.X_train_)
# Posterior mean
mu = K_s @ np.linalg.solve(B, pi)
# Posterior variance
v = np.linalg.solve(B, K_s.T)
var = np.diag(K_ss) - np.sum(v * (W @ v), axis=0)
# Probit approximation for classification
prob = self._logistic(mu / np.sqrt(1 + np.pi * var / 8))
return prob
def predict(self, X):
"""Predict class labels."""
prob = self.predict_proba(X)
return (prob > 0.5).astype(int)
class SparseGaussianProcess:
"""Sparse GP using inducing points (Variational Free Energy)."""
def __init__(self, n_inducing=50, kernel='rbf', length_scale=1.0,
signal_variance=1.0, noise=1e-2):
self.n_inducing = n_inducing
self.kernel_type = kernel
self.length_scale = length_scale
self.signal_variance = signal_variance
self.noise = noise
self.Z_ = None # Inducing points
self.m_ = None # Variational mean
self.S_ = None # Variational covariance
def _compute_kernel(self, X1, X2):
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return self.signal_variance * np.exp(-0.5 * sq_dists / self.length_scale**2)
def _variational_free_energy(self, Z, m, S, X, y):
"""Compute variational free energy."""
n = len(y)
m_ind = len(Z)
# Kernel matrices
K_mm = self._compute_kernel(Z, Z)
K_mn = self._compute_kernel(Z, X)
K_nn = self._compute_kernel(X, X)
# Cholesky of K_mm
L = np.linalg.cholesky(K_mm + 1e-6 * np.eye(m_ind))
# Trace term
trace_term = np.trace(np.linalg.solve(K_mm, S))
# KL term
L_S = np.linalg.cholesky(S + 1e-6 * np.eye(m_ind))
kl = 0.5 * (np.trace(np.linalg.solve(K_mm, S)) +
m @ np.linalg.solve(K_mm, m) - m_ind +
2 * np.sum(np.log(np.diag(L))) -
2 * np.sum(np.log(np.diag(L_S))))
# Data fit (approximate)
v = np.linalg.solve(L, K_mn)
Q_ff = v.T @ v
# Variational expectation
K_nn_diag = np.diag(K_nn)
trace_term_2 = np.sum(K_nn_diag - np.diag(Q_ff))
# Expected log likelihood
mu = K_mn.T @ np.linalg.solve(K_mm + S, m)
var = K_nn_diag - np.diag(K_mn.T @ np.linalg.solve(K_mm, K_mn)) + \
np.diag(K_mn.T @ np.linalg.solve(K_mm, S @ np.linalg.solve(K_mm, K_mn)))
# Gaussian approximation to log likelihood
elbo = -0.5 * n * np.log(2 * np.pi * self.noise) - \
0.5 * np.sum((y - mu)**2 + var) / self.noise - \
kl
return -elbo # Negative for minimization
def fit(self, X, y, n_optim=100):
"""Fit sparse GP."""
n, d = X.shape
# Initialize inducing points
idx = np.random.choice(n, self.n_inducing, replace=False)
self.Z_ = X[idx].copy()
# Initialize variational parameters
self.m_ = np.zeros(self.n_inducing)
self.S_ = np.eye(self.n_inducing)
# Optimize
def objective(params):
Z = params[:self.n_inducing * d].reshape(self.n_inducing, d)
m = params[self.n_inducing * d:self.n_inducing * (d + 1)]
S_flat = params[self.n_inducing * (d + 1):]
S = S_flat.reshape(self.n_inducing, self.n_inducing)
S = (S + S.T) / 2 # Ensure symmetry
return self._variational_free_energy(Z, m, S, X, y)
# Pack parameters
x0 = np.concatenate([self.Z_.flatten(), self.m_, self.S_.flatten()])
# Optimize (simplified - in practice use more sophisticated methods)
from scipy.optimize import minimize
result = minimize(objective, x0, method='L-BFGS-B',
options={'maxiter': n_optim})
# Unpack results
params = result.x
self.Z_ = params[:self.n_inducing * d].reshape(self.n_inducing, d)
self.m_ = params[self.n_inducing * d:self.n_inducing * (d + 1)]
S_flat = params[self.n_inducing * (d + 1):]
self.S_ = S_flat.reshape(self.n_inducing, self.n_inducing)
return self
def predict(self, X):
"""Predict at new points."""
K_mm = self._compute_kernel(self.Z_, self.Z_)
K_mn = self._compute_kernel(self.Z_, X)
K_nn = self._compute_kernel(X, X)
# Mean
mu = K_mn.T @ np.linalg.solve(K_mm + self.S_, self.m_)
# Variance
v = np.linalg.solve(K_mm, K_mn)
var = np.diag(K_nn) - np.trace(v.T @ np.linalg.solve(K_mm, K_mn)) + \
np.diag(K_mn.T @ np.linalg.solve(K_mm, self.S_ @ np.linalg.solve(K_mm, K_mn)))
var = np.maximum(var, 0)
return mu, np.sqrt(var)
def plot_gp_prediction(gp, X_train, y_train, X_test, title="GP Prediction"):
"""Plot GP prediction with uncertainty bands."""
import matplotlib.pyplot as plt
mu, std = gp.predict(X_test, return_std=True)
plt.figure(figsize=(10, 6))
plt.scatter(X_train, y_train, c='red', s=20, label='Training data')
plt.plot(X_test, mu, 'b-', label='Mean')
plt.fill_between(X_test.flatten(), mu - 2*std, mu + 2*std,
alpha=0.3, color='blue', label='95% confidence')
plt.fill_between(X_test.flatten(), mu - std, mu + std,
alpha=0.5, color='blue', label='68% confidence')
plt.xlabel('x')
plt.ylabel('y')
plt.title(title)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
10. Summary
Gaussian processes provide a powerful Bayesian nonparametric framework:
- Full predictive distributions with uncertainty quantification
- Automatic Occam's razor via the marginal likelihood
- Flexible modeling through diverse kernel functions
- Theoretical guarantees on generalization
- Scalable approximations for large datasets (sparse GPs, variational methods)
The key insight is that GPs model distributions over functions, enabling principled uncertainty quantification while maintaining the flexibility of nonparametric methods.