Fast likelihood evaluation for Gaussian distribution with diagonal plus low rank covariance Let's assume the likelihood
$$
y \sim\mathcal N_p(0, \Sigma + \Lambda\Lambda^\top)
$$
where $\Sigma$ is diagonal and $\Lambda$ is a $p \times d$ matrix with $d \ll p$.
What is the fastest way to compute the Gaussian density at one specific value $y$?
Naive implementations would still scale cubically in $p$ since one needs to invert and compute the determinant of $\Sigma + \Lambda \Lambda^\top)$. Are there some more efficient ways to go?

$\newcommand{\L}{\Lambda}$We need to be able to efficiently evaluate $y^T(\Sigma + \L\L^T)^{-1} y$ and $\det (\Sigma + \L\L^T)$. The standard approach here is to use the Woodbury matrix identity and matrix determinant lemma to change our hard operations from being on $p\times p$ matrices to being on $d\times d$ matrices.
The Woodbury matrix identity gives us
$$
(\Sigma + \L\L^T)^{-1} = \Sigma^{-1} - \Sigma^{-1}\L(I + \L^T \Sigma^{-1} \L)^{-1}\L^T\Sigma^{-1}.
$$
$\Sigma$ is diagonal so $\Sigma^{-1}$ takes $O(p)$ operations. $\Sigma^{-1} \L$ then is scaling each row of $\Lambda$ by the corresponding element of $\Sigma^{-1}$ so this takes $O(dp)$ time. $\L^T (\Sigma^{-1} \L)$ then naively takes $O(d^2 p)$, so the time taken to get to $I + \L^T \Sigma^{-1} \L$ is $O(d^2p)$. This is a $d\times d$ matrix now so inversion is naively $O(d^3)$ so we're at $O(d^3 + d^2p)$.
$\Sigma^{-1}\L(I + \L^T \Sigma^{-1} \L)^{-1}$ is $O(d^2p)$ as well, but then we end up with a $(p\times d) (d\times p)$ matrix multiplication which is $O(p^2d)$, so in all we're at $O(d^3 + d^2p + dp^2)$ and $d \ll p$ means this is effectively $O(dp^2)$ which is better than the $O(p^3)$ we'd have with a direct naive inversion.
For the determinant, the matrix determinant lemma gives us
$$
\det(\Sigma + \L\L^T) = \det(I + \L^T \Sigma^{-1} \L)\det \Sigma.
$$
$\det \Sigma$ is just $O(p)$, and we can do $I + \L^T \Sigma^{-1} \L$ in $O(d^2p)$ as we saw above. We now need the determinant of this matrix which is $d\times d$ so we're at $O(d^3 + d^2 p)$ which is effectively $O(d^2 p)$ given $d \ll p$.
All together we can do this in $O(dp^2)$ operations (and it would be a little faster if we did this for real since matrix multiplication and inversion can be done a little faster than the naive way).
