Bayesian Updating for Gaussians: A Gentle Derivation of the Posterior Formula

Posterior derivations for multivariate Gaussians are often presented as a “known formula” without explanation. In this post, we will derive the posterior for Gaussian Processes step by step, starting from the scalar case and building up to the multivariate case.

1-minute read: In essence, everything about a Gaussian distribution is contained in its exponent—specifically, the quadratic form inside the exponential. In the univariate case, the coefficient of the squared term \((x-\mu)^2\) is \(-\frac{1}{2\sigma^2}\). This means the mean \(\mu\) shifts the centre, while the coefficient determines the precision, i.e., the inverse of the variance \(\sigma^2\).

The elegance of the Gaussian process is that both the prior over functions and the likelihood of the noisy observations are Gaussian. Because of this, the posterior is also Gaussian. Deriving the posterior mean and covariance reduces to a simple algebraic operation: multiply the two Gaussian exponents (which means adding their quadratic forms) and then complete the square to read off the new mean and covariance.

In the multivariate setting of a GP, this same principle applies to vectors and matrices. It is just that the quadratic form changes: the scalar squared term \((x-\mu)^2\) becomes a bilinear form \((f - \mu)^T \Sigma^{-1} (f - \mu)\). Division by a scalar variance becomes multiplication by an inverse covariance matrix (the precision matrix). When we combine the prior and the likelihood, we add these precision matrices—not the covariances themselves. Completing the square now requires a matrix inversion, but the conceptual core remains identical: the posterior mean is a precision‑weighted average of the prior mean and the observed data, and the posterior covariance is simply the inverse of the summed precision matrices. So, while the algebra grows heavier, the Gaussian magic stays the same: add the exponents, complete the square, and read off the updated Gaussian.

1. The Intuition: Estimating a Hidden Value from a Noisy Measurement

Imagine you have a hidden scalar value \(f\), and you believe before seeing anything that it is drawn from:

\[f \sim \mathcal{N}(0, k)\]

So your prior says: \(f\) is probably near 0, with variance \(k\). And the prior density is:

\[ p(f) \propto \exp\!\left(-\frac{\left(f-0\right)^2}{2k}\right) \]

Likelihood: now you observe a noisy measurement

\[ y = f + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \sigma^2) \]

so

\[ p(y \mid f) \propto \exp\!\left(-\frac{(y - f)^2}{2\sigma^2}\right) \]

Posterior: By Bayes

\[ p(f \mid y) \propto p(y \mid f)\, p(f) \]

\[ p(f \mid y) \propto \mathcal{N}(y \mid f, \sigma^2) \, \mathcal{N}(f \mid 0, k) \]

The product of two Gaussians is Gaussian. You can complete the square in the exponent, and you get:

\[f \mid y \sim \mathcal{N}\!\left( \frac{k}{k + \sigma^2} \, y, \;\; k - \frac{k^2}{k + \sigma^2} \right)\]

Let’s walk through step by step to see how this works.

1.1 Combine the Exponents

The posterior exponent, ignoring constant factors, is the sum of the two exponents (from the standard Gaussian density):

\[ -\frac{f^2}{2k} - \frac{(y - f)^2}{2\sigma^2} \]

Multiply by -2 to make the algebra cleaner:

\[ \mathcal{Q}(f) = \frac{f^2}{k} + \frac{(y - f)^2}{\sigma^2} \]

1.2 Expand the Quadratic

Expand the second term:

\[ \frac{(y-f)^2}{\sigma^2} = \frac{y^2 - 2yf + f^2}{\sigma^2} \]

So:

\[ \mathcal{Q}(f) = \frac{f^2}{k} + \frac{f^2}{\sigma^2} - \frac{2y}{\sigma^2}f + \frac{y^2}{\sigma^2} \]

Group the quadratic and linear terms in \(f\):

\[ \mathcal{Q}(f) = \left(\frac{1}{k} + \frac{1}{\sigma^2}\right) f^2 - \frac{2y}{\sigma^2} f + \frac{y^2}{\sigma^2} \]

The last term is constant with respect to \(f\), so we can ignore it for now.

Define:

\[ A = \frac{1}{k} + \frac{1}{\sigma^2} \]

Then:

\[ \mathcal{Q}(f) = A f^2 - \frac{2y}{\sigma^2} f + \text{const} \]

1.3 Complete the Square

We want to write this as:

\[ A(f - \mu)^2 + \text{const} \]

