Special considerations when using the Woodbury matrix identity numerically.

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.

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.

• 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? – fedja Nov 8 '11 at 3:19
• 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. – fedja Nov 8 '11 at 4:45
• Oops, that should be $A^{−1}[I−UB(I+U^TA^{−1}UB)^{−1}U^TA^{−1}]$, of course :). – fedja Nov 8 '11 at 12:46
• Thanks for the tips! I'll try this later today... stay tuned. – Kiyo Nov 8 '11 at 13:37
• By the way, what inversion algorithm are you using? – fedja Nov 8 '11 at 19:46

Alternative idea: switch to solving with the extended matrix $\begin{bmatrix}A & U\\\\ U^* & B\end{bmatrix}$ --- this could be stabler.
• 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. – Federico Poloni Nov 8 '11 at 13:54
• 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. – Federico Poloni Nov 8 '11 at 22:30