# Efficient computation of marginalized multivariate normal likelihood

In general,if we know that the marginal Gaussian distribution for some variable $$\textbf{x}$$ and a conditional Gaussian distribution for some $$\textbf{y}|\textbf{x}$$ of the forms: $$p(\textbf{x}) = \mathcal{N}(\textbf{x}|\boldsymbol{\mu},\Lambda^{-1})$$

$$p(\textbf{y}|\textbf{x}) = \mathcal{N}(\textbf{y}|A\textbf{x}+\textbf{b},L^{-1})$$ then the marginal distribution of $$\textbf{y}$$ and the conditional distribution of $$\textbf{x}|\textbf{y}$$ are given by:

$$p(\textbf{y}) = \mathcal{N}(\textbf{y}|A\boldsymbol{\mu}+\textbf{b},L^{-1}+A\Lambda^{-1}A^{T})$$ $$p(\textbf{x}|\textbf{y}) = \mathcal{N}(\textbf{x}| \Sigma \left\{ A^{T} L ( \textbf{y} - \textbf{b} ) + \Lambda \boldsymbol{\mu} \right\} , \Sigma)$$

where : $$\Sigma = (\Lambda + A^{T}LA)^{-1}$$

In my situation, I have a prior and likelihood: $$\boldsymbol{\beta} \sim \mathcal{N}(0,I_p\sigma^2_\beta)$$ $$\hat{\boldsymbol{\beta}} | \boldsymbol{\beta} \sim \mathcal{N}(\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}}^{-1}\boldsymbol{\beta},\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}}).$$ Where $$\hat{\textbf{S}}$$ is a diagonal matrix and $$\hat{\textbf{R}}$$ is a symmetric correlation matrix. By making substitutions I arrive at the following marginalized likelihood :

$$\hat{\boldsymbol{\beta}}|\sigma_\beta^2 \sim \mathcal{N}(0,\sigma_\beta^2\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}}^{-2}\hat{\textbf{R}}\hat{\textbf{S}}+\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}})$$

My question is this, if I want to find the value of $$\sigma_\beta^2$$ that maximizes the marginalized likelihood, do I have to recompute the cholesky of that entire covariance matrix from scratch every time I try out a new value of $$\sigma_\beta^2$$, or do I have other options?

If $$\hat{\textbf{S}}=\text{diag}(1)$$, then we can simply take the eigenvalue decomposition of $$\textbf{R}$$: $$\hat{\boldsymbol{\beta}}|\sigma_\beta^2 \sim \mathcal{N}(0,\sigma_\beta^2\textbf{R}^2+\textbf{R}) = \mathcal{N}(0,\sigma_\beta^2\textbf{Q}\textbf{D}^2\textbf{Q}+\textbf{Q}\textbf{D}\textbf{Q}) \\ \rightarrow \textbf{Q}^T\hat{\boldsymbol{\beta}}|\sigma_\beta^2 \sim \mathcal{N}(0,\sigma_\beta^2\textbf{D}^2+\textbf{D})$$ Where $$\textbf{Q}$$ are the eigenvectors of $$\textbf{R}$$ and $$\textbf{D}$$ is the diagonal matrix of eigenvalues. That's really easy to compute for different values of $$\sigma_\beta^2$$. With $$\hat{\textbf{S}}\neq\text{diag}(1)$$ None of that works. My intuition is that I'm trying to get a cheap full rank update of an eigenvalue decomposition, which is in general a no-go, but I would love to hear thoughts.

• Could you clarify in linear algebra terms what it is you are trying to achieve? For example, do you want to solve linear systems of the form A + \sigma^2 B, for positive definite matrices A and B? That would be easy, so I guess that’s not it. – Amit Hochman Nov 13 '20 at 17:43