because a Gaussian posterior has exponent:

\[ -\frac{(f-\mu)^2}{2v} \]

which is proportional to:

\[ \frac{(f-\mu)^2}{v} \]

Expanding \(A(f-\mu)^2\):

\[ A(f-\mu)^2 = A f^2 - 2A\mu f + A\mu^2 \]

Compare with \(A f^2 - \frac{2y}{\sigma^2} f\).

We need:

\[ 2A\mu = \frac{2y}{\sigma^2} \]

so

\[ \mu = \frac{y}{\sigma^2 A} \]

But \(A = \frac{1}{k} + \frac{1}{\sigma^2}\), hence:

\[ \mu = \frac{y}{\sigma^2} \cdot \frac{1}{\frac{1}{k} + \frac{1}{\sigma^2}} \]

Simplify the denominator:

\[ \frac{1}{k} + \frac{1}{\sigma^2} = \frac{\sigma^2 + k}{k\sigma^2} \]

Thus:

\[ \mu = \frac{y}{\sigma^2} \cdot \frac{k\sigma^2}{k + \sigma^2} = \frac{k}{k + \sigma^2}\, y \]

So the posterior mean is a weighted combination of the prior mean 0 and the observation \(y\), with weight \(\frac{k}{k+\sigma^2}\) on the data.

1.4 Find the Posterior Variance

The posterior precision is \(A\), because the exponent is:

\[ -\frac{1}{2} A (f - \mu)^2 \]

and in a Gaussian density the exponent is:

\[ -\frac{(f-\mu)^2}{2v} \]

so \(v = 1/A\).

Therefore:

\[ v = \frac{1}{\frac{1}{k} + \frac{1}{\sigma^2}} = \frac{k\sigma^2}{k + \sigma^2} \]

This can also be written as:

\[ v = k - \frac{k^2}{k + \sigma^2} \]

which is the prior variance \(k\) minus the uncertainty removed by the observation.

1.5 Final Posterior

We have:

\[ f \mid y \sim \mathcal{N}\!\left( \frac{k}{k + \sigma^2}\, y,\; \frac{k\sigma^2}{k + \sigma^2} \right) \]

Often written as:

\[f \mid y \sim \mathcal{N}\!\left( \frac{k}{k + \sigma^2} \, y, \;\; k - \frac{k^2}{k + \sigma^2} \right)\]

1.6 Intuition

• If the prior variance \(k\) is large compared with the noise \(\sigma^2\), then the posterior mean is close to \(y\); you trust the observation.
• If the noise \(\sigma^2\) is large compared with \(k\), then the posterior mean stays near 0; you trust the prior.
• The posterior variance is always smaller than both \(k\) and the noise, because information from both sources is combined.

This scalar computation is exactly the pattern behind the full GP posterior. The only difference is that in the GP case, \(f\) becomes a vector \(\mathbf{f}\), the scalar variance \(k\) becomes the covariance matrix \(K\), and the scalar division becomes matrix multiplication by an inverse.

2. The GP Setting: Multiple Latent Values

Now suppose we have \(N\) latent function values at inputs:

\[\mathbf{f} = [f(x_1), \dots, f(x_N)]^\top\]

Our GP prior says:

\[\mathbf{f} \sim \mathcal{N}(\mathbf{0}, K)\]

where \(K\) is the \(N \times N\) covariance matrix with \(K_{ij} = k(x_i, x_j)\).

We observe noisy values:

\[\mathbf{y} = \mathbf{f} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \sigma_n^2 I)\]

So the likelihood is:

\[p(\mathbf{y} \mid \mathbf{f}) = \mathcal{N}(\mathbf{y} \mid \mathbf{f}, \sigma_n^2 I)\]

We want:

\[p(\mathbf{f} \mid \mathbf{y}) \propto p(\mathbf{y} \mid \mathbf{f})\, p(\mathbf{f})\]

Because both are Gaussian, the posterior is Gaussian. We need its mean and covariance.

2.1 Gaussian Density for Scalars and Vectors

For a scalar random variable \(x \in \mathbb{R}\), with mean \(\mu\) and variance \(\sigma^2\),

\[ \mathcal{N}(x \mid \mu, \sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \]

where \(\sigma > 0\).

For a \(N\)-dimensional vector \(\mathbf{x} \in \mathbb{R}^N\), with mean vector \(\boldsymbol{\mu} \in \mathbb{R}^N\) and covariance matrix \(\boldsymbol{\Sigma} \in \mathbb{R}^{N \times N}\),

\[ \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{N/2} |\boldsymbol{\Sigma}|^{1/2}} \exp\left( -\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}) \right) \]

