9
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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.) $\endgroup$ – Jed Brown Jun 2 '12 at 1:55
  • $\begingroup$ $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. $\endgroup$ – Costis Jun 2 '12 at 1:58
  • 2
    $\begingroup$ 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? $\endgroup$ – Jed Brown Jun 2 '12 at 2:36
  • $\begingroup$ Currently I am only solving for one RHS ($\textbf{b}$ vector) with each $A$ matrix before modifying $A$. $\endgroup$ – Costis Jun 2 '12 at 4:31
4
$\begingroup$
  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.

$\endgroup$
  • $\begingroup$ Thanks, Prof. Neumaier. I will try this out. Could you perhaps give me a brief explanation as to how this works? $\endgroup$ – Costis Jun 5 '12 at 1:41
  • $\begingroup$ 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? $\endgroup$ – Costis Jun 5 '12 at 1:44
  • $\begingroup$ @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. $\endgroup$ – Arnold Neumaier Jun 5 '12 at 8:00
  • $\begingroup$ @Costis: I added a little more detail to step 2. - If you write up the application, please send ma a preprint. $\endgroup$ – Arnold Neumaier Jun 5 '12 at 14:34
  • $\begingroup$ 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! $\endgroup$ – Costis Jun 6 '12 at 2:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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