Deriving Log Marginal Likelihood for Gaussian Process I am trying to evaluate the following integral marginalized across all possible functions.
$$\mathbb{P}(y|X,\theta) = \int \mathbb{P}(y|f)\  \mathbb{P}(f|X,\theta) \ df$$
In G.P. we assume prior to be of zero mean and a covariance matrix i.e. $\mathbb{P}(f|X,\theta) \sim \mathcal{N}(f;0, K)$. The likelihood is modeled as $\mathbb{P}(y|f) \sim \mathcal{N}(y;f,\sigma_n^2I)$ with a mean $f$ and covariance $\sigma_n^2I$.
So we have;
$$\begin{align}
\mathbb{P}(y|X,\theta) &= \int \mathbb{P}(y|f)\  \mathbb{P}(f|X,\theta)\\
&= \int \mathcal{N}(y;f,\sigma_n^2I)\ \ \mathcal{N}(f;0, K) \tag{1} 
\end{align}$$
As per http://gaussianprocess.org/gpml/chapters/RW.pdf equation (2.30); the log marginal likelihood is $$\text{log}\ \mathbb{P}(y|X,\theta) = -\frac{1}{2}y^TK_y^{-1}y - \frac{1}{2}\text{log}|K_y| - \frac{n}{2}\text{log}\ 2\pi \tag{2}$$
This means the above integral should give me $\mathcal{N}(y;0,K_y)$ where $K_y = K + \sigma_n^2I$. Furthermore, Appendix (A.7) in the same book gives the product of two Gaussians as: $$\mathcal{N}(x;a,A)\ \mathcal{N}(x;b,B) = z^{-1}\mathcal{N}(x;c,C) \tag{3}$$
where $z^{-1} = \frac{1}{(2\pi)^{n/2}}|A+B|^{-\frac{1}{2}}\text{exp}\big\{-\frac{1}{2}(a-b)^T(A+B)^{-1}(a-b)\big\}$
i.e. $z^{-1} = \frac{1}{(2\pi)^{n/2}}|K_y|^{-\frac{1}{2}}\text{exp}\big\{-\frac{1}{2}f^T(K_y)^{-1}f\big\}$
$c = C(A^{-1}a+B^{-1}b) = (\sigma_n^{-2}I + K^{-1})\sigma_n^{-2}f$
$C = (A^{-1} + B^{-1})^{-1} = \sigma_n^{-2}I + K^{-1}$
With this in mind, if I plug (3) in (1) and carry out the integration I should simply get $z^{-1}$ with $f$ replaced by $y$ (not quite sure if I can do this) and since $\mathcal{N}(c,C)$ is a pdf it should give me a value of 1. If I do this I end up with (2) after taking log. But I am not 100% confident on my approach.
Any help in solving the integral would be appreciated.

I was able to get something.
Let $\Sigma = \sigma_n^2I$. Now, (1) can be fleshed out as follows;
\begin{equation*}
\begin{split}
    \mathbb{P}(y|X,\theta) & = \int \frac{1}{(2\pi)^{n/2}} |\Sigma|^{-1/2} \text{exp}\ (-\frac{1}{2} (f-y)^T\Sigma^{-1}(f-y)) \times \frac{1}{(2\pi)^{n/2}} |K|^{-1/2} \text{exp}\ (-\frac{1}{2} (f)^TK^{-1}(f))\ df \\
    & = \frac{1}{(2\pi)^n}\frac{1}{\sqrt{|\Sigma||K|}} \int \text{exp}\ \bigg(-\frac{1}{2}\big[(f-y)^T\Sigma^{-1}(f-y) + f^TK^{-1}f\big]\bigg)\ df \\
    \end{split}
\end{equation*}
Looking at the exponent term:
\begin{equation*}
    \begin{split}
        &= f^T(\Sigma^{-1}+K^{-1})f - 2f^T\Sigma^{-1}y + y^T\Sigma^{-1}y \\
        &= f^T\Pi^{-1}f - 2f^T\Pi^{-1}\nu + y^T\Sigma^{-1}y\\
        &= (f-\nu)^T\Pi^{-1}(f-\nu) - \nu^T\Pi^{-1}\nu + y^T\Sigma^{-1}y
    \end{split}
\end{equation*}
where $\Pi = (\Sigma^{-1}+K^{-1})^{-1}$ and $\nu = \Pi\Sigma^{-1}y$.
By definition we have;
$$\frac{1}{\sqrt{|2\pi\Pi|}} \int \text{exp}\bigg[  -\frac{1}{2} (f-\nu)^T\Pi^{-1}(f-\nu)  \bigg] \ df = 1$$
Plugging this back to the integral gives the following expression;
$$\frac{\sqrt{(2\pi)^n|\Pi|}}{(2\pi)^n\sqrt{|\Sigma||K|}}\ \text{exp}\ \bigg[ \frac{1}{2}(\nu^T\Pi^{-1}\nu - y^T\Sigma^{-1}y)   \bigg]$$
Substitute values for $\Pi$ and $\nu$ we get;
\begin{equation*}
    \begin{split}
        \mathbb{P}(y|X,\theta) &= \frac{1}{(2\pi)^{n/2}} \bigg(\big|\Sigma\big|\big|K\big|\big|\Sigma^{-1}+K^{-1}\big|\bigg)^{-1/2} \text{exp}\ \bigg[-\frac{1}{2}\big(y^T\Sigma^{-1}(\Sigma^{-1}+K^{-1})^{-1}K^{-1}y\big) \bigg] \\
        &= \frac{1}{(2\pi)^{n/2}} \bigg(\big|\Sigma\big|\big|K\big|\big|\frac{\Sigma + K}{\Sigma K}\big|\bigg)^{-1/2} \text{exp}\ \bigg[-\frac{1}{2}\big(y^T\Sigma^{-1}(  \frac{K + \Sigma}{\Sigma K}  )^{-1}K^{-1}y\big) \bigg] \\
        &= \frac{1}{(2\pi)^{n/2}} \bigg(\big|\Sigma + K \big|\bigg)^{-1/2} \text{exp}\ \bigg[-\frac{1}{2}\big(y^T(  K + \Sigma)^{-1}y\big) \bigg] \\
        &= \frac{1}{(2\pi)^{n/2}} \bigg(\big|\sigma_n^2I + K \big|\bigg)^{-1/2} \text{exp}\ \bigg[-\frac{1}{2}\big(y^T(  K + \sigma_n^2I)^{-1}y\big) \bigg] \\
        &= \frac{1}{(2\pi)^{n/2}} \big|K_y \big|^{-1/2} \text{exp}\ \bigg[-\frac{1}{2}y^TK_y^{-1}y \bigg] \\
    \end{split}
\end{equation*}
Taking log on both sides yields (2)
