3

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?

DanielRch
  • 472
  • 1
  • 4
  • 12
  • 1
    Are the nonzero blocks in $C_{i}$ much smaller than $B$ and always in the same position in the matrix? – Brian Borchers Oct 26 '18 at 21:25
  • Is your intent to solve $A_i x = f_i$ for multiple $A_i$ and (possibly) multiple $f_i$? – user7440 Oct 26 '18 at 22:48
  • The blocks are always much smaller than $B$; @user7440 yes, as part of time stepping – DanielRch Oct 27 '18 at 10:14
  • 1
    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 – Richard Zhang Oct 31 '18 at 17:19
  • 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 :-) – Richard Zhang Oct 31 '18 at 17:20
  • I think $B$ is PD, but $A_i$ is PSD in that case, would the Cholesky update work here? – DanielRch Oct 31 '18 at 18:58

2 Answers2

4

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.

Anton Menshov
  • 8,672
  • 7
  • 38
  • 94
2
  • 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$.
user7440
  • 617
  • 5
  • 15