At one step of an algorithm (i=1: something big) I have to draw a vector from the normal distribution $X_p^{(i)}\sim N((A^{(i)} + D^{(i)}_p)^{-1}m, (A^{(i)} + D^{(i)}_p)^{-1})$ for $p=1:P$, with a large $P$ (say 5,000).
$A^{(i)}$ is symmetric positive definite and rank $K$, $A^{(i)}$ is $K\times K$ symmetric positive definite with $K$ is much smaller than $P$. $(A^{(i)})^{-1}$ is cheap relative to the rest of the calculation.
So I'd like to efficiently get $(A + D^{(i)}_p)^{-1}$ (or a matrix square root).
I've tried using the matrix inversion lemma (aka Woodbury) to get a quick update from $A^{-1}$ or $(A+D_p)^{-1}$ to $(A+D_{p+1})^{-1}$, but it takes me in circles because it's a full rank update.
It feels like this should have a simple solution since the $D^{(i)}_p$ are all diagonal, but if there is it escapes me...