When to stop Gauss-Seidel-iterations?

Question

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);
    }
}

Answer 1

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.

Answer 2

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$.