Gaussian Processes (GPs)

A Gaussian process is a distribution over functions. Where a Gaussian distribution is specified by a mean and a variance for a number, a GP is specified by a mean function and a covariance function for a whole curve:

\[ f(x) \sim \mathcal{GP}\big(m(x),\; k(x, x’)\big) \]

The defining property is that any finite set of points drawn from it is jointly Gaussian, so conditioning on observed data is ordinary multivariate Gaussian arithmetic. Instead of fitting parameters and returning one curve, a GP returns a mean prediction and a variance at every point.

That variance is the reason to use one, and it behaves in a specific way.

The uncertainty tracks where the data is

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

rng = np.random.default_rng(0)
train_x = np.concatenate([rng.uniform(0, 3, 20), rng.uniform(7, 10, 20)]).reshape(-1, 1)
train_y = np.sin(train_x).ravel() + rng.normal(scale=0.1, size=len(train_x))
"""Observations cluster in [0, 3] and [7, 10]. Nothing was measured in between.
   The kernel is fixed rather than fitted, so the numbers below are reproducible."""

gp = GaussianProcessRegressor(kernel=1.0 * RBF(length_scale=1.0) + WhiteKernel(0.01),
                              optimizer=None).fit(train_x, train_y)

print(f"{'x':>6}  {'prediction':>11} {'std dev':>9}  {'nearest observation':>21}")
for x in (1.5, 3.0, 4.5, 5.0, 5.5, 7.0, 9.0, 13.0):
    mean, sd = gp.predict([[x]], return_std=True)
    print(f"{x:>6.1f}  {mean[0]:>11.4f} {sd[0]:>9.4f}  {np.abs(train_x - x).min():>21.4f}")
print(f"\n  the latent-function prior standard deviation is {np.sqrt(1.0):.4f}")
print(f"  the observation prior standard deviation, including noise, is {np.sqrt(1.0 + 0.01):.4f}")
#      x   prediction   std dev    nearest observation
#    1.5       1.0387    0.1123                 0.1244
#    3.0       0.2755    0.1647                 0.1948
#    4.5       0.0398    0.9106                 1.6948
#    5.0       0.0155    0.9620                 2.0850
#    5.5       0.0213    0.8903                 1.5850
#    7.0       0.6372    0.1498                 0.0850
#    9.0       0.4124    0.1064                 0.0119
#   13.0      -0.0056    1.0046                 3.0084
#
#   the latent-function prior standard deviation is 1.0000
#   the observation prior standard deviation, including noise, is 1.0050

Observations cluster in two ranges with a gap between them. Where they are dense, the predictive standard deviation is around 0.11. In the middle of the gap it is 0.9620. Beyond every observation, at \(x = 13\), it is 1.0046 — close to the observation prior standard deviation of 1.0050, which includes the WhiteKernel noise variance of 0.01. The latent-function prior standard deviation alone is 1.0000. Under this fixed RBF kernel, the distant observations provide little reduction in predictive uncertainty at that point.

A point-prediction interface for linear regression or boosting does not itself provide uncertainty, but the model families are not incapable of doing so. Linear-regression prediction intervals can widen with leverage away from the observed inputs, and Bayesian linear models also return predictive distributions. A GP supplies a distribution through its kernel and likelihood assumptions; those assumptions, rather than the existence of an interval alone, determine how it extrapolates.

So a GP is the natural choice when the cost of a confidently wrong extrapolation is high, and when the question includes “how sure are you” rather than only “what is the answer.” That is why they underpin Bayesian optimization — the search decides where to sample next by trading off the predicted mean against exactly this variance.

The kernel is the model

“Flexible — you just need to choose a kernel” undersells the choice. The kernel specifies what functions are considered plausible before any data arrives, and the data adjusts that belief far less than it looks.

• RBF — \(k(x,x’) = \exp\left(-\frac{\|x-x’\|^2}{2\ell^2}\right)\) — asserts infinitely smooth functions. Nearby points are strongly correlated, and the length scale \(\ell\) decides how fast that correlation decays. It is why the variance in the gap above grew as it did.
• Linear gives straight-line behavior, periodic gives repetition, and Matérn gives a tunable roughness that is usually a better default than RBF for real measurements.
• Kernels add and multiply. A trend plus a seasonal cycle plus noise is a sum of three kernels, and that compositionality is where most of a GP’s modeling power lives.

A GP with the wrong kernel is confidently wrong in exactly the way the variance is supposed to prevent — an RBF kernel on a series with jumps will interpolate smoothly across them and report narrow intervals while doing it.

What it costs

Conditioning requires inverting an \(n \times n\) covariance matrix, which is \(O(n^3)\) in time and \(O(n^2)\) in memory. A single dense 10,000-by-10,000 float64 covariance matrix occupies 800 MB before factorization workspace and other arrays. Runtime depends on hardware, kernel evaluation, implementation and numerical settings; larger datasets often require approximation. Sparse and inducing-point methods exist for that, at the cost of some of the exactness that was the attraction.

Which places GPs where the data is small, each observation is expensive, and the uncertainty matters — physical experiments, hyperparameter search, spatial interpolation, small clinical datasets. On a million rows of tabular data, compare scalable alternatives such as gradient-boosted trees and approximate GPs on the task, runtime budget and uncertainty requirement. Faster fitting and better predictive performance are separate claims to measure.

References

• scikit-learn User Guide, Gaussian Processes
• scikit-learn API reference, GaussianProcessRegressor