Francis double-shift QR step in Python
One QR step of a real Hessenberg matrix with a shift pair, done in real arithmetic by chasing a 3 by 3 bulge. Repeated with deflation, it is the standard dense real eigensolver.
- File
francis_double_shift_qr.py, 164 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/francis-double-shift-qr/. The authors are listed in prov/ADBYNC.jsonl.
"""Francis implicit double-shift QR algorithm for a real matrix, in real arithmetic.
francis_step(h, l, hi, s, t, ...) applies one double-shift step to the window l..hi of an upper
Hessenberg matrix, with shifts the roots of z^2 - s z + t. schur(a) reduces a real square matrix
to Hessenberg form and iterates to real Schur form: it returns (t, q, eigenvalues, iterations)
with a = q t q^T, t quasi upper triangular with 1x1 and 2x2 diagonal blocks, and the number of
double-shift steps used. Standard library only; matrices are lists of lists of floats.
"""
import math
U = 2.0 ** -53 # unit roundoff of IEEE double precision
def house3(x):
"""Householder vector v and beta with (I - beta v v^T) x = -sign(x_1) ||x|| e_1."""
norm = math.sqrt(sum(c * c for c in x))
if norm == 0.0:
return None, 0.0
v = list(x)
v[0] += norm if x[0] >= 0 else -norm
vv = sum(c * c for c in v)
return v, 2.0 / vv
def francis_step(h, l, hi, s, t, q=None, full=True):
"""One implicit double-shift step on rows and columns l..hi, shifts the roots of z^2 - s z + t."""
n = len(h)
lo = 0 if full else l
top = n if full else hi + 1
x = h[l][l] * h[l][l] + h[l][l + 1] * h[l + 1][l] - s * h[l][l] + t # @line L1
y = h[l + 1][l] * (h[l][l] + h[l + 1][l + 1] - s)
z = h[l + 1][l] * h[l + 2][l + 1]
for k in range(l, hi - 1): # @line L2
v, beta = house3([x, y, z]) # @line L3
if v is not None:
# @line L4
reflect_window(h, q, v, beta, k, max(l, k - 1), min(k + 3, hi), lo, top, n)
x = h[k + 1][k] # @line L5
y = h[k + 2][k]
z = h[k + 3][k] if k < hi - 2 else 0.0
v, beta = house3([x, y]) # @line L6
if v is not None:
reflect_window(h, q, v, beta, hi - 1, hi - 2, hi, lo, top, n) # @line L7
for k in range(l, hi - 1): # entries below the subdiagonal are zero in exact arithmetic
for i in range(k + 2, hi + 1):
h[i][k] = 0.0
def reflect_window(h, q, v, beta, k, col_from, row_to, lo, top, n):
"""Apply P from the left to rows k.. (columns col_from..top-1) and from the right (rows lo..row_to)."""
r = len(v)
for j in range(col_from, top):
d = beta * sum(v[i] * h[k + i][j] for i in range(r))
if d != 0.0:
for i in range(r):
h[k + i][j] -= d * v[i]
for i in range(lo, row_to + 1):
d = beta * sum(h[i][k + c] * v[c] for c in range(r))
if d != 0.0:
for c in range(r):
h[i][k + c] -= d * v[c]
if q is not None:
for i in range(n):
d = beta * sum(q[i][k + c] * v[c] for c in range(r))
if d != 0.0:
for c in range(r):
q[i][k + c] -= d * v[c]
def hessenberg(a, q):
"""Householder reduction to upper Hessenberg form, accumulating into q."""
n = len(a)
h = [row[:] for row in a]
for k in range(n - 2):
v, beta = house3([h[i][k] for i in range(k + 1, n)])
if v is None:
continue
reflect_window(h, q, v, beta, k + 1, k, n - 1, 0, n, n)
for i in range(k + 2, n):
h[i][k] = 0.0
return h
def standardize(h, q, i, full=True):
"""Split a 2x2 block at rows i, i+1 with real eigenvalues by a rotation; return its eigenvalues."""
n = len(h)
a, b, c, d = h[i][i], h[i][i + 1], h[i + 1][i], h[i + 1][i + 1]
p = 0.5 * (a - d)
disc = p * p + b * c
if disc < 0.0:
root = math.sqrt(-disc)
mid = 0.5 * (a + d)
return [complex(mid, root), complex(mid, -root)]
lam = 0.5 * (a + d) + math.copysign(math.sqrt(disc), p)
x1, x2 = b, lam - a
if abs(lam - d) + abs(c) > abs(x1) + abs(x2):
x1, x2 = lam - d, c
r = math.hypot(x1, x2)
if r == 0.0:
return [a, d]
cs, sn = x1 / r, x2 / r
for j in range(i, n if full else i + 2):
u, w = h[i][j], h[i + 1][j]
h[i][j], h[i + 1][j] = cs * u + sn * w, -sn * u + cs * w
for row in range(0 if full else i, i + 2):
u, w = h[row][i], h[row][i + 1]
h[row][i], h[row][i + 1] = cs * u + sn * w, -sn * u + cs * w
if q is not None:
for row in range(n):
u, w = q[row][i], q[row][i + 1]
q[row][i], q[row][i + 1] = cs * u + sn * w, -sn * u + cs * w
h[i + 1][i] = 0.0
return [h[i][i], h[i + 1][i + 1]]
def schur(a, want_q=True, max_its_per_eig=30, trace=None):
"""Real Schur form a = q t q^T by Hessenberg reduction and Francis double-shift steps."""
n = len(a)
q = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)] if want_q else None
h = hessenberg(a, q)
eig = [None] * n
its_total = 0
hi = n - 1
its = 0
while hi >= 0:
l = hi
while l > 0: # @line L8
tol = U * (abs(h[l - 1][l - 1]) + abs(h[l][l]))
if tol == 0.0:
tol = U * max(abs(h[l][l - 1]), U)
if abs(h[l][l - 1]) <= tol:
h[l][l - 1] = 0.0
break
l -= 1
if l == hi: # @line L9
eig[hi] = complex(h[hi][hi])
hi -= 1
its = 0
continue
if l == hi - 1:
pair = standardize(h, q, hi - 1, full=want_q)
eig[hi - 1], eig[hi] = complex(pair[0]), complex(pair[1])
hi -= 2
its = 0
continue
its += 1
its_total += 1
if its > max_its_per_eig:
raise RuntimeError("QR iteration did not converge")
if its % 10 == 0: # exceptional shift, as in LAPACK's xLAHQR
w = abs(h[hi][hi - 1]) + abs(h[hi - 1][hi - 2])
h11 = 0.75 * w + h[hi][hi]
s, t = 2 * h11, h11 * h11 + 0.4375 * w * w
else:
s = h[hi - 1][hi - 1] + h[hi][hi] # @line L10
t = h[hi - 1][hi - 1] * h[hi][hi] - h[hi - 1][hi] * h[hi][hi - 1]
if trace is not None:
trace.append((hi, abs(h[hi][hi - 1]), abs(h[hi - 1][hi - 2])))
francis_step(h, l, hi, s, t, q, full=want_q) # @line L11
return h, q, eig, its_total
test_francis_double_shift_qr.py 208 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/francis-double-shift-qr/. The authors are listed in prov/ADBYNC.jsonl.
import math
import random
from fractions import Fraction
import pytest
from francis_double_shift_qr import U, francis_step, hessenberg, schur
def mm(a, b):
return [[sum(a[i][k] * b[k][j] for k in range(len(b))) for j in range(len(b[0]))] for i in range(len(a))]
def tr(a):
return [list(r) for r in zip(*a)]
def fro(a):
return math.sqrt(sum(x * x for r in a for x in r))
def close_to(found, expected, tol):
pool = list(found)
for lam in expected:
best = min(range(len(pool)), key=lambda i: abs(pool[i] - lam))
assert abs(pool[best] - lam) <= tol, (lam, pool[best])
pool.pop(best)
def check_schur(a, t, q, c=20):
n = len(a)
r = mm(mm(q, t), tr(q))
assert fro([[r[i][j] - a[i][j] for j in range(n)] for i in range(n)]) <= c * n * U * max(fro(a), 1.0)
g = mm(tr(q), q)
assert max(abs(g[i][j] - (i == j)) for i in range(n) for j in range(n)) <= c * n * U
for i in range(n):
for j in range(i - 1):
assert t[i][j] == 0.0
for i in range(n - 2):
assert t[i + 1][i] == 0.0 or t[i + 2][i + 1] == 0.0 # 2x2 blocks do not overlap
for i in range(n - 1):
if t[i + 1][i] != 0.0: # a 2x2 block holds a complex pair
a2, b2, c2, d2 = t[i][i], t[i][i + 1], t[i + 1][i], t[i + 1][i + 1]
assert (a2 - d2) ** 2 / 4 + b2 * c2 < 0
def test_page_example():
a = [[5.0, -7.0, 5.0, -6.0], [1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]
trace = []
t, q, eig, its = schur(a, trace=trace)
assert its == 6
close_to(eig, [3, 2, 1j, -1j], 1e-14)
check_schur(a, t, q)
h32 = [x[2] for x in trace]
expected = [1.0, 0.263987, 0.0728072, 7.34319e-4, 5.10336e-8]
for got, want in zip(h32, expected):
assert got == pytest.approx(want, rel=1e-5)
for e1, e2 in zip(h32[2:], h32[3:]):
assert e2 <= 0.2 * e1 * e1 / min(h32[1:]) + 1e-15 # quadratic in the last steps
def test_random_schur_and_iteration_count():
rng = random.Random(3)
n = 40
a = [[rng.gauss(0, 1) for _ in range(n)] for _ in range(n)]
t, q, eig, its = schur(a)
check_schur(a, t, q)
assert abs(sum(eig).real - sum(a[i][i] for i in range(n))) <= 1e-12 * n
assert all(abs(e.imag) == 0 or any(abs(f - e.conjugate()) <= 1e-12 for f in eig) for e in eig)
assert its / n <= 4 # observed about 2 double steps per eigenvalue
t2, q2, eig2, its2 = schur(a, want_q=False)
assert q2 is None
close_to(eig2, eig, 1e-10)
def det_frac(m):
m = [row[:] for row in m]
n = len(m)
d = Fraction(1)
for c in range(n):
p = next((r for r in range(c, n) if m[r][c] != 0), None)
if p is None:
return Fraction(0)
if p != c:
m[c], m[p] = m[p], m[c]
d = -d
d *= m[c][c]
for r in range(c + 1, n):
f = m[r][c] / m[c][c]
m[r] = [x - f * y for x, y in zip(m[r], m[c])]
return d
def test_exact_trace_and_determinant():
ai = [[2, -1, 3, 0, 1], [4, 1, -2, 5, 0], [0, 3, 1, -1, 2], [1, 0, 2, 3, -4], [-2, 1, 0, 1, 1]]
exact_det = det_frac([[Fraction(x) for x in r] for r in ai])
a = [[float(x) for x in r] for r in ai]
t, q, eig, its = schur(a)
check_schur(a, t, q)
prod = 1.0 + 0j
for e in eig:
prod *= e
assert abs(prod - float(exact_det)) <= 1e-11 * abs(float(exact_det))
assert abs(sum(eig) - sum(ai[i][i] for i in range(5))) <= 1e-12
def test_known_spectra():
n = 15
toe = [[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)]
t, q, eig, _ = schur(toe)
check_schur(toe, t, q)
close_to(eig, [2 - 2 * math.cos(k * math.pi / (n + 1)) for k in range(1, n + 1)], 1e-13)
s = [[1.0, 1.0, 0.0], [0.0, 1.0, 1.0], [1.0, 0.0, 1.0]]
a = [[3.0, 1.0, -1.0], [-1.0, 5.0, 1.0], [-2.0, 2.0, 4.0]] # S diag(2,4,6) S^-1
assert mm(a, s) == [[2.0, 4.0, 0.0], [0.0, 4.0, 6.0], [2.0, 0.0, 6.0]]
t, q, eig, _ = schur(a)
close_to(eig, [2, 4, 6], 1e-13)
coeffs = [1.0, -21.0, 175.0, -735.0, 1624.0, -1764.0, 720.0] # (z-1)...(z-6)
comp = [[0.0] * 6 for _ in range(6)]
comp[0] = [-c for c in coeffs[1:]]
for i in range(1, 6):
comp[i][i - 1] = 1.0
t, q, eig, _ = schur(comp)
check_schur(comp, t, q, c=200)
close_to(eig, [1, 2, 3, 4, 5, 6], 1e-8)
def test_trivial_and_structured_cases():
for a, expected in (([[4.0]], [4]), ([[0.0] * 3 for _ in range(3)], [0, 0, 0]),
([[3.0, 0, 0], [0, 1.0, 0], [0, 0, 2.0]], [3, 1, 2]),
([[1.0, 5.0, -2.0], [0, 4.0, 7.0], [0, 0, -3.0]], [1, 4, -3])):
t, q, eig, its = schur(a)
assert its == 0
close_to(eig, expected, 0.0)
t, q, eig, _ = schur([[1.0, -2.0], [3.0, 1.0]])
close_to(eig, [1 + math.sqrt(6) * 1j, 1 - math.sqrt(6) * 1j], 1e-15)
a = [[2.0, 1.0], [1.0, 2.0]]
t, q, eig, _ = schur(a)
check_schur(a, t, q)
assert t[1][0] == 0.0
close_to(eig, [3, 1], 1e-15)
def test_complex_pairs_repeated_jordan_nearly_equal():
rng = random.Random(8)
n = 8
v = [rng.uniform(-1, 1) for _ in range(n)]
nv = math.sqrt(sum(x * x for x in v))
p = [[(i == j) - 2 * v[i] * v[j] / nv ** 2 for j in range(n)] for i in range(n)]
def hide(b):
return mm(mm(p, b), p)
base = [[0.0] * n for _ in range(n)]
base[0][0], base[0][1], base[1][0], base[1][1] = 1.0, 2.0, -2.0, 1.0
base[2][2], base[2][3], base[3][2], base[3][3] = -3.0, 0.5, -0.5, -3.0
for i, d in zip(range(4, 8), (2.0, 2.0, 2.0, 0.5)):
base[i][i] = d
a = hide(base)
t, q, eig, _ = schur(a)
check_schur(a, t, q)
close_to(eig, [1 + 2j, 1 - 2j, -3 + 0.5j, -3 - 0.5j, 2, 2, 2, 0.5], 1e-7)
jb = [[0.0] * n for _ in range(n)]
for i in range(4):
jb[i][i] = 1.0
if i < 3:
jb[i][i + 1] = 1.0
for i, d in zip(range(4, 8), (3.0, 1.0 + 1e-10, 1.0 - 1e-10, -1.0)):
jb[i][i] = d
a = hide(jb)
t, q, eig, _ = schur(a)
check_schur(a, t, q)
near_one = [e for e in eig if abs(e - 1) < 0.1]
assert len(near_one) == 6
assert all(abs(e - 1) <= 1e-3 for e in near_one) # about u^(1/4) for the Jordan block
close_to([e for e in eig if abs(e - 1) >= 0.1], [3, -1], 1e-12)
def test_exceptional_shift_on_cyclic_permutation():
n = 6
c = [[0.0] * n for _ in range(n)]
for i in range(1, n):
c[i][i - 1] = 1.0
c[0][n - 1] = 1.0
h = [r[:] for r in c]
francis_step(h, 0, n - 1, 0.0, 0.0) # standard shifts of the trailing block are 0, 0
# the step only changes signs (a diagonal similarity): no progress toward deflation
assert max(abs(abs(h[i][j]) - c[i][j]) for i in range(n) for j in range(n)) <= 1e-15
t, q, eig, its = schur(c)
assert its > 10 # converges only after an exceptional shift
check_schur(c, t, q)
close_to(eig, [complex(math.cos(2 * math.pi * k / n), math.sin(2 * math.pi * k / n)) for k in range(n)], 1e-13)
def test_hessenberg_reduction_keeps_spectrum():
rng = random.Random(1)
n = 10
a = [[rng.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
q = [[float(i == j) for j in range(n)] for i in range(n)]
h = hessenberg(a, q)
for i in range(n):
for j in range(i - 1):
assert h[i][j] == 0.0
r = mm(mm(q, h), tr(q))
assert fro([[r[i][j] - a[i][j] for j in range(n)] for i in range(n)]) <= 20 * n * U * fro(a)
What the next level needs
- the tests did not run in this build