Prev
Unit 4 · Linear Systems: Iterative Methods
Successive Over-Relaxation
55m100 XP
Next
🔁Linear Systems: Iterative Methods

Successive Over-Relaxation

SOR, Young's optimal ω, cost versus grid size, red–black ordering and stopping criteria

Cambridge, Massachusetts, 1950
David M. Young is a graduate student in mathematics at Harvard, supervised by Garrett Birkhoff. Electronic computers are just beginning to solve the partial differential equations of physics on grids, and the bottleneck is obvious: Gauss–Seidel on a fine grid needs tens of thousands of sweeps. Young studies a variant engineers had tried by hand: compute the Gauss–Seidel correction, then multiply it by a factor ω\omega a little larger than 1 — over-relax. Engineers chose ω\omega by feel. Young's thesis, Iterative Methods for Solving Partial Difference Equations of Elliptic Type, does something different: for a large class of matrices it finds the exact optimal ω\omega, and proves how fast the iteration then converges. The result is startling. On a fine grid, the right ω\omega makes the method not twice as fast, or ten times, but hundreds of times faster. Stanley Frankel at Caltech found the same idea independently that year. SOR became the workhorse of reactor physics and fluid dynamics codes for the next two decades.
Let x~i\tilde x_i be the value Gauss–Seidel would assign to xix_i. SOR moves xix_i further in the same direction:
xi(k+1)=xi(k)+ω(x~i−xi(k)),0<ω<2x_i^{(k+1)} = x_i^{(k)} + \omega\left(\tilde x_i - x_i^{(k)}\right),\qquad 0 < \omega < 2

