# solving series of linear systems with diagonal perturbations

I would like to solve a series of linear systems Ax=b as quickly as quickly as possible. However, the systems are related. Specifically, each matrix A is given by:

cI + E

where E is a fixed sparse, symmetric positive definite real matrix (unchanged in all the linear systems), I is the identity matrix, and c is a varying complex number.

In other words, I am wondering how to quickly solve a series of complex linear systems which are all identical except for complex perturbations along the diagonal. I should say that the resulting matrices are not necessarily Hermitian, so currently I compute the LU decomposition. This works, but given the large number of rather closely related systems to be solved, I wonder if there is a better way to solve the problem, perhaps by using a more expensive (e.g. QR) decomposition up front.

(Edit for Jiahao: Yes, the bs are all the same.) (Edit for J. Mangaldan: The matrices are of order n=10^5 ~ 10^6, with about 10 times that many nonzeros.)

Update:

I'd like to thank everyone here for their suggestions. My implementation is ugly, but in the end interpolation was the key to a reasonable (10x) speedup. Since the c are quite close (imagine a small region of the complex plane, small in the sense that the spectrum of the matrix E is much larger) I could get away with computing solutions for a subset of the values of c and interpolating a solution for a given value of c using the precomputed values. It isn't elegant at all but it's something.

-
I assume the bs are all the same then? – Jiahao Chen Aug 26 '10 at 11:55
Exactly how big are your matrices anyway? This can crucially affect the practicality of proposed solutions. – J. M. Aug 26 '10 at 12:33
$10^6$... you've essentially said none of the proposals thus far are practical. ;) So what you in fact want is the vector $(E+cI)^{-1}b$ with varying c; I'll check my notes on special methods for sparse matrices and report back. Only one thing: definitely QR decomposition will not be a practical solution to your problem either! – J. M. Aug 26 '10 at 12:51

You want the resolvent of $E$ (at $z:=-c$). Recall it's an analytic function of $z$ defined on the resolvent set, $\mathbb{C}\setminus\operatorname{spec}(E)$. According to the complexity of the matrix $E$, and with the number and the location of the $c$ you need to consider, it may be worth computing a power series expansion at various centers so as to cover the set $\{c\}$ of the data. For $|z|$ larger than the spectral radius you have of course the Laurent expansion $(z-E)^{-1}=1/z+ E/z^2+ E^2/z^3+\dots$

-
Well, (E+cI)x=b can indeed be turned into something like (E/c+I)x=b/c, and thus you can use the geometric series technique. Of course, if the spectral radius of E/c is bigger than 1, this idea is shot. – J. M. Aug 26 '10 at 12:29
Is there a method to use the resolvent without computing it explicitly? It seems to use the resolvent would have to be recalculated for all cs, and explicit computation could result in massive fill-in (loss of sparsity) as @J. Mangaldan pointed out above once products like E^2 are computed. – Jiahao Chen Aug 26 '10 at 12:50
Well, one can do it Krylov-style, assuming there is a nice black box for matrix-vector multiplications. Maintain a vector v initialized to the right-hand side, and at every iteration multiply this by E. But again, this is only feasible if it can be assured that the spectral radius of E/c never exceeds unity. – J. M. Aug 26 '10 at 13:00
Yes, as I wrote. Nevertheless, finitely many expansions suffice to cover the set of {c}; whether this approach is efficient depends on the details. – Pietro Majer Aug 26 '10 at 14:22

a) There are formulae such as the Woodbury identity that allow for rank k updates to a previously solved problem, which I think fits your problem nicely.

b) In addition, using a reasonably smart iterative algorithm such as conjugate gradients (or whatever is appropriate for your problem) will also be helpful since you can feed it the solution from your previous problem, and for small perturbations the new solution can be computed very quickly.

In practice I have found it sufficient to use just (b), but it might be worth trying both separately or together.

