I am trying to derive a distribution shown in the paper Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models.
We are given \begin{align} p(x) &= \mathcal{N}(x | 0, I) \\ p(y | x , W , \beta) &= \mathcal{N}(y|Wx, \beta^{-1}I), \end{align} where $\beta$ is a scalar.
I want to show
$$ \int p(y|x,W,\beta)p(x)dx = \mathcal{N}(y| 0, WW^\top + \beta^{-1}I). $$
To do this I thought I would use the method decribed here, which utilize "completion of squares", which states
$$ \frac{1}{2}z^\top A z + b^\top z + c = \frac{1}{2}(z + A^{-1}b)A(z + A^{-1 }b) + c - \frac{1}{2}b^\top A^{-1}b $$
Then I have that
\begin{align} p(y| x, W, b)p(x) &\propto \exp(-\frac{1}{2}((y-Wx)^\top(\beta I)(y-Wx) + x^\top I x)) \\ &= \exp(-\frac{1}{2}(\beta y^\top y - \beta y^\top Wx - \beta(Wx)^\top y + \beta(Wx)^\top Wx + x^\top x)) \\ &= \exp(-\frac{1}{2}(\beta y^\top y + 2(- \beta y^\top W) x + x^\top(\beta W^\top W + I )x)) \\ &= \exp(-\frac{1}{2}\beta y^\top y +\beta y^\top W x -\frac{1}{2}x^\top(\beta W^\top W + I )x)\\ &= \exp(-\frac{1}{2}x^\top(\beta W^\top W + I )x +(\beta y^\top W) x-\frac{1}{2} \beta y^\top y) \text{, apply completion of squares} \\ &=\exp(-\frac{1}{2} (x + (\beta W^\top W + I)^{-1} (\beta y^\top W)^\top)^\top(\beta W^\top W + I)(x + (\beta W^\top W + I)^{-1} (\beta y^\top W)^\top) \\ &-\frac{1}{2}\beta y^\top y -\frac{1}{2}(\beta y^\top W)(\beta W^\top W + I)^{-1}(\beta y^\top W)^\top)\\ &= \exp(-\frac{1}{2} (x + (\beta W^\top W + I)^{-1} (\beta y^\top W)^\top)^\top(\beta W^\top W + I)(x + (\beta W^\top W + I)^{-1} (\beta y^\top W)^\top) \\ & \cdot\exp(-\frac{1}{2}y^\top(\beta I + \beta W(\beta W^\top W + I)^{-1}\beta W^\top) y) \\ &\implies \text{covariance of marginal is } (\beta I + \beta W(\beta W^\top W + I)^{-1}\beta W^\top)^{-1} \end{align}
from which I don't see how to get $WW^\top + \beta^{-1}I$.
First, you have the covariance of the marginal slightly wrong - the plus sign should be a minus, as you need to add the extra term to complete the square and subtract the same term to compensate. Thus the marginal covariance is $$(\beta I \color{red}{-} \beta W(\beta W^\top W + I)^{-1}\beta W^\top)^{-1}$$
You can then use the Woodbury matrix identity where $$(A+UCV)^{-1}=A^{-1}-A^{-1}U\;\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}$$ setting $A=\beta^{-1}I$, $U=W$, $C=I$ and $V=W^\top$ results in $$(\beta^{-1}I+WW^\top)^{-1}=\beta I-\beta W(I+\beta W^\top W)^{-1}\beta W^\top$$ Inverting both sides leads to the desired result $$(\beta I-\beta W(\beta W^\top W+I)^{-1}\beta W^\top)^{-1}=WW^\top+\beta^{-1}I$$