The Infinite Gaussian Vector
7 min read
Suppose a handful of measurements are all we have of a smooth function, such as several runs of an expensive simulator. We want to predict what another run would produce and express how uncertain that prediction is. A fitted curve supplies a prediction, but we also need a model for how the observed and unobserved values relate.
A Gaussian process (GP) treats the function values at any finite collection of inputs as entries of one joint Gaussian vector. Observing some entries tells us about the others through the formula for conditioning a Gaussian. That calculation gives both the prediction and its uncertainty under the model [RW, ’06].
One Gaussian Over Every Input
What does one Gaussian over every input mean concretely? 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. For any finite set of inputs, the corresponding values have a multivariate Gaussian distribution:
\[ (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.
Training Is Conditioning
How do noisy observations turn into a prediction and an error bar at a new input? 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.
Claim (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
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}\). 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.
Learning Is Maximizing One Likelihood
Where do the lengthscale and the noise level come from? 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. 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. 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.
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.
Where the Band Gets Used
Where does the variance, rather than the mean, earn its keep? 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. The same posterior variance directs search when a Gaussian process models a costly black-box function. Examples include simulations, physical experiments, and full model retraining.
A Gaussian process predicts by conditioning one joint Gaussian and fits its kernel by maximizing the marginal likelihood.