Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|cite|improve this question
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$.

share|cite|improve this answer
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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.