Iterative Methods

Why Iterative Methods?

Direct methods — Gaussian elimination, LU decomposition, Cholesky — solve Ax = b exactly (up to rounding) in O(n³) operations. For n = 10,000, that is 10¹² floating-point operations: expensive, but feasible on a workstation. For n = 10⁶ (the scale of finite-element models or discretized PDEs), O(n³) requires 10¹⁸ operations — entirely out of reach. Worse, storing a dense 10⁶ × 10⁶ matrix requires 8 petabytes of memory.

The saving grace is that large-scale problems are almost always sparse: each row of A has only a handful of nonzero entries (typically O(1) or O(√n) in PDE problems). Iterative methods never form the full matrix; they only require matrix-vector products Av, which cost O(nnz) where nnz is the number of nonzeros — often O(n) for sparse problems. The price is that iterative methods return an approximate solution after a finite number of iterations, and convergence is not always guaranteed.

Jacobi Iteration

The simplest iterative method decomposes A = D + L + U, where D is the diagonal, L is strictly lower triangular, and U is strictly upper triangular. The Jacobi method rewrites each equation i so that xᵢ is expressed in terms of all other variables, then updates all components simultaneously using the previous iterate:

Jacobi Iteration

x_i^{(k+1)} = \frac{1}{a_{ii}}\!\left(b_i - \sum_{j \neq i} a_{ij}\, x_j^{(k)}\right)

Each component of x is updated using only the previous iterate x⁽ᵏ⁾. The iteration matrix is T_J = −D⁻¹(L + U). The iteration converges if and only if the spectral radius ρ(T_J) < 1. Diagonal dominance — |aᵢᵢ| > Σⱼ≠ᵢ |aᵢⱼ| for all i — guarantees ρ(T_J) < 1 and therefore convergence.

In matrix form, the Jacobi update is x⁽ᵏ⁺¹⁾ = D⁻¹(b − (L + U)x⁽ᵏ⁾) = D⁻¹b − D⁻¹(L + U)x⁽ᵏ⁾. The method requires no factorization: each iteration costs one sparse matrix-vector product and n diagonal divisions. It is trivially parallelizable because all components are updated independently — a significant advantage on modern multi-core and GPU hardware. The downside is slow convergence: the error at step k is proportional to ρ(T_J)^k, and for realistic problems ρ(T_J) may be close to 1, requiring thousands of iterations.

Convergence Analysis

Define the error at step k as e⁽ᵏ⁾ = x⁽ᵏ⁾ − x*, where x* is the true solution. Subtracting the fixed-point equation x* = T_J x* + c from the iteration, we get e⁽ᵏ⁺¹⁾ = T_J e⁽ᵏ⁾. Therefore ‖e⁽ᵏ⁾‖ ≤ ‖T_J‖^k ‖e⁽⁰⁾‖. The iteration matrix norm ‖T_J‖ must be less than 1 for convergence, but the asymptotic convergence rate is governed by the spectral radius ρ(T_J) = max|λᵢ(T_J)|. The number of iterations needed to reduce the error by a factor of ε is approximately log(ε) / log(ρ(T_J)). When ρ(T_J) = 0.99, reducing the error by 10⁶ requires about 1,380 iterations; when ρ(T_J) = 0.5, it only takes about 20.

Gauss-Seidel Iteration

A simple improvement over Jacobi: when updating component i, immediately use the already-updated components x₁⁽ᵏ⁺¹⁾, …, xᵢ₋₁⁽ᵏ⁺¹⁾ rather than the old values. This is the Gauss-Seidel method:

Gauss-Seidel Iteration

x_i^{(k+1)} = \frac{1}{a_{ii}}\!\left(b_i - \sum_{j < i} a_{ij}\, x_j^{(k+1)} - \sum_{j > i} a_{ij}\, x_j^{(k)}\right)

The Gauss-Seidel iteration matrix is T_GS = −(D + L)⁻¹U. For symmetric positive definite matrices, Gauss-Seidel always converges; for strictly diagonally dominant matrices, it converges at least as fast as Jacobi (often about twice as fast). Unlike Jacobi, Gauss-Seidel updates sequentially, so it is less amenable to parallelization.

In matrix form, the Gauss-Seidel update solves (D + L)x⁽ᵏ⁺¹⁾ = b − Ux⁽ᵏ⁾. The lower-triangular system (D + L)x⁽ᵏ⁺¹⁾ = b − Ux⁽ᵏ⁾ can be solved by forward substitution in O(nnz) operations. For the Poisson equation on an n×n grid (a standard PDE benchmark), the Gauss-Seidel spectral radius is ρ(T_GS) ≈ 1 − π²/n², so the number of iterations to converge scales as O(n²) — prohibitively slow for large n, which motivates multigrid and preconditioned conjugate gradient methods.

