Solve $Ax=b$ repeatedly where $A$ is a sparse weighted Laplacian matrix with changing weights

Question

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.

Answer 1

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.

Answer 2

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.