Gaussian Processes: Non-Parametric Bayesian Regression

Processus Gaussiens : Guide Complet — Régression Bayésienne Non-Paramétrique

Gaussian Process: Complete Guide — Non-Parametric Bayesian Regression

Summary — Gaussian Processes (GP) are a non-parametric Bayesian approach to regression that directly models a distribution over functions. Unlike parametric models (such as linear regression) that learn a fixed set of parameters, a GP learns a complete distribution of possible functions compatible with the observed data. The result: predictions with quantifiable uncertainty bars, which is crucial for sensitive applications like industrial control, medicine, or decision-making under uncertainty.


Mathematical Principle

1. Definition of a Gaussian Process

A Gaussian Process is a distribution over functions f(x) completely specified by its mean function m(x) and its covariance function (kernel) k(x,x’):

f(x) ~ GP(m(x), k(x,x'))

For any finite set of points X = {x₁, …, xₙ}, the values f(X) follow a multivariate Gaussian distribution:

f(X) ~ N(m(X), K(X,X))

Where K is the Gram matrix (or covariance matrix) with K_{ij} = k(x_i, x_j).

2. Bayesian prediction

Given observations y = f(X) + ε where ε ~ N(0, sigma_n² · I), the predictive distribution at new points X* is:

f* | X, y, X* ~ N(μ*, Σ*)

μ* = k(X*, X) · [K(X, X) + sigma_n² · I]^{-1} · y

Σ* = k(X*, X*) - k(X*, X) · [K(X, X) + sigma_n² · I]^{-1} · k(X, X*)

The key point: Σ gives epistemic uncertainty* — how uncertain the model is due to lack of data in that region of the input space.

3. Covariance functions (kernels)

RBF (Radial Basis Function / Squared Exponential):

k(x, x') = sigma² · exp(-||x - x'||² / (2 · l²))

This is the most commonly used kernel. l (correlation length) controls smoothness: small l = rapid variations, large l = smooth variations.

Matérn:

k_Matern(x, x') = sigma² · 2^(1-nu) / Gamma(nu) · (sqrt(2·nu) · r / l)^nu · K_nu(sqrt(2·nu) · r / l)

Where r = ||x – x’||, K_nu is the modified Bessel function, and nu controls regularity. nu = 3/2 or 5/2 are common choices (smoother than RBF for nu > 0.5).

Periodic:

k_periodic(x, x') = sigma² · exp(-2 · sin²(pi · ||x - x'|| / p) / l²)

For data with periodic patterns (seasonality, cycles).

4. Hyperparameter optimization

Kernel hyperparameters (l, sigma, etc.) are optimized by maximizing the log marginal likelihood:

log p(y | X, theta) = -1/2 · y^T · K_y^{-1} · y - 1/2 · log|K_y| - n/2 · log(2·pi)

Where K_y = K(X, X) + sigma_n² · I. The first term penalizes poor fit, the second controls model complexity (automatic regularization known as Occam’s razor).


Intuition

Imagine a draftsman who is shown a few points on a graph and asked to draw the curve connecting them.

The draftsman doesn’t draw a single curve — he imagines a bundle of thousands of curves that all pass near the observed points. Where the points are numerous and close together, all the curves look alike (high confidence). Where there are no points, the curves diverge in all directions (high uncertainty).

A Gaussian Process is exactly that: it doesn’t predict a single value, it predicts a distribution of possible values with their probability. It’s the only machine learning model that honestly says “I don’t know” — and that’s its greatest strength.


Python Implementation

1. GP Regression with scikit-learn

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel
from sklearn.gaussian_process.kernels import Matern
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Noisy 1D synthetic data
np.random.seed(42)
X = np.array([[0.0], [0.3], [0.6], [0.8], [1.0]])
y = np.sin(3 * X.ravel() + 0.5) + 0.1 * np.random.randn(5)

# Kernel: Constant * RBF + White (noise)
kernel = ConstantKernel(1.0) * RBF(length_scale=0.5) + WhiteKernel(noise_level=0.1)

# Create and train the GP
gp = GaussianProcessRegressor(
    kernel=kernel,
    alpha=0.0,  # noise already in kernel
    n_restarts_optimizer=10,  # avoid local minima
    normalize_y=True
)
gp.fit(X, y)

print(f"Optimized kernel: {gp.kernel_}")
print(f"Log-marginal-likelihood: {gp.log_marginal_likelihood(gp.kernel_.theta):.3f}")

# Predictions with confidence interval
X_plot = np.linspace(-0.3, 1.5, 300).reshape(-1, 1)
y_pred, sigma = gp.predict(X_plot, return_std=True)

# Visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Prediction + uncertainty
ax1.fill_between(X_plot.ravel(),
                 y_pred - 1.96 * sigma,
                 y_pred + 1.96 * sigma,
                 alpha=0.2, label='95% confidence interval')
ax1.plot(X_plot, y_pred, 'b-', label='GP prediction')
ax1.plot(X, y, 'ro', markersize=8, label='Observations')
ax1.fill_between(X_plot.ravel(),
                 np.sin(3*X_plot.ravel() + 0.5) - 0.2,
                 np.sin(3*X_plot.ravel() + 0.5) + 0.2,
                 alpha=0.1, color='green', label='True function +/- noise')
ax1.set_title('Gaussian Process: Prediction + Uncertainty')
ax1.legend()
ax1.set_xlabel('x')
ax1.set_ylabel('y')

# Posterior process samples
y_samples = gp.sample_y(X_plot, n_samples=8)
ax2.plot(X_plot, y_samples, alpha=0.6)
ax2.plot(X, y, 'ro', markersize=8, label='Observations')
ax2.set_title('8 posterior process samples')
ax2.legend()

plt.tight_layout()
plt.savefig('gaussian_process.png', dpi=150)
print('Visualization saved!')

# Comparison with different kernels
kernels_to_test = [
    RBF(length_scale=1.0),
    Matern(length_scale=1.0, nu=1.5),
    Matern(length_scale=1.0, nu=2.5),
    ConstantKernel(1.0) * RBF(length_scale=1.0),
]

for i, k in enumerate(kernels_to_test):
    gp_test = GaussianProcessRegressor(kernel=k, n_restarts_optimizer=15)
    gp_test.fit(X, y)
    print(f'Kernel {i+1}: {gp_test.kernel_} | LML: {gp_test.log_marginal_likelihood():.3f}')

Hyperparameters

Hyperparameter Typical value Description
kernel (RBF/Matern) RBF(1.0) or Matern(1.0, nu=2.5) Covariance function
alpha 1e-10 to 1.0 Measurement noise variance (if not in kernel)
n_restarts_optimizer 10-50 Kernel optimization restarts
normalize_y True/False Target normalization (recommended=True)

Advantages

  1. Quantified uncertainty: The GP provides exact confidence intervals, not approximations. Essential for critical decision-making.
  2. No overfitting: Marginal likelihood maximization acts as automatic regularization. The GP chooses the complexity level adapted to the data.
  3. Little data needed: Unlike deep networks, GPs work well with as few as 10-50 points. Ideal for expensive experiments.
  4. No iterative training: The “fit” is an analytical optimization of kernel hyperparameters. No learning rate, no epochs, no risk of non-convergence.

Limitations

  1. Cubic complexity O(n³): Inverting the K matrix (of size n×n) costs O(n³). Beyond ~10,000 points, it’s impractical without approximations (Sparse GPs, inducing points).
  2. Kernel choice: The result depends heavily on the chosen kernel. A bad kernel gives bad predictions. Practitioner expertise is necessary.
  3. Gaussian assumption: The GP assumes data follows a Gaussian distribution. For strongly skewed or multi-modal data, it’s inadequate.
  4. Non-stationarity: Standard kernels assume covariance depends only on the distance between points. If function regularity varies by region, non-stationary kernels are needed (more complex).

4 concrete use cases

1. Bayesian hyperparameter optimization

The GP is the surrogate model of choice in Bayesian optimization (see article 133). It models the unknown cost function f(learning_rate, batch_size, …) and guides the search toward promising regions.

2. Industrial sensor calibration

When a sensor is expensive to calibrate (each measurement takes hours), the GP interpolates between measured calibration points and provides uncertainty for each prediction. The engineer then knows when to recalibrate.

3. Energy consumption forecasting

With seasonal data (day/night, weekly cycles), a GP with periodic + RBF kernel captures consumption patterns and quantifies uncertainty, enabling appropriate sizing of energy reserves.

4. Pharmacokinetic modeling

In the pharmaceutical industry, the GP models the concentration of a drug in the blood as a function of time and dose, with rigorous confidence intervals for clinical trials.


Conclusion

Gaussian Processes represent one of the most elegant tools in machine learning: a rigorously Bayesian approach that provides predictions with uncertainty guarantees. While deep learning dominates large-scale applications, GPs remain unbeatable for small-scale problems where uncertainty matters as much as the prediction itself.

Modern approximations (Sparse GPs with inducing points, Deep Kernel Learning combining GPs and deep networks) extend the applicability of GPs to ever larger scales, preserving their fundamental advantage: quantifying uncertainty.


See also


About Salah YAHIAOUI

Salah YAHIAOUI publie sur The Pingouin des tutoriels Python et des décryptages sur l’intelligence artificielle. Les cas pratiques récents précisent leurs données, leurs tests et leurs limites. À propos, méthode éditoriale et contact.

View all posts by Salah YAHIAOUI →