The Infinite Gaussian Vector
6 min read
A curve can fit every measurement and still be badly wrong between them. If each measurement requires another run of an expensive simulator, we want to know how much confidence to place in the predictions between runs.
A Gaussian process gives us both predictions and a measure of their uncertainty [RW, ’06]. The beautiful simplification is that we only need a finite Gaussian vector, even though the unknown function has infinitely many values.
One Gaussian Over Every Input
Let \(f\) be the unknown function. A Gaussian process (GP) models its values at any finite collection of inputs \(x_1, \dots, x_m\) as a Gaussian vector. We specify a mean for each value and a covariance for each pair of values. The covariance describes how the values vary together.
Take the mean to be zero, centering the data first, and let \(k(x,x')\) give the covariance between \(f(x)\) and \(f(x')\). This function \(k\) is called the kernel. The resulting distribution is:
\[ (f(x_1), \dots, f(x_m)) \sim \mathcal{N}(\mathbf{0}, \mathbf{K}) \]
with \([\mathbf{K}]_{ij} = k(x_i, x_j)\). Sampling a GP on a grid therefore means sampling this Gaussian vector.
The kernel encodes how similar the function values should be at nearby inputs. The default choice here 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 nearly independent. This decay with distance specifies the prior relation between function values.
In the figure, each new point is one more coordinate of a larger Gaussian vector, consistent with the coordinates already drawn. “Infinite-dimensional” means that this construction works for every finite set of inputs and agrees on overlaps. Ordinary Gaussian conditioning of this vector gives the prediction and its uncertainty, and maximizing its likelihood fits the kernel; the next two sections do each in turn.
Conditioning on Observations
Suppose we observe the function with noise at \(n\) training points \(x_1, \dots, x_n\):
\[ y_i = f(x_i) + \epsilon_i, \qquad \epsilon_i \stackrel{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2). \]
Fix a separate query point \(x_*\) and write \(f_* = f(x_*)\) for its unknown function value.
With one observation \(y_1\), this is ordinary bivariate Gaussian conditioning. The observation has variance \(k(x_1,x_1)+\sigma^2\) and covariance \(k(x_1,x_*)\) with \(f_*\), so:
\[ f_* \mid y_1 \sim \mathcal{N}\!\left( \frac{k(x_1,x_*)}{k(x_1,x_1)+\sigma^2}\,y_1,\; k(x_*,x_*)-\frac{k(x_1,x_*)^2}{k(x_1,x_1)+\sigma^2} \right). \]
The posterior mean weights the observation by its covariance with the query, and the posterior variance subtracts the variance explained by that observation. For \(n\) observations, the same calculation replaces scalar division with a matrix solve.
Write \(\mathbf{y} \in \mathbb{R}^n\) for the observations and \(\mathbf{K} \in \mathbb{R}^{n \times n}\) for the training kernel matrix. The vector \(\mathbf{k}_* \in \mathbb{R}^n\) contains the cross-covariances \([\mathbf{k}_*]_i = k(x_i, x_*)\).
Claim (Conditioning formulas): The posterior at the query 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
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. Because \(\mathrm{Cov}(r,\mathbf{y})=\mathbf{0}\), \(r\) and \(\mathbf{y}\) are independent, so conditioning on \(\mathbf{y}\) leaves \(r\) with its unconditional mean zero and variance \(\mathrm{Var}(r)\). In the identity \(f_* = r + \mathbf{b}^\top \mathbf{y}\), conditioning makes the second term constant. Thus \(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\).
\(\mu_*\) weights the observations by their relation to the query after accounting for correlations among the training points through \((\mathbf{K}+\sigma^2\mathbf{I})^{-1}\). Under the chosen kernel and noise model, the variance \(\sigma_*^2\) starts at the prior variance \(k(x_*,x_*)\) and subtracts the part explained by the training data. It shrinks near observations and returns toward the prior far from them. The mean \(\mathbf{k}_*^\top(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y}\) is also the kernel ridge regression predictor with ridge parameter \(\sigma^2\). Kernel ridge regression derives it by minimizing a penalized squared loss.
In the plot, the posterior remains Gaussian, but conditioning removes variance in directions explained by the observations. The draws agree near the data and spread out farther away.
Choosing the Kernel
Everything above assumed the kernel’s hyperparameters \(\theta = (\ell, \sigma^2)\) were known. In practice we estimate them by maximizing the marginal likelihood.
The observation \(y_i = f(x_i) + \epsilon_i\) is a sum of two independent Gaussians: the GP’s value at \(x_i\) and the observation noise. Sums of independent Gaussians stay Gaussian, with their covariances adding (see Why Sums Become Gaussian). For the training vector, \(\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. Their sum has distribution:
\[ \mathbf{y} \sim \mathcal{N}\big(\mathbf{0},\, \mathbf{K} + \sigma^2\mathbf{I}\big). \]
The density of this distribution at the observed data is the marginal likelihood, viewed 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). \]
Fitting the kernel means choosing \(\theta\) to maximize this expression, typically by gradient-based search because there is no closed-form maximizer.
How the likelihood balances fit and flexibility
The gradient is available in closed form from \(\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}\), the identities implemented by GP libraries.
The last term is constant in \(\theta\). The quadratic term measures the size of \(\mathbf{y}\) in the metric \((\mathbf{K}+\sigma^2\mathbf{I})^{-1}\), so it is small when \(\theta\)’s correlation structure matches the pattern in the data. To isolate the effect of the lengthscale, hold the noise variance fixed. The log-determinant term is a complexity penalty whose behavior follows from the diagonal of \(\mathbf{K}+\sigma^2\mathbf{I}\). 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 (the arithmetic mean of nonnegative numbers is at least their geometric mean), 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 the training points and pushes \(\mathbf{K}+\sigma^2\mathbf{I}\) toward a multiple of the identity, where its determinant is close to the maximum possible value. A long lengthscale concentrates the same total variance in fewer directions, spreading the eigenvalues and shrinking the determinant. The negative log-determinant term therefore penalizes the flexible short-lengthscale kernel whether or not the data require that flexibility.
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. Examples include simulations, physical experiments, and full model retraining.