10
$\begingroup$

Are there any special considerations when using the Woodbury matrix identity numerically? What is the best metric for numerical stability in this case? Can anyone point me to a good reference?

The special case of the identity that I'm using is: $ (A + UBU^T)^{-1} = A^{-1} - A^{-1}U(B^{-1} + U^TA^{-1}U)^{-1}U^TA^{-1}$

I'm using the identity to speed up a matrix inverse. In my case $A$ is diagonal, and $U$ is rectangular by a factor that varies between 10 and 200. I'm using a sparse form for $U$.

This is failing catastrophically for me. Normal inversion seems stable for condition number of the LHS well beyond $10^{10}$, while the identity is failing around condition number of $10^7$ (everything in double precision). However, condition number of the LHS doesn't seem to be the best metric for when things fail, as I can arrange cases where the failure happens at much higher or lower condition number.

Also, if there is another way to compute this inverse fast, I'm open to suggestions.

Thanks in advance!

Update:

To answer a few questions: yes the dimensions of $B$ are much less than the dimensions of $A$ and both are symmetric and positive definite.

I've checked and the matrices $A$, $B$, and $(B^{-1} + U^TA^{-1}U)$ all have relatively low condition number.

I've also done the following: I've factored the identity to: $(A + UBU^T)^{-1} = A^{-1}[I - U(B^{-1} + U^TA^{-1}U)^{-1}U^TA^{-1}]$

Then one would expect this fail when the second term in square brackets has eigenvalues smaller or equal to $-1$, because adding the identity just translates the eigenvalues by $+1$. This is indeed the case. The translation of the eigenvalues seems to be working fine, which tells me that there is no problem with taking the difference of the two terms. This issue must be in computing the second term.

Further Update:

This reference appears to address the problem. I haven't gotten all the way through it yet but it looks like you can precondition the problem by making sure that the columns of $U$ are all mutually orthogonal. After reading it, it turns out this doesn't help, since this only helps with the condition number of $(B^{-1} + U^TA^{-1}U)$, which I don't think is the problem.

Improve this question
$\endgroup$
7
  • $\begingroup$ Well, one obvious thing is that the left hand side may be well-defined and having low condition number without $A$ or $B$ being individually invertible at all, which would be a computational disaster for the right hand side. It is hard to say more without additional information. Anyway, as far as I can see, you just want a quick inversion of the LHS and do not really care whether it comes from this formula or from something else. Maybe it is a good idea to ask if someone knows a good alternative. Do I understand it right that the size of B is less than that of A and that B is symmetric? $\endgroup$
    – fedja
    Commented Nov 8, 2011 at 3:19
  • $\begingroup$ What happens if you use this formula in a slightly different form: $LHS=A^{-1}[I-UB(I+U^TA^{-1}UB)^{-1}UA^{-1}]$? You'll have to choose between $(U^TA^{-1}U)B$ and $(U^TA^{-1})(UB)$ to see which one works faster but it takes only marginally more time than your original formula anyway and eliminates the inversion of B completely. You still have to invert A, but, since it is diagonal, you can detect any problems there very easily and act accordingly. $\endgroup$
    – fedja
    Commented Nov 8, 2011 at 4:45
  • $\begingroup$ Oops, that should be $A^{−1}[I−UB(I+U^TA^{−1}UB)^{−1}U^TA^{−1}]$, of course :). $\endgroup$
    – fedja
    Commented Nov 8, 2011 at 12:46
  • $\begingroup$ Thanks for the tips! I'll try this later today... stay tuned. $\endgroup$
    – Kiyo
    Commented Nov 8, 2011 at 13:37
  • $\begingroup$ By the way, what inversion algorithm are you using? $\endgroup$
    – fedja
    Commented Nov 8, 2011 at 19:46

1 Answer 1

Reset to default
1
$\begingroup$

Have you tried checking the condition number of the three matrices that you invert in the RHS, instead of the LHS? I am afraid that's the only sure way to tell.

I've heard in many conferences that many experts in matrix computations (led by the late Gene Golub) consider the SMW identity to be "numerically dangerous"; I am not sure that I can find a paper where this is stated black on white, though (sorry for that, I was always curious about that myself).

Alternative idea: switch to solving with the extended matrix $\begin{bmatrix}A & U\\\\ U^* & B\end{bmatrix}$ --- this could be stabler.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Since Woodbury is derived from block wise matrix inversion, shouldn't your alternative be equivalent? $\endgroup$
    – Kiyo
    Commented Nov 8, 2011 at 13:45
  • $\begingroup$ No. When you do SMW, you constrain the inversion to be done "starting from a certain block", which is not necessarily the best thing to do. Think for instance to the case in which $B$ is ill-conditioned. When you invert the big matrix, say, with LU, suitable pivoting is done so that you never have trouble with ill-conditioned subblocks; on the other hand, if you go with SMW, you start by inverting $B$, which is troublesome. $\endgroup$
    – Federico Poloni
    Commented Nov 8, 2011 at 13:54
  • $\begingroup$ Can this be done fast the way SMW lets you? Or do I have to do the whole rank(A) + rank(B) inversion? I guess I'm a little confused about what you are proposing. $\endgroup$
    – Kiyo
    Commented Nov 8, 2011 at 15:09
  • $\begingroup$ Do you really need an inverse, or do you only need to solve linear system with it (i.e., compute $A^{-1}b$ or $c^T A^{-1}$)? Chances are you only need the second; in this case, your choice for a general linear system is between sparse LU and iterative solvers (like GMRES), depending on the size and sparseness of the matrix. Both can exploit the structure of the augmented matrix I am suggesting. In any case, don't use inv() or compute an inverse explicitly unless you really have to. $\endgroup$
    – Federico Poloni
    Commented Nov 8, 2011 at 22:30
  • $\begingroup$ I do need the inverse. It's a statistical code map making code, but you can think of it as adding up a bunch of Fisher matrices. $\endgroup$
    – Kiyo
    Commented Nov 9, 2011 at 12:56

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .