Successive over-relaxation in Python
Gauss-Seidel with a relaxation factor omega that can speed convergence.
- File
successive_over_relaxation.py, 89 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/successive-over-relaxation/. The authors are listed in prov/MJSE1P.jsonl.
"""Successive over-relaxation (SOR) for A x = b, with A stored by sparse rows.
A matrix is a list of rows; row i is a list of (j, a_ij) pairs with a_ij != 0.
from_dense converts a list of lists. Entries may be floats or Fractions.
"""
import math
from fractions import Fraction
from typing import TypeVar
T = TypeVar("T", float, Fraction)
SparseRows = list[list[tuple[int, T]]]
class ZeroDiagonal(ValueError):
"""Raised when some a_ii is zero, so the SOR sweep is undefined."""
def from_dense(a: list[list[T]]) -> SparseRows:
"""Keep the nonzero entries of each row."""
return [[(j, v) for j, v in enumerate(row) if v != 0] for row in a]
def residual(rows: SparseRows, b: list[T], x: list[T]) -> list[T]:
"""r = b - A x, in 2 nnz(A) flops."""
return [b[i] - sum(v * x[j] for j, v in row) for i, row in enumerate(rows)]
def _diag(rows: SparseRows) -> list:
d = []
for i, row in enumerate(rows):
aii = next((v for j, v in row if j == i), 0)
if aii == 0:
raise ZeroDiagonal("a[%d][%d] is zero" % (i, i))
d.append(aii)
return d
def sor_sweep(rows: SparseRows, b: list[T], x: list[T], omega: T, diag: list[T]) -> None:
"""One forward SOR sweep, in place: x_i <- (1 - omega) x_i + omega (b_i - sum_{j != i} a_ij x_j) / a_ii.
Cost: 2 nnz(A) + 2n flops.
"""
for i, row in enumerate(rows): # @line L3
s = b[i]
for j, v in row: # @line L4
if j != i:
s -= v * x[j]
x[i] = (1 - omega) * x[i] + omega * s / diag[i] # @line L5
def sor(rows: SparseRows, b: list[float], omega: float, x0: list[float] | None = None,
tol: float = 1e-10, maxiter: int = 10000) -> tuple[list[float], int, bool]:
"""Run SOR sweeps until ||b - A x||_2 <= tol ||b||_2 or maxiter sweeps.
Returns (x, sweeps, converged). Precondition: a_ii != 0 for all i (raises
ZeroDiagonal). Converges for every x0 iff rho(G_omega) < 1; for symmetric
positive definite A this holds exactly when 0 < omega < 2.
Cost per sweep: 2 nnz(A) + 2n flops for the sweep and 2 nnz(A) + 2n for the residual test.
"""
n = len(rows)
diag = _diag(rows)
x = list(x0) if x0 is not None else [0.0] * n # @line L1
nb = math.sqrt(sum(v * v for v in b))
target = tol * (nb if nb > 0 else 1.0)
for k in range(1, maxiter + 1): # @line L2
sor_sweep(rows, b, x, omega, diag)
r = residual(rows, b, x)
if math.sqrt(sum(v * v for v in r)) <= target: # @line L6
return x, k, True
return x, maxiter, False # @line L7
def optimal_omega(mu: float) -> float:
"""Young's optimal parameter 2 / (1 + sqrt(1 - mu^2)) for a consistently ordered matrix with rho(Jacobi) = mu < 1."""
return 2.0 / (1.0 + math.sqrt(1.0 - mu * mu))
def sor_radius(mu: float, omega: float) -> float:
"""rho(G_omega) for a consistently ordered matrix with real Jacobi spectrum and rho(Jacobi) = mu < 1, 0 < omega < 2."""
wopt = optimal_omega(mu)
if omega >= wopt:
return omega - 1.0
t = omega * mu + math.sqrt(omega * omega * mu * mu - 4.0 * (omega - 1.0))
return (t / 2.0) ** 2
test_successive_over_relaxation.py 216 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/successive-over-relaxation/. The authors are listed in prov/MJSE1P.jsonl.
import math
import random
from fractions import Fraction
import pytest
from successive_over_relaxation import (
ZeroDiagonal, from_dense, optimal_omega, residual, sor, sor_radius, sor_sweep,
)
U = 2.0 ** -53
def model(n):
return [[2.0 if i == j else (-1.0 if abs(i - j) == 1 else 0.0) for j in range(n)] for i in range(n)]
def exact_solve(a, b):
n = len(b)
g = [[Fraction(v) for v in row] + [Fraction(bi)] for row, bi in zip(a, b)]
for k in range(n):
p = next(i for i in range(k, n) if g[i][k] != 0)
g[k], g[p] = g[p], g[k]
piv = g[k][k]
g[k] = [v / piv for v in g[k]]
for i in range(n):
if i != k and g[i][k] != 0:
f = g[i][k]
g[i] = [vi - f * vk for vi, vk in zip(g[i], g[k])]
return [row[n] for row in g]
def inf_norm_inverse(a):
n = len(a)
cols = [exact_solve(a, [1 if i == j else 0 for i in range(n)]) for j in range(n)]
return max(sum(abs(cols[j][i]) for j in range(n)) for i in range(n))
def random_spd(rng, n):
b = [[Fraction(rng.randint(-3, 3)) for _ in range(n)] for _ in range(n)]
return [[sum(b[k][i] * b[k][j] for k in range(n)) + (n if i == j else 0) for j in range(n)] for i in range(n)]
def gauss_seidel_sweep(a, b, x):
n = len(b)
for i in range(n):
x[i] = (b[i] - sum(a[i][j] * x[j] for j in range(n) if j != i)) / a[i][i]
def test_one_by_one():
x, k, ok = sor(from_dense([[5.0]]), [10.0], 1.0)
assert ok and k == 1 and x == [2.0]
def test_diagonal_exact_in_one_sweep():
d = [3.0, -2.0, 0.5, 8.0]
a = [[d[i] if i == j else 0.0 for j in range(4)] for i in range(4)]
x, k, ok = sor(from_dense(a), [6.0, 4.0, 1.0, 16.0], 1.0)
assert ok and k == 1 and x == [2.0, -2.0, 2.0, 2.0]
def test_hand_sweep():
a = [[Fraction(v) for v in r] for r in [[4, -1, 0], [-1, 4, -1], [0, -1, 4]]]
x = [Fraction(0)] * 3
sor_sweep(from_dense(a), [Fraction(3), Fraction(2), Fraction(3)], x, Fraction(6, 5), [Fraction(4)] * 3)
assert x == [Fraction(9, 10), Fraction(87, 100), Fraction(1161, 1000)]
def test_omega_one_is_gauss_seidel_exactly():
rng = random.Random(1)
for _ in range(30):
n = rng.randint(2, 6)
a = [[Fraction(rng.randint(-4, 4)) for _ in range(n)] for _ in range(n)]
for i in range(n):
a[i][i] = Fraction(rng.choice([-1, 1]) * rng.randint(1, 9))
b = [Fraction(rng.randint(-9, 9)) for _ in range(n)]
x = [Fraction(rng.randint(-3, 3)) for _ in range(n)]
y = list(x)
diag = [a[i][i] for i in range(n)]
for _ in range(3):
sor_sweep(from_dense(a), b, x, Fraction(1), diag)
gauss_seidel_sweep(a, b, y)
assert x == y
def test_solution_is_a_fixed_point():
rng = random.Random(2)
for _ in range(30):
n = rng.randint(1, 6)
a = random_spd(rng, n)
b = [Fraction(rng.randint(-9, 9)) for _ in range(n)]
xs = exact_solve(a, b)
x = list(xs)
omega = Fraction(rng.randint(1, 19), 10)
sor_sweep(from_dense(a), b, x, omega, [a[i][i] for i in range(n)])
assert x == xs
def test_random_spd_converges_and_error_is_bounded_by_kappa_times_residual():
rng = random.Random(3)
for _ in range(40):
n = rng.randint(2, 7)
a = random_spd(rng, n)
b = [Fraction(rng.randint(-9, 9)) for _ in range(n)]
if all(v == 0 for v in b):
b[0] = Fraction(1)
omega = rng.choice([0.4, 1.0, 1.3, 1.8])
af = [[float(v) for v in row] for row in a]
x, k, ok = sor(from_dense(af), [float(v) for v in b], omega, tol=1e-12, maxiter=20000)
assert ok, (n, omega)
xs = exact_solve(a, b)
r = [b[i] - sum(a[i][j] * Fraction(x[j]) for j in range(n)) for i in range(n)]
norm_a = max(sum(abs(v) for v in row) for row in a)
kappa = norm_a * inf_norm_inverse(a)
err = max(abs(Fraction(x[i]) - xs[i]) for i in range(n)) / max(abs(v) for v in xs)
relres = max(abs(v) for v in r) / max(abs(v) for v in b)
assert err <= kappa * relres # exact: ||e||/||x|| <= kappa ||r||/||b||
def predicted(rho, tol):
return math.log(tol) / math.log(rho)
def test_iteration_counts_match_young_theory():
n, tol = 20, 1e-8
rows, b = from_dense(model(n)), [1.0] * n
mu = math.cos(math.pi / (n + 1))
counts = {}
for omega in (1.0, 1.5, 1.9):
_, k, ok = sor(rows, b, omega, tol=tol)
assert ok
ratio = k / predicted(sor_radius(mu, omega), tol)
assert 0.85 <= ratio <= 1.2, (omega, k, ratio)
counts[omega] = k
wopt = optimal_omega(mu)
_, kopt, ok = sor(rows, b, wopt, tol=tol)
assert ok and kopt < min(counts.values())
assert kopt <= 1.5 * predicted(wopt - 1, tol) # a Jordan block at omega_opt slows it a little
def test_optimal_sor_scales_like_n_gauss_seidel_like_n_squared():
tol = 1e-8
gs, opt = [], []
for n in (10, 20, 40):
rows, b = from_dense(model(n)), [1.0] * n
mu = math.cos(math.pi / (n + 1))
opt.append(sor(rows, b, optimal_omega(mu), tol=tol)[1])
if n < 40:
gs.append(sor(rows, b, 1.0, tol=tol)[1])
assert 3.3 <= gs[1] / gs[0] <= 4.2
for p, q in zip(opt, opt[1:]):
assert 1.7 <= q / p <= 2.3
def test_radius_formula_against_observed_contraction():
n = 10
a = model(n)
rows = from_dense(a)
mu = math.cos(math.pi / (n + 1))
assert sor_radius(mu, 1.0) == pytest.approx(mu * mu)
w = optimal_omega(mu)
assert sor_radius(mu, w) == pytest.approx(w - 1)
rng = random.Random(4)
for omega in (0.8, 1.0, 1.3):
x = [rng.uniform(-1, 1) for _ in range(n)]
zero = [0.0] * n
for _ in range(150):
sor_sweep(rows, zero, x, omega, [2.0] * n)
before = math.sqrt(sum(v * v for v in x))
sor_sweep(rows, zero, x, omega, [2.0] * n)
after = math.sqrt(sum(v * v for v in x))
assert after / before == pytest.approx(sor_radius(mu, omega), rel=0.02)
def test_omega_outside_zero_two_diverges():
n = 20
rows, b = from_dense(model(n)), [1.0] * n
for omega in (2.0, 2.2, -0.5):
x, k, ok = sor(rows, b, omega, maxiter=300)
assert not ok
x, _, _ = sor(rows, b, 2.2, maxiter=60)
assert max(abs(v) for v in residual(rows, b, x)) > 1e3
def test_symmetric_indefinite_diverges_for_every_omega():
rows = from_dense([[1.0, 2.0], [2.0, 1.0]])
for omega in (0.5, 1.0, 1.5):
x, k, ok = sor(rows, [1.0, 1.0], omega, maxiter=100)
assert not ok
assert not max(abs(v) for v in x) <= 1e3 # grows without bound (or overflows to inf, nan)
def test_nonsymmetric_gauss_seidel_divergence():
rows = from_dense([[1.0, 2.0, -2.0], [1.0, 1.0, 1.0], [2.0, 2.0, 1.0]]) # rho(G_1) = 2
x, k, ok = sor(rows, [1.0, 3.0, 5.0], 1.0, maxiter=100)
assert not ok
def test_zero_diagonal_raises():
with pytest.raises(ZeroDiagonal):
sor(from_dense([[0.0, 1.0], [1.0, 0.0]]), [1.0, 1.0], 1.0)
def test_stagnation_at_rounding_level():
n = 20
rows, b = from_dense(model(n)), [1.0] * n
x, k, ok = sor(rows, b, 1.74, tol=1e-18, maxiter=1500)
assert not ok # the tolerance is below what rounding allows
r = max(abs(v) for v in residual(rows, b, x))
xn = max(abs(v) for v in x)
assert r <= n * U * 4.0 * xn # observed: about 0.4 u ||A|| ||x||, not a proved bound
exact = [(i + 1) * (n - i) / 2 for i in range(n)]
assert max(abs(p - q) for p, q in zip(x, exact)) <= 1e-12
What the next level needs
- the tests did not run in this build