solving series of linear systems with diagonal perturbations

Question

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.

Answer 1

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$

Answer 2

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.

Answer 3

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

Answer 4

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?)

Answer 5

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