Successive Over-Relaxation (SOR)

SOR accelerates Gauss-Seidel by introducing a relaxation parameter ω. After computing the Gauss-Seidel update x̃ᵢ, the SOR update takes a weighted combination of the old value and the Gauss-Seidel result:

SOR Update

x_i^{(k+1)} = (1-\omega)\,x_i^{(k)} + \frac{\omega}{a_{ii}}\!\left(b_i - \sum_{j < i} a_{ij}\,x_j^{(k+1)} - \sum_{j > i} a_{ij}\,x_j^{(k)}\right)

When ω = 1, SOR reduces to Gauss-Seidel. For ω ∈ (1, 2), the method over-relaxes (extrapolates past the Gauss-Seidel update), which can dramatically accelerate convergence. For ω ∈ (0, 1), the method under-relaxes, which can stabilize an otherwise divergent iteration. The optimal ω for the Poisson equation on an n×n grid is ω* = 2/(1 + sin(π/n)), giving ρ(T_SOR) ≈ 1 − 2π/n — an O(n) reduction in iterations compared to Gauss-Seidel's O(n²).

Finding the optimal ω analytically requires knowledge of the spectral radius of the Gauss-Seidel iteration matrix, which is rarely available in general. In practice, ω is tuned experimentally or estimated via an adaptive procedure. SOR with the optimal ω reduces the iteration count from O(n²) to O(n) for Poisson problems — a factor of n improvement, often making the difference between hours and minutes. For the 100×100 grid, this means ~3,000 Gauss-Seidel iterations versus ~100 SOR iterations with optimal ω.

Convergence Conditions

A unified framework for analyzing stationary iterative methods (Jacobi, Gauss-Seidel, SOR) considers the splitting A = M − N, where x⁽ᵏ⁺¹⁾ = M⁻¹(b + Nx⁽ᵏ⁾). The iteration matrix is T = M⁻¹N = M⁻¹(M − A) = I − M⁻¹A. Key convergence results:

Necessary and sufficient: ρ(T) < 1 (spectral radius of iteration matrix strictly less than 1).
Sufficient — diagonal dominance: If |aᵢᵢ| > Σⱼ≠ᵢ |aᵢⱼ| for all i, both Jacobi and Gauss-Seidel converge.
Sufficient — symmetric positive definite: If A is SPD and 0 < ω < 2, then SOR converges. Gauss-Seidel (ω = 1) always converges for SPD matrices.
Divergence guaranteed: For SOR, ρ(T_SOR) ≥ |ω − 1| for any matrix, so ω must be in (0, 2) for any hope of convergence.
Comparison theorem: For M-matrices (special class including diagonally dominant matrices with nonpositive off-diagonal entries), ρ(T_GS) = ρ(T_J)² — Gauss-Seidel's spectral radius is the square of Jacobi's.

The Conjugate Gradient Method

For symmetric positive definite (SPD) matrices — which arise in Poisson equations, least-squares normal equations, covariance matrices, and many physical systems — the conjugate gradient (CG) method is far superior to the stationary methods above. CG is a Krylov subspace method: it builds an optimal approximation from the subspace spanned by {b, Ab, A²b, …, Aᵏb}, which captures the most important directions for the solution.

Conjugate Gradient Algorithm

\alpha_k = \frac{r_k^T r_k}{p_k^T A p_k},\quad x_{k+1} = x_k + \alpha_k p_k,\quad r_{k+1} = r_k - \alpha_k A p_k,\quad \beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k},\quad p_{k+1} = r_{k+1} + \beta_k p_k

CG maintains three vectors: x (solution), r = b − Ax (residual), and p (search direction). At each step, α is chosen to minimize the error along p; then β updates p to be A-conjugate (A-orthogonal) to the previous direction. The method terminates in at most n steps in exact arithmetic — it is a direct method in principle — but with floating-point arithmetic, one typically stops when ‖r‖/‖b‖ < tolerance.

The CG method has two remarkable optimality properties. First, xᵏ minimizes the A-norm of the error ‖e‖_A = (eᵀAe)^{1/2} over the Krylov subspace 𝒦ₖ(A, b). Second, the residuals rₖ are mutually orthogonal and the search directions pₖ are A-conjugate: pᵢᵀApⱼ = 0 for i ≠ j. These properties together ensure that CG never repeats work: each step adds genuinely new information. The convergence rate depends on the condition number:

CG Convergence Bound

\frac{\|e_k\|_A}{\|e_0\|_A} \leq 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^k

Here κ = κ₂(A) is the condition number. The error decays geometrically at rate (√κ − 1)/(√κ + 1). For κ = 100, the rate is 9/11 ≈ 0.82 — each iteration reduces the error by 18%. For κ = 10,000, the rate is 99/101 ≈ 0.98 — only 2% error reduction per iteration, requiring ~300 iterations for 6-digit accuracy. Preconditioning replaces A with M⁻¹A to reduce the effective condition number.

