Stopping criteria for iterative linear solvers applied to nearly singular systems

Question

Consider $Ax=b$ with $A$ nearly singular which means there is an eigenvalue $\lambda_0$ of $A$ that is very small. The usual stop criterion of an iterative method is based on the residual $r_n:=b-Ax_n$ and regards the iterations can stop when $\|r_n\|/\|r_0\|<tol$ with $n$ the iteration number. But in the case we are considering, there could be large error $v$ living in the eigenspace associated with the small eigenvalue $\lambda_0$ which gives small residual $Av=\lambda_0v$. Suppose the initial residual $r_0$ is large, then it might happen we stop at $\|r_n\|/\|r_0\|<tol$ but the error $x_n-x$ is still big. What is a better error indicator in this case? Is $\|x_{n}-x_{n-1}\|$ a good candidate?

Answer 1

Please never use the difference between successive iterates to define a stopping criteria. This misdiagnoses stagnation for convergence. Most nonsymmetric matrix iterations are not monotone, and even GMRES in exact arithmetic with no restarts may stagnate for an arbitrary number of iterations (up to the dimension of the matrix) before converging suddenly. See examples in Nachtigal, Reddy, and Trefethen (1993).

A better way to define convergence

We are usually interested in the accuracy of our solution more than the size of the residual. Specifically, we might like to guarantee that the difference between an approximate solution $x_n$ and the exact solution $x$ satisfies $$|x_n - x| < c$$ for some user-specified $c$. It turns out that can achieve this by finding an $x_n$ such that $$|A x_n - b| < c\epsilon$$ where $\epsilon$ is the smallest singular value of $A$, due to

$$\begin{align} |x_n - x| &= |A^{-1} A (x_n - x)| \\ & \le \frac 1 \epsilon |A x_n - A x| \\ & = \frac 1 \epsilon |A x_n - b| \\ & < \frac 1 \epsilon \cdot c \epsilon = c \end{align}$$

where we have used that $1/\epsilon$ is the largest singular value of $A^{-1}$ (second line) and that $x$ exactly solves $A x = b$ (third line).

Estimating the smallest singular value $\epsilon$

An accurate estimate of the smallest singular value is usually not directly available from the problem, but it can be estimated as a byproduct of a conjugate gradient or GMRES iteration. Note that although estimates of the largest eigenvalues and singular values is usually quite good after only a few iterations, an accurate estimate of the smallest eigen/singular value $\epsilon$ is usually only obtained once convergence is reached. Before convergence, the estimate will generally be significantly larger than the true value. This suggests that you must actually solve the equations before you can define the correct tolerance $c\epsilon$. An automatic convergence tolerance that takes a user-provided accuracy $c$ for the solution and estimates the smallest singular value $\epsilon$ with the current state of the Krylov method might converge too early because the estimate of $\epsilon$ was much larger than the true value.

Notes

1. The above discussion also works with $A$ replaced by the left-preconditioned operator $P^{-1}A$ and the preconditioned residual $P^{-1} (A x^n - b)$ or with the right-preconditioned operator $A P^{-1}$ and the error $P (x_n - x)$. If $P^{-1}$ is a good preconditioner, the preconditioned operator will be well-conditioned. For left-preconditioning, this means the preconditioned residual can be made small, but the true residual may not be. For right preconditioning, $|P(x_n - x)|$ is easily made small, but the true error $|x_n-x|$ may not be. This explains why left-preconditioning is better for making error small while right-preconditioning is better for making the residual small (and for debugging unstable preconditioners).
2. See this answer for more discussion of norms minimized by GMRES and CG.
3. The estimates of extremal singular values can be monitored using -ksp_monitor_singular_value with any PETSc program. See KSPComputeExtremeSingularValues() to compute singular values from code.
4. When using GMRES to estimate singular values, it is crucial that restarts not be used (e.g. -ksp_gmres_restart 1000 in PETSc).

Answer 2

Another way of looking at this problem is to consider the tools from discrete inverse problems, that is, problems which involve solving $Ax=b$ or $\min ||Ax-b||_2$ where $A$ is very ill-conditioned (i.e. the ratio between the first and last singular value $\sigma_1/\sigma_n$ is large).

Here, we have several methods for choosing the stopping criterion, and for an iterative method, I would recommend the L-curve criterion since it only involves quantities that are available already (DISCLAIMER: My advisor pioneered this method, so I am definitely biased towards it). I have used this with success in an iterative method.

The idea is to monitor the residual norm $\rho_k=||Ax_k-b||_2$ and the solution norm $\eta_k=||x_k||_2$, where $x_k$ is the $k$'th iterate. As you iterate, this begins to draw the shape of an L in a loglog(rho,eta) plot, and the point at the corner of that L is the optimal choice.

This allows you to implement a criterion where you keep an eye on when you have passed the corner (i.e. looking at the gradient of $(\rho_k,\eta_k)$), and then choose the iterate that was located at the corner.

The way I did it involved storing the last 20 iterates, and if the gradient $abs(\frac{\log(\eta_k)-\log(\eta_{k-1})}{\log(\rho_k)-\log(\rho_{k-1})})$ was larger than some threshold for 20 successive iterations, I knew that I was on the vertical part of the curve and that I had passed the corner. I then took the first iterate in my array (i.e. the one 20 iterations ago) as my solution.

There are also more detailed methods for finding the corner, and these work better but require storing a significant number of iterates. Play around with it a bit. If you are in matlab, you can use the toolbox Regularization Tools, which implements some of this (specifically the "corner" function is applicable).

Note that this approach is particularly suitable for large-scale problems, since the extra computing time involved is minuscule.