4
$\begingroup$

In the problem I am dealing with, I require to repeatedly solve $Ax=b$ where $A$ is a weighted Laplacian matrix of a sparse graph. The right-hand side remains constant. However each time I solve the system, only one weight of the graph changes, effectively changing 4 coefficients in the Laplacian matrix (2 diagonal and 2 off-diagonal entries).

I am currently using the GMRES solver in Petsc, and using the previous solution as the initial guess. However the coefficient changes by 1 or 2 orders of magnitude, and it is not as fast as if the coefficient were changing slightly.

I was wondering if there is anyway to solve this problem any faster than what I have been currently doing. Perhaps something which takes advantage of the linear nature of the problem and involves a direct method.

$\endgroup$
  • 3
    $\begingroup$ Doesn't changing one single weight lead to a rank-1 update of $A$? If so, why not use Sherman-Morrison? $\endgroup$ – Rodrigo de Azevedo Jan 26 '18 at 13:18
  • 3
    $\begingroup$ How big and dense is your matrix? Is using a sparse direct factorization possible? If so, then as Rodrigo commented, using Sherman-Morrison-Woodbury would be a good approach. $\endgroup$ – Brian Borchers Jan 26 '18 at 15:09
  • 1
    $\begingroup$ Have you considered Krylov subspace recycling? This idea is suitable for sequences of large linear systems where the matrix and/or the right hand side change "slowly". I'll refer you to this paper by Parks and related work by Eric de Sturler and collaborators. $\endgroup$ – GoHokies Jan 26 '18 at 18:41
  • 1
    $\begingroup$ You don't want or need to store the inverse. You might be able to store an LU factorization. $\endgroup$ – Brian Borchers Jan 27 '18 at 1:44
  • 1
    $\begingroup$ 1e6 entries in a sparse matrix is quite small, with storage requirements measured in megabytes rather than gigabytes. If you're saying that the size of the matrix is 1e6 by 1e6, then depending on the sparsity it may or may not be reasonable to compute a sparse factorization of the matrix. $\endgroup$ – Brian Borchers Jan 27 '18 at 5:46
4
$\begingroup$

It appears to me that you are dealing with a sequence of linear systems $A_j x_j = b$, where $A_{j+1}$ is a low rank modification of $A_{j}$.

In your case, I would investigate if the Krylov subspace $K = \text{range $V$}$ which you have built to solve one linear system, say, $A x = b$ is relevant for the solution of the next problem, i.e., $$(A + \Delta A) (x + \Delta x) =b.$$ Specifically, I would first solve the small dense problem $$V^T(A + \Delta A) V z = V^T b$$ and use $x_0 = Vz$ as my initial guess for $x + \Delta x$ when applying GMRES to the next problem.

I am gambling that this initial guess will serve you better than the solution of the original problem.

In your case, the small dense problem is a low rank modification of the Hessenberg system $$ H w = V^T b = e_1, \quad H = V^T A V$$ which is solved internally by GMRES. Here you may be well served with an application of the Sherman-Morrison-Woodbury formula mentioned in the comments.

$\endgroup$
4
$\begingroup$

Michael Saunders wrote a sparse LU package that can do rank-1 updates, LUSOL. You could try to use that, since you write that direct solvers are viable for your problem.

$\endgroup$

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.