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

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.

• Since Woodbury is derived from block wise matrix inversion, shouldn't your alternative be equivalent?
– Kiyo
Nov 8, 2011 at 13:45
• 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. Nov 8, 2011 at 13:54
• 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.
– Kiyo
Nov 8, 2011 at 15:09
• 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. Nov 8, 2011 at 22:30
• 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.
– Kiyo
Nov 9, 2011 at 12:56