Preconditioning

Preconditioning is the most powerful tool for accelerating Krylov methods. Instead of solving Ax = b directly, we solve the preconditioned system M⁻¹Ax = M⁻¹b (or the symmetric variant M⁻¹/²AM⁻¹/²y = M⁻¹/²b), where M is chosen so that M⁻¹A has a much smaller condition number than A, and M⁻¹ can be applied cheaply. Common preconditioners include:

Diagonal (Jacobi) preconditioning: M = D. Scales rows so diagonal entries are all 1. Cheap but often insufficient for highly ill-conditioned problems.
Incomplete LU (ILU): Compute an approximate LU factorization that keeps only the sparsity pattern of A (or a slightly denser pattern). Costs O(nnz) per apply. Very effective for PDE-based systems.
Incomplete Cholesky (IC): The SPD variant of ILU. Drop small entries from the Cholesky factor to control fill-in.
Algebraic multigrid (AMG): The gold standard for elliptic PDEs. Exploits the multi-scale structure of the problem to achieve nearly O(n) convergence regardless of mesh size.
SSOR preconditioning: Symmetric SOR applied as a preconditioner; often reduces κ(M⁻¹A) from O(n²) to O(n) for Laplacian-type problems.

The goal of preconditioning is to cluster the eigenvalues of M⁻¹A near 1, since CG converges in exactly k iterations when the preconditioned matrix has only k distinct eigenvalues. An ideal preconditioner is cheap to apply (O(n) operations), cheap to construct (O(nnz) preprocessing), and produces M⁻¹A with condition number close to 1. In practice, one balances these trade-offs: more aggressive preconditioning reduces iterations but increases per-iteration cost.

When Iterative Beats Direct Methods

The choice between direct and iterative solvers depends on problem size, structure, and the number of right-hand sides:

Large sparse systems (n > 10⁴–10⁵): Direct methods fill in the sparse structure during factorization, producing a dense triangular factor that may require O(n^{3/2}) or more storage for 2D problems. Iterative methods maintain sparsity throughout.
Many right-hand sides with the same A: Direct methods win — factor A once (O(n³) or O(n^{3/2}) for sparse), then solve each new b in O(n) or O(n log n) via forward/backward substitution.
Poorly conditioned matrices: Neither direct nor iterative methods can overcome fundamental ill-conditioning. For nearly singular matrices, direct methods detect this via small pivots; iterative methods diverge or converge extremely slowly.
Matrix-free problems: If A is defined implicitly (e.g., as a differential operator or a function evaluation), only iterative methods are possible — they only need matrix-vector products Av.
3D PDE problems: Direct methods fill even more aggressively in 3D (O(n²) fill-in), making iterative methods with multigrid preconditioning the only practical option.

A rule of thumb: for 1D and small 2D problems (n ≲ 10⁴), use direct methods for robustness. For large 2D problems (n ∼ 10⁵–10⁶), preconditioned CG or GMRES is standard. For 3D problems and any problem where n > 10⁶, multigrid or domain-decomposition methods with Krylov acceleration are the only viable approaches.

Engineering Implication

In scipy: scipy.sparse.linalg.cg(A, b) runs preconditioned CG; pass M=preconditioner for dramatic speedup. scipy.sparse.linalg.spsolve uses a sparse direct solver (UMFPACK). For PDE-based problems in Python, the pyamg library provides algebraic multigrid preconditioners. In MATLAB, the backslash operator automatically selects a sparse direct solver; pcg() is the preconditioned CG. When solving a 10⁶ × 10⁶ sparse Poisson system, CG with an incomplete Cholesky preconditioner typically converges in ~100 iterations — each costing O(n) — versus a direct sparse Cholesky factorization costing O(n^{3/2}) time and storage in 2D.

Key Takeaways

Iterative methods avoid O(n³) factorization by refining an approximate solution through matrix-vector products. Jacobi updates all components simultaneously using old values; Gauss-Seidel immediately uses updated components for faster convergence. SOR adds a relaxation parameter ω ∈ (0,2) that, at its optimal value, cuts iteration counts from O(n²) to O(n) for Poisson problems. Convergence requires the spectral radius of the iteration matrix to be less than 1; diagonal dominance and positive definiteness are key sufficient conditions. The conjugate gradient method is optimal for SPD matrices: it minimizes the A-norm error over the Krylov subspace in k ≤ n steps, converging at rate ((√κ−1)/(√κ+1))^k. Preconditioning — transforming A to M⁻¹A with smaller κ — is the primary tool for accelerating CG. Iterative methods dominate direct methods for large sparse systems, 3D PDE problems, and matrix-free settings.