# Inverting $(A + D_p)$ efficiently for many diagonal $D_p$ (and $A$ SPD)

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...

-
possible duplicate of Diagonal update of the inverse of $XDX^T$ – J. M. Apr 30 '11 at 1:27
Not a duplicate; that question has $A+XDX'$ which is not full rank, the update there is on $D$ not $A$ ($A=XDX'$ in the notation of that question), and further the update there isn't full rank either since it's a diagonal matrix with mostly 0 on the diagonal. So yeah... they're pretty different :) – JMS May 5 '11 at 23:35
Yes, sorry about that; I somehow misread your condition for $D_p^{(i)}$. Anyway, the "close vote" seems to be gone now... – J. M. May 5 '11 at 23:39
No problem. I'll confess to a less-than-clear problem statement as well! – JMS May 6 '11 at 16:09

I assume you mean $A$ is symmetric positive semidefinite with rank $k$. If you can do an efficient low rank factorization of $A = U U^T$, then Woodbury update is the way to go for this. Further if you know the rank $k$ before hand then you could actually get a low-rank factorization in $\mathcal{O}(k^2n)$ where $n$ is the size of $A$.
Sorry, I stated the conditions poorly: $A$ is $K\times K$ symmetric positive definite. So like I mentioned the Woodbury identity doesn't buy me much; taking $A = LL^T$ gives me $D_p^{-1}-D_p^{-1}L(I+L^T D_p^{-1}L)^{-1} L^TD_p^{-1}$ – JMS Jan 29 '11 at 21:58