# Evaluate the multivariate normal using variance matrices $\boldsymbol{\Lambda}+\alpha_{i}\mathbf{a}\mathbf{a}^{T}+\beta_{j}\mathbf{b}\mathbf{b}^{T}$

I need to calculate a huge amount of inverses and determinants to evaluate the pdf of the multivariate Gaussian.

Specifically I need to compute the inverses and determinants of the following variance matrices:

$\boldsymbol{\Sigma}_{ij}=\boldsymbol{\Lambda}+\alpha_{i}\mathbf{a}\mathbf{a}^{T}+\beta_{j}\mathbf{b}\mathbf{b}^{T}\,\,\,\,\,\,\, i=1,...,N\,\,\,\,\,\,\,j=1,...,M$

where

$\mathbf{A}$ is $K\times K$ (Constant)

$\mathbf{a},\mathbf{b}$ are $K\times 1$ (Constants)

$\alpha_{i}, \beta_{j}$ are $1\times 1$ and varies with $i$ and $j$ respectively.

In my application $N=10000$, $M=50$ and $K=10$.

Currently I am using the Woodbury formula twice to compute the inverse

$\boldsymbol{\Sigma}_{ij}^{-1}=\left(\Lambda+\alpha_{i}\mathbf{a}\mathbf{a}^{T}+\beta_{j}\mathbf{b}\mathbf{b}^{T}\right)^{-1}=\tilde{\Lambda}_{i}^{-1}-\tilde{\Lambda}_{i}^{-1}\mathbf{b}\left(\beta_{j}^{-1}+\mathbf{b}^{T}\tilde{\Lambda}_{i}^{-1}\mathbf{b}\right)^{-1}\mathbf{c}^{T}\tilde{\Lambda}_{i}^{-1}$

where

$\tilde{\Lambda}_{i}^{-1}=\left(\Lambda+\alpha_{i}\mathbf{a}\mathbf{a}^{T}\right)^{-1}=\Lambda^{-1}-\Lambda^{-1}\mathbf{a}\left(\alpha_{i}^{-1}+\mathbf{a}^{T}\Lambda^{-1}\mathbf{a}\right)^{-1}\mathbf{a}^{T}\Lambda^{-1}$

I need to do this in a faster way.

Intuitively I think this should be possible since I also need to compute the determinant. I am thinking a smart way of doing cholesky/QR decomposition etc.?

Btw I am using MATLAB as the programming language.

• The Woodbury formula (correctly executed) is quadratic in $K$. It's hard to do better than that, given that usually all $K^2$ entries of $\Lambda$ have to be modified! You ought to consider whether you need to compute the determinants--those could be expensive. They are used only for normalization, which often can be skipped or evaluated in a different way. But since you haven't given any statistical context--why are you doing this?--you have closed off all opportunities for this community to point you in creative directions. Consider adding an explanation of what all this is about. – whuber Mar 24 '15 at 17:54
• Do you really need the inverses? I have serious doubt that you need exactly the inverse and not a product of it with something else. I think you can get a long way with Cholesky. Can you please add a few code of exactly what are you trying to do? For example give a positive definite $A$, instead of using det(A) you could use prod(diag(chol(A)))^2. – usεr11852 says Reinstate Monic May 9 '15 at 1:19