# Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation?

Is there an $O(n^3+n^2 k)$ method to solve $k$ linear systems of the form $(D_i + A) x_i = b_i$ where $A$ is a fixed SPD matrix and $D_i$ are positive diagonal matrices?

For example, if each $D_i$ is scalar, it suffices to compute the SVD of $A$. However, this breaks down for general $D$ due to lack of commutativity.

Update: The answers so far are "no". Does anyone have any interesting intuition as to why? A no answer means that there's no nontrivial way to compress the information between two noncommuting operators. It isn't terribly surprisingly, but it'd be great to understand it better.

-
SPD = semi-positive definite? –  rcollyer Jan 2 '12 at 15:50
Yes, though the problem is essentially the same without SPD. I added that constraint only to ensure that the systems are never singular. –  Geoffrey Irving Jan 2 '12 at 19:40

The closest positive answers to your question that I could find is for sparse diagonal perturbations (see below).

With that said, I do not know of any algorithms for the general case, though there is a generalization of the technique you mentioned for scalar shifts from SPD matrices to all square matrices:

Given any square matrix $A$, there exists a Schur decomposition $A=U T U^H$, where $U$ is unitary and $T$ is upper triangular, and $A+\sigma I = U (T + \sigma I) U^H$ provides a Schur decomposition of $A + \sigma I$. Thus, your precomputation idea extends to all square matrices through the algorithm:

• Compute $[U,T]=\mathrm{schur}(A)$ in at most $\mathcal{O}(n^3)$ work.
• Solve each $(A+\sigma I) x = b$ via $x := U (T +\sigma I)^{-1} U^H b$ in $\mathcal{O}(n^2)$ work (the middle inversion is simply back substitution).

This line of reasoning reduces to the approach you mentioned when $A$ is SPD since the Schur decomposition reduces to an EVD for normal matrices, and the EVD coincides with the SVD for Hermitian positive-definite matrices.

Response to update: Until I have a proof, which I do not, I refuse to claim that the answer is "no". However, I can give some insights as to why it's hard, as well as another subcase where the answer is yes.

The essential difficulty is that, even though the update is diagonal, it is still in general full rank, so the primary tool for updating an inverse, the Sherman-Morrison-Woodbury formula, does not appear to help. Even though the scalar shift case is also full rank, it is an extremely special case since it commutes with every matrix, as you mentioned.

With that said, if each $D$ was sparse, i.e., they each had $\mathcal{O}(1)$ nonzeros, then the Sherman-Morrison-Woodbury formula yields an $\mathcal{O}(n^2)$ solve with each pair $\{D,b\}$. For example, with a single nonzero at the $j$th diagonal entry, so that $D=\delta e_j e_j^H$:

$$[A^{-1}+\delta e_j e_j^H]^{-1} = A^{-1} - \frac{\delta A^{-1} e_j e_j^H A^{-1}}{1+\delta (e_j^H A^{-1} e_j)},$$

where $e_j$ is the $j$th standard basis vector.

Another update: I should mention that I tried the $A^{-1}$ preconditioner that @GeoffOxberry suggested on a few random SPD $1000 \times 1000$ matrices using PCG and, perhaps not surprisingly, it seems to greatly reduce the number of iterations when $||D||_2/||A||_2$ is small, but not when it is $\mathcal{O}(1)$ or greater.

-

If $(D_{i} + A)$ is diagonally dominant for each $i$, then recent work by Koutis, Miller, and Peng (see Koutis' website for work on symmetric diagonally dominant matrices) could be used to solve each system in $\mathcal{O}(n^2 \log(n))$ time (actually $\mathcal{O}(m\log(n))$ time, where $m$ is the maximum number of nonzero entries in $(D_{i} + A)$ over all $i$, so you could take advantage of sparsity as well). Then, the total running time would be $\mathcal{O}(n^2 \log(n) k)$, which is better than the $\mathcal{O}(n^3 k)$ approach of solving each system naïvely using dense linear algebra, but slightly worse than the quadratic run time you're asking for.

Significant sparsity in $(D_{i} + A)$ for all $i$ could be exploited by sparse solvers to yield an $\mathcal{O}(n^2 k)$ algorithm, but I'm guessing that if you had significant sparsity, then you would have mentioned it.

You could also use $A^{-1}$ as a preconditioner to solve each system using iterative methods, and see how that works out.

Response to update: @JackPaulson makes a great point from the standpoint of numerical linear algebra and algorithms. I will focus on computational complexity arguments instead.

The computational complexity of the solution of linear systems and the computational complexity of matrix multiplication are essentially equal. (See Algebraic Complexity Theory.) If you could find an algorithm that could compress the information between two non-commuting operators (ignoring the positive semidefinite part) and directly solve the collection of systems you're proposing in quadratic time in $n$, then it's likely that you could use such an algorithm to make inferences about faster matrix multiplication. It's difficult to see how positive semidefinite structure could be used in a dense, direct method for linear systems to decrease its computational complexity.

Like @JackPaulson, I'm unwilling to say that the answer is "no" without a proof, but given the connections above, the problem is very difficult and of current research interest. The best you could do from an asymptotic standpoint without leveraging special structure is an improvement on the Coppersmith and Winograd algorithm, yielding a $\mathcal{O}(n^{\alpha}k)$ algorithm, where $\alpha \approx 2.375$. That algorithm would be difficult to code, and would likely be slow for small matrices, because the constant factor preceding the asymptotic estimate is probably huge relative to Gaussian elimination.

-
I have yet to see a concrete statement of where the crossover might be, but several reputable sources have stated that (implementation issues aside), Coppersmith-Winograd cannot beat standard methods for matrix sizes that will be able to fit in memory in the near future (a few decades). Given that the Linpack benchmark takes more than a day to run on the current top machines, it does not seem likely that Coppersmith-Winograd will ever be used in practice. Strassen is actually practical for large problems, though it is somewhat less numerical stable. –  Jed Brown Dec 25 '11 at 5:16
That does not surprise me. +1 for the implementation details. –  Geoff Oxberry Dec 25 '11 at 13:42

A first order Taylor expansion can be used to improve convergence over simple lagging. Suppose we have a preconditioner (or factors for a direct solve) available for $A+D$, and we want to use it for preconditioning $A$. We can compute

\begin{align} A^{-1} &= (A+D-D)^{-1} (A+D) (A+D)^{-1} \\ &= [(A+D)^{-1} (A+D-D)]^{-1} (A+D)^{-1} \\ &= [I - (A+D)^{-1} D]^{-1} (A+D)^{-1} \\ &\approx [I + (A+D)^{-1} D] (A+D)^{-1} \end{align}

where the Taylor expansion was used to write the last line. Application of this preconditioner requires two solves with $A+D$.

It works fairly well when the preconditioner is shifted away from 0 by a similar or larger amount than the operator we are trying to solve with (e.g. $D\gtrsim 0$). If the shift in the preconditioner is smaller ($D \lesssim \min \sigma(A)$), the preconditioned operator becomes indefinite.

If the shift in the preconditioner is much larger than in the operator, this method tends to produce a condition number about half that of preconditioning by the lagged operator (in the random tests I ran, it could be better or worse for a specific class of matrices). That factor of 2 in condition number gives a factor of $\sqrt 2$ in iteration count. If the iteration cost is dominated by the solves with $A+D$, then this is not a sufficient factor to justify the first order Taylor expansion. If matrix application is proportionately expensive (e.g. you only have an inexpensive-to-apply preconditioner for $A+D$), then this first order method may make sense.

-