ω = 1 is Gauss–Seidel. ω > 1 over-relaxes, ω < 1 under-relaxes (useful to damp oscillating nonlinear iterations). Outside 0 < ω < 2 SOR always diverges (Kahan's theorem).

Derivation
SOR as a splitting, and why ω\omega must lie between 0 and 2
Write the update for component ii with the Gauss–Seidel value spelled out, and multiply by aiia_{ii}:
aiixinew=(1−ω) aiixi+ω(bi−∑j<iaijxjnew−∑j>iaijxj).a_{ii}x_i^{\text{new}} = (1 - \omega)\,a_{ii}x_i + \omega\Big(b_i - \sum_{j<i}a_{ij}x_j^{\text{new}} - \sum_{j>i}a_{ij}x_j\Big).
Move the new values to the left. With A=D+L+UA = D + L + U as in lesson 4-1, all rows together read
(D+ωL) xnew=((1−ω)D−ωU)x+ωb,Tω=(D+ωL)−1((1−ω)D−ωU).(D + \omega L)\,\mathbf x^{\text{new}} = \big((1 - \omega)D - \omega U\big)\mathbf x + \omega\mathbf b, \qquad T_\omega = (D + \omega L)^{-1}\big((1 - \omega)D - \omega U\big).
Kahan's bound. Both bracketed matrices are triangular, so their determinants are the products of their diagonals: det⁡(D+ωL)=det⁡D\det(D + \omega L) = \det D and det⁡((1−ω)D−ωU)=(1−ω)ndet⁡D\det\big((1 - \omega)D - \omega U\big) = (1 - \omega)^n\det D. Hence det⁡Tω=(1−ω)n\det T_\omega = (1 - \omega)^n. A determinant is the product of the eigenvalues, so at least one eigenvalue has ∣λ∣≥∣1−ω∣|\lambda| \ge |1 - \omega|, and
ρ(Tω)≥∣ω−1∣.\rho(T_\omega) \ge |\omega - 1|.
For ω≤0\omega \le 0 or ω≥2\omega \ge 2 this is at least 1, and SOR cannot converge, whatever the matrix.
💡Young's optimal ω
For matrices with 'Property A' — which includes every five-point grid matrix — Young showed that the SOR spectral radius depends only on ω\omega and the Jacobi spectral radius μ=ρJ\mu = \rho_J, and is smallest at ωopt=21+1−μ2,ρSOR(ωopt)=ωopt−1.\omega_{opt} = \dfrac{2}{1 + \sqrt{1 - \mu^2}}, \qquad \rho_{SOR}(\omega_{opt}) = \omega_{opt} - 1. For the model problem μ=cos⁡πh\mu = \cos\pi h, so ωopt=2/(1+sin⁡πh)\omega_{opt} = 2/(1 + \sin\pi h) and 1−ρSOR≈2πh1 - \rho_{SOR} \approx 2\pi h. Compare Gauss–Seidel's 1−ρ≈π2h21 - \rho \approx \pi^2 h^2: SOR's convergence rate degrades like hh, not h2h^2.
Derivation
From Young's relation to ωopt\omega_{opt}
Young's thesis proved that, for matrices with Property A, every eigenvalue λ\lambda of the SOR matrix is linked to an eigenvalue μ\mu of the Jacobi matrix by
(λ+ω−1)2=λ ω2μ2.(\lambda + \omega - 1)^2 = \lambda\,\omega^2\mu^2.
We take that relation as given and find the best ω\omega from it. Multiplied out, it is a quadratic in λ\lambda:
λ2−(ω2μ2−2(ω−1))λ+(ω−1)2=0.\lambda^2 - \big(\omega^2\mu^2 - 2(\omega - 1)\big)\lambda + (\omega - 1)^2 = 0.
Its two roots multiply to (ω−1)2(\omega - 1)^2. When the discriminant is negative the roots are complex conjugates of equal size, so each has ∣λ∣=ω−1|\lambda| = \omega - 1, which grows as ω\omega grows. When it is positive the roots are real, and the larger one shrinks as ω\omega grows. The best ω\omega is where the two regimes meet. The discriminant of the quadratic works out to ω2μ2(ω2μ2−4(ω−1))\omega^2\mu^2\big(\omega^2\mu^2 - 4(\omega - 1)\big), so the meeting point is where the bracket is zero:
ω2μ2−4(ω−1)=0⟹ω=2(1−1−μ2)μ2=21+1−μ2,\omega^2\mu^2 - 4(\omega - 1) = 0 \quad\Longrightarrow\quad \omega = \frac{2\big(1 - \sqrt{1 - \mu^2}\big)}{\mu^2} = \frac{2}{1 + \sqrt{1 - \mu^2}},
the last step multiplying top and bottom by 1+1−μ21 + \sqrt{1 - \mu^2}. There ρSOR=ωopt−1\rho_{SOR} = \omega_{opt} - 1. The worst case is the largest Jacobi eigenvalue, μ=ρJ\mu = \rho_J.
The model problem. With ρJ=cos⁡πh\rho_J = \cos\pi h, 1−μ2=sin⁡πh\sqrt{1 - \mu^2} = \sin\pi h, so
ωopt=21+sin⁡πh,ρSOR=1−sin⁡πh1+sin⁡πh≈1−2πh.\omega_{opt} = \frac{2}{1 + \sin\pi h}, \qquad \rho_{SOR} = \frac{1 - \sin\pi h}{1 + \sin\pi h} \approx 1 - 2\pi h.
One digit costs ln⁡10/(2πh)\ln 10/(2\pi h) sweeps: proportional to 1/h1/h instead of Gauss–Seidel's 1/h21/h^2. For N=320N = 320 that is the 24,039 against 117.6 of the table.
NJacobiGauss–Seideloptimal SOR
1056284.0
4078439215.0
16012,0946,04759.0
32048,07824,039117.6

Sweeps per decimal digit for the model Laplace problem on an N × N grid (from the spectral radii)

📐Choosing ω on the plate
On the 19×1919 \times 19 plate grid, ωopt=2/(1+sin⁡(π/20))=1.7295\omega_{opt} = 2/(1 + \sin(\pi/20)) = 1.7295. Sweeps to a change below 10−810^{-8}: ω=1.0\omega = 1.0: 740 · 1.51.5: 248 · 1.71.7: 108 · 1.731.73: 80 · 1.81.8: 105 · 1.91.9: 213 · 1.951.95: 442. The curve has a sharp minimum, and it is asymmetric: overshooting ωopt\omega_{opt} a little costs far less than undershooting. When ωopt\omega_{opt} is not known exactly, err on the high side.
The first principle is that you must not fool yourself — and you are the easiest person to fool.
– Richard Feynman, Cargo Cult Science (1974)
When to stop. The cheapest test — stop when no value changes by more than some τ\tau — is also the most dangerous. For an iteration with rate ρ\rho, the remaining error is roughly
∥e(k)∥≈ρ1−ρ ∥x(k)−x(k−1)∥.\|\mathbf e^{(k)}\| \approx \dfrac{\rho}{1-\rho}\,\|\mathbf x^{(k)} - \mathbf x^{(k-1)}\|.
When ρ=0.999\rho = 0.999, a change of 10−610^{-6} hides an error of about 10−310^{-3}. Slow iterations look converged long before they are (the derivation below shows where the factor comes from). On the 63×6363 \times 63 plate, Gauss–Seidel stopped by a 10−610^{-6} update test is still 4×10−44 \times 10^{-4} °C wrong; SOR, with a much smaller ρ/(1−ρ)\rho/(1-\rho), stops at 6×10−66 \times 10^{-6}. Prefer a test on the relative residual ∥b−Ax(k)∥/∥b∥\|\mathbf b - A\mathbf x^{(k)}\|/\|\mathbf b\|, and remember lesson 3-2: a small residual means a small error only when AA is well conditioned.
Derivation
How much error a small update hides
Once one mode dominates, each sweep multiplies the error by about ρ\rho: e(k)≈ρ e(k−1)\mathbf e^{(k)} \approx \rho\,\mathbf e^{(k-1)}. The update you can see is the difference of two consecutive iterates, which is also the difference of their errors (the exact solution cancels):
x(k)−x(k−1)=e(k)−e(k−1)≈e(k)−e(k)ρ=−1−ρρ e(k).\mathbf x^{(k)} - \mathbf x^{(k-1)} = \mathbf e^{(k)} - \mathbf e^{(k-1)} \approx \mathbf e^{(k)} - \frac{\mathbf e^{(k)}}{\rho} = -\frac{1 - \rho}{\rho}\,\mathbf e^{(k)}.
Turn it round:
∥e(k)∥≈ρ1−ρ ∥x(k)−x(k−1)∥.\|\mathbf e^{(k)}\| \approx \frac{\rho}{1 - \rho}\,\|\mathbf x^{(k)} - \mathbf x^{(k-1)}\|.
For ρ=0.5\rho = 0.5 the error is about the size of the update; for ρ=0.999\rho = 0.999 it is a thousand times larger. It is the same bound as for a scalar contraction in lesson 2-3, with LL replaced by ρ\rho.
💡Red–black ordering
Colour the grid like a chessboard. Every red point's four neighbours are black, and vice versa. So all red points can be updated simultaneously using the current black values, then all black points using the new red values. Each half-sweep is a single vectorised (or parallel) operation, yet the result is a genuine Gauss–Seidel/SOR sweep in a different ordering — with the same optimal ω\omega and convergence rate for the model problem. It is how SOR is written in NumPy and on GPUs.
💡Where you will use this
CFD: pressure-correction methods (SIMPLE and its descendants) under-relax velocity and pressure updates with factors like 0.7 and 0.3 — the same ω\omega, used for stability rather than speed. Geotechnics and Hydraulics: SOR remains a quick, memory-light solver for seepage and groundwater grids. Power Systems: accelerated Gauss–Seidel load flow used acceleration factors around 1.6.
📚Further reading
  • David M. Young, Iterative Solution of Large Linear Systems (Academic Press, 1971; Dover reprint 2003) — Young's own book, with the full theory of the optimal ω\omega.
  • Yousef Saad, Iterative Methods for Sparse Linear Systems (2nd ed., SIAM, 2003), chapter 4 — SOR, red–black ordering and stopping criteria.
  • William Kahan, Gauss–Seidel Methods of Solving Large Systems of Linear Equations, PhD thesis, University of Toronto (1958) — where the bound ρ≥∣ω−1∣\rho \ge |\omega - 1| first appeared.
↓But there is a problem
SOR is a big improvement, but it shares a weakness with Jacobi and Gauss–Seidel: it is stationary. Every sweep applies the same fixed rule, and the method remembers nothing but its latest iterate. And it needs ωopt\omega_{opt}, which is known exactly only for model problems. In 1952 Magnus Hestenes and Eduard Stiefel published a method that remembers. Each step chooses its direction so that it never undoes the progress of any earlier step. In exact arithmetic it finds the solution of an n×nn \times n system in at most nn steps. Then, for almost twenty years, almost nobody used it.