5
$\begingroup$

Say if I already have the LU factorization of a square matrix $A$, is there an efficient way to get the LU factorization of $I+A$? (We may assume all the matrix I mentioned is invertible.)

I know from SMW formula that a low-rank update of an invertible matrix can be computed efficiently given the original matrix can be inverted efficiently, but I am not sure whether a similar method can be applied here.

By the way, I am not sticking to LU factorization. We may also assume $Ax=b$ can be solved in another way, such as an iterative method, or even $A^{-1}$ is available explicitly. We may also assume $A$ is symmetric if this helps.

$\endgroup$
3
  • 2
    $\begingroup$ Not an answer, but implicit integrators often have to repeatedly solve systems with matrices of the form $I - \gamma A$ and these typically only use the sparsity pattern of $A$ instead of factorizations of $A$ so I suspect that means the answer is no in general. If there were an efficient way to do this, this would be the primary application. $\endgroup$
    – whpowell96
    Commented May 19 at 1:33
  • $\begingroup$ Related math.stackexchange.com/q/4350639/532409 $\endgroup$
    – Quillo
    Commented May 19 at 12:04
  • 1
    $\begingroup$ Related scicomp.stackexchange.com/q/32801/19136 $\endgroup$ Commented May 19 at 16:10

1 Answer 1

2
$\begingroup$

It is doubtful in the extreme that the ability to solve linear systems of the form $$Ax=b$$ efficiently can generally be used to accelerate the solution of linear systems of the form $$(I + A)x = b.$$ Why do I say so? There is a variety of applications that would benefit tremendously if such a device existed. Examples include the implicit solvers for which @whpowell96 has already mentioned, the solution of parabolic differential equation Crank-Nicolson method, the construction and use of rational Krylov subspaces to approximate functions of matrices or the solution of Sylvester matrix equations using the low-rank ADI method. These methods all hinge on our ability to solve linear systems of the form $$(s I - A)x=b$$ for different values of the scalar $s$. If such as scheme had been discovered, say, within the last few months, then the entire community for numerical linear algebra would know about it by now.

That said, there at least three very special cases where information about $A$ could be utilized.

  1. If $A$ is symmetric positive definite and $s \ge 0$. Then $sI+A$ is also symmetric positive definite. If $A$ is sparse and $PAP^T = LL^T$ has a Cholesky factor $L$ that is suitably sparse, then the Cholesky factor of the matrix $$P(sI+A)P^T$$ will of course have the same sparsity pattern. In other words, the data structures needed to carry out the factorization of $A$ can be utilized to carry out the factorization of $sI+A$. Admittedly, this is a small bonus, but if we are doing 50 large solves pr. second to make a real-time simulation, then there is no time to spare.

  2. If $\|A\|<1$, then the Neumann series imply that $I-A$ is nonsingular and $$(I - A)^{-1} = \sum_{j=0}^{\infty} A^j.$$

  3. If $A$ is nonsingular and $\|A^{-1}\| < 1$, then $A-I$ is nonsingular, because $$A-I = A(I - A^{-1})$$ and we can apply the second point.

  4. As noted by @whpowell96, the Krylov subspace $$K_j(A,b) = \text{span}\{A^ib \: : \: i = \in \{0,1,2,\dots,j-1\}\}$$ is identical to $K_j(I+A,b)$. Hence, it is possible to recycle a Krylov subspace from the problem $Ax=b$ to the problem $(I+A)x=b$ with the same right-hand side. If $I$ can be regarded as small relative to $A$, then one would expect the convergence to be similar in both cases.

Unfortunately, the Neumann series tends to converge quite slowly, so it is very rare that we gain enough accuracy/speed from utilizing a truncated Neumann series.

$\endgroup$
1
  • 1
    $\begingroup$ If $Ax=b$ is solved via a Krylov method, then Krylov subspace recycling methods could potentially solve $(A+I)x=b$ quicker because the Krylov subspaces of $A$ and $A+I$ with the same starting vector are identical. $\endgroup$
    – whpowell96
    Commented May 22 at 20:12

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .