# When to stop Gauss-Seidel-iterations?

I want to have an estimation, that my solution has an error, let's say less than 1e-8.

Usually, I stop the Gauss-Seidel algorithm, when the residual is "small enough" and this is already the problem. How do I know when the residual is small enough, because even when the residual is small, the solution may still have too much error. So this is no good method.

What do you use as stopping criterion?

On another website (math-linux.com) I found a stopping criterion: $$\|r\|/\|b\|\leq \epsilon$$ But again, what theory is behind that?

This, by the way, is the code I used in my last project, just for information, how I did it:

void relax(double epsilon, vector<double> &x, SparseMatrix &A, const vector<double> &f) {
int maxIter = 100;
int iter = 0;
double residual = 1.0;
double minResidual 0.000001; //I also tried 1e-14;

while (iter < maxIter && residual >= minResidual) {
for (int i = 0; i < A.dim; ++i) {
double ls = A.lineScalar(i, x);
x[i] = (1.0/A(i,i)) * (f[i] - ls);
}
vector<double> temp = A.multiply(x);
residual = L2(temp - f);
}
}

-
Seems to me pretty inefficient to test the residual every iteration this way. – Vladimir F Jul 25 '12 at 8:29
Yes, that's true. – vanCompute Jul 25 '12 at 8:35
I would like to pose another question for this topic: Can we estimate the lowest singular value of $A$ from below in a reasonable way provided the matrix has additional structure, like strict diagonal dominance? As far as I know, Kaasschieter has described such a procedure for the CGM. – Martin Jul 25 '12 at 13:13
@Martin: You should pose this as a separate question rather than a comment. – Arnold Neumaier Jul 25 '12 at 13:52

One cannot conclude from a residual how accurate the solution is. Between the best and the worst case in norm, there is a factor of exactly the condition number. More precisely, if the residual norm is r and the error norm is e then $\|A\|^{-1}\le e/r \le \|A^{-1}\|$, and both bounds are attainable. Taking the quotient of the bounds proves the claim.

The order of magnitude of the norm of the inverse can often be estimated from a few steps of the Lanczos/Arnoldi method. Thus one can apply the bounds.

However, generally, model equations are not really exact anyway, and getting the residual smaller than the accuracy with which uncertain data determine the residual is considered the scientifically corrrect procedure. Depending on the data, this may be far earlier that stagnation.

The theory behind this is called backward error analysis. For example, if the only source of uncertainty are measurement errors in the right hand side of 1e-4 in some norm, and your residuum norm is smaller than 1e-4, it cannot be distinguished from the residuum the exact solution of your solved problem would have when used as approximation in a problem with a different rhs of the same accuracy.

-
As one needs $\|A^{-1}\|$ this approach is not feasible, as one has to know the solution for calculation of the error. It may be better to look, wether the residuum stagnates? – vanCompute Jul 25 '12 at 12:45
The order of magnitude of the inverse norm can often be estimated from a few steps of the Lanczos/Arnoldi method. - The fact that one cannot conclude from a residual how accurate the solution is is independent of the knowledge of the value of the inverse norm. - Generally, model equations are not really exact anyway, and getting the residual smaller than the accuracy with which uncertain data determine the residual is considered the scientifically corrrect procedure. Depending on the data, this may be far earlier that stagnation. – Arnold Neumaier Jul 25 '12 at 13:03
Ok, if my matrix A comes from a PDE and the right hand side is a measurement with accuracy 1e-4 (for example), then I would stop ich the residuum is smaller then 1e-4? But what theory is behind that? – vanCompute Jul 25 '12 at 13:26
The theory behind it is called backward error analysis. If your residuum is smaller than 1e-4 it cannot be distinguished from the residuum the the exact solution of your particular problem would have when used as approximation in a problem with a different rhs of the same accuracy. – Arnold Neumaier Jul 25 '12 at 13:51
Ok, I will look it up that evening. Thank you for the hint. – vanCompute Jul 25 '12 at 14:24
Usually, when I'm estimating a solution of a system of linear equations, I save the approximation $\vec x^{n-1}$ and use it to compute $x_{err}=max|x_i^n-x_i^{n-1}|$ over each component i. This is faster than calculating the residual, since it doesn't require a matrix-vector product. Although, if your system is ill-conditioned, it still may not give you the correct solution for a given tolerance $\epsilon$.