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?


1 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 )$$

  • $\begingroup$ I don't follow your mathematics. $\endgroup$ Commented Jun 7, 2018 at 2:26
  • $\begingroup$ Read more carefully maybe? Its pretty clear. $\endgroup$
    – debo
    Commented Jun 7, 2018 at 15:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.