I found in a paper that the multi-dimensional integral of a squared exponential correlation function times a multivariate normal has an analytical expression as follows.
\begin{equation}
\int\exp{\left\{-(\theta-c)^T\Omega_t (\theta-c)\right\}}(2\pi)^{-k/2}|V_\theta|^{-\frac{1}{2}}\exp{\left\{-\frac{1}{2}(\theta-m_\theta)^TV_\theta^{-1} (\theta-m_\theta)\right\}} d\theta=
\end{equation}
\begin{equation}
=|I+2V_\theta\Omega_t|^{-1/2}\exp{\left\{-(m_\theta-c)^T(2V_\theta+\Omega_t^{-1})^{-1} (m_\theta-c)\right\}}
\end{equation}
where $c$ is a constant, $\Omega_t$ is a $k$-dimensional diagonal matrix with the inverse of the length scales of the correlation function and $m_\theta$ and $V_\theta$ are the mean and covariance of the MVN distribution, respectively. I would like to see a proof or have a hint as to how the analytic expression is as shown above. I've searched and found the following related posts
Integral of product of two normal distribution densities
https://en.wikipedia.org/wiki/Gaussian_integral#n-dimensional_and_functional_generalization
Thank you in advance.
Best Answer
I was able to find a breakthrough to crack it. Essentially I've used the equation shown in this answer https://math.stackexchange.com/a/2996439/608359, and realised that the matrix in the numerator can be written much more simply as $(\Gamma + \Sigma)^{-1}$. Similarly, the determinants in the denominator can also be simplified in a similar manner. Therefore, it is necessary to convert the correlation function to a MVN distribution, and then to apply the simplified formula
$$\int_{\mathbb{R}^D} \rho_{\mu, \Sigma}(x)\cdot\rho_{\mathbf{0},\Gamma}(x)\,dx=\frac{\exp\left(-\frac 12 (\mu^T(\Gamma+\Sigma)^{-1}\mu) \right)} {\sqrt{(2\pi)^D |\Sigma||\Gamma\Sigma^{-1}+I|}}$$