Given \begin{equation} A_i=B+C_i \end{equation} where $A_i$,$B$ and $C_i$, $i=1,\dotsc,N$ are large square matrices, $B$ is symmetric, $C_i$ are zero matrices aside for a square block on the diagonal.

I would like to factor the $A_i$, preferably but not necessarily by $LU$ decomposition, but I don't want to do this repeatedly $N$ times, as $N$ may be large. Is there a method of doing this?

  • 1
    $\begingroup$ Are the nonzero blocks in $C_{i}$ much smaller than $B$ and always in the same position in the matrix? $\endgroup$ – Brian Borchers Oct 26 '18 at 21:25
  • $\begingroup$ Is your intent to solve $A_i x = f_i$ for multiple $A_i$ and (possibly) multiple $f_i$? $\endgroup$ – user7440 Oct 26 '18 at 22:48
  • $\begingroup$ The blocks are always much smaller than $B$; @user7440 yes, as part of time stepping $\endgroup$ – DanielRch Oct 27 '18 at 10:14
  • 1
    $\begingroup$ If you can guarantee $A_i$ and $B$ to be positive definite somehow, then you can use the product form Cholesky update of Goldfarb and Scheinberg. It's similar to the SMW formula, but is far more numerically stable. It is implemented in interior-point solvers like SeDuMi for example $\endgroup$ – Richard Zhang Oct 31 '18 at 17:19
  • $\begingroup$ Conversely, if you can extend the Goldfarb and Scheinberg technique to indefinite symmetric matrices while guaranteeing stability, then that would be a nice contribution worth of publication in Math. Prog. Series A :-) $\endgroup$ – Richard Zhang Oct 31 '18 at 17:20

The structure of $C_i$ gives a hope that $C_i$ can be described as a low-rank update. If that is the case (size of the square block is much smaller than the size of the matrix $B$), one can apply classical methods for low-rank updates based on Sherman-Morrisson-Woodbury identity.

While updating the LU factorization is not the easiest task (the required pivoting might be sensitive to the performed modification) it is certainly possible. To begin with, you might want to read the classic low-rank update section of QR in Golub, Van Loan (Section 12.6 in 2nd edition). Then, proceed with the intricacies of updating the $PA=LU$ factorization.

Note, since $B$ is symmetric, one can take advantage of Cholesky factorization whenever $C_i$ is also symmetric.

However, if $C_i$ is not (enough) low-rank, I don't think there is much hope on any savings here.

  • $\begingroup$ The SMW identity often leads to instabilities, do I have reason to believe it would be stable in this case? $\endgroup$ – DanielRch Oct 29 '18 at 19:05
  • 1
    $\begingroup$ @DanielRch you do not, in general. There are just not a lot of other options. $\endgroup$ – Anton Menshov Oct 29 '18 at 21:09
  • If you plan on solving $A_i x = f_i$, one alternative option is to use a Krylov iterative solver (like GMRES) with $B$ as preconditioner and $B^{-1} f_i$ as starting vector. Then you can prove that all the residuals will belong to the range of $C_i$. The iterations will be performed in the range of $C_i$ (which reduces the dimension of vectors from $N$ down to the rank of $C_i$).
  • For a reference with GMRES, you could check Section 4.2 of Heikkola, Rossi, Toivanen, "A parallel fictitious domain method for three-dimensional Helmholtz equation", SISC, vol. 24, pp. 1567-1588 (2003).
  • Finally, you could recycle the Krylov subspace at every step $i$, as long as $C_{i+1}$ and $f_{i+1}$ are relatively "close" to $C_{i}$ and $f_i$.

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.