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?

  • 3
    $\begingroup$ 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. $\endgroup$ Jan 6, 2012 at 11:05
  • 1
    $\begingroup$ Also, $||r_n/r_0||$ seems like the wrong notation. $||r_n||/||r_0||$ is more typical, no? $\endgroup$
    – Bill Barth
    Jan 6, 2012 at 13:17
  • $\begingroup$ Yes, you are right, Bill! I will correct this mistake. $\endgroup$
    – Hui Zhang
    Jan 6, 2012 at 13:57
  • 1
    $\begingroup$ What about $\| b - Ax \| / \| b \|$? and what is your algorithm, exactly? $\endgroup$
    – shuhalo
    Jan 6, 2012 at 14:08
  • 2
    $\begingroup$ 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 $\endgroup$
    – shuhalo
    Jan 6, 2012 at 14:12

2 Answers 2


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.


  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).
  • 1
    $\begingroup$ ''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. $\endgroup$ Jul 25, 2012 at 13:08
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jed Brown
    Jul 25, 2012 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.

  • 1
    $\begingroup$ 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? $\endgroup$
    – Hui Zhang
    Jan 11, 2012 at 19:15
  • $\begingroup$ 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. $\endgroup$
    – OscarB
    Jan 12, 2012 at 12:13
  • 4
    $\begingroup$ @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.) $\endgroup$ Jul 25, 2012 at 13:10
  • $\begingroup$ 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? $\endgroup$
    – nukeguy
    Jun 6, 2016 at 14:40
  • $\begingroup$ 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. $\endgroup$
    – OscarB
    Jun 7, 2016 at 7:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.