# Fast computation/estimation of a low-rank linear system

Linear systems of equations are pervasive in computational statistics. One special system I have encountered (e.g., in factor analysis) is the system

$$Ax=b$$

where $$A=D+ B \Omega B^T$$ Here $D$ is a $n\times n$ diagonal matrix with a strictly positive diagonal, $\Omega$ is an $m\times m$ (with $m\ll n$) symmetric positive semi-definite matrix, and $B$ is an arbitrary $n\times m$ matrix. We are asked to solve a diagonal linear system (easy) that has been perturbed by a low-rank matrix. The naive way to solve the problem above is to invert $A$ using Woodbury's formula. However, that doesn't feel right, since Cholesky and QR factorizations can usually speed up solution of linear systems (and normal equations) dramatically. I recently came up on the following paper, that seems to take the Cholesky approach, and mentions numerical instability of Woodbury's inversion. However, the paper seems in draft form, and I could not find numerical experiments or supporting research. What is the state of the art to solve the problem I described?

• @gappy, did you consider using QR (or Cholesky) decomposition for matrix $\Omega^{-1}+BD^{-1}B^T$ (the middle term in Woodburry formula)? The remaining operations are simple matrix multiplications. The main source of instability is then the calculation of $\Omega^{-1}$. Since $m<<n$ I suspect this application of QR or Cholesky combined with Woodbury will be faster than QR on all matrix $A$. This of course is no state of art, just general observations. Jun 23 '11 at 13:25
• I suspect that what Matthias Seeger advocates is within $\epsilon$ of the state of the art, he is a very bright bloke and these sorts of issues crop up repeatedly in the kind of models he investigates. I use Cholesky based methods for the same reasons. I suspect there is discussion in "Matrix Computations" by Golub and Van Loan, which is the standard reference for this sort of thing (although I don't have my copy to hand). Jun 23 '11 at 15:07
• Note that by taking $\bar{B} = D^{-1/2} B$ your problem is equivalent to solving the system $(I + \bar{B}\Omega \bar{B}^T)x = \bar{b}$ where $\bar{b} = D^{-1/2} b$. So, that simplifies the problem a bit right there. Now, letting $\Sigma = \bar{B} \Omega \bar{B}^T$, we know that $\Sigma$ is positive semidefinite with at most $m$ positive eigenvalues. Since $m \ll n$, finding the $m$ largest eigenvalues and corresponding eigenvectors can be done in various ways. The solution is then $x = Q(I + \Lambda)^{-1} Q^T \bar{b}$ where $\Sigma = Q \Lambda Q^T$ gives the eigendecomposition of $\Sigma$. Jun 23 '11 at 15:38
• Small corrections: (1) Equivalent system is $(I + \bar{B} \Omega \bar{B}^T) D^{1/2} x = \bar{b}$ and (2) Final solution is $x = D^{-1/2} Q (I + \Lambda)^{-1} Q^T D^{-1/2} b$. (I had dropped a $D^{1/2}$ in front of $x$ in both cases.) Notice that all inverses are of diagonal matrices and so are trivial. Jun 23 '11 at 15:51
• @mpiktas: I think you meant $\Omega^{-1} + B^T D^{-1} B$ since in the version you wrote the matrix product is not well-defined due a dimension mismatch. :) Jun 23 '11 at 17:43