# Updating matrix diagonal with Woodbury matrix identity and maintaining numerical accuracy

I have a dense matrix A and its corresponding inverse $A^{-1}$. The Woodbury matrix identity states:

$$(A + UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}$$

I wish to perform small updates to the diagonal of $A$ and compute corresponding $A^{-1}$. The updates can actually be described as a constant $k$ multiplying an identity matrix: $kI$ and in my case $U=V^T$ where $V$ is a tall matrix with just ones and zeros that selects the desired diagonal positions of $A$ to modify.

Everything is fine so far, except $k$ needs to be quite large (e.g., $10^6$). I am noticing that for large $k$, changes made are not "reversible" without significant error. What I mean by reversible is performing two update steps, one with $kI$ on the original $A$ and another with $-kI$ on the updated matrix to recover the original $A^{-1}$. If $k$ is small enough, it works fine, but otherwise I run into serious accuracy issues and the resulting matrix does not match the initial one closely. For my target application, I need to be able to apply $k$ to diagonal elements or $-k$ at will, so my question is: Is there a clever way to do this which preserves numerical accuracy even for large values of $k$?

Note: I am using double precision for all values stored and for computation. Furthermore, $k$ is a real scalar, but the matrix $A$ is complex, although that should not really matter in this case.

• You were missing an inverse in your formula as per Wikipedia's rendition - I'm not sure if there are other forms. And if $C$ is the identity matrix it wouldn't matter of course. – horchler Dec 7 '13 at 19:18
• No ideas or suggestions? Still unsure of the right way to go about this (besides using really high precision-- I am already using double precision arithmetic for everything currently.) – Costis Jan 19 '14 at 10:47
• It might make the problem more approachable if you provided a bit of code (Matlab -best for me, Python or C) that succinctly demonstrates the issue. It would make your matrix multiplication and inversion methods, order of operations, and even values of A and k concrete as they might be responsible for the issue at hand rather than some inherent problem with Woodbury. – horchler Jan 19 '14 at 15:43
• This sounds like a conditioning problem. Have you considered scaling your domain, or moving your precision? The finest granularity of the IEEE 784 happens when values are around 1.0, and the farther you get the worse they are. Eventually roundoff kills you. If typical values are order-of $10^0$ and you want to interact with $10^6$ there can be pain. It might be less if you divide both by 1e3. They span the same scale - 6 decades - but they both have consistent discretization. You could also use quadruple precision (or better). IEEE 754-2008 binary128 form is quaruple precision. – EngrStudent Mar 8 '15 at 14:28