Why are systems with clustered eigenvalues easy to solve?

Question

I came across the following slide by Theo Diamandis & Zachary Frangella on what makes the linear system $Ax=b$ easy to solve using the conjugate gradient method.

---

---

Transcription:

CG converges quickly when

1. The condition number of $A$ is small
2. The eigenvalues of $A$ are clustered

Why is this the case?

I suspect this may be a well-known phenomenon, so any textbook pointers are appreciated.

Answer 1

A good explanation of this phenomena with many examples is given in Iterative Methods for Linear and Nonlinear Equations by Tim Kelley. The crux of it comes down to the fact that each step of a Krylov method is actually the minimizer of some subproblem involving approximating $b$ (in $Ax=b$) via applications of matrix polynomials of $A$ applied to $b$ and that it is easy to bound the application of matrix polynomials with clustered roots on the spectrum of $A$.

Specifically, for CG (this analysis is easiest due to normality), we have the following result

Let $A$ be spd and let $\{x_k\}$ be the conjugate gradient iterates. Let $k$ be fixed and given and let $p_k$ be any degree $k$ polynomial such that $p_k(0)=1$. Then $$\frac{\|x_k - x^*\|_A}{\|x_0 - x^*\|_A} \leq \max_{z\in\sigma(A)}|p_k(z)|.$$ This is true for any such polynomial. This means that if we can find a polynomial that makes this bound small, then CG iterations will converge at least as fast as that polynomial bound.

To illustrate with an example, suppose that $A$ contains only eigenvalues in the intervals $[1-\epsilon, 1+\epsilon]$ and $[10^6 - \epsilon, 10^6 + \epsilon]$. The condition number of this matrix is relatively large, but CG will converge quickly by noticing that the polynomial $$p_{2k}(z) = \frac{(1-z)^k(10^6 - z)^k}{10^{6k}}$$ satisfies $p_{2k}(0) = 1$ and satisfies $$\max_{z\in\sigma(A)} |p_{2k}(z)| \leq \epsilon^k\left(\frac{10^6-1 + \epsilon}{10^6}\right)^k \leq \epsilon^k,$$ which is obtained by substituting $z = 1-\epsilon$. This bounds the relative reduction in error in $2k$ steps, so eigenvalue clustering (small $\epsilon$) will result in rapid convergence.

Kelley, C. T., Iterative methods for linear and nonlinear equations, Frontiers in Applied Mathematics. 16. Philadelphia, PA: SIAM, Society for Industrial and Applied Mathematics. xiii, 165 p. (1995). ZBL0832.65046.

Answer 2

At the $k$th iteration, typical Krylov methods for solving $Ax=b$ (such as CG, MINRES, and GMRES) implicitly construct a $k$th order polynomial $Q(x)$ such that:

• $Q(0) = 1$.
• $|Q(\lambda_i)|$ is as small as possible for all eigenvalues $\lambda_i$ of $A$.

The meaning of "as small as possible" depends on the exact method, and also on the right hand side $b$. However, in general smaller values of $|Q(\lambda_i)|$ mean less error for the $k$th iterate, and larger values of mean more error.

Picture the eigenvalues as dots on the $x$-axis. If you can imagine a low order polynomial that is zero or very small at each of these dots, then the Krylov method will perform well. But if making the polynomial small on some dots forces it to be large at other dots, then the Krylov method will perform poorly.

From this picture it is clear why Krylov methods perform well when the eigenvalues are clustered: if the polynomial is small for one eigenvalue in the cluster, then it will also be small for all other points in the cluster for free.

Here is a picture illustrating the idea: This picture is from Section 3.3 of my dissertation (https://repositories.lib.utexas.edu/bitstream/handle/2152/75559/ALGER-DISSERTATION-2019.pdf), where I discuss this in more detail.

A good reference is Chapter 6 of the following book:

Saad, Yousef. Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics, 2003. https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=1e76c1586c460787b06cae3247a34c94600975e0

Another good reference is Section 9 of the following paper:

Shewchuk, J. R. (1994). An introduction to the conjugate gradient method without the agonizing pain. https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf

---

Edit 1: (comment on non-symmetry)

If the matrix $A$ is non-symmetric, then then the eigenvalues do still matter, but they do not tell the whole story; the conditioning of the eigenvector matrix matters as well. In particular, suppose that $A$ has eigenvalue decomposition $A=X \Lambda X^{-1}$. Convergence bounds are typically worsened by a factor of $\operatorname {cond} \left(X\right) := ||X||_2 ||X^{-1}||_2$.

For example, in Proposition 6.32 of Saad's book the following bound for GMRES convergence is proven: $$||r_k||_2 \le \operatorname{cond}\left(X\right) ||r_0||_2 \min_Q \max_i |Q(\lambda_i)|$$ where $r_k$ is the residual at the $k$th iteration, the minimization is taken over all $k$th order polynomials $Q$ satisfying $Q(0)=1$, and the maximization is taken over all eigenvalues.

---

Edit 2: (proof sketch)

Here is a sketch of the proof for the convergence bound in Edit 1, adapted from appendix D in my dissertation. Similar proofs can be found in all the references mentioned.

To see why the convergence bound in Edit 1 holds, notice that elements of the Krylov subspace $$\mathcal{K}_k = \operatorname{span}\{b, Ab, A^2 b, \dots, A^{k-1} b\}$$ are in bijective correspondence with the set of all $(k-1)$th order polynomials, $\mathcal{P}_{k-1}$ (including polynomials $P$ with $P(0) \neq 1$) via the map $$P \leftrightarrow P(A) b.$$ Starting from the minimal residual condition that defines GMRES, and using this bijection, we have: \begin{align*} ||b - A x_k|| &= \min_{x \in \mathcal{K}_k} ||b - A x|| \\ &= \min_{P \in \mathcal{P}_{k-1}} ||b - A P(A) b|| \\ &= \min_{Q \in \mathcal{Q}_j} ||Q(A)b||, \end{align*} where $\mathcal{Q}_k$ is the set of all $k$th order polynomials satisfying $Q(0)=1$.

The desired convergence bound now follows from expressing $Q(A)$ in the basis of eigenvectors of $A$ so that it becomes diagonal, i.e., $$Q(A)=X Q(\Lambda) X^{-1},$$ then using the submultiplicative property of the induced 2-norm.