# Stopping criteria for iterative linear solvers applied to nearly singular systems

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?

• You may want to think about your definition of "nearly singular". The matrix $I \cdot \epsilon$ (with $\epsilon\ll 1$ and $I$ the identity matrix) has a very small eigenvalue, but is as far from singular as any matrix could be. – David Ketcheson Jan 6 '12 at 11:05
• Also, $||r_n/r_0||$ seems like the wrong notation. $||r_n||/||r_0||$ is more typical, no? – Bill Barth Jan 6 '12 at 13:17
• Yes, you are right, Bill! I will correct this mistake. – Hui Zhang Jan 6 '12 at 13:57
• What about $\| b - Ax \| / \| b \|$? and what is your algorithm, exactly? – shuhalo Jan 6 '12 at 14:08
• Addendum: I think the following paper pretty much adresses the ill-conditioned systems you worry about, at least if you are using CG: Axelson, Kaporin: Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. DOI: 10.1002/nla.244 – shuhalo Jan 6 '12 at 14:12

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).
• ''also works with A replaced by a preconditioned operator'' - However, it then applies only to the preconditioned residual $P^{-1}r$ if $P^{-1}A$ is used, resp. to the preconditioned error $P^{-1}\delta x$ if $AP^{-1}$ is used. – Arnold Neumaier Jul 25 '12 at 13:08
• Good point, I edited my answer. Note that the right-preconditioned case gives you control of $P\delta x$, unwinding the preconditioner (applying $P^{-1}$) typically amplifies low-energy modes in the error. – Jed Brown Jul 25 '12 at 16:19

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.

• Thanks a lot! So in loglog(rho,eta) plot we begin from the right of the L curve and end at the top of L, is it? I just do not know the principle behind this criterion. Can you explain why it always behave like an L curve and why we choose the corner? – Hui Zhang Jan 11 '12 at 19:15
• You're welcome :-D. For an iterative method, we begin from right and end at top always. It behaves as an L due to the noise in the problem - the vertical part happens at $||Ax-b||_2=||e||_2$, where $e$ is the noise vector $b_{exact}=b+e$. For more analysis see Hansen, P. C., & O'Leary, D. P. (1993). The use of the L-curve in the regularization of discrete ill-posed problems. SIAM Journal on Scientific Computing, 14. Note that I just made a slight update to the post. – OscarB Jan 12 '12 at 12:13
• @HuiZhang: it isn't always an L. If the regularization is ambiguous it may be a double L, leading to two candidates for the solution, one with gross featurse better resolved, the other with certain details better resolved. (And of course, mor ecomplex shapes may appear.) – Arnold Neumaier Jul 25 '12 at 13:10
• Does the L-curve apply to ill-conditioned problems where there should be a unique solution? That is, I'm interested in problems Ax = b where b is known "exactly" and A is nearly singular but still technically invertible. It would seem to me that if you use something like GMRES the norm of your current guess of x doesn't change too much over time, especially after the first however many iterations. It seems to me that the vertical part of the L-curve occurs because there is no unique/valid solution in an ill-posed problem; would this vertical feature be present in all ill-conditioned problems? – nukeguy Jun 6 '16 at 14:40
• At one point, you will reach such a vertical line, typically because the numerical errors in your solution method result in ||Ax-b|| not decreasing. However, you are right that in such noise-free problems the curve does not always look like an L, meaning that you typically have a few corners to choose from and choosing one over the other can be hard. I believe that the paper I referenced in my comment above discusses noise-free scenarios briefly. – OscarB Jun 7 '16 at 7:19