where \(\boldsymbol{\Sigma}\) is symmetric positive definite and \(|\boldsymbol{\Sigma}|\) is its determinant.

The scalar case is the special case when \(N=1\), \(\boldsymbol{\mu} = \mu\), and \(\boldsymbol{\Sigma} = \sigma^2\).

Within the exponent, \((\mathbf{x} - \boldsymbol{\mu})^\top (\mathbf{x} - \boldsymbol{\mu})\) is a scalar (a single number), because:

• Each \(\mathbf{x}\) is from one of the \(N\) dimensions, having its own \(\mathbf{\mu}\).
• \(\mathbf{x}\) and \(\boldsymbol{\mu}\) are \(N\)-dimensional column vectors, so \(\mathbf{x} - \boldsymbol{\mu}\) is also a \(N \times 1\) column vector.
• Its transpose \((\mathbf{x} - \boldsymbol{\mu})^\top\) is a \(1 \times N\) row vector.
• Multiplying a \(1 \times N\) row vector by a \(N \times 1\) column vector gives a \(1 \times 1\) matrix, which is effectively a scalar.

Explicitly:

\[ (\mathbf{x} - \boldsymbol{\mu})^\top (\mathbf{x} - \boldsymbol{\mu}) = \sum_{i=1}^N (x_i - \mu_i)^2 \]

which is the squared Euclidean distance between \(\mathbf{x}\) and \(\boldsymbol{\mu}\).

In the Gaussian exponent, we actually have:

\[ (\mathbf{x}-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}) \]

which is also a scalar, because \(\boldsymbol{\Sigma}^{-1}\) is a \(N \times N\) matrix, and the whole product is \(1 \times 1\). This scalar is called the Mahalanobis distance squared between \(\mathbf{x}\) and \(\boldsymbol{\mu}\).

2.2 The Quadratic Form for Multivariate Gaussians

As shown above, the Gaussian density for a vector \(\mathbf{x}\) has the form:

\[ p(\mathbf{x}) \propto \exp\left(-\frac{1}{2}(\mathbf{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})\right) \]

The expression inside the exponent is a quadratic form — a generalisation of a square. When we multiply two Gaussians (prior and likelihood), their exponents add, producing a new quadratic form. To read off the resulting Gaussian’s mean and covariance, we must complete the square: rewrite that combined quadratic as a single perfect square plus a constant.

2.3 Completing the Square for Vectors

We will first repeat the scalar case, explained in section 1, but with a generalised notation that will make it easier to see how to generalise to vectors.

For a single variable \(x\), suppose we have:

\[ \text{exponent} = -\frac{(x-a)^2}{2s_1} - \frac{(x-b)^2}{2s_2} \]

We combine terms in \(x^2\) and \(x\):

\[ = -\frac{1}{2}\left(\frac{1}{s_1} + \frac{1}{s_2}\right)x^2 + \left(\frac{a}{s_1} + \frac{b}{s_2}\right)x + \text{const} \]

Completing the square means writing:

\[ = -\frac{1}{2}P(x - m)^2 + \text{const} \]

where \(P = \frac{1}{s_1} + \frac{1}{s_2}\) is the precision (inverse variance) and \(m\) is the mean. Solving gives:

\[ m = \frac{\frac{a}{s_1} + \frac{b}{s_2}}{\frac{1}{s_1} + \frac{1}{s_2}} \]

which is exactly the posterior mean of a Gaussian–Gaussian conjugate update.

the same process works with matrices. We start with two quadratic terms:

• Prior: \(\frac{1}{2}\mathbf{f}^\top K^{-1} \mathbf{f}\)
• Likelihood: \(\frac{1}{2}(\mathbf{y} - \mathbf{f})^\top (\sigma_n^2 I)^{-1} (\mathbf{y} - \mathbf{f})\)
• Because \((\sigma_n^2 I)^{-1} = \frac{1}{\sigma_n^2} I\), the likelihood: \((\mathbf{y}-\mathbf{f})^\top (\sigma_n^2 I)^{-1}(\mathbf{y}-\mathbf{f}) = \frac{1}{\sigma_n^2}(\mathbf{y}-\mathbf{f})^\top I(\mathbf{y}-\mathbf{f})\)

