Solving a system with a small rank diagonal update

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?

• Are you starting with a good preconditioner for $A$ and want to know how to update it? What rank is the update? (A rank $1000$ update is "small" compared to a matrix of size $10^9$ but not small in terms of iteration count.) – Jed Brown Jun 2 '12 at 1:55
• $A$ is around size $10^6$ to $10^7$, and the update is <1000 (likely <100) elements. I am using a diagonal type of preconditioner for A which works really well, so updating that would be trivial, but I was wondering if there is anything better that I can do rather than resolving the new system from scratch. – Costis Jun 2 '12 at 1:58
• The solution of one system doesn't tell you much about it. If you solve the same system multiple times, the inverse map on those vectors (and/or associated Krylov spaces) give you some information that can be used to accelerate convergence. How many systems are you solving in each case? – Jed Brown Jun 2 '12 at 2:36
• Currently I am only solving for one RHS ($\textbf{b}$ vector) with each $A$ matrix before modifying $A$. – Costis Jun 2 '12 at 4:31

1. 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$.

2. 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!

3. 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.

• Thanks, Prof. Neumaier. I will try this out. Could you perhaps give me a brief explanation as to how this works? – Costis Jun 5 '12 at 1:41
• Also, what if I want to solve the same system for many different RHS vectors? i.e. $A\textbf{x}_0=\textbf{b}_0$, $A\textbf{x}_1=\textbf{b}_1$, $A\textbf{x}_2=\textbf{b}_2$, etc. Is there any information that I can use from the previous solves to speed up the subsequent ones? – Costis Jun 5 '12 at 1:44
• @Costis: A resolve with the same matrix is just the special case $D=0$ of the general problem. For your first question see the edit. – Arnold Neumaier Jun 5 '12 at 8:00
• @Costis: I added a little more detail to step 2. - If you write up the application, please send ma a preprint. – Arnold Neumaier Jun 5 '12 at 14:34
• Thanks for the explanation! Why can't I just solve the overdetermined system $(C + DB)y \approx b'$ by using a QR factorization based approach and using all the rows instead of just a subset? I guess as the number of columns of C and B increase, I may have to get rid of some rows to make the operation run faster. Sure, I will write up a description of the system and e-mail it to you. I actually think that it is possible to come up with an application-specific scheme that might work better than the most general case. Thanks! – Costis Jun 6 '12 at 2:44