Marginal distribution of Gaussians

Question

I want to calculate the log marginal likelihood for a Gaussian Process regression, for that and by GP definition I have the prior: $$ p(\textbf{f} \mid X) = \mathcal{N}(\textbf{0} , K)$$

Where $ K $ is the covariance matrix given by the kernel. And the likelihood is (a factorized gaussian):

$$ p(\textbf{y} \mid \textbf{f}, X) = \mathcal{N}(\textbf{f} ,\sigma_n²I)$$

where $ \textbf{f} $ are the training outputs (the values of the function) with some random gaussian noise term with variance $ \sigma_n²I $ So the log marginal likelihood is calculated as follows: $$ p(\textbf{y}\mid X) = \int p(\textbf{y} \mid \textbf{f}, X) p(\textbf{f} \mid X) d\textbf{f}$$

which ends up in the solution

$$ y \sim \mathcal{N}(\textbf{0} , K + \sigma_n²I) $$ so the log marginal likehood is given by: $$ \log p(\textbf{y} \mid X) = -\frac{1}{2} \textbf{y}^T (K + \sigma_n²I)^{-1}\textbf{y} -\frac{1}{2} \log |K+\sigma_n²I| -\frac{n}{2}\log(2\pi)$$ ($n$ is the number of training points).

The book I'm currently following (Rasmussen, C. E., & Williams, C. K. (2006). Gaussian processes for machine learning, Page 19 ) uses Gaussian identities to get to the solution, specifically the the multiplication of two gaussian functions, but I'm a bit loss on how to apply them to get to the solution, ¿how do I calculate the marginal distribution?

Answer 1

There is an easy way to do this kind of problem and a hard way.

The easy way is to notice that

\begin{align} y & = f + \varepsilon \\ & \varepsilon \sim \mathcal{N}(0, \sigma^2I) \end{align}

The expectation and variance of $ y $ give

\begin{align} E[y] & = E[f] + E[\varepsilon] = 0 \\ V[y] & = V[f] + V[\varepsilon] = K + \sigma^2 I \end{align}

Sums of normals are normal, so $ y \sim \mathcal{N}(0, K + \sigma^2 I) $

The hard way is to multiply out the pdf's and refactor it. Inside the exponential term,

\begin{align} y^T I/\sigma^2 y - 2 f^T I/\sigma^2 y + f^T I/\sigma^2 f + f^T K^{-1} f \\ \end{align}

Move terms outside of the exponential that do not have to do with $f$ and complete the quadratic form,

\begin{align} & f^T A^{-1} f - 2f^Tb \\ A &= [ K^{-1} + I/\sigma^2]^{-1} \\ b & = I/\sigma^2 y \\ \end{align}

The completed quadratic form will have mean $ h $ defined as

\begin{align} h & = A b \\ c & = - b^T A b \\ (x & - h) A^{-1} (x-h) + c \\ \end{align}

Pull out the $ c $ term and notice that this is Gaussian, so it has a known normalizing constant. The $ c $ combines with the terms involving $ y $ that we moved outside of the integral in the first step, and find another Gaussian:

\begin{align} y ^T[ I/\sigma^2 + I/\sigma^2[I/\sigma^2 + K^{-1} ]^{-1} I/\sigma^2] y \end{align}

The Woodbury matrix inversion formula gives another form of the covariance matrix,

$$ I/\sigma^2 + I/\sigma^2[I/\sigma^2 + K^{-1} ]^{-1} I/\sigma^2 = \sigma^2I + K $$

See https://en.wikipedia.org/wiki/Woodbury_matrix_identity and let \begin{align} A & = \sigma^2I \\ U & = I \\ V & = I \\ C & = K^{-1} \end{align}

Therefore, $$ y \sim \mathcal{N} (0, \sigma^2I + K )$$