Hager's condition estimator in Python
Estimate the 1-norm condition number cheaply from a few solves with the factorization.
- File
hagers_condition_estimator.py, 153 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/hagers-condition-estimator/. The authors are listed in prov/TXN1K8.jsonl.
"""Hager's estimator of the 1-norm of a matrix B that is known only through
products B x and B^T x, and its use for the condition number
kappa_1(A) = ||A||_1 ||A^{-1}||_1 with B = A^{-1} applied by LU solves.
variant="hager" is Hager's 1984 iteration run to its stopping test.
variant="higham" is Higham's refinement used in LAPACK (xLACON/xLACN2):
at most itmax iterations, a test for a repeated sign vector, and an extra
estimate from the vector b_i = (-1)^i (1 + i/(n-1)).
Works with floats or fractions.Fraction: the probe vectors are exact
Fractions, so the arithmetic follows the type of B. Standard library only.
"""
from fractions import Fraction
class Singular(ArithmeticError):
"""The matrix is singular to working precision (a zero pivot)."""
def norm1(a):
"""The 1-norm of a dense matrix: the largest column sum of absolute values."""
n = len(a[0]) if a else 0
return max((sum(abs(row[j]) for row in a) for j in range(n)), default=0)
def _sign(v):
return -1 if v < 0 else 1
def norm1_estimate(apply, apply_t, n, variant="higham", itmax=5, stats=None):
"""Return (est, v) with est = ||B v||_1 / ||v||_1 <= ||B||_1.
apply(x) returns B x and apply_t(x) returns B^T x. stats (a dict) receives
"products" (the number of calls to apply and apply_t) and "iterations".
"""
if n == 0:
return 0, []
one = Fraction(1)
prods = 0
x = [one / n] * n if n > 1 else [one] # @line L1
est, v = None, None
xi_old = None
it = 0
limit = itmax if variant == "higham" else n + 1
while it < limit: # @line L2
it += 1
y = apply(x) # @line L3
prods += 1
new = sum(abs(t) for t in y) # @line L4
if est is not None and new <= est:
break # no increase: keep the previous estimate
est, v = new, x
xi = [_sign(t) for t in y] # @line L5
if variant == "higham" and xi == xi_old:
break # repeated sign vector: the next step would repeat too
z = apply_t(xi) # @line L6
prods += 1
j = max(range(n), key=lambda i: abs(z[i])) # @line L7
ztx = sum(p * q for p, q in zip(z, x))
if (variant == "hager" or it > 1) and abs(z[j]) <= ztx: # @line L8
break # Hager's test: x is a local maximum
x = [0 * one] * n # @line L9
x[j] = one
xi_old = xi
if variant == "higham" and n > 1: # @line L10
b = [(-1) ** i * (one + one * i / (n - 1)) for i in range(n)]
y = apply(b)
prods += 1
alt = 2 * sum(abs(t) for t in y) / (3 * n)
if alt > est:
est, v = alt, b
if stats is not None:
stats["products"] = prods
stats["iterations"] = it
return est, v # @line L11
def lu_factor(a):
"""LU with partial pivoting, P A = L U, stored in one matrix. Raises Singular."""
n = len(a)
m = [row[:] for row in a]
piv = list(range(n))
for k in range(n):
p = max(range(k, n), key=lambda i: abs(m[i][k]))
if m[p][k] == 0:
raise Singular("zero pivot in column %d" % k)
if p != k:
m[k], m[p] = m[p], m[k]
piv[k], piv[p] = piv[p], piv[k]
for i in range(k + 1, n):
f = m[i][k] / m[k][k]
m[i][k] = f
if f:
rk, ri = m[k], m[i]
for j in range(k + 1, n):
ri[j] -= f * rk[j]
return m, piv
def lu_solve(f, b):
"""Solve A x = b with the factors of lu_factor."""
m, piv = f
n = len(m)
y = [b[p] for p in piv]
for i in range(n):
s = y[i]
row = m[i]
for j in range(i):
s -= row[j] * y[j]
y[i] = s
for i in range(n - 1, -1, -1):
s = y[i]
row = m[i]
for j in range(i + 1, n):
s -= row[j] * y[j]
y[i] = s / row[i]
return y
def lu_solve_t(f, b):
"""Solve A^T x = b with the factors of lu_factor (A^T = U^T L^T P)."""
m, piv = f
n = len(m)
w = list(b)
for i in range(n):
s = w[i]
for j in range(i):
s -= m[j][i] * w[j]
w[i] = s / m[i][i]
for i in range(n - 1, -1, -1):
s = w[i]
for j in range(i + 1, n):
s -= m[j][i] * w[j]
w[i] = s
x = [0 * b[0]] * n if n else []
for k, p in enumerate(piv):
x[p] = w[k]
return x
def cond1_estimate(a, variant="higham", itmax=5, stats=None):
"""Estimate kappa_1(A) = ||A||_1 ||A^{-1}||_1 from one LU factorization.
The estimate is a lower bound of kappa_1(A) in exact arithmetic.
Raises Singular for a matrix with an exactly zero pivot.
"""
f = lu_factor(a)
n = len(a)
est, _ = norm1_estimate(lambda x: lu_solve(f, x), lambda x: lu_solve_t(f, x), n, variant, itmax, stats)
return norm1(a) * est
test_hagers_condition_estimator.py 186 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/hagers-condition-estimator/. The authors are listed in prov/TXN1K8.jsonl.
import random
from fractions import Fraction as F
import pytest
from hagers_condition_estimator import (
Singular,
cond1_estimate,
lu_factor,
lu_solve,
lu_solve_t,
norm1,
norm1_estimate,
)
def ops(b):
n = len(b)
return (lambda x: [sum(r[j] * x[j] for j in range(n)) for r in b],
lambda x: [sum(b[i][j] * x[i] for i in range(n)) for j in range(n)])
def inverse(a):
a = [[F(v) for v in r] for r in a]
f = lu_factor(a)
n = len(a)
cols = [lu_solve(f, [F(int(i == j)) for i in range(n)]) for j in range(n)]
return [[cols[j][i] for j in range(n)] for i in range(n)]
def test_page_example():
a = [[F(2), F(-1), F(0)], [F(-1), F(2), F(-1)], [F(0), F(-1), F(2)]]
st = {}
assert cond1_estimate(a, stats=st) == 8
assert st == {"products": 4, "iterations": 2} # stops on the repeated sign vector
seen = []
f = lu_factor(a)
def ap(x):
y = lu_solve(f, x)
seen.append(sum(abs(t) for t in y))
return y
est, v = norm1_estimate(ap, lambda x: lu_solve_t(f, x), 3, "hager")
assert seen == [F(5, 3), 2] and est == 2 and v == [0, 1, 0]
def test_lower_bound_and_attained_exact_500():
rng = random.Random(1)
exact = total = 0
for t in range(500):
n = rng.randint(1, 9)
b = [[F(rng.randint(-9, 9)) for _ in range(n)] for _ in range(n)]
ap, at = ops(b)
true = norm1(b)
for var in ("hager", "higham"):
st = {}
est, v = norm1_estimate(ap, at, n, var, stats=st)
assert est <= true
assert est == sum(abs(t) for t in ap(v)) / sum(abs(t) for t in v)
if var == "hager":
assert st["iterations"] <= n + 1
else:
assert st["products"] <= 2 * 5 + 1 and st["iterations"] <= 5
total += 1
exact += est == true
assert exact >= 0.5 * total # observed: 57 percent exact for random integer B
def test_estimates_increase_strictly():
rng = random.Random(2)
for _ in range(500):
n = rng.randint(2, 9)
b = [[F(rng.randint(-9, 9)) for _ in range(n)] for _ in range(n)]
ap, at = ops(b)
seen = []
def rec(x):
y = ap(x)
seen.append(sum(abs(t) for t in y))
return y
norm1_estimate(rec, at, n, "hager")
# All but possibly the last product strictly increase the estimate.
assert all(p < q for p, q in zip(seen[:-2], seen[1:-1]))
def test_kappa_estimate_floats_lower_bound_and_quality():
rng = random.Random(3)
worst = 1.0
for t in range(150):
n = rng.randint(2, 14)
a = [[rng.gauss(0, 1) for _ in range(n)] for _ in range(n)]
if t % 4 == 0:
for i in range(n):
a[i][i] *= 1e-8
k = float(norm1(a) * norm1(inverse(a)))
e = cond1_estimate(a)
assert e <= k * (1 + 1e-6)
worst = min(worst, e / k)
assert worst > 0.2 # observed, not a theorem
def test_hilbert():
for n in range(2, 11):
h = [[1.0 / (i + j + 1) for j in range(n)] for i in range(n)]
k = float(norm1(h) * norm1(inverse(h)))
st = {}
e = cond1_estimate(h, stats=st)
assert abs(e / k - 1) < 1e-4
assert st["products"] <= 5
h = [[1.0 / (i + j + 1) for j in range(10)] for i in range(10)]
assert 3.4e13 < cond1_estimate(h) < 3.6e13
def test_pure_hager_can_fail_arbitrarily():
"""B = I - (1 - d) e e^T / n has B e = d e: Hager stops at once with d."""
for n in (2, 5, 9):
for d in (F(1, 10), F(1, 1000), F(1, 10 ** 9)):
b = [[(1 if i == j else 0) - (1 - d) / n for j in range(n)] for i in range(n)]
ap, at = ops(b)
assert norm1_estimate(ap, at, n, "hager")[0] == d
assert norm1_estimate(ap, at, n, "higham")[0] == norm1(b)
b = [[F(1), F(-1)], [F(-1), F(1)]]
ap, at = ops(b)
assert norm1_estimate(ap, at, 2, "hager")[0] == 0
assert norm1_estimate(ap, at, 2, "higham")[0] == 2
def test_alternative_vector_helps():
# Every column has the same norm, signs alternate: the iteration stalls
# but the alternative vector catches the norm.
rng = random.Random(4)
helped = 0
for _ in range(500):
n = rng.randint(2, 7)
b = [[F(rng.randint(-3, 3)) for _ in range(n)] for _ in range(n)]
ap, at = ops(b)
e1 = norm1_estimate(ap, at, n, "higham", itmax=1)[0]
if n > 1:
bb = [(-1) ** i * (1 + F(i, n - 1)) for i in range(n)]
alt = 2 * sum(abs(t) for t in ap(bb)) / (3 * n)
assert e1 >= alt
helped += alt == e1 and alt > 0
assert helped > 0
def test_singular_reports_failure():
with pytest.raises(Singular):
cond1_estimate([[1.0, 2.0], [2.0, 4.0]])
with pytest.raises(Singular):
cond1_estimate([[F(1), F(2), F(3)], [F(4), F(5), F(6)], [F(7), F(8), F(9)]])
with pytest.raises(Singular):
cond1_estimate([[0.0]])
def test_edge_cases():
assert cond1_estimate([[F(5)]]) == 1
d = [F(3), F(-1, 2), F(7), F(1, 4)]
a = [[d[i] if i == j else F(0) for j in range(4)] for i in range(4)]
assert cond1_estimate(a) == 7 / F(1, 4)
# Pivot swap needed: a permutation matrix has kappa 1.
assert cond1_estimate([[0.0, 1.0], [1.0, 0.0]]) == 1.0
# Triangular, with known inverse: ||A^-1||_1 = 2^(n-1) for the unit
# upper triangular matrix with -1 above the diagonal.
n = 12
a = [[F(1) if i == j else F(-1) if j > i else F(0) for j in range(n)] for i in range(n)]
assert cond1_estimate(a) == norm1(a) * 2 ** (n - 1)
# SPD and indefinite symmetric.
assert cond1_estimate([[F(4), F(1)], [F(1), F(3)]]) == norm1(inverse([[4, 1], [1, 3]])) * 5
assert cond1_estimate([[F(1), F(3)], [F(3), F(1)]]) == 4 * norm1(inverse([[1, 3], [3, 1]]))
# Transposed solves are right.
rng = random.Random(5)
for _ in range(50):
n = rng.randint(1, 8)
a = [[F(rng.randint(-5, 5)) for _ in range(n)] for _ in range(n)]
try:
f = lu_factor(a)
except Singular:
continue
b = [F(rng.randint(-5, 5)) for _ in range(n)]
x = lu_solve_t(f, b)
assert [sum(a[i][j] * x[i] for i in range(n)) for j in range(n)] == b
What the next level needs
- the tests did not run in this build