Add them and keep only terms involving \(\mathbf{f}\) in posterior likelihood \(\mathcal{L}(\mathbf{f})\):

\[ \mathcal{L}(\mathbf{f}) = \frac{1}{2}\mathbf{f}^\top K^{-1} \mathbf{f} + \frac{1}{2\sigma_n^2}(\mathbf{y} - \mathbf{f})^\top (\mathbf{y} - \mathbf{f}) \]

Expand the likelihood term:

\[ (\mathbf{y}-\mathbf{f})^\top(\mathbf{y}-\mathbf{f}) = \mathbf{y}^\top\mathbf{y} - \mathbf{y}^\top\mathbf{f} - \mathbf{f}^\top\mathbf{y} + \mathbf{f}^\top\mathbf{f} \]

Now note:

• \(\mathbf{y}\) and \(\mathbf{f}\) are both \(N \times 1\) column vectors.
• \(\mathbf{y}^\top\mathbf{f}\) is a scalar.
• \(\mathbf{f}^\top\mathbf{y}\) is the transpose of that scalar.
• Therefore \(\mathbf{y}^\top\mathbf{f} = \mathbf{f}^\top\mathbf{y}\).

So:

\[ -\mathbf{y}^\top\mathbf{f} - \mathbf{f}^\top\mathbf{y} = -2\mathbf{y}^\top\mathbf{f} \]

Thus:

\[ (\mathbf{y}-\mathbf{f})^\top(\mathbf{y}-\mathbf{f}) = \mathbf{f}^\top\mathbf{f} - 2\mathbf{y}^\top\mathbf{f} + \mathbf{y}^\top\mathbf{y} \]

\[ \mathcal{L}(\mathbf{f})= \frac{1}{2}\mathbf{f}^\top K^{-1} \mathbf{f} + \frac{1}{2\sigma_n^2}\mathbf{f}^\top \mathbf{f} - \frac{1}{\sigma_n^2}\mathbf{y}^\top \mathbf{f} + \text{const} \]

Combine the quadratic coefficients:

\[ = \frac{1}{2}\mathbf{f}^\top \left(K^{-1} + \frac{1}{\sigma_n^2}I\right)\mathbf{f} - \frac{1}{\sigma_n^2}\mathbf{y}^\top \mathbf{f} + \text{const} \]

Now this has the form:

\[ \frac{1}{2}\mathbf{f}^\top A \mathbf{f} - \mathbf{b}^\top \mathbf{f} \]

where \(A = K^{-1} + \sigma_n^{-2} I\) and \(\mathbf{b} = \sigma_n^{-2}\mathbf{y}\).

We want to write this as:

\[ \frac{1}{2}(\mathbf{f} - \boldsymbol{\mu})^\top A (\mathbf{f} - \boldsymbol{\mu}) + \text{const} \]

Expanding that perfect square gives:

\[ \frac{1}{2}\mathbf{f}^\top A \mathbf{f} - \boldsymbol{\mu}^\top A \mathbf{f} + \frac{1}{2}\boldsymbol{\mu}^\top A \boldsymbol{\mu} \]

Comparing with our expression, we need:

\[ \boldsymbol{\mu}^\top A = \mathbf{b}^\top \quad \Rightarrow \quad \left(\boldsymbol{\mu}^\top A\right)^\top = \left(\mathbf{b}^\top\right)^\top \quad \Rightarrow \quad A^\top \boldsymbol{\mu} = \mathbf{b} \]

In our posterior derivation, \(A\) is the posterior precision:

\[ A = K^{-1} + \frac{1}{\sigma_n^2}I \]

Both \(K^{-1}\) and \(I\) are symmetric, so \(A\) is symmetric:

\[ A^\top = A \]

Therefore:

\[ \boldsymbol{\mu}^\top A = \mathbf{b}^\top \quad \Rightarrow \quad A \boldsymbol{\mu} = \mathbf{b} \quad \Rightarrow \quad \boldsymbol{\mu} = A^{-1} \mathbf{b} \]

Thus:

\[ \boldsymbol{\mu} = \left(K^{-1} + \frac{1}{\sigma_n^2}I\right)^{-1} \frac{1}{\sigma_n^2}\mathbf{y} \]

and the posterior covariance is \(\Sigma_{\text{post}} = A^{-1}\), i.e.

\[ \Sigma_{\text{post}} = \left(K^{-1} + \frac{1}{\sigma_n^2}I\right)^{-1} \]

2.4 From Precision Form to the Familiar Covariance Form

