# Why are systems with clustered eigenvalues easy to solve?

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.

• If you write your error in an eigenvector basis of $A$ then you can study how descent methods reduce the error w.r.t. different eigenvectors. For instance you can choose a descent direction along an eigenvector, and then after one (optimal) step you would have eliminated the error w.r.t. this eigenvector. It turns out that if you have that all eigenvalues are equal you can reduce the error in all of those to zero in a single step. The motivation in the clustered case is similar. See Shewchuk's write up on the conjugate gradient method - he discusses this with examples. Nov 2 at 13:00
• A related resource is Chapter 5 of Nocedal and Wright's Numerical optimization, in particular Figure 5.3 and Theorem 5.5 ("If A has only r distinct eigenvalues, then the CG iteration will terminate at the solution in at most r iterations.") Nov 3 at 10:27
• While not specific to conjugate gradient methods, consideration of condition number seems relevant.
– hardmath
Nov 5 at 2:04
• @hardmath would it be correct to say that "well-conditioned" implies eigenvalues are tightly clustered with one cluster, whereas two clusters implies the corresponding data matrix is a sum of two well-conditioned data matrices? Nov 5 at 17:04
• As the Wikipedia article I linked explains, the condition number is defined as a ratio of largest singular value to smallest singular value. For a real symmetric positive definite matrix, this amounts to dividing the largest eigenvalue by the smallest eigenvalue. So a low condition number in that case does imply eigenvalues are clustered. Unfortunately we often cannot expect matrices that arise in applications to be so nice.
– hardmath
Nov 5 at 18:04

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.

• I will leave a comment to note that I don't think this really explains "why" such phenomena happens in an intuitive way (@lightxbulb's comment does that better), but this is something concrete that you can point to that can be derived without much effort as to why such phenomena happens in a quantitative way Nov 2 at 20:29
• Interesting, thanks for the pointers! Trying to figure out how how small $\epsilon$ must be for your two-cluster bound to be reasonably close to reality. Also, does "epsilon is small" definition depend on the number of clusters? It seems it must, otherwise we could grow the number of clusters without limit and the problem would remain easy Nov 3 at 13:48
• The above inequality breaks when $\epsilon \geq 1$ but could be modified slightly to incorporate differently sized clusters. I highly recommend checking out the cited book. Notice that this bounds the $2k$-th iteration by $\epsilon^k$, so we have a per-step convergence rate of $\epsilon^{1/2}$. With 3 clusters, you can construct an order $3k$ polynomial with maximum $\approx \epsilon^k$, yielding a per-step convergence rate of $\epsilon^{1/3}$. Nov 3 at 15:07

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.

• It is worth mentioning that for nonsymmetric matrices, Krylov iterations can be slow even with clustered eigenvalues if the condition number of the eigenvector matrix is large. Nov 3 at 22:12
• @whpowell96 Yes of course; thank you for pointing that out. In that case the performance is typically made worse by a factor given by the condition number of the eigenvector matrix. For GMRES, this is proven in Proposition 6.32 in Saad's book. Nov 3 at 22:27
• I don't know why Anne Greenbaum's, Liesen's and Strakos' works are not well known as some others, but GMRES convergence behaviour is highly nonlinear and you can have a terrible convergence behaviour even for problems with "acceptable" condition numbers. ( e.g. Any nonincreasing convergence curve is possible for GMRES ). As a shameless plug, you can find an example in my thesis where the convergence is achieved at 501st iteration, with no drop in residual until that point (Preconditioning of Hybridizable Discontinuous Galerkin Discretizations of the Navier-Stokes Equations, Abdullah Ali Sivas) Nov 4 at 5:38
• How Fast are Nonsymmetric Matrix Iterations? by Nachtigal, Reddy, and Trefethen is also a good work on this topic that actually determines which characteristics of the matrices govern the convergence behaviors of different Krylov algorithms. Nov 5 at 16:16
• $\kappa_2(A^TA) = \kappa_2(A)^2$, so solving the normal equations can exascerbate conditioning issues. Additionally, Kyrlov methods are often used in applications where one doesn't explicitly store $A$ in memory, but instead just has a routime to compute the map $x\mapsto Ax$. Using a Krylov method on the normal equations requires one to have an additional routine to compute the action of $A^T$, which may be unavailable or expensive. Nov 5 at 22:19