Mathematical foundations of Gauss-Seidel & SOR.
The Successive Over-Relaxation (SOR) method—developed independently in 1950 by the American mathematician David M. Young Jr. in his doctoral thesis and by Stanley P. Frankel for nuclear weapon calculations on the ENIAC—is one of the most celebrated algorithmic breakthroughs in the history of computational linear algebra. In the 1940s, solving the finite-difference equations of heat transfer and neutron diffusion was severely bottlenecked by the slow convergence of the Jacobi and Gauss-Seidel iterations. Gauss-Seidel improved upon Jacobi by immediately utilizing newly computed components x_j^{(k+1)} within the same iteration: x_i^{GS} = (b_i − ∑_{j<i} a_{ij}x_j^{(k+1)} − ∑_{j>i} a_{ij}x_j^{(k)}) / a_{ii}. However, on fine computational meshes, Gauss-Seidel error reduction is plagued by sluggishness: the updates Δx_i = x_i^{GS} − x_i^{(k)} crawl forward monotonically. Young and Frankel introduced a bold modification: why not accelerate the step by multiplying the correction by a relaxation weight ω > 1? The resulting SOR iteration formula takes the form: x_i^{(k+1)} = (1 − ω)x_i^{(k)} + (ω/a_{ii})(b_i − ∑_{j<i} a_{ij}x_j^{(k+1)} − ∑_{j>i} a_{ij}x_j^{(k)}). In matrix notation, this corresponds to the matrix splitting (D − ωL)x^{(k+1)} = [(1 − ω)D + ωU]x^{(k)} + ωb. The Ostrowski-Reich Theorem establishes that for any symmetric positive-definite matrix, SOR converges if and only if 0 < ω < 2. In a triumph of applied analysis, David Young derived an exact closed-form expression for the optimal relaxation parameter minimizing the spectral radius on 2-cyclic consistently ordered matrices: ω_opt = 2 / [1 + √(1 − ρ(T_J)²)]. By tuning ω to ω_opt, the spectral radius drops from ρ_{GS} = 1 − O(h²) to ρ_{SOR} = 1 − O(h). On a 1000×1000 grid, this reduces the required iterations from hundreds of thousands down to a few hundred, cutting computation time by several orders of magnitude. In this computational laboratory, students evaluate the SOR method applied to a linear system with Jacobi spectral radius ρ = 0.500, investigate parameter sensitivity across the stability domain ω ∈ (0, 2), and extract the theoretical optimal parameter ω_opt = 1.268 via linear regression.
Viva voce questions & analytical theory
What is the physical and mathematical intuition behind Successive Over-Relaxation (SOR)?
In standard Gauss-Seidel iteration (ω = 1), the newly calculated value x_i^{GS} is immediately substituted into the remaining equations. However, on smooth elliptic problems, Gauss-Seidel iterates approach the exact solution monotonically from one direction, meaning the step Δx_i = x_i^{GS} − x_i^{(k)} consistently points toward the root but falls short. Successive Over-Relaxation accelerates this by extrapolating past the Gauss-Seidel estimate: x_i^{(k+1)} = x_i^{(k)} + ω(x_i^{GS} − x_i^{(k)}) = (1 − ω)x_i^{(k)} + ω x_i^{GS}. Choosing ω > 1 ('over-relaxation') overshoots the sluggish Gauss-Seidel step, leaping further into the basin of attraction and reducing the required iterations by orders of magnitude.
What is the Ostrowski-Reich Theorem regarding SOR convergence?
Formulated by Alexander Ostrowski (1954) and Edgar Reich (1949), the theorem states that for any real symmetric matrix A with positive diagonal entries, the SOR iteration converges if and only if: (1) A is positive-definite, and (2) the relaxation parameter satisfies 0 < ω < 2. If ω ≤ 0 or ω ≥ 2, the spectral radius ρ(T_ω) ≥ 1, causing the method to diverge. This establishes that the open interval (0, 2) is the universal stability domain for SOR on SPD systems.
What is David M. Young's formula for the optimal relaxation parameter ω_opt?
In his 1950 Harvard doctoral dissertation, David M. Young proved that if A possesses 'Property A' and is consistently ordered (which holds for 5-point Laplacian stencils on regular grids), the eigenvalues of the SOR iteration matrix μ satisfy: (μ + ω − 1)² = ω² μ λ², where λ are eigenvalues of the Jacobi matrix. Differentiating with respect to ω reveals that the spectral radius ρ(T_ω) is minimized at the unique optimal parameter: ω_opt = 2 / [1 + √(1 − ρ(T_J)²)]. At this optimal value, the spectral radius drops to ρ(T_{ω_opt}) = ω_opt − 1.
How much faster is optimal SOR than Jacobi and Gauss-Seidel on fine discretization grids?
For a 2D Poisson equation on an N×N grid with grid spacing h = 1/N, the Jacobi spectral radius is ρ_J ≈ 1 − (1/2)π²h². Gauss-Seidel has ρ_{GS} ≈ 1 − π²h², requiring O(N²) iterations to converge. In contrast, optimal SOR has ω_opt ≈ 2 − 2πh, giving a spectral radius ρ_{SOR} = ω_opt − 1 ≈ 1 − 2πh. Because 1 − 2πh is much smaller than 1 − π²h², the number of iterations required by SOR scales as O(N) rather than O(N²). On a 1000×1000 grid, SOR converges roughly 500 times faster than Gauss-Seidel!
What is under-relaxation (ω < 1) and when is it required?
Under-relaxation chooses 0 < ω < 1, taking a damped step that falls between the old iterate and the Gauss-Seidel prediction. While under-relaxation is slower on linear SPD systems, it is essential in highly nonlinear problems (such as the Navier-Stokes equations in CFD codes like OpenFOAM), where full Gauss-Seidel steps cause wild oscillations, pressure-velocity divergence, and numerical explosion.
How does Gauss-Seidel differ from Jacobi in memory access and caching?
Gauss-Seidel overwrites the solution vector in-place: as soon as x_1^{(k+1)} is computed, it is stored in array element x[1], and the very next line computing x_2^{(k+1)} immediately reads the fresh value x[1]. This halves memory storage (requiring only one vector instead of two) and dramatically improves CPU cache locality. However, this sequential data dependency makes standard Gauss-Seidel harder to parallelize on GPUs than Jacobi, requiring Red-Black (checkerboard) coloring.
What is the Symmetric Successive Over-Relaxation (SSOR) method?
Standard SOR uses an unsymmetric iteration matrix, preventing it from being used as a preconditioner for the Conjugate Gradient (CG) method. SSOR resolves this by executing a forward SOR sweep followed immediately by a backward SOR sweep. The resulting two-step iteration matrix is self-adjoint, making SSOR one of the most effective classical preconditioners for CG.
Algorithm precautions and stability rules
- Never select ω ≥ 2; the Ostrowski-Reich theorem proves that SOR is unconditionally divergent for ω ≥ 2.
- Accurately estimate the Jacobi spectral radius ρ(T_J) before choosing ω; an inaccurate ω degrades performance back to Gauss-Seidel.
- If using Red-Black coloring for parallel GPU execution, verify that the matrix graph is bipartite.
- Do not use standard unsymmetric SOR as a preconditioner for Conjugate Gradient; use SSOR instead.
- For non-symmetric systems, verify that all eigenvalues of the iteration matrix satisfy |λ| < 1.
- Maintain 64-bit precision to prevent accumulated round-off in backward-forward sweeps.
Theoretical limitations
- Optimal parameter formula ω_opt = 2/[1+√(1−ρ²)] requires matrix Property A and consistent ordering; general matrices require empirical tuning.
- In-place dependencies make naive SOR sequential and resistant to GPU vectorization without multi-color ordering.
- Overshooting with too large an ω (ω > ω_opt) increases the spectral radius faster than undershooting.