Suppose I have the original large, sparse linear system: $A\textbf{x}_0=\textbf{b}_0$. Now, I do not have $A^{-1}$ as A is too large to factor or any sort of decomposition of $A$, but assume that I have the solution $\textbf{x}_0$ found with an iterative solve.Now, I wish to apply a small rank update to the diagonal of A (change a few of the diagonal entries): $(A+D)\textbf{x}_1=\textbf{b}_0$ where $D$ is a diagonal matrix with mostly 0's on its diagonal and a few nonzero values. If I had $A^{-1}$ I would be able to take advantage of the Woodbury formula to apply an update to the inverse. However, I do not have this available. Is there anything I can do short of just resolving the whole system all over again? Is there some way perhaps that I can come up with a preconditioner $M$ which is easy\easier to invert, such that $MA_1 \approx A_0$, so that all I would have to do if I have $\textbf{x}_0$ is apply $M^{-1}$ and an iterative method would converge in a couple/few iterations?
Solving a system with a small rank diagonal update
linear algebra;iterative method;sparse
Save in the columns of two matrices $B$ and $C$ all vectors $b_j$ to which you applied the matrix in the previous iterations and the results $c_j=Ab_j$.For each new system $(A+D)x'=b'$ (or $Ax=b'$, which is the special case $D=0$), approximately solve the overdetermined linear system $(C+DB)y\approx b'$, e.g., by selecting a subset of the rows (possibly all) and using a dense least square method. Note that only the selected part of $C+DB$ needs to be assembled; so this is a fast operation!Put $x_0=By$. This is a good initial approximation with which to start the iteration for solving $(A+D)x'=b'$. In case further systems must be processed, use the matrix vector products in this new iteration to extend the matrices $B$ and $C$ on the resulting subsystem.If the matrices $B$ and $C$ do not fit into main memory, store $B$ on disk, and select the subset of rows in advance. This allows you to keep in core the relevant part of $B$ and $C$ needed to form the least squares system, and the next $x_0$ can be computed by one pass through $B$ with little use of core memory.The rows should be selected in such a way that they approximately correspond to a coarse discretization of the full problem. Taking five times more rows than the total number of expected matrix vector multiplies should be enough.Edit: Why does this work? By construction, the matrices $B$ and $C$ are related by $C=AB$. If the subspace spanned by the columns of $B$ contains the exact solution vector $x'$ (a rare but simple situation) then $x'$ has the form $x'=By$ for some $y$. Substituting this into the equation defining $x'$ gives the equation $(C+DB)y= b'$. Thus in this case, the above process gives as starting point $x_0=By=x'$, which is the exact solution. In general, one cannot expect $x'$ to lie in the column space of $B$, but the starting point generated will be the point in this cloumn space closest to $x'$, in a metric determined by the selected rows. Thus it is likely to be a sensible approximation. As more systems are processed, the column space grows and the approximation will be likely to improve a lot, so that one can hope to converge in fewer and fewer iterations.Edit2: About the subspace generated: If one solves each system with a Krylov method, the vectors used to get the starting point for the second system span the Krylov subspace of the first right hand side. Thus one gets a good approximation whenever this Krylov subspace contains a vector close to the solution of your second system.In general, the vectors used to get the starting point for the $(k+1)$st system span a space containing the Krylov subspace of the first $k$ right hand sides.