The Woodbury Matrix Identity:

\[ (A + U C V)^{-1} = A^{-1} - A^{-1} U \left(C^{-1} + V A^{-1} U\right)^{-1} V A^{-1} \]

The matrix \(A = K^{-1} + \sigma_n^{-2} I\) is the posterior precision. Using the Woodbury matrix identity and assuming \(A = K^{-1}, \qquad U = I, \qquad V = I, \qquad C = \sigma^{-2} I\), we can rewrite \(\Sigma_{\text{post}}\) into a standard and familiar form:

\[ \Sigma_{\text{post}} = K - K(K + \sigma_n^2 I)^{-1} K \]

and the posterior mean as:

\[ \boldsymbol{\mu} = K(K + \sigma_n^2 I)^{-1} \mathbf{y} \]

which matches the standard GP formula.

2.5 Summary of the “Square in Exponent” Step

1. Start with prior and likelihood exponents.
2. Expand and collect terms quadratic in \(\mathbf{f}\) and linear in \(\mathbf{f}\).
3. Combine quadratic coefficients into one matrix \(A\) (the precision).
4. Complete the square to read off the posterior mean \(A^{-1}\mathbf{b}\) and covariance \(A^{-1}\).
5. Use the Woodbury identity to express everything in terms of the original covariance \(K\) and noise \(\sigma_n^2 I\).

The inverse \((K + \sigma_n^2 I)^{-1}\) appears precisely because we first worked with the precision form and then transformed back to covariance form. This is the algebraic mechanism behind the GP posterior.

3. The Shortcut: Deriving the Posterior by conditional Gaussian formula

The above sections show a precision route to derive the posterior of a Gaussian process, by multiplying prior and likelihood, completing the squre then using the Woodbury identity. There is a more direct route by applying the conditional Gaussian formula to the joint distribution of \((\mathbf{y}, \mathbf{f})\). This is a standard result in multivariate statistics.

Since \(\mathbf{y} = \mathbf{f} + \boldsymbol{\varepsilon}\) and \(\mathbf{f}, \boldsymbol{\varepsilon}\) are independent, the joint vector is Gaussian with mean zero and covariance:

• \(\text{Cov}(\mathbf{f}, \mathbf{f}) = K\)
• \(\text{Cov}(\mathbf{y}, \mathbf{y}) = \text{Cov}(\mathbf{f}+\boldsymbol{\varepsilon}, \mathbf{f}+\boldsymbol{\varepsilon}) = K + \sigma_n^2 I\)
• \(\text{Cov}(\mathbf{y}, \mathbf{f}) = \text{Cov}(\mathbf{f}+\boldsymbol{\varepsilon}, \mathbf{f}) = K\)

So the joint distribution is:

\[ \begin{bmatrix} \mathbf{y} \\ \mathbf{f} \end{bmatrix} \sim \mathcal{N}\!\left( \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix}, \begin{bmatrix} K + \sigma_n^2 I & K \\ K & K \end{bmatrix} \right) \]

This is the key object. It tells us how the observations and the latent values are related.

We now use the standard formula for conditioning a joint Gaussian. If:

\[ \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \end{bmatrix} \sim \mathcal{N}\!\left( \begin{bmatrix} \boldsymbol{\mu}_a \\ \boldsymbol{\mu}_b \end{bmatrix}, \begin{bmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ba} & \Sigma_{bb} \end{bmatrix} \right) \]

then the conditional distribution of \(\mathbf{b}\) given \(\mathbf{a}\) is:

\[\mathbf{b} \mid \mathbf{a} \sim \mathcal{N}\Big( \boldsymbol{\mu}_b + \Sigma_{ba} \Sigma_{aa}^{-1} (\mathbf{a} - \boldsymbol{\mu}_a),\; \Sigma_{bb} - \Sigma_{ba} \Sigma_{aa}^{-1} \Sigma_{ab} \Big)\]

This formula is the vector analogue of the scalar result above.

Now set:

• \(\mathbf{a} = \mathbf{y}\)
• \(\mathbf{b} = \mathbf{f}\)
• \(\boldsymbol{\mu}_a = \mathbf{0}\), \(\boldsymbol{\mu}_b = \mathbf{0}\)
  
• \(\Sigma_{aa} = K + \sigma_n^2 I\)
  
• \(\Sigma_{ab} = K\)
• \(\Sigma_{ba} = K\)
• \(\Sigma_{bb} = K\)

Plug these into the formula.

Posterior mean:

\[\boldsymbol{\mu}_{\mathbf{f}|\mathbf{y}} = \mathbf{0} + K (K + \sigma_n^2 I)^{-1} (\mathbf{y} - \mathbf{0}) = K (K + \sigma_n^2 I)^{-1} \mathbf{y}\]

Posterior covariance:

\[\Sigma_{\mathbf{f}|\mathbf{y}} = K - K (K + \sigma_n^2 I)^{-1} K\]

Therefore:

\[\mathbf{f} \mid \mathbf{y} \sim \mathcal{N}\Big( K (K+\sigma_n^2 I)^{-1}\mathbf{y}, \; K - K (K+\sigma_n^2 I)^{-1} K \Big)\]

That is exactly the posterior we derived earlier.

4. Why Does This Match the Scalar Case?

Compare with the scalar result:

• Scalar prior variance \(k\) → matrix \(K\)
• Scalar noise variance \(\sigma^2\) → matrix \(\sigma_n^2 I\)
• Scalar division \(\frac{k}{k+\sigma^2}\) → matrix multiplication \(K(K+\sigma_n^2 I)^{-1}\) (because matrix division is multiplication by the inverse)
• Scalar posterior variance \(k - \frac{k^2}{k+\sigma^2}\) → matrix expression \(K - K(K+\sigma_n^2 I)^{-1}K\)

So the matrix formula is nothing but the scalar formula with matrix operations.

5. The Posterior Process for New Points

Once we have the posterior over \(\mathbf{f}\) at the observed inputs, we can extend it to any new input \(x_*\). When we want the distribution of the function value at a new input \(x_*\), we are not observing a new \(y_*\). We are asking:

“Given the posterior over the function after seeing the old data, what do I believe about \(f(x_*)\) at a new location?”

This is a prediction, not a new observation.

If we take our old \(N\) observations \(\mathbf{y}\) and one new latent value \(f_* = f(x_*)\), their joint prior is an \((N+1)\)-dimensional Gaussian:

\[ \begin{bmatrix} \mathbf{y} \\ f_* \end{bmatrix} \sim \mathcal{N}\!\left( \mathbf{0}, \begin{bmatrix} K + \sigma_n^2 I & \mathbf{k}_* \\ \mathbf{k}_*^\top & k(x_*, x_*) \end{bmatrix} \right) \]

where \(\mathbf{k}_*\) is the column vector of covariances between \(x_*\) and the \(N\) training inputs:

\[ \mathbf{k}_* = \begin{bmatrix} k(x_*, x_1) \\ k(x_*, x_2) \\ \vdots \\ k(x_*, x_N) \end{bmatrix} \]

and \(\mathbf{k}_*^\top\) is its transpose, giving a \(1 \times N\) row vector. This makes the full covariance matrix symmetric and \((N+1) \times (N+1)\).

As we already know \(\mathbf{y}\). So we condition this \((N+1)\)-dimensional joint distribution on the observed \(\mathbf{y}\). The result is a 1-dimensional posterior for \(f_*\):

\[ f_* \mid \mathbf{y}, x_* \sim \mathcal{N}\big(m_{\text{post}}(x_*),\, k_{\text{post}}(x_*, x_*)\big) \]

So the final predictive distribution is univariate for a single new point, even though we formed a larger joint distribution to get there.

This is exactly the same Gaussian conditioning formula we used before. Here:

• \(\mathbf{a} = \mathbf{y}\)
• \(\mathbf{b} = f_*\)
• \(\Sigma_{aa} = K + \sigma_n^2 I\)
• \(\Sigma_{ab} = \mathbf{k}(x_*)\)
• \(\Sigma_{ba} = \mathbf{k}(x_*)^\top\)
• \(\Sigma_{bb} = k(x_*, x_*)\)

Plugging these into the conditioning formula:

\[m_{\text{post}}(x_*) = 0 + \Sigma_{ba} \Sigma_{aa}^{-1} (\mathbf{y} - 0) = \mathbf{k}(x_*)^\top (K + \sigma_n^2 I)^{-1} \mathbf{y}\]

\[k_{\text{post}}(x_*, x_*) = \Sigma_{bb} - \Sigma_{ba} \Sigma_{aa}^{-1} \Sigma_{ab} = k(x_*, x_*) - \mathbf{k}(x_*)^\top (K + \sigma_n^2 I)^{-1} \mathbf{k}(x_*)\]