Iterative solvers
Gauss-Seidel and Jacobi Iteration: Splitting a Matrix to Solve It
Instead of inverting a matrix, what if you just guessed a solution and let it improve itself, sweep after sweep, until it stopped changing? That is the entire idea behind Jacobi and Gauss-Seidel iteration: split the coefficient matrix A into an easy-to-invert piece and a leftover, then turn the equation Ax = b into a fixed-point map x ↦ Mx + c that you iterate. Convergence is governed by a single number — the spectral radius ρ(M) of the iteration matrix — and the iteration converges to the true solution for every starting vector if and only if ρ(M) < 1.
For a strictly diagonally dominant matrix both methods are guaranteed to converge, and Gauss-Seidel — which reuses updated components the instant they are computed — typically halves the number of iterations Jacobi needs. These were among the first practical algorithms for large linear systems (Gauss, 1823; Jacobi, 1845; Seidel, 1874), and their descendants still power finite-element and PDE solvers today.
- FieldNumerical linear algebra / iterative solvers
- First appearedGauss 1823 (letter to Gerling); Jacobi 1845; Seidel 1874
- Convergence criterionρ(M) < 1 (spectral radius of iteration matrix) — necessary and sufficient
- Key sufficient hypothesisA strictly diagonally dominant, or A symmetric positive definite (Gauss-Seidel)
- Iteration matrixJacobi: M = −D⁻¹(L+U); Gauss-Seidel: M = −(D+L)⁻¹U
- Generalizes toSOR, block iterations, and the general splitting A = M − N framework
Interactive visualization
Press play, or step through manually. The visualization is yours to drive — try it before reading on.
Watch the 60-second explainer
A condensed visual walkthrough — narrated, captioned, under a minute.
The precise setup and what it claims
Write the n×n system Ax = b with the matrix decomposed by position as A = D + L + U, where D is the diagonal, L the strictly lower-triangular part, and U the strictly upper-triangular part (assume every diagonal entry aii ≠ 0). A stationary iterative method picks a splitting A = M − N with M invertible, and iterates
- x(k+1) = M⁻¹N x(k) + M⁻¹b = G x(k) + c,
where G = M⁻¹N is the iteration matrix. Jacobi takes M = D, so componentwise xi(k+1) = (bi − ∑j≠i aij xj(k))/aii. Gauss-Seidel takes M = D + L, using the just-updated components: xi(k+1) = (bi − ∑j<i aij xj(k+1) − ∑j>i aij xj(k))/aii. The central theorem: the iteration converges to the unique solution x* = A⁻¹b for every initial guess x(0) if and only if the spectral radius ρ(G) = max|λ| over eigenvalues λ of G satisfies ρ(G) < 1.
The picture: error that shrinks in the worst direction
Subtract the fixed-point relation x* = Gx* + c from the update to get the clean error recursion e(k+1) = G e(k), where e(k) = x(k) − x*. So after k steps the error is exactly Gke(0). The question of convergence is therefore the question: does Gk → 0?
Geometrically, expand e(0) in eigenvectors of G. Each component is scaled by λk every sweep. A component along an eigenvector with |λ| < 1 decays geometrically; a component with |λ| ≥ 1 does not. The slowest-decaying mode — the eigenvalue of largest modulus, ρ(G) — controls the asymptotic rate: the error norm shrinks roughly like ρ(G)k, so the number of iterations to gain one decimal digit is about 1/(−log₁₀ ρ(G)). Splitting the matrix well means engineering G to have a small spectral radius.
The mechanism: why ρ(G) < 1 is exactly the right threshold
The proof has two halves, and both hinge on the spectral radius. Sufficiency (ρ < 1 ⟹ convergence): Gelfand's formula says ρ(G) = limk→∞ ‖Gk‖1/k. If ρ(G) < 1, pick ε with ρ(G) + ε < 1; then for large k, ‖Gk‖ ≤ (ρ(G)+ε)k → 0, so Gke(0) → 0 for every e(0). (Equivalently: there exists a matrix norm with ‖G‖ < 1, making G a contraction and the Banach fixed-point theorem applies.) Necessity (convergence ⟹ ρ < 1): if some eigenvalue λ has |λ| ≥ 1 with eigenvector v, then choosing e(0) = v gives Gkv = λkv, whose norm |λ|k‖v‖ does not tend to 0 — the iteration fails for that start. So ρ(G) < 1 is genuinely necessary and sufficient. A cheap sufficient shortcut: if any induced matrix norm gives ‖G‖ < 1 then ρ(G) ≤ ‖G‖ < 1; strict diagonal dominance makes the ∞-norm of the Jacobi G less than 1 in one line.
A worked 2×2 example (and why Gauss-Seidel wins)
Solve A = [[2, 1], [1, 2]], b = [3, 3]; the true solution is x* = (1, 1). Here D = diag(2,2). Jacobi iteration matrix GJ = −D⁻¹(L+U) = [[0, −½], [−½, 0]], eigenvalues ±½, so ρ(GJ) = ½. Starting from x(0) = (0, 0): (0,0) → (1.5, 1.5) → (0.75, 0.75) → (1.125, 1.125) → …, each error exactly half the last.
Gauss-Seidel: M = D+L = [[2,0],[1,2]], and GGS = −M⁻¹U = [[0, −½], [0, ¼]], eigenvalues 0 and ¼, so ρ(GGS) = ¼ = ρ(GJ)². From (0,0): x₁ = 3/2 = 1.5, then x₂ = (3 − 1.5)/2 = 0.75; next sweep x₁ = (3 − 0.75)/2 = 1.125, x₂ = (3 − 1.125)/2 = 0.9375. The error contracts by ¼ per sweep — twice the digits-per-step of Jacobi. This ρGS = ρJ² relation holds exactly for consistently ordered / tridiagonal matrices (Young, 1950).
Why the hypotheses matter — and where it breaks
Convergence is not automatic; the splitting must control ρ(G). Two clean sufficient conditions:
- Strict diagonal dominance (|aii| > ∑j≠i|aij| for all i) ⟹ both Jacobi and Gauss-Seidel converge. Drop it and both can diverge.
- Symmetric positive definite A ⟹ Gauss-Seidel (and SOR for 0 < ω < 2) converges — the Householder-John / Ostrowski-Reich theorem. Note Jacobi need not converge for SPD matrices.
The subtlety: the two methods are not comparable in general. There exist matrices where Jacobi converges but Gauss-Seidel diverges, and others the reverse — a fact that surprises students who assume Gauss-Seidel is strictly better. Counterexample flavor: take A SPD but not diagonally dominant; Gauss-Seidel converges (Reich), yet the Jacobi iteration matrix D⁻¹(L+U) can have spectral radius ≥ 1 and diverge. Neither diagonal dominance nor positive-definiteness is necessary; ρ(G) < 1 alone is. This ties the topic to the general convergence theory of matrix iterations and to the spectral theorem, since for symmetric A the analysis reduces to eigenvalues.
Why it matters: the ancestor of modern solvers
Direct methods like Gaussian elimination cost O(n³) and destroy the sparsity of a matrix — fatal for the sparse systems from discretizing PDEs, where n can be 10⁶ or more but each row has a handful of nonzeros. Jacobi and Gauss-Seidel cost only O(nnz) per sweep and never store more than the sparse A, so they were the first practical way to attack such systems and remain foundational.
Their real modern role is as building blocks: (1) Successive over-relaxation (SOR) accelerates Gauss-Seidel by extrapolating with a relaxation factor ω, and optimally chosen ω can turn ρ ≈ 1 − ch² into ρ ≈ 1 − c√h; (2) multigrid methods use Gauss-Seidel as a smoother that kills high-frequency error while a coarse grid handles the rest, achieving optimal O(n) complexity; (3) they serve as preconditioners for Krylov methods like conjugate gradients and GMRES. Understanding when ρ(G) < 1 — and how the splitting shapes it — is the gateway to all of numerical PDE solving.
| Aspect | Jacobi | Gauss-Seidel | SOR (relaxation ω) |
|---|---|---|---|
| Splitting M (invert this) | D (diagonal) | D + L (lower triangular) | (1/ω)D + L |
| Iteration matrix | −D⁻¹(L+U) | −(D+L)⁻¹U | (D+ωL)⁻¹((1−ω)D − ωU) |
| Uses new values same sweep? | No (fully parallel) | Yes (sequential) | Yes (sequential) |
| Guaranteed if A strictly diag. dominant | Yes | Yes | Yes for 0 < ω ≤ 1 |
| Guaranteed if A symmetric positive definite | Not in general | Yes | Yes for 0 < ω < 2 |
| Typical speed | Baseline | ≈ 2× faster (ρ_GS ≈ ρ_J² for many matrices) | Optimal ω can give order-of-magnitude speedup |
Frequently asked questions
What is the exact convergence condition for Jacobi and Gauss-Seidel?
The iteration x⁽ᵏ⁺¹⁾ = Gx⁽ᵏ⁾ + c converges to A⁻¹b for every starting vector if and only if the spectral radius ρ(G) < 1, where G is the iteration matrix (−D⁻¹(L+U) for Jacobi, −(D+L)⁻¹U for Gauss-Seidel). This condition is both necessary and sufficient. The asymptotic error contraction rate per step is ρ(G).
Why does strict diagonal dominance guarantee convergence?
If |aᵢᵢ| > ∑ⱼ≠ᵢ|aᵢⱼ| for every row, then the Jacobi iteration matrix G = −D⁻¹(L+U) has row-sums of absolute values less than 1, so its induced ∞-norm ‖G‖∞ < 1. Since ρ(G) ≤ ‖G‖∞ for any induced norm, ρ(G) < 1 and convergence follows. A parallel argument works for Gauss-Seidel. Diagonal dominance is sufficient but not necessary.
Is Gauss-Seidel always faster than Jacobi?
Not always. For an important class — symmetric consistently-ordered or tridiagonal matrices — Young's theory (1950) gives ρ_GS = ρ_J², so Gauss-Seidel roughly doubles the convergence rate. (For nonnegative Jacobi iteration matrices, the separate Stein-Rosenberg theorem gives a convergence dichotomy — Jacobi and Gauss-Seidel converge or diverge together, and when both converge Gauss-Seidel is faster — but it is not the squaring identity.) But in general the two are not comparable: there exist matrices where Jacobi converges and Gauss-Seidel diverges, and vice versa. Only ρ(G) < 1 decides convergence for either.
Does convergence depend on the initial guess x⁽⁰⁾?
No. When ρ(G) < 1 the iteration converges to the unique solution for every x⁽⁰⁾, because the error obeys e⁽ᵏ⁾ = Gᵏe⁽⁰⁾ and Gᵏ → 0 regardless of e⁽⁰⁾. A good initial guess only reduces the number of iterations, not whether it converges. Conversely, if ρ(G) ≥ 1, choosing x⁽⁰⁾ to excite the offending eigenvector makes the iteration fail.
What guarantees Gauss-Seidel convergence for symmetric positive definite matrices?
The Householder-John theorem (and the Ostrowski-Reich theorem for SOR) states that if A is symmetric positive definite, then Gauss-Seidel converges, and SOR converges for any relaxation factor 0 < ω < 2. The mechanism is a Lyapunov/energy argument: the A-norm of the error strictly decreases each sweep. Notably, Jacobi is not guaranteed to converge for SPD matrices.
How do these relate to SOR, multigrid, and Krylov methods?
All arise from the same splitting framework A = M − N. SOR generalizes Gauss-Seidel with M = (1/ω)D + L and an optimal ω can dramatically lower ρ(G). Multigrid uses Gauss-Seidel as a high-frequency error smoother to reach O(n) complexity. And both methods serve as preconditioners M⁻¹ that improve the conditioning of A for conjugate gradients or GMRES.