Why are ill-conditioned systems of equations hard to solve iteratively?

Question

Is there some intuition as to why ill conditioned system of equations hard to solve iteratively ( i.e. the convergence is slow) ? I've read convergence proofs of several methods, but still don't have any real intuition.

Answer 1

Assume for simplicity that you want to solve a system $Ax=b$ where $A$ is Hermitian and positive definite (HDP) and hence there is a unitary $U$ and real diagonal $D$ (with the positive eigenvalues $\lambda_1\geq\ldots\geq\lambda_n>0$ on the diagonal) such that $A=VDV^*$.

Many methods including the Richardson's method, the conjugate gradient method, GMRES, etc. start from a given initial guess $x_0$ and seek for the approximation $x_k$ at step $k$ in the shifted Krylov subspace $$x_0+\mathcal{K}_k(A,b):=x_0+\mathrm{span}\{r_0,Ar_0,\ldots,A^{k-1}r_0\},$$ where $r_0:=b-Ax_0$ is the residual vector of $x_0$. Consequently, $x_k$ can be written in the form $$ x_k=x_0+q(A)r_0, $$ where $q$ is a polynomial of the degree (at most) $k-1$. Let's have a look at the error $e_k:=x-x_k$. We have $$\tag{1} e_k=x-x_k=x-x_0-q(A)r_0=e_0-q(A)Ae_0=[I-q(A)A]e_0=p(A)e_0, $$ where the polynomial $p$ of the degree (at most) $k$ satisfies $p(0)=1$. This characterizes all the so-called polynomial methods. The error at step $k$ is given by a polynomial in $A$ of the degree at most $k$ with normalized at origin (having the unit constant coefficient). For example, in the Richardson's method $x_k=x_{k-1}+\omega(b-Ax_{k-1})$, the polynomial $p$ is given by $p(t)=(1-\omega t)^k$.

The convergence analysis of polynomial methods is usually based on analysing the worst case convergence rate which is obtained by eliminating the effects of the right-hand side (or, equivalently, the initial error). Say, we are interested in the energy norm of the error $\|e_k\|_A:=(e_k^*Ae_k)^{1/2}$. We have $$ \|e_k\|_A=\|p(A)e_0\|_A=\|A^{1/2}p(A)e_0\|_2=\|p(A)A^{1/2}e_0\|_2\leq\|p(A)\|_2\|A^{1/2}e_0\|_2=\|p(A)\|_2\|e_0\|_A. $$ So at step $k$, the relative error can be bounded from above by $$\tag{2} \frac{\|e_k\|_A}{\|e_0\|_A}\leq \|p(A)\|_2=\|p(D)\|_2=\max_{1\leq i\leq n}|p(\lambda_i)|. $$

Now comes the issue with the ill-conditioned problems:

As indicated by (2), you have to find a polynomial $p$ of a certain degree fixated such that $p(0)=1$ which is "small" on the whole spectrum of $A$!

This is certainly not the case for the Richardson method. If the parameter $\omega$ is chosen optimally, that is, $\omega=2/(\lambda_1+\lambda_n)$, the polynomial $p$ is given by $$ p(t)=\left(1-\frac{2t}{\lambda_1+\lambda_n}\right)^k. $$ You can easily verify, that $p$ is small in the middle part of the spectrum but fairly large on its boundaries, in fact, $$\tag{3} p(\lambda_1)=p(\lambda_n)=\left(\frac{\lambda_1-\lambda_n}{\lambda_1+\lambda_n}\right)^k =\left(\frac{\kappa-1}{\kappa+1}\right)^k, $$ where $\kappa:=\lambda_1/\lambda_n$ is the spectral condition number of $A$. That is, the Richardson method eliminates quickly the components of $e_0$ lying in the middle of the spectrum but slowly the components combined from the eigenvectors corresponding to small and large eigenvalues (see the attached figure, where $\lambda_1=10$, $\lambda_n=0.1$). The bound on the convergence rate of the Richardson's method following from (2) and (3) is of course the worst-case bound. If by any luck your initial error were combined from the eigenvectors lying in the middle part of the spectrum, the real convergence rate would be much better.

The Richardson's method is not very optimal one. The best polynomial method for HPD problems which minimizes the $A$-norm of the error at each step is the conjugate gradient method (CG). If we denote by $\Pi_k$ the set of $k$ degree polynomials normalized at origin, the error $e_k$ of CG satisfies $$ \|e_k\|_A=\min_{p\in\Pi_k}\|p(A)e_0\|_A. $$ Again, we can obtain the worst-case convergence bound $$ \frac{\|e_k\|_A}{\|e_0\|_A}\leq\min_{p\in\Pi_k}\|p(A)\|_2=\min_{p\in\Pi_k}\max_{1\leq i\leq n}|p(\lambda_i)|. $$ Consequently, CG faces the real polynomial interpolation problem: find $p\in\Pi_k$ which attains the smallest possible values at all the eigenvalues of $A$.

This does not necessarily must be (at least theoretically) be a problem, e.g., when the spectrum is somewhat clustered. Imagine an "ideal" case, where $A$ has only few (say, $m\ll n$) distinct eigenvalues so that the minimum polynomial of $A$ has degree $m$. Since we can normalize it so that it attains the value $1$ at origin, we have then that CG converges in at most $m$ steps as $p_m(A)=0$. A useful, but sometimes misleading, idea is hence that if the spectrum of $A$ consists of a small number of tight clusters of eigenvalues, then CG converges fast. The worst convergence, on the other hand, can be expected when the spectrum is equidistant.

I hope this somewhat helps with what you were asking :-) I would suggest you to have a look on this book for some interesting overview on the topic.

---