Iterative solvers
Successive Over-Relaxation: Tuning ω for Faster Convergence
Take Gauss–Seidel, the classic iterative solver, and multiply each correction by a single number ω > 1 before applying it — a change so small it fits on one line of code. On the model Poisson problem with n interior points, that one tweak drops the number of iterations from O(n²) to O(n): a 5-point grid that needed thousands of Gauss–Seidel sweeps converges in tens. Successive Over-Relaxation (SOR) is the art of choosing that ω optimally.
Precisely: SOR splits the matrix A = D − L − U (diagonal, strict lower, strict upper parts) and iterates x⁽ᵏ⁺¹⁾ = (D − ωL)⁻¹[(1−ω)D + ωU] x⁽ᵏ⁾ + ω(D − ωL)⁻¹ b. Convergence is governed by the spectral radius ρ(L_ω) of that iteration matrix; the whole subject is about pushing ρ(L_ω) as far below 1 as the relaxation parameter ω ∈ (0, 2) allows.
- FieldNumerical linear algebra / iterative solvers
- IntroducedDavid Young and Stanley Frankel, independently, 1950
- Optimal ω (model problem)ω_opt = 2 / (1 + √(1 − μ²)), μ = ρ(Jacobi)
- Convergence conditionρ(L_ω) < 1; requires ω ∈ (0, 2) (Kahan)
- Key sufficient hypothesisA symmetric positive definite ⇒ converges for all ω ∈ (0,2) (Ostrowski–Reich)
- SpeedupO(n²) → O(n) iterations on the 2D Poisson model problem
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 statement: the SOR iteration and what optimality means
Let A ∈ ℝⁿˣⁿ have nonzero diagonal and write A = D − L − U, where D is the diagonal, −L the strictly lower-triangular part, and −U the strictly upper-triangular part. For a fixed relaxation parameter ω, the SOR iteration updates each component in place:
- x_i⁽ᵏ⁺¹⁾ = (1 − ω)·x_i⁽ᵏ⁾ + (ω / a_ii)·(b_i − ∑_{j
In matrix form x⁽ᵏ⁺¹⁾ = L_ω x⁽ᵏ⁾ + c_ω with the SOR iteration matrix L_ω = (D − ωL)⁻¹[(1 − ω)D + ωU]. The iteration converges for every starting vector iff the spectral radius ρ(L_ω) < 1. ω = 1 recovers Gauss–Seidel; ω < 1 is under-relaxation, ω > 1 over-relaxation. The optimization problem is: find ω_opt = argmin_{ω} ρ(L_ω), because the number of iterations to reach tolerance ε scales like log ε / log ρ(L_ω_opt).
The picture: nudging the correction past the target
Gauss–Seidel computes, for each unknown, the value that would exactly satisfy its equation given the current neighbors, then moves there. On smooth elliptic problems those corrections are systematically biased in one direction — the error is dominated by a slowly-decaying low-frequency mode, and each sweep only chips away at it. Over-relaxation says: since the correction keeps pointing the same way, overshoot it. Take the Gauss–Seidel step and extrapolate past it by a factor ω.
Geometrically, if Δ = x_GS − x_old is the Gauss–Seidel displacement, SOR moves to x_old + ωΔ. Too timid (ω near 1) and the low-frequency error crawls; too aggressive (ω near 2) and you overshoot into oscillation and eigenvalues climb back toward 1. There is a sweet spot where the largest eigenvalue of L_ω is pushed to its minimum modulus — and at exactly that ω the previously real, slowly-decaying eigenvalues collide and become a complex-conjugate pair of equal, minimal modulus. That collision is the mechanism behind the dramatic speedup.
The key idea: Young's eigenvalue relation for consistently ordered matrices
The clean theory (David Young, 1950) needs a structural hypothesis: A has property A — a red/black 2-coloring so the graph is bipartite, as finite-difference stencils do — and is ordered consistently (property A guarantees a consistent ordering exists, e.g. red/black; the two are not equivalent in general though they coincide for the red/black-ordered finite-difference Laplacian here). Young's theorem gives an exact algebraic link between the Jacobi eigenvalues λ (of D⁻¹(L+U)) and the SOR eigenvalues φ (of L_ω):
- (φ + ω − 1)² = ω² λ² φ.
This quadratic in φ is the whole engine. For fixed ω you solve for φ in terms of λ; you then choose ω to minimize the maximum |φ|. The maximum is minimized precisely when the discriminant vanishes so the two roots coincide (the eigenvalue collision), giving ω_opt = 2 / (1 + √(1 − μ²)) where μ = ρ(Jacobi) = max|λ|, and at that value ρ(L_ω_opt) = ω_opt − 1. Since ρ(Gauss–Seidel) = μ², when μ = 1 − c·h² the Gauss–Seidel rate is 1 − 2c·h² but ρ(L_ω_opt) ≈ 1 − 2√(2c)·h — square-root better, hence O(1/h) fewer iterations.
Worked example: the 1D discrete Laplacian
Discretize −u″ = f on (0,1) with n interior points, spacing h = 1/(n+1). The matrix A = tridiag(−1, 2, −1) is symmetric positive definite and consistently ordered, so the theory applies exactly. Its Jacobi spectral radius is
- μ = cos(π/(n+1)).
Take n = 20. Then μ = cos(π/21) ≈ 0.98883, so Gauss–Seidel has ρ = μ² ≈ 0.97779 — each sweep removes barely 2.2% of the error. Plug μ into the formula: ω_opt = 2/(1 + √(1 − μ²)) ≈ 1.7406, and ρ(L_ω_opt) = ω_opt − 1 ≈ 0.7406. Now each sweep removes 26% of the error. Reaching a 10⁻⁶ reduction takes ⌈−6·ln 10 / ln 0.9778⌉ ≈ 615 Gauss–Seidel sweeps versus ⌈−6·ln 10 / ln 0.7406⌉ ≈ 46 SOR sweeps — a 13× reduction, growing without bound as n → ∞. (Direct eigenvalue computation of L_ω confirms these radii to five digits.)
Why the hypotheses matter: Kahan, Ostrowski–Reich, and where it breaks
Two limits are non-negotiable. Kahan's theorem (1958): for any matrix, ρ(L_ω) ≥ |ω − 1|. The proof is a one-liner: det(L_ω) = det[(1−ω)D + ωU] / det(D − ωL) = (1 − ω)ⁿ (both matrices are triangular), and the spectral radius is at least the geometric mean of eigenvalue moduli = |1 − ω|. Hence ω ∈ (0, 2) is necessary for convergence — pick ω ≥ 2 and it always diverges.
Sufficiency needs structure. Ostrowski–Reich theorem (1949–1954): if A is symmetric positive definite, then ρ(L_ω) < 1 for every ω ∈ (0, 2). Drop positive-definiteness and it can fail: SOR may diverge for all ω, or converge only in a strict sub-interval of (0, 2). The exact formula ω_opt = 2/(1+√(1−μ²)) additionally needs consistent ordering; without it the elegant φ–λ relation is false and ω_opt must be estimated numerically or by adaptive schemes. Non-normal A can also make ρ < 1 while transient error grows for many iterations.
Significance: from the 1950s workhorse to a preconditioner
SOR (David Young's 1950 Harvard thesis under Garrett Birkhoff; Stanley Frankel independently the same year) was the dominant solver for elliptic PDEs through the 1960s–70s, the method that made large finite-difference computations feasible on early machines. Its intellectual legacy is larger than its current standalone use:
- Symmetric SOR (SSOR) — a forward sweep followed by a backward sweep — yields a symmetric positive-definite splitting used as a preconditioner for conjugate gradients and other Krylov methods, where it is very much alive today.
- Multigrid superseded SOR as a standalone solver (it achieves the optimal O(N) work), but uses relaxation sweeps — often weighted Jacobi or SOR — as its smoother, exploiting exactly the low-/high-frequency structure that motivates over-relaxation.
- The eigenvalue-collision analysis and the O(h²)→O(h) rate improvement are a template for understanding acceleration by extrapolation, echoed in Chebyshev semi-iteration and momentum methods.
Knowing when the exact ω_opt formula applies — SPD and consistently ordered — separates a textbook speedup from a divergent mess.
| Method | Iteration matrix | Spectral radius (model) | Asymptotic rate |
|---|---|---|---|
| Jacobi | D⁻¹(L + U) | μ = cos(π/(n+1)) | O(h²), h = 1/(n+1) |
| Gauss–Seidel | (D − L)⁻¹U | μ² = cos²(π/(n+1)) | 2× faster than Jacobi |
| SOR, generic ω∈(0,2) | (D − ωL)⁻¹[(1−ω)D + ωU] | depends on ω; ≥ |ω − 1| | varies |
| SOR, optimal ω_opt | same, ω = ω_opt | ω_opt − 1 = (1−√(1−μ²))/(1+√(1−μ²)) | O(h): an order-of-magnitude jump |
Frequently asked questions
What is the optimal relaxation parameter ω, and when is the formula valid?
For a consistently ordered matrix with Jacobi spectral radius μ = ρ(Jacobi) < 1, the optimum is ω_opt = 2 / (1 + √(1 − μ²)), and then ρ(L_ω_opt) = ω_opt − 1. The formula requires consistent ordering (property A) — true for standard 5-point/7-point finite-difference Laplacians. Without that structure you must estimate ω_opt numerically or adaptively, since the underlying Young eigenvalue relation no longer holds.
Why must ω lie strictly between 0 and 2?
Kahan's theorem: for any matrix, det(L_ω) = (1 − ω)ⁿ because both triangular factors have determinants (1−ω)ⁿ·(stuff) that cancel, forcing the product of eigenvalues to have modulus |1−ω|ⁿ. So ρ(L_ω) ≥ |ω − 1|. If |ω − 1| ≥ 1, i.e. ω ≤ 0 or ω ≥ 2, the spectral radius is at least 1 and the iteration cannot converge. Thus ω ∈ (0, 2) is necessary (though not always sufficient).
Does symmetric positive definiteness guarantee convergence?
Yes. The Ostrowski–Reich theorem states that if A is symmetric positive definite, then ρ(L_ω) < 1 for every ω ∈ (0, 2), so SOR converges for any such ω and any starting vector. This is a sufficient condition on A alone; it does not tell you the optimal ω, only that the whole open interval (0,2) is safe. Optimality still requires the consistent-ordering hypothesis.
How much faster is SOR than Gauss–Seidel, concretely?
On the 2D model Poisson problem with mesh size h, Gauss–Seidel has ρ ≈ 1 − O(h²), needing O(h⁻²) = O(n²) iterations, while optimal SOR has ρ(L_ω_opt) ≈ 1 − O(h), needing only O(h⁻¹) = O(n). That is an order-of-magnitude reduction in iteration count that grows without bound as the grid refines — for the n=20 1D example it is already a 13× reduction (615 vs 46 sweeps for 10⁻⁶).
What breaks if the matrix is not positive definite or not consistently ordered?
Two different failures. Lose positive-definiteness (or diagonal dominance) and Ostrowski–Reich no longer applies: SOR may converge only on a strict sub-interval of (0,2), or diverge for all ω. Lose consistent ordering and the exact φ–λ relation (φ+ω−1)² = ω²λ²φ fails, so the clean ω_opt formula is wrong — the true optimum can be quite different and must be found numerically. Non-normal A can also give ρ < 1 yet large transient error growth.
How does SOR relate to modern solvers like multigrid and conjugate gradient?
SOR is rarely a standalone solver now, but its ideas persist. Symmetric SOR (a forward then backward sweep) is an SPD splitting used as a preconditioner for conjugate gradient and other Krylov methods. Multigrid, which achieves optimal O(N) work, uses relaxation sweeps (weighted Jacobi or SOR) as smoothers — leveraging the very high-/low-frequency error structure that over-relaxation was designed to exploit.