3
$\begingroup$

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

  • 0
    possible duplicate of [Diagonal update of the inverse of $XDX^T$](http://math.stackexchange.com/questions/3627/diagonal-update-of-the-inverse-of-xdxt)2011-04-30
  • 0
    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 :)2011-05-05
  • 0
    Yes, sorry about that; I somehow misread your condition for $D_p^{(i)}$. Anyway, the "close vote" seems to be gone now...2011-05-05
  • 0
    No problem. I'll confess to a less-than-clear problem statement as well!2011-05-06

1 Answers 1