Biconjugate gradient method in Python
Extends conjugate gradients to nonsymmetric systems using a shadow sequence.
- File
biconjugate_gradient.py, 115 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/biconjugate-gradient/. The authors are listed in prov/4SSPDC.jsonl.
"""The biconjugate gradient method (BiCG) for A x = b, A nonsymmetric.
Works with floats or with fractions.Fraction (then the run is exact).
A is a dense list of rows or a CSR object; BiCG needs products with A
and with A^T. Standard library only.
"""
class CSR:
"""Compressed sparse row matrix: row i holds data[ptr[i]:ptr[i+1]]."""
def __init__(self, n, ptr, ind, data):
self.n, self.ptr, self.ind, self.data = n, ptr, ind, data
@classmethod
def from_dense(cls, a):
ptr, ind, data = [0], [], []
for row in a:
for j, v in enumerate(row):
if v != 0:
ind.append(j)
data.append(v)
ptr.append(len(ind))
return cls(len(a), ptr, ind, data)
def matvec(self, x):
out = []
for i in range(self.n):
acc = 0 * x[0] if x else 0
for p in range(self.ptr[i], self.ptr[i + 1]):
acc += self.data[p] * x[self.ind[p]]
out.append(acc)
return out
def rmatvec(self, y):
"""A^T y, scattering row i of A times y_i."""
out = [0 * y[0] for _ in range(self.n)] if y else []
for i in range(self.n):
yi = y[i]
for p in range(self.ptr[i], self.ptr[i + 1]):
out[self.ind[p]] += self.data[p] * yi
return out
def _mv(a, x):
if isinstance(a, CSR):
return a.matvec(x)
return [sum((aij * xj for aij, xj in zip(row, x)), 0 * x[0]) for row in a]
def _rmv(a, y):
if isinstance(a, CSR):
return a.rmatvec(y)
n = len(y)
return [sum((a[i][j] * y[i] for i in range(n)), 0 * y[0]) for j in range(n)]
def _dot(u, v):
return sum((p * q for p, q in zip(u, v)), 0 * u[0])
class BiCGResult:
def __init__(self, x, status, iterations, history):
self.x = x
self.status = status # "converged", "maxit", "rho-breakdown" or "pivot-breakdown"
self.converged = status == "converged"
self.iterations = iterations
self.history = history # ||r_k||_2^2 for k = 0, 1, ...
def bicg(a, b, x0=None, shadow=None, tol=1e-10, maxit=None, trace=None):
"""Run BiCG until ||r_k||_2 <= tol ||b||_2.
shadow is the shadow residual r~_0 (default r_0). trace, if a list,
receives (r_k, r~_k, p_k, p~_k) for every k, for tests of
biorthogonality.
"""
n = len(b)
maxit = 2 * n if maxit is None else maxit
zero = 0 * b[0]
x = [zero] * n if x0 is None else list(x0)
r = [bi - t for bi, t in zip(b, _mv(a, x))] # @line L1
rt = list(r) if shadow is None else list(shadow)
p, pt = list(r), list(rt) # @line L2
rho = _dot(rt, r)
target = tol * tol * _dot(b, b)
history = [_dot(r, r)]
for k in range(maxit + 1): # @line L3
if trace is not None:
trace.append((list(r), list(rt), list(p), list(pt)))
if history[-1] <= target: # @line L4
return BiCGResult(x, "converged", k, history)
if k == maxit:
break
if rho == 0:
return BiCGResult(x, "rho-breakdown", k, history)
q = _mv(a, p) # @line L5
qt = _rmv(a, pt)
sigma = _dot(pt, q)
if sigma == 0:
return BiCGResult(x, "pivot-breakdown", k, history)
alpha = rho / sigma # @line L6
x = [xi + alpha * pi for xi, pi in zip(x, p)] # @line L7
r = [ri - alpha * qi for ri, qi in zip(r, q)]
rt = [ri - alpha * qi for ri, qi in zip(rt, qt)] # @line L8
rho_new = _dot(rt, r) # @line L9
beta = rho_new / rho
rho = rho_new
p = [ri + beta * pi for ri, pi in zip(r, p)] # @line L10
pt = [ri + beta * pi for ri, pi in zip(rt, pt)]
history.append(_dot(r, r))
return BiCGResult(x, "maxit", maxit, history)
test_biconjugate_gradient.py 147 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/biconjugate-gradient/. The authors are listed in prov/4SSPDC.jsonl.
import math
import random
from fractions import Fraction as F
from biconjugate_gradient import CSR, bicg
U = 2.0**-53
def fr(m):
return [[F(t) for t in row] for row in m]
def mv(a, x):
return [sum(p * q for p, q in zip(row, x)) for row in a]
def dot(u, v):
return sum(p * q for p, q in zip(u, v))
def test_one_by_one():
res = bicg([[F(5)]], [F(10)], tol=0)
assert res.converged and res.iterations == 1 and res.x == [F(2)]
def test_page_example_exact():
a = fr([[2, 1, 0], [-1, 2, 1], [0, -1, 2]])
b = [F(1)] * 3
tr = []
res = bicg(a, b, tol=0, trace=tr)
assert res.converged and res.iterations == 3
assert res.x == [F(1, 3), F(1, 3), F(2, 3)]
assert tr[1][0] == [F(-1, 2), 0, F(1, 2)]
assert tr[2][0] == [F(1, 7), F(-2, 7), F(1, 7)]
assert res.history == [3, F(1, 2), F(6, 49), 0]
def test_random_exact_termination_and_biorthogonality():
rng = random.Random(4242)
breakdowns = 0
for trial in range(500):
n = rng.randint(1, 5)
a = fr([[rng.randint(-3, 3) for _ in range(n)] for _ in range(n)])
for i in range(n):
a[i][i] += 6
b = [F(rng.randint(-4, 4)) for _ in range(n)]
if all(t == 0 for t in b):
b[0] = F(1)
shadow = [F(rng.randint(-4, 4)) for _ in range(n)] if trial % 2 else None
if shadow is not None and all(t == 0 for t in shadow):
shadow[0] = F(1)
tr = []
res = bicg(a, b, shadow=shadow, tol=0, trace=tr)
if not res.converged:
assert res.status in ("rho-breakdown", "pivot-breakdown")
breakdowns += 1
continue
assert res.iterations <= n
assert mv(a, res.x) == b
k = len(tr)
for i in range(k):
for j in range(k):
if i != j:
assert dot(tr[i][1], tr[j][0]) == 0 # r~_i . r_j
assert dot(tr[i][3], mv(a, tr[j][2])) == 0 # p~_i . A p_j
assert breakdowns < 50
def test_symmetric_shadow_equals_residual():
# For symmetric A and r~_0 = r_0, r~_k = r_k: BiCG is CG done twice.
a = fr([[4, 1, 0, 0], [1, 3, 1, 0], [0, 1, 2, 1], [0, 0, 1, 5]])
b = [F(1), F(2), F(3), F(4)]
tr = []
res = bicg(a, b, tol=0, trace=tr)
assert res.converged and res.iterations == 4
for r, rt, p, pt in tr:
assert r == rt and p == pt
def test_pivot_and_rho_breakdown():
res = bicg(fr([[0, 1], [1, 0]]), [F(1), F(0)], tol=0)
assert res.status == "pivot-breakdown"
res = bicg(fr([[1, 0], [0, 1]]), [F(1), F(0)], shadow=[F(0), F(1)], tol=0)
assert res.status == "rho-breakdown"
def test_singular_reports_failure():
res = bicg(fr([[1, 0], [0, 0]]), [F(0), F(1)], tol=0)
assert not res.converged
res = bicg([[1.0, 0.0], [0.0, 0.0]], [1.0, 1.0])
assert not res.converged
def convection_diffusion(n, c):
rows = []
for i in range(n):
row = [0.0] * n
row[i] = 2.0
if i > 0:
row[i - 1] = -1.0 - c
if i < n - 1:
row[i + 1] = -1.0 + c
rows.append(row)
return rows
def test_float_convection_diffusion_and_nonmonotone():
n = 150
rows = convection_diffusion(n, 0.3)
a = CSR.from_dense(rows)
rng = random.Random(7)
b = [rng.uniform(-1, 1) for _ in range(n)]
res = bicg(a, b, tol=1e-10, maxit=2000)
assert res.converged
r = [bi - t for bi, t in zip(b, mv(rows, res.x))]
assert math.sqrt(dot(r, r)) <= 1e-9 * math.sqrt(dot(b, b))
h = res.history
assert any(h[k + 1] > h[k] for k in range(len(h) - 1)) # not monotone
def test_spd_matches_cg_iteration_count_bound():
n, kappa = 200, 50.0
d = [1.0 + (kappa - 1.0) * i / (n - 1) for i in range(n)]
a = CSR(n, list(range(n + 1)), list(range(n)), d)
rng = random.Random(9)
b = [rng.uniform(-1, 1) for _ in range(n)]
res = bicg(a, b, tol=1e-8, maxit=500)
q = (math.sqrt(kappa) - 1) / (math.sqrt(kappa) + 1)
# CG bound in the A-norm turns into a residual bound with sqrt(kappa).
kbound = math.ceil(math.log(2 * math.sqrt(kappa) / 1e-8) / math.log(1 / q))
assert res.converged and res.iterations <= kbound
def test_hilbert_backward_error():
n = 6
h = [[1.0 / (i + j + 1) for j in range(n)] for i in range(n)]
b = [sum(row) for row in h]
res = bicg(h, b, tol=1e-13, maxit=300)
assert res.converged
r = [bi - t for bi, t in zip(b, mv(h, res.x))]
hn = max(sum(map(abs, row)) for row in h)
eta = max(map(abs, r)) / (hn * max(map(abs, res.x)) + max(map(abs, b)))
assert eta < 1e-12
What the next level needs
- the tests did not run in this build