I would like to solve the following integral: $$ \int{\exp{\left(\left(y_{ij}-\sum_{k=1}^M{A_{ik}r_{jk}}\right)^2\right)}dr_{11}\ldots dr_{PM}}, $$ where index $i$ ranges from $1$ to $N$, $j$ from $1$ to $P$. I know it is closely related to the marginalisation of Gaussian process regression coefficients: $$ \int{\exp{\left(\frac{1}{2}\left(\overline{y}-\overline{\overline{F}}\overline{\beta}\right)^T\overline{\overline{\Sigma}}^{-1}\left(\overline{y}-\overline{\overline{F}}\overline{\beta}\right)\right)}}d\overline{\beta}, $$ see e.g. Andrianakis and Challenor, section 3.1 however, I can't seem to find any derivation of this either.
Thanks in advance.