QZ algorithm in Python
Computes the generalized Schur form Q* A Z = S, Q* B Z = T of a pencil with unitary Q, Z and triangular S, T, by QR-type steps that never invert B. Backward stable; handles singular B.
- File
qz_algorithm.py, 220 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/qz-algorithm/. The authors are listed in prov/WCHJ4Y.jsonl.
"""The QZ algorithm for the generalized eigenproblem A x = lambda B x.
qz(a, b) computes unitary Q and Z and upper triangular S and T with
Q* A Z = S and Q* B Z = T, the complex generalized Schur form. The
generalized eigenvalues are the pairs (alpha_i, beta_i) = (s_ii, t_ii):
lambda_i = alpha_i / beta_i when beta_i != 0, and an infinite eigenvalue when
beta_i = 0. The method first reduces (A, B) to Hessenberg-triangular form by
Givens rotations, then runs single-shift complex QZ steps with a shift from
the trailing 2 by 2 pencil, and deflates infinite eigenvalues by chasing a
zero diagonal entry of T to the top of the active block.
Matrices are lists of rows of Python numbers (real or complex); all work is in
complex arithmetic. Standard library only.
"""
import cmath
import math
from typing import NamedTuple
EPS = 2.0**-52
class QZResult(NamedTuple):
s: list # upper triangular, Q* A Z
t: list # upper triangular, Q* B Z
q: list # unitary
z: list # unitary
alpha: list # s_ii
beta: list # t_ii
iterations: int # the number of QZ steps
converged: bool
history: list # (hi, |h_(hi,hi-1)|) after each QZ step, for rate measurements
def identity(n):
return [[complex(i == j) for j in range(n)] for i in range(n)]
def frob(m):
return math.sqrt(sum(abs(x) ** 2 for row in m for x in row))
def givens(f, g):
"""(c, s, r) with c real, [[c, s], [-conj(s), c]] (f, g)^T = (r, 0)^T."""
if g == 0:
return 1.0, 0j, f
if f == 0:
return 0.0, 1 + 0j, g
af = abs(f)
nrm = math.hypot(af, abs(g))
ph = f / af
return af / nrm, ph * g.conjugate() / nrm, ph * nrm
def rot_rows(m, i, j, c, s, cols):
"""Rows i, j of m <- G (rows i, j), G = [[c, s], [-conj(s), c]]."""
sc = s.conjugate()
ri, rj = m[i], m[j]
for k in cols:
x, y = ri[k], rj[k]
ri[k] = c * x + s * y
rj[k] = -sc * x + c * y
def rot_cols(m, i, j, c, s, rows):
"""Columns i, j of m <- (columns i, j) G, G = [[c, s], [-conj(s), c]]."""
sc = s.conjugate()
for k in rows:
row = m[k]
x, y = row[i], row[j]
row[i] = c * x - sc * y
row[j] = s * x + c * y
class _Pencil:
"""H, T, Q, Z with Q* A Z = H and Q* B Z = T; every rotation keeps this invariant."""
def __init__(self, a, b):
n = len(a)
self.n = n
self.h = [[complex(x) for x in row] for row in a]
self.t = [[complex(x) for x in row] for row in b]
self.q = identity(n)
self.z = identity(n)
self.all = range(n)
def left(self, i, j, f, g):
"""Rotate rows i, j of H and T so that the pair (f, g) becomes (r, 0); Q <- Q G*."""
c, s, _ = givens(f, g)
rot_rows(self.h, i, j, c, s, self.all)
rot_rows(self.t, i, j, c, s, self.all)
# Q G*: columns i, j of Q times [[c, -s], [conj(s), c]]
sc = s.conjugate()
for row in self.q:
x, y = row[i], row[j]
row[i] = c * x + sc * y
row[j] = -s * x + c * y
def right(self, i, j, x, y):
"""Rotate columns i < j of H and T so that a row with entries (x, y) in them becomes (0, r); Z <- Z R."""
c, s, _ = givens(y, x)
rot_cols(self.h, i, j, c, s, self.all)
rot_cols(self.t, i, j, c, s, self.all)
rot_cols(self.z, i, j, c, s, self.all)
def hessenberg_triangular(p):
"""Reduce H to upper Hessenberg and T to upper triangular form (GVL Algorithm 7.7.1, with Givens)."""
n, h, t = p.n, p.h, p.t
for j in range(n - 1): # T upper triangular by left rotations (a Givens QR of B)
for i in range(n - 1, j, -1):
if t[i][j] != 0:
p.left(i - 1, i, t[i - 1][j], t[i][j])
t[i][j] = 0j
for j in range(n - 2):
for i in range(n - 1, j + 1, -1):
if h[i][j] == 0:
continue
p.left(i - 1, i, h[i - 1][j], h[i][j]) # zero h_ij; fills t_(i, i-1)
h[i][j] = 0j
p.right(i - 1, i, t[i][i - 1], t[i][i]) # zero t_(i, i-1); keeps h_ij = 0
t[i][i - 1] = 0j
def pencil_shift(h, t, k):
"""The eigenvalue of the 2 by 2 pencil in rows and columns k-1, k nearer to h_kk / t_kk."""
h11, h12, h21, h22 = h[k - 1][k - 1], h[k - 1][k], h[k][k - 1], h[k][k]
t11, t12, t22 = t[k - 1][k - 1], t[k - 1][k], t[k][k]
# det(H2 - x T2) = qa x^2 + qb x + qc
qa = t11 * t22
qb = -(h11 * t22 + h22 * t11 - h21 * t12)
qc = h11 * h22 - h12 * h21
target = h22 / t22
if qa == 0:
return target
disc = cmath.sqrt(qb * qb - 4 * qa * qc)
den = -qb + disc if abs(-qb + disc) >= abs(-qb - disc) else -qb - disc
r1 = den / (2 * qa)
r2 = (2 * qc) / den if den != 0 else r1
return r1 if abs(r1 - target) <= abs(r2 - target) else r2
def qz(a, b, maxit_per_eig=30):
"""Return a QZResult for the pencil (A, B) of square matrices of equal order."""
n = len(a)
if len(b) != n or any(len(r) != n for r in a) or any(len(r) != n for r in b):
raise ValueError("A and B must be square of the same order")
p = _Pencil(a, b)
hessenberg_triangular(p) # @line L1
h, t = p.h, p.t
hnorm = max(frob(h), 1e-300)
tnorm = max(frob(t), 1e-300)
iters, since, history = 0, 0, []
maxit = maxit_per_eig * max(n, 1)
hi = n - 1 # @line L2
converged = True
while hi > 0: # @line L3
for k in range(1, hi + 1): # @line L4
scale = abs(h[k][k]) + abs(h[k - 1][k - 1]) or hnorm
if abs(h[k][k - 1]) <= EPS * scale:
h[k][k - 1] = 0j
if h[hi][hi - 1] == 0: # @line L5
hi -= 1
since = 0
continue
lo = hi - 1 # @line L6
while lo > 0 and h[lo][lo - 1] != 0:
lo -= 1
for k in range(lo, hi + 1): # negligible diagonal entries of T in the active block
if abs(t[k][k]) <= EPS * tnorm:
t[k][k] = 0j
zero = next((k for k in range(lo, hi + 1) if t[k][k] == 0), None)
if zero is not None: # @line L7
for j in range(zero, lo, -1): # chase the zero of T up to t_(lo, lo)
p.right(j - 1, j, t[j - 1][j - 1], t[j - 1][j])
t[j - 1][j - 1] = 0j
if j + 1 <= hi:
p.left(j, j + 1, h[j][j - 1], h[j + 1][j - 1])
h[j + 1][j - 1] = 0j
t[j + 1][j] = 0j
p.left(lo, lo + 1, h[lo][lo], h[lo + 1][lo]) # deflate an infinite eigenvalue at the top
h[lo + 1][lo] = 0j
t[lo + 1][lo] = 0j
since = 0
continue
if iters >= maxit:
converged = False
break
iters += 1
since += 1
if since % 11 == 0: # an exceptional shift breaks rare cycles
mu = h[hi][hi] / t[hi][hi] + 0.75 * abs(h[hi][hi - 1]) / abs(t[hi][hi])
else:
mu = pencil_shift(h, t, hi) # @line L8
p.left(lo, lo + 1, h[lo][lo] - mu * t[lo][lo], h[lo + 1][lo]) # @line L9
for k in range(lo, hi): # @line L10
p.right(k, k + 1, t[k + 1][k], t[k + 1][k + 1]) # zero the bulge t_(k+1, k)
t[k + 1][k] = 0j
if k + 2 <= hi:
p.left(k + 1, k + 2, h[k + 1][k], h[k + 2][k]) # zero the bulge h_(k+2, k)
h[k + 2][k] = 0j
history.append((hi, abs(h[hi][hi - 1])))
for i in range(n): # clean the strictly lower parts, which are zero up to rounding
for j in range(i):
h[i][j] = 0j
t[i][j] = 0j
alpha = [h[i][i] for i in range(n)]
beta = [t[i][i] for i in range(n)]
return QZResult(h, t, p.q, p.z, alpha, beta, iters, converged, history) # @line L11
def eigenvalues(res):
"""The generalized eigenvalues alpha / beta, with math.inf for beta = 0."""
out = []
for a, b in zip(res.alpha, res.beta):
out.append(a / b if b != 0 else complex(math.inf))
return out
test_qz_algorithm.py 249 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/qz-algorithm/. The authors are listed in prov/WCHJ4Y.jsonl.
import math
import random
from fractions import Fraction
import pytest
from qz_algorithm import eigenvalues, frob, qz
U = 2.0**-53
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 ctrans(a):
return [[a[j][i].conjugate() for j in range(len(a))] for i in range(len(a[0]))]
def rel_residual(a, q, s, z):
r = mm(mm(q, s), ctrans(z))
n = len(a)
return math.sqrt(sum(abs(r[i][j] - a[i][j]) ** 2 for i in range(n) for j in range(n))) / max(frob(a), 1e-300)
def unitarity(q):
g = mm(ctrans(q), q)
n = len(q)
return max(abs(g[i][j] - (i == j)) for i in range(n) for j in range(n))
def det_fraction(m):
"""Exact determinant by fraction Gaussian elimination."""
a = [[Fraction(x) for x in row] for row in m]
n = len(a)
det = Fraction(1)
for j in range(n):
piv = next((i for i in range(j, n) if a[i][j] != 0), None)
if piv is None:
return Fraction(0)
if piv != j:
a[j], a[piv] = a[piv], a[j]
det = -det
det *= a[j][j]
for i in range(j + 1, n):
f = a[i][j] / a[j][j]
for k in range(j, n):
a[i][k] -= f * a[j][k]
return det
def pencil_poly(a, b):
"""Exact coefficients c_0..c_n of p(x) = det(A - x B), by interpolation at x = 0..n."""
n = len(a)
af = [[Fraction(x) for x in row] for row in a]
bf = [[Fraction(x) for x in row] for row in b]
xs = list(range(n + 1))
ys = [det_fraction([[af[i][j] - x * bf[i][j] for j in range(n)] for i in range(n)]) for x in xs]
coeffs = [Fraction(0)] * (n + 1)
for k, xk in enumerate(xs): # Lagrange basis polynomials, expanded
basis = [Fraction(1)]
denom = Fraction(1)
for m, xm in enumerate(xs):
if m == k:
continue
basis = [Fraction(0)] + basis
for i in range(len(basis) - 1):
basis[i] -= xm * basis[i + 1]
denom *= xk - xm
for i in range(n + 1):
coeffs[i] += ys[k] * basis[i] / denom
return coeffs
def poly_from_roots(roots):
"""Coefficients c_0..c_d of prod (x - r), lowest first."""
c = [1 + 0j]
for r in roots:
c = [0j] + c
for i in range(len(c) - 1):
c[i] -= r * c[i + 1]
return c
def orthogonal_rational(n, rng, k=3):
q = [[Fraction(int(i == j)) for j in range(n)] for i in range(n)]
for _ in range(k):
v = [Fraction(rng.randint(-4, 4)) for _ in range(n)]
if all(t == 0 for t in v):
v[0] = Fraction(1)
vv = sum(t * t for t in v)
for row in q:
s = sum(r * t for r, t in zip(row, v)) * 2 / vv
for j in range(n):
row[j] -= s * v[j]
return q
def check_schur(a, b, r, c=10):
n = len(a)
assert r.converged
assert rel_residual(a, r.q, r.s, r.z) <= c * max(n, 1) * U
if frob(b) > 0:
assert rel_residual(b, r.q, r.t, r.z) <= c * max(n, 1) * U
assert unitarity(r.q) <= c * max(n, 1) * U
assert unitarity(r.z) <= c * max(n, 1) * U
for i in range(n):
for j in range(i):
assert r.s[i][j] == 0 and r.t[i][j] == 0
def test_backward_stability_many_seeds():
for seed in range(40):
rng = random.Random(seed)
n = rng.randint(1, 9)
a = [[rng.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
b = [[rng.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
if seed % 4 == 0:
for i in range(n):
b[i][0] = 0.0 # B singular: at least one infinite eigenvalue
check_schur(a, b, qz(a, b))
def test_exact_pencil_polynomial_and_determinants():
for seed in range(16):
rng = random.Random(500 + seed)
n = rng.randint(2, 6)
a = [[float(rng.randint(-5, 5)) for _ in range(n)] for _ in range(n)]
b = [[float(rng.randint(-5, 5)) for _ in range(n)] for _ in range(n)]
if seed % 3 == 0:
b[n - 1] = [0.0] * n # rank-deficient B
c = pencil_poly(a, b)
if all(x == 0 for x in c):
continue # a singular pencil; covered separately
deg = max(i for i, x in enumerate(c) if x != 0)
r = qz(a, b)
check_schur(a, b, r)
scale_b = max(abs(x) for x in r.beta)
finite = [x for x, beta in zip(eigenvalues(r), r.beta) if abs(beta) > 1e-12 * scale_b]
assert len(finite) == deg # n - deg infinite eigenvalues
mine = poly_from_roots(finite)
lead = c[deg]
exact = [complex(float(x / lead)) for x in c[: deg + 1]]
size = max(abs(x) for x in exact)
for u_, v_ in zip(mine, exact):
assert abs(u_ - v_) <= 1e-8 * size
# |det A| = |prod alpha| and |det B| = |prod beta|, since Q and Z are unitary
da, db = det_fraction(a), det_fraction(b)
pa = abs(math.prod(r.alpha))
pb = abs(math.prod(r.beta))
assert abs(pa - abs(float(da))) <= 1e-10 * max(1.0, abs(float(da)))
assert abs(pb - abs(float(db))) <= 1e-10 * max(1.0, abs(float(db)))
def test_known_spectra_with_infinite_eigenvalues():
for seed in range(10):
rng = random.Random(900 + seed)
n = rng.randint(2, 7)
uq, vq = orthogonal_rational(n, rng), orthogonal_rational(n, rng)
da = [Fraction(rng.randint(1, 9)) * rng.choice((-1, 1)) for _ in range(n)]
db = [Fraction(rng.randint(1, 5)) for _ in range(n)]
ninf = rng.randint(0, 2)
for i in range(min(ninf, n)):
db[i] = Fraction(0)
sa = [[da[i] if i == j else (Fraction(rng.randint(-3, 3)) if j > i else Fraction(0)) for j in range(n)] for i in range(n)]
sb = [[db[i] if i == j else (Fraction(rng.randint(-3, 3)) if j > i else Fraction(0)) for j in range(n)] for i in range(n)]
a = [[float(x) for x in row] for row in mm(mm(uq, sa), vq)]
b = [[float(x) for x in row] for row in mm(mm(uq, sb), vq)]
r = qz(a, b)
check_schur(a, b, r)
scale = max(abs(x) for x in r.beta)
got_inf = sum(1 for x in r.beta if abs(x) <= 1e-12 * scale)
assert got_inf == min(ninf, n)
want = sorted(float(x / y) for x, y in zip(da, db) if y != 0)
got = sorted(x.real for x, beta in zip(eigenvalues(r), r.beta) if abs(beta) > 1e-12 * scale)
assert got == pytest.approx(want, abs=1e-8 * max(abs(w) for w in want))
def test_identity_b_gives_ordinary_eigenvalues():
a = [[0.0, -1.0], [1.0, 0.0]]
r = qz(a, [[1.0, 0.0], [0.0, 1.0]])
ev = sorted(eigenvalues(r), key=lambda z: z.imag)
assert ev == [pytest.approx(-1j, abs=1e-14), pytest.approx(1j, abs=1e-14)]
rng = random.Random(3)
n = 6
q = orthogonal_rational(n, rng)
lam = [Fraction(k) for k in (5, 3, 2, -1, -4, 7)]
a = [[float(sum(q[i][k] * lam[k] * q[j][k] for k in range(n))) for j in range(n)] for i in range(n)]
r = qz(a, [[float(i == j) for j in range(n)] for i in range(n)])
assert sorted(x.real for x in eigenvalues(r)) == pytest.approx(sorted(float(x) for x in lam), abs=1e-12)
def test_quadratic_convergence_of_the_subdiagonal():
ratios = []
steps = []
for seed in range(30):
rng = random.Random(seed)
n = rng.randint(4, 9)
a = [[rng.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
b = [[rng.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
r = qz(a, b)
steps.append(r.iterations / n)
h = r.history
for (k1, e1), (k2, e2) in zip(h, h[1:]):
if k1 == k2 and 1e-12 <= e1 <= 1e-4:
ratios.append(e2 / e1**2)
assert ratios
ratios.sort()
assert ratios[len(ratios) // 2] <= 10.0 # e_(k+1) = O(e_k^2): quadratic in the typical case
assert sum(steps) / len(steps) <= 4.0 # typically 2 to 3 QZ steps per eigenvalue
def test_page_example():
a = [[1.0, 2, 3], [4, 5, 6], [7, 8, 10]]
b = [[1.0, 0, 0], [0, 1.0, 0], [0, 0, 0.0]]
assert pencil_poly(a, b) == [-3, 9, 10, 0]
r = qz(a, b)
check_schur(a, b, r)
ev = eigenvalues(r)
assert math.isinf(ev[0].real)
finite = sorted(x.real for x in ev[1:])
assert finite == pytest.approx([(-9 - 201**0.5) / 20, (-9 + 201**0.5) / 20], abs=1e-14)
assert r.iterations == 1
rng = random.Random(1)
n = 6
a = [[rng.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
b = [[rng.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
r = qz(a, b)
assert r.iterations == 16
tail = [e for k, e in r.history if k == 5][3:]
assert [float("%.1e" % e) for e in tail] == [4.1e-04, 9.7e-09, 1.5e-18]
def test_edge_cases():
r = qz([[3.0]], [[2.0]])
assert r.alpha[0] / r.beta[0] == 1.5 and r.iterations == 0
tri_a = [[1.0, 2.0, 3.0], [0.0, 4.0, 5.0], [0.0, 0.0, 6.0]]
tri_b = [[1.0, 1.0, 1.0], [0.0, 2.0, 1.0], [0.0, 0.0, 3.0]]
r = qz(tri_a, tri_b)
assert r.iterations == 0 and [x.real for x in eigenvalues(r)] == [1.0, 2.0, 2.0]
r = qz([[1.0, 2.0], [3.0, 4.0]], [[0.0, 0.0], [0.0, 0.0]])
assert all(x == 0 for x in r.beta) # B = 0 and det A != 0: two infinite eigenvalues
r = qz([[1.0, 0.0], [0.0, 0.0]], [[1.0, 0.0], [0.0, 0.0]])
assert any(abs(x) <= 1e-15 and abs(y) <= 1e-15 for x, y in zip(r.alpha, r.beta)) # singular pencil
with pytest.raises(ValueError):
qz([[1.0, 2.0]], [[1.0]])
What the next level needs
- the tests did not run in this build