-
(a) does not look useful since the updates are full-rank (b) seems useful only if the $c$s are close. Or am I missing something? – Federico Poloni Aug 26 '10 at 12:14
Sherman-Morrison-Woodbury, as already stated, is of no help since it is intended for low-rank corrections; e.g., corrections expressible in the form $UV^T$, where $U$ and $V$ are rectangular. CG is of no use since he said in the problem statement "I should say that the resulting matrices are not necessarily Hermitian", and computing $(E+cI)^T(E+cI)$ so that you can apply CG can result in a much denser matrix. – J. M. Aug 26 '10 at 12:24
@Frederico: That's true, I forgot about fill-in, since using Woodbury would require the computation of E^-1 . E^-1. I am assuming that the c's are small since the OP did say diagonal perturbations. – Jiahao Chen Aug 26 '10 at 12:54
@J. Mangaldan: that is true, CG is not guaranteed to work. Sometimes it does anyway though! But really without knowing more about the problem it is difficult to recommend any particular iterative solver. Tricks like level shifting (regularizing the matrix A to make it always posdef) can sometimes be useful in using CG. BCG certainly doesn't seem feasible due to fillin. My first thought was regularized CG or GMRES. – Jiahao Chen Aug 26 '10 at 12:54
A method that sometimes works and sometimes doesn't, IMHO, cannot be recommended in good conscience to somebody whose problem's properties still have to be explored fully. Besides, he said c can be complex, and that can play havoc with regularization. – J. M. Aug 26 '10 at 12:57

If you're doing a full LU decomposition and ignoring sparsity, then you could switch to a Schur decomposition (costs $25n^3$ instead of $2/3n^3$, but allows you to solve any of the resulting systems within $O(n^2)$). If you're using sparsity, as far as I know it is an open research problem how to exploit fully this property (see e.g. the rational Krylov method).

-
E might be sparse, the Schur decomposition is definitely dense, and computational time and storage can be prohibitive. I'm assuming that E is a large enough matrix that the "LU decomposition" alluded to in the original post is in fact set to do something like the Cuthill-McKee ordering to maintain sparsity in the triangular factors. – J. M. Aug 26 '10 at 12:32
Unfortunately, I do need to maintain sparsity or the problem becomes too big to handle. (I am using the KLU package, but I expect that most LU solvers would be able to handle the linear systems I have.) – Fumiyo Eda Aug 26 '10 at 12:47
As I said, the reason why those LU solvers can manage your large matrices is that they do a preliminary analysis of sparsity pattern before performing the LU decomposition. Blindly triangularizing or pivoting can result in disastrous fill-in, and thus pattern analysis is a crucial step for these LU solvers. – J. M. Aug 26 '10 at 12:54

Here's an idea, perhaps naive? Introduce auxiliary variable $y = x \sqrt{c}$. $$\left( \begin{array}{cc} E & \sqrt{c} I \\ \sqrt{c} I & -I \\ \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} b \\ 0 \end{array} \right).$$ Now an LU decomposition of the larger matrix is $$\left( \begin{array}{cc} L_E & 0 \\ L_{10} & L_{11} \end{array} \right) \left( \begin{array}{cc} U_E & U_{01} \\ 0 & U_{11} \end{array} \right),$$ where $E = L_E U_E$. When $c$ changes to $d$, the new LU decomposition is given by $$\left( \begin{array}{cc} L_E & 0 \\ \sqrt{\frac{d}{c}} L_{10} & \tilde L_{11} \end{array} \right) \left( \begin{array}{cc} U_E & \sqrt{\frac{d}{c}} U_{01} \\ 0 & \tilde U_{11} \end{array} \right)$$ where $$\frac{d}{c} L_{10} U_{01} + \tilde L_{11} \tilde U_{11} = -I.$$ But we know $L_{10} U_{01} + L_{11} U_{11} = -I$, so substituting $$\tilde L_{11} \tilde U_{11} = -I + \frac{d}{c} (I + L_{11} U_{11}).$$ Find a solution to that and then LU solve (oops ... maybe this is as hard as the original problem?)

-

Here's another silly idea. Given the scale of the problem, this is perhaps a ridiculous suggestion, but: if you have the SVD of $E = U \Sigma V^\top$, then the SVD of $E + c I = U (\Sigma + c I) V^\top$. So $x (c) = V (\Sigma + c I)^{-1} U^\top b$. (This is what is done for computing ridge regression regularization paths.)

-