Iterative linear solvers and preconditioning

A stationary residual correction applies a fixed approximate inverse M inverse to b minus A times the current guess. Subtract the exact solution to find a linear recurrence for the error. Convergence from every starting vector in finite dimensions occurs when the error propagation matrix has spectral radius below one. A bound in a convenient induced norm gives an easier sufficient condition, but failing that norm test does not automatically prove divergence. The chosen diagnostic must match the theorem being used.

Jacobi iteration chooses M as the diagonal of A and requires nonzero diagonal entries. Strict row diagonal dominance is one useful sufficient condition for convergence, but arbitrary nonsingular matrices need not satisfy it. Gauss–Seidel uses newly computed coordinates within each sweep, producing a different propagation matrix. Rearranging or damping updates can change convergence. None of these choices makes a small residual synonymous with small forward error: the original system's conditioning remains relevant at termination.

A Krylov method uses vectors generated by repeatedly applying A to the initial residual. The resulting subspace can capture useful directions without storing a dense factorization. Conjugate gradients exploits real symmetry and positive definiteness, minimizing the error in the A-energy norm over its affine Krylov space in exact arithmetic. Its textbook finite-termination property is an ideal exact-arithmetic fact, not a promise that a large floating-point run finishes with zero error after exactly n iterations.

GMRES instead minimizes the Euclidean residual over a growing Krylov subspace and applies to general nonsymmetric square systems. Building and orthogonalizing that basis consumes increasing memory and communication, so restarted variants discard part of the history and may stagnate. An algorithm choice follows from operator structure and available resources. Using conjugate gradients on an arbitrary nonsymmetric or indefinite problem invalidates its usual guarantees; changing the solver name in code does not change the mathematics of the matrix.

A preconditioner approximates the action of A well enough that applying its inverse makes the transformed problem easier, while remaining cheaper than solving the original system directly. Diagonal scaling is inexpensive; incomplete factorizations and multilevel methods can be stronger but need setup work and storage. For conjugate gradients, a symmetric positive-definite preconditioner is used with a compatible symmetric formulation. Merely left-multiplying an SPD matrix by an arbitrary inverse does not preserve ordinary Euclidean symmetry.

Measure total time: setup, sparse products, preconditioner applications, orthogonalization, communication, and the number of right-hand sides. A preconditioner that halves iterations but triples each step can lose overall. Monitor a meaningful residual in the original system, since a transformed residual has a different scale. In floating-point arithmetic, recursively updated residuals may drift from b−Ax, so occasional explicit recomputation can be necessary. Stop at an error level justified by the data and the outer algorithm instead of solving an inner system far more accurately than needed.