MINRES in Python
Krylov method for symmetric indefinite systems that minimizes the residual norm.
- File
minres.py, 101 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/minres/. The authors are listed in prov/NAYDFS.jsonl.
"""MINRES (Paige and Saunders, 1975) for A x = b with A symmetric, possibly indefinite.
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. The method needs only products A v, so
any symmetric operator given as a function can be passed instead (see minres).
"""
import math
from typing import Callable
SparseRows = list[list[tuple[int, float]]]
MatVec = Callable[[list[float]], list[float]]
def from_dense(a: list[list[float]]) -> 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 matvec(rows: SparseRows, x: list[float]) -> list[float]:
"""y = A x, in 2 nnz(A) flops."""
return [sum(v * x[j] for j, v in row) for row in rows]
def _dot(x: list[float], y: list[float]) -> float:
return math.fsum(a * b for a, b in zip(x, y))
def minres(a, b: list[float], x0: list[float] | None = None, tol: float = 1e-10,
maxiter: int | None = None, history: list[float] | None = None
) -> tuple[list[float], int, bool]:
"""Solve A x = b for symmetric A by minimizing ||b - A x_k||_2 over x0 + K_k(A, r0).
a is SparseRows or a function v -> A v. Stops when the recurrence residual
|eta| <= tol ||b||_2. Returns (x, k, converged); if history is a list, the
recurrence residual norm after each step is appended to it.
Cost per step: one product with A plus about 16n flops; 7 vectors of storage.
A must be symmetric; this is not checked.
"""
av: MatVec = a if callable(a) else (lambda v: matvec(a, v))
n = len(b)
if maxiter is None:
maxiter = 10 * n
x = list(x0) if x0 is not None else [0.0] * n
ax = av(x)
r = [bi - axi for bi, axi in zip(b, ax)] # @line L1
gamma = math.sqrt(_dot(r, r))
nb = math.sqrt(_dot(b, b))
target = tol * (nb if nb > 0 else 1.0)
if gamma <= target:
return x, 0, True
v_old = [0.0] * n # @line L2
v = r
w_old = [0.0] * n
w = [0.0] * n
eta = gamma
c_old, c, s_old, s = 1.0, 1.0, 0.0, 0.0
anorm = 0.0 # running estimate of ||A||_2 from the Lanczos coefficients
small = 10.0 * n * 2.0 ** -53
for k in range(1, maxiter + 1): # @line L3
v = [vi / gamma for vi in v] # @line L4
p = av(v)
delta = _dot(p, v)
v_new = [pi - delta * vi - gamma * voi for pi, vi, voi in zip(p, v, v_old)] # @line L5
gamma_new = math.sqrt(_dot(v_new, v_new))
anorm = max(anorm, math.sqrt(gamma * gamma + delta * delta + gamma_new * gamma_new))
a0 = c * delta - c_old * s * gamma # @line L6
a2 = s * delta + c_old * c * gamma
a3 = s_old * gamma
a1 = math.hypot(a0, gamma_new)
if a1 <= small * anorm: # T_k singular and K_k invariant: b is not in range(A)
return x, k - 1, False
c_new, s_new = a0 / a1, gamma_new / a1 # @line L7
w_new = [(vi - a3 * woi - a2 * wi) / a1 for vi, woi, wi in zip(v, w_old, w)] # @line L8
x = [xi + c_new * eta * wi for xi, wi in zip(x, w_new)] # @line L9
eta = -s_new * eta
if history is not None:
history.append(abs(eta))
if abs(eta) <= target: # @line L10
return x, k, True
if gamma_new <= small * anorm: # invariant subspace found: x is the exact minimizer
return x, k, False
v_old, v, gamma = v, v_new, gamma_new # @line L11
w_old, w = w, w_new
c_old, c, s_old, s = c, c_new, s, s_new
return x, maxiter, False # @line L12
def residual_norm(a, b: list[float], x: list[float]) -> float:
"""||b - A x||_2 computed directly."""
av: MatVec = a if callable(a) else (lambda v: matvec(a, v))
return math.sqrt(_dot(*(2 * [[bi - y for bi, y in zip(b, av(x))]])))
def interval_bound(c_lo: float, d_hi: float, k: int) -> float:
"""Bound on ||r_k|| / ||r_0|| for spectrum in [-d, -c] U [c, d]: 2 ((d - c)/(d + c))^floor(k/2)."""
q = (d_hi - c_lo) / (d_hi + c_lo)
return 2.0 * q ** (k // 2)
test_minres.py 179 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/minres/. The authors are listed in prov/NAYDFS.jsonl.
import math
import random
from fractions import Fraction
from minres import from_dense, interval_bound, matvec, minres, residual_norm
U = 2.0 ** -53
def diagop(d):
return lambda v: [a * b for a, b in zip(d, v)]
def random_symmetric(rng, n, indefinite=True):
"""Q diag(d) Q^T with Q a product of random Givens rotations, so the spectrum is known."""
d = [rng.uniform(1.0, 10.0) * (rng.choice([-1, 1]) if indefinite else 1) for _ in range(n)]
a = [[d[i] if i == j else 0.0 for j in range(n)] for i in range(n)]
for _ in range(3 * n):
p, q = rng.sample(range(n), 2)
t = rng.uniform(0, 2 * math.pi)
c, s = math.cos(t), math.sin(t)
for row in a: # a <- a G
row[p], row[q] = c * row[p] - s * row[q], s * row[p] + c * row[q]
a[p], a[q] = ([c * x - s * y for x, y in zip(a[p], a[q])],
[s * x + c * y for x, y in zip(a[p], a[q])])
for i in range(n): # symmetrize the rounding
for j in range(i):
a[i][j] = a[j][i] = (a[i][j] + a[j][i]) / 2
return a, d
def exact_min_residual(a, b, k):
"""min ||b - A y||_2 over y in K_k(A, b), in exact rational arithmetic (x0 = 0)."""
A = [[Fraction(v) for v in row] for row in a]
r0 = [Fraction(v) for v in b]
mv = lambda v: [sum(x * y for x, y in zip(row, v)) for row in A]
K, v = [], r0
for _ in range(k):
K.append(v)
v = mv(v)
M = [mv(col) for col in K] # columns A K
g = [[sum(x * y for x, y in zip(M[i], M[j])) for j in range(k)] + [sum(x * y for x, y in zip(M[i], r0))]
for i in range(k)]
for c in range(k): # Gauss-Jordan on the normal equations (M^T M is SPD here)
piv = g[c][c]
g[c] = [v / piv for v in g[c]]
for i in range(k):
if i != c and g[i][c] != 0:
f = g[i][c]
g[i] = [x - f * y for x, y in zip(g[i], g[c])]
y = [g[i][k] for i in range(k)]
rr = sum(v * v for v in r0) - sum(y[i] * sum(x * z for x, z in zip(M[i], r0)) for i in range(k))
return math.sqrt(rr)
def test_one_by_one():
x, k, conv = minres(from_dense([[-2.0]]), [4.0])
assert conv and k == 1 and abs(x[0] + 2.0) <= 4 * U
def test_zero_rhs():
x, k, conv = minres(from_dense([[1.0, 2.0], [2.0, -1.0]]), [0.0, 0.0])
assert conv and k == 0 and x == [0.0, 0.0]
def test_hand_example_stagnates_then_solves():
# A = [[0,1],[1,0]], b = e1: A b is orthogonal to b, so step 1 makes no progress.
h = []
x, k, conv = minres(from_dense([[0.0, 1.0], [1.0, 0.0]]), [1.0, 0.0], history=h)
assert conv and k == 2
assert h[0] == 1.0 and h[1] == 0.0
assert x == [0.0, 1.0]
def test_matches_exact_minimal_residual():
rng = random.Random(11)
for trial in range(25):
n = rng.randint(3, 7)
a, _ = random_symmetric(rng, n)
b = [rng.uniform(-1, 1) for _ in range(n)]
h = []
minres(from_dense(a), b, tol=0.0, maxiter=n - 1, history=h)
nb = math.sqrt(sum(v * v for v in b))
for k in range(1, len(h) + 1):
ref = exact_min_residual(a, b, k)
assert abs(h[k - 1] - ref) <= 1e-8 * nb, (trial, k)
def test_random_converges_and_error_bound():
rng = random.Random(5)
for trial in range(30):
n = rng.randint(2, 25)
a, d = random_symmetric(rng, n)
xs = [rng.uniform(-1, 1) for _ in range(n)]
b = matvec(from_dense(a), xs)
h = []
x, k, conv = minres(from_dense(a), b, tol=1e-10, history=h)
assert conv
assert all(h[i + 1] <= h[i] * (1 + 1e-12) for i in range(len(h) - 1)) # monotone
nb = math.sqrt(sum(v * v for v in b))
rel = residual_norm(from_dense(a), b, x) / nb
assert rel <= 1e-9
kappa = max(map(abs, d)) / min(map(abs, d))
err = math.sqrt(sum((p - q) ** 2 for p, q in zip(x, xs))) / math.sqrt(sum(q * q for q in xs))
assert err <= kappa * rel * 1.01 + 1e-14
def test_finite_termination_with_few_distinct_eigenvalues():
vals = [-3.0, -1.0, 0.5, 2.0, 7.0]
d = [vals[i % 5] for i in range(60)]
rng = random.Random(2)
b = [rng.uniform(0.5, 1.5) for _ in range(60)]
x, k, conv = minres(diagop(d), b, tol=1e-12)
assert conv and k <= 6
def test_iteration_counts_against_interval_bounds():
# spectrum [c, d] (SPD) or [-d, -c] U [c, d], uniformly filled; b = ones
for kap in (10.0, 100.0):
N = 300
spd = [1 + (kap - 1) * i / (N - 1) for i in range(N)]
ind = [-v for v in spd] + spd
for d, bound_k in ((spd, None), (ind, None)):
b = [1.0] * len(d)
h = []
minres(diagop(d), b, tol=1e-8, history=h, maxiter=5000)
r0 = math.sqrt(len(d))
for k, rk in enumerate(h, start=1):
if d is spd:
q = (math.sqrt(kap) - 1) / (math.sqrt(kap) + 1)
assert rk / r0 <= 2 * q ** k + 1e-12
else:
assert rk / r0 <= interval_bound(1.0, kap, k) + 1e-12
# the indefinite problem takes many more steps than the definite one
k_spd = minres(diagop(spd), [1.0] * N, tol=1e-8)[1]
k_ind = minres(diagop(ind), [1.0] * (2 * N), tol=1e-8, maxiter=5000)[1]
assert k_ind > 3 * k_spd
def test_odd_steps_stagnate_for_symmetric_spectrum():
d = [-(i + 1.0) for i in range(20)] + [i + 1.0 for i in range(20)]
h = []
minres(diagop(d), [1.0] * 40, tol=1e-10, history=h)
r0 = math.sqrt(40)
assert abs(h[0] - r0) <= 1e-12 * r0
for j in range(1, 8):
assert abs(h[2 * j] - h[2 * j - 1]) <= 1e-10 * r0
def test_singular_inconsistent_never_below_least_squares_residual():
h = []
x, k, conv = minres(from_dense([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 0.0]]),
[1.0, 1.0, 1.0], maxiter=30, history=h)
assert not conv and k == 2
assert min(h) >= 1.0 - 1e-12 # dist(b, range A) = 1
assert abs(residual_norm(from_dense([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 0.0]]),
[1.0, 1.0, 1.0], x) - 1.0) <= 1e-12
def test_nonsymmetric_is_not_minimal():
# GMRES solves this in 2 steps; MINRES assumes symmetry and does not.
h = []
x, k, conv = minres(from_dense([[1.0, 1.0], [0.0, 1.0]]), [0.0, 1.0], maxiter=10, history=h)
assert not conv and h[1] > 0.5
def test_recurrence_residual_runs_below_true_residual():
rng = random.Random(9)
a, _ = random_symmetric(rng, 30)
b = [rng.uniform(-1, 1) for _ in range(30)]
h = []
x, k, conv = minres(from_dense(a), b, tol=1e-22, maxiter=200, history=h)
nb = math.sqrt(sum(v * v for v in b))
true = residual_norm(from_dense(a), b, x) / nb
assert min(h) / nb < 1e-18
assert 1e-18 < true < 1e-11
What the next level needs
- the tests did not run in this build