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

Question

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.