The Infinite Gaussian Vector
9 min read
Suppose a handful of measurements are all you have of some smooth function: a few temperature readings over a day, three or four runs of an expensive simulator, a curve you cannot afford to sample everywhere. A least-squares fit through those points commits to one function and offers no sense of how much to trust it between them. A Gaussian process (GP) instead returns every function consistent with the data at once, each one weighted by how well it fits, so the width of the prediction band is not a diagnostic bolted on afterward but a number the model computes for free. Long the default surrogate in geostatistics and formalized for machine learning by Rasmussen & Williams [RW, ’06], a GP earns that width honestly: the whole family of consistent functions comes from ordinary Gaussian conditioning, and the one part that is actually fit to data, the kernel’s hyperparameters, is fit by maximizing a single likelihood rather than minimizing an ad hoc loss.
One Gaussian Over Every Input
A Gaussian process on an input domain is a collection of random variables \(f(x)\), one for every input \(x\), such that for any finite set of inputs \(x_1, \dots, x_m\) the vector \((f(x_1), \dots, f(x_m))\) is jointly Gaussian. Two functions specify the whole object: a mean function, which we take to be zero throughout by centering the data first, and a kernel \(k(x, x')\), the covariance \(\mathrm{Cov}(f(x), f(x'))\) between any two inputs. Any finite draw is then literally
\[ (f(x_1), \dots, f(x_m)) \sim \mathcal{N}(\mathbf{0}, \mathbf{K}) \]
with \([\mathbf{K}]_{ij} = k(x_i, x_j)\), so sampling a GP on a grid of query points is nothing but sampling one multivariate Gaussian.
The kernel is where the modeling lives, and the belief it has to encode is simple: inputs that are close together should have similar function values. The default choice, used throughout this page, is the squared-exponential kernel
\[ k(x, x') = \exp\!\left(-\frac{(x - x')^2}{2\ell^2}\right), \qquad \ell > 0, \]
with lengthscale \(\ell\) setting how fast that similarity decays. At \(\ell = 1\), two points a unit apart have covariance \(k(0,1) = e^{-1/2} \approx 0.61\), while two points three units apart have covariance \(k(0,3) = e^{-4.5} \approx 0.011\): nearby inputs are still strongly linked, distant ones essentially independent. That single number, decaying with distance, is the entire prior belief a GP encodes before it sees any data.
Nothing stops you from asking for a next point: it arrives as one more coordinate of a slightly bigger Gaussian vector, consistent with the points you already had, and there is no largest \(m\) where the process runs out of coordinates to hand you. “Infinite-dimensional” names exactly that: not a new kind of randomness beyond the ordinary multivariate Gaussian, but the promise that every finite question gets a Gaussian answer of matching size, and all those answers agree where they overlap.
Training Is Conditioning
Now specialize the arbitrary index set above to a concrete one: suppose \(n\) of the inputs, the training points \(x_1, \dots, x_n\), have been observed with noise, giving \(y_i = f(x_i) + \epsilon_i\) and \(\epsilon_i \sim \mathcal{N}(0, \sigma^2)\) i.i.d., and fix a separate query point \(x_*\) where we want to know \(f_* := f(x_*)\).
Conditioning Formulas: Writing \(\mathbf{y} \in \mathbb{R}^n\) for the observations, \(\mathbf{K} \in \mathbb{R}^{n \times n}\) for the training kernel matrix, and \(\mathbf{k}_* \in \mathbb{R}^n\) for the cross-covariances \([\mathbf{k}_*]_i = k(x_i, x_*)\), the general posterior is
\[ f_* \mid \mathbf{y} \sim \mathcal{N}(\mu_*, \sigma_*^2), \qquad \mu_* = \mathbf{k}_*^\top (\mathbf{K} + \sigma^2 \mathbf{I})^{-1} \mathbf{y}, \qquad \sigma_*^2 = k(x_*, x_*) - \mathbf{k}_*^\top (\mathbf{K} + \sigma^2 \mathbf{I})^{-1} \mathbf{k}_*. \]
Proof of the Conditioning Formulas
The pair \((f_*, \mathbf{y})\) is jointly Gaussian, since both are linear in the same underlying Gaussian vector \((f_*, f(x_1), \dots, f(x_n), \epsilon_1, \dots, \epsilon_n)\). Its covariance blocks are:
\(\mathrm{Var}(f_*) = k(x_*, x_*)\),
\(\mathrm{Cov}(f_*, \mathbf{y}) = \mathbf{k}_*^\top\), and
\(\mathrm{Cov}(\mathbf{y}, \mathbf{y}) = \mathbf{K} + \sigma^2 \mathbf{I}\),
the last because \(\mathrm{Cov}([\mathbf{y}]_i, [\mathbf{y}]_j) = k(x_i, x_j) + \sigma^2 \mathbf{1}[i=j]\) with the noise independent of everything else.
Write \(\mathbf{b} = (\mathbf{K} + \sigma^2\mathbf{I})^{-1} \mathbf{k}_* \in \mathbb{R}^n\), so the claimed mean is \(\mu_* = \mathbf{b}^\top \mathbf{y}\). Define the residual \(r = f_* - \mathbf{b}^\top \mathbf{y}\), what is left of \(f_*\) after subtracting that candidate mean. The vector \(\mathbf{b}\) is exactly the choice that decorrelates the residual from the data:
\[ \begin{aligned} \mathrm{Cov}(r, \mathbf{y}) &= \mathrm{Cov}(f_*, \mathbf{y}) - \mathbf{b}^\top \mathrm{Cov}(\mathbf{y}, \mathbf{y}) \\ &= \mathbf{k}_*^\top - \mathbf{k}_*^\top (\mathbf{K}+\sigma^2\mathbf{I})^{-1} (\mathbf{K}+\sigma^2\mathbf{I}) \\ &= \mathbf{0}. \end{aligned} \]
The pair \((r, \mathbf{y})\) is jointly Gaussian, and jointly Gaussian variables with zero covariance are independent, so conditioning on \(\mathbf{y}\) leaves \(r\)’s distribution untouched: still mean zero, with whatever variance it had unconditionally. Now condition the identity \(f_* = r + \mathbf{b}^\top \mathbf{y}\) on the data: the second term is a known constant and the first is an independent mean-zero Gaussian, so \(f_* \mid \mathbf{y}\) is Gaussian with mean \(\mathbf{b}^\top \mathbf{y} = \mu_*\) and variance \(\mathrm{Var}(r)\). Expanding that variance by bilinearity and substituting \(\mathbf{b}\),
\[ \begin{aligned} \mathrm{Var}(r) &= \mathrm{Var}(f_*) - 2\, \mathbf{k}_*^\top \mathbf{b} + \mathbf{b}^\top (\mathbf{K}+\sigma^2\mathbf{I}) \mathbf{b} \\ &= k(x_*,x_*) - 2\,\mathbf{k}_*^\top \mathbf{b} + \mathbf{k}_*^\top \mathbf{b} \\ &= k(x_*,x_*) - \mathbf{k}_*^\top(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{k}_*, \end{aligned} \]
where the middle step used \(\mathbf{b}^\top(\mathbf{K}+\sigma^2\mathbf{I}) = \mathbf{k}_*^\top\) so that \(\mathbf{b}^\top(\mathbf{K}+\sigma^2\mathbf{I})\mathbf{b} = \mathbf{k}_*^\top \mathbf{b}\), which is exactly \(\sigma_*^2\).
Interpretation: \(\mu_*\) pushes each observed value out through the kernel, weighted by how much the query point resembles it, after untangling the training points’ own correlations with each other through \((\mathbf{K}+\sigma^2\mathbf{I})^{-1}\); \(\sigma_*^2\) starts at the prior variance \(k(x_*,x_*)\) and subtracts exactly the variance the training data explains, so it shrinks to (nearly) zero at an observed point and relaxes back toward the prior far from all of them. The mean formula is also a familiar predictor under a different name: \(\mathbf{k}_*^\top(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y}\) is precisely the kernel ridge regression predictor with ridge parameter \(\sigma^2\), the same function an optimization derives by minimizing a penalized squared loss. Conditioning a Gaussian prior and regularizing a least-squares fit land on the same formula from two different directions.
Learning Is Maximizing One Likelihood
Everything above assumed the kernel’s hyperparameters \(\theta = (\ell, \sigma^2)\) were already known. In practice they are fit to the very data being conditioned on, and the mechanism is a second, completely different piece of math: not conditioning, but maximizing a likelihood.
The observation \(y_i = f(x_i) + \epsilon_i\) is a sum of two independent Gaussians, the GP’s own value at \(x_i\) and the observation noise, and sums of independent Gaussians stay Gaussian with the covariances adding (see Why Sums Become Gaussian). Applied to the whole training vector at once, \(\mathbf{y} = \mathbf{f} + \boldsymbol{\epsilon}\) with \(\mathbf{f} \sim \mathcal{N}(\mathbf{0}, \mathbf{K})\) and \(\boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2\mathbf{I})\) independent, so
\[ \mathbf{y} \sim \mathcal{N}\big(\mathbf{0},\, \mathbf{K} + \sigma^2\mathbf{I}\big). \]
That is the marginal likelihood: the density of the data the GP would assign before seeing it, as a function of \(\theta\) through \(\mathbf{K}\). Taking the log of the multivariate Gaussian density at \(\mathbf{y}\) gives
\[ \ln p(\mathbf{y} \mid \theta) = -\frac{1}{2}\, \mathbf{y}^\top (\mathbf{K}+\sigma^2\mathbf{I})^{-1} \mathbf{y} \;-\; \frac{1}{2} \ln \big|\mathbf{K}+\sigma^2\mathbf{I}\big| \;-\; \frac{n}{2}\ln(2\pi), \]
and fitting the kernel means choosing \(\theta\) to maximize this, typically by gradient-based search since no closed form solves it directly (the gradient itself is available in closed form, via the identities \(\partial_\theta \ln|\mathbf{A}| = \mathrm{tr}(\mathbf{A}^{-1}\partial_\theta \mathbf{A})\) and \(\partial_\theta \mathbf{A}^{-1} = -\mathbf{A}^{-1}(\partial_\theta \mathbf{A})\mathbf{A}^{-1}\), which is what libraries actually call).
Interpretation: the last term is a constant. The quadratic term measures how surprised the model is by \(\mathbf{y}\), in the metric \((\mathbf{K}+\sigma^2\mathbf{I})^{-1}\) the kernel itself defines, so it is small exactly when \(\theta\)’s correlation structure matches the pattern actually in the data. The log-determinant term is a complexity penalty, and the diagonal of \(\mathbf{K}+\sigma^2\mathbf{I}\) shows exactly why. It is \(k(x,x)+\sigma^2\) everywhere, a constant that does not depend on the lengthscale, so the trace of \(\mathbf{K}+\sigma^2\mathbf{I}\), the sum of its eigenvalues, is fixed as \(\ell\) varies. By the AM–GM inequality, a determinant (a product of nonnegative eigenvalues summing to a fixed trace) is largest when the eigenvalues are equal and shrinks as they spread apart. A short lengthscale decorrelates every pair of training points, pushing \(\mathbf{K}+\sigma^2\mathbf{I}\) toward a multiple of the identity, equal eigenvalues, and close to the maximum possible determinant; a long lengthscale concentrates the same total variance onto a few directions, spreading the eigenvalues and shrinking the determinant. So the penalty is largest exactly for the most flexible kernel, the one willing to explain any pattern in the data, and it is paid whether or not that flexibility was needed. Maximizing the sum of the two terms is an automatic Occam’s razor: a lengthscale too short fits the data well but pays the maximal penalty regardless, a lengthscale too long pays little penalty but fits poorly, and the optimum sits where the two derivatives cancel.
In the plot, conditioning is visibly a narrowing, not a re-fitting: the posterior draws on the right are still draws from a Gaussian, just one whose covariance has had the observed directions removed, so they agree with the data and with each other near it, and spread back out into ordinary prior-like variation between points, and beyond them.
Where the Band Gets Used
Bayesian optimization uses exactly \(\sigma_*^2\) to decide where to spend the next expensive evaluation, trading a high posterior mean against a wide band rather than exploring blindly. Any surrogate model for a black-box function that is costly to query, a simulation, a physical experiment, a full model retraining, inherits the same posterior variance as its native measure of where to look next.
A Gaussian process turns “fit a function” and “know how much to trust the fit” into one object, and both halves of getting there, the prediction and the kernel, reduce to Gaussian identities: condition a joint Gaussian to predict, and maximize its density to learn.