dqds algorithm in Python
Computes all singular values of a bidiagonal matrix to high relative accuracy by shifted Cholesky LR steps on the squared entries, with no subtraction of large terms.
- File
dqds_algorithm.py, 85 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/dqds-algorithm/. The authors are listed in prov/BE4Y2X.jsonl.
"""The dqds algorithm (differential quotient-difference with shifts) for the
singular values of an upper bidiagonal matrix, in the positive qd form
q_k = a_k^2, e_k = b_k^2. Python standard library only.
singular_values(a, b) returns (sigma, stats) with sigma[0] >= ... >= sigma[n-1] >= 0.
"""
import math
U = 2.0 ** -53
def dqds_step(q, e, n, tau):
"""One dqds transform of the leading n entries with shift tau.
Returns (qh, eh, dmin) on success, or None if a new q would be <= 0 or a d < 0,
which means tau is not below the smallest eigenvalue."""
qh, eh = [0.0] * n, [0.0] * (n - 1)
d = q[0] - tau # @line L4
dmin = d
for k in range(n - 1): # @line L5
if d < 0:
return None
qh[k] = d + e[k] # @line L6
if qh[k] == 0.0: # zero block boundary: e_k = 0 and d = 0
eh[k] = 0.0
d = q[k + 1] - tau
else:
t = q[k + 1] / qh[k]
eh[k] = e[k] * t # @line L7
d = d * t - tau # @line L8
dmin = min(dmin, d)
if d < 0:
return None
qh[n - 1] = d # @line L9
return qh, eh, dmin
def newton_shift(q, e, n):
"""Newton step for det(B B^T - x I) from x = 0: 1 / trace((B B^T)^-1) = 1 / ||B^-1||_F^2.
Since every eigenvalue is real and positive, this is a lower bound for the
smallest one. The column norms c_j of B^-1 satisfy c_1 = 1 / q_1 and
c_j = (1 + e_j-1 c_j-1) / q_j, all with positive terms."""
c = 1.0 / q[0]
total = c
for j in range(1, n):
c = (1.0 + e[j - 1] * c) / q[j]
total += c
return 1.0 / total
def singular_values(a, b, max_steps_per_value=60):
"""Singular values of the bidiagonal matrix with diagonal a and superdiagonal b."""
q = [float(x) * float(x) for x in a] # @line L1
e = [float(x) * float(x) for x in b]
n = len(q)
stats = {"steps": 0, "failures": 0}
sigma2 = []
shift = 0.0
limit = max_steps_per_value * max(n, 1)
while n > 0:
if n == 1:
sigma2.append(shift + q[0])
break
if e[n - 2] <= U * U * (shift + q[n - 1]) or e[n - 2] == 0.0: # @line L2
sigma2.append(shift + q[n - 1])
n -= 1
continue
if stats["steps"] >= limit:
raise ArithmeticError("dqds did not converge")
# Newton's lower bound for the smallest eigenvalue; tau = 0 if B is singular
tau = newton_shift(q, e, n) if min(q[:n]) > 0.0 else 0.0 # @line L3
out = dqds_step(q, e, n, tau)
if out is None: # tau was not below the smallest eigenvalue after rounding
stats["failures"] += 1
tau = 0.0
out = dqds_step(q, e, n, tau)
qh, eh, _ = out
q[:n], e[:n - 1] = qh, eh
shift += tau # @line L10
stats["steps"] += 1
sigma = sorted((math.sqrt(max(s, 0.0)) for s in sigma2), reverse=True) # @line L11
return sigma, stats
test_dqds_algorithm.py 125 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/dqds-algorithm/. The authors are listed in prov/BE4Y2X.jsonl.
import math
import random
from fractions import Fraction
from dqds_algorithm import U, dqds_step, newton_shift, singular_values
def exact_singular_values(a, b):
"""Singular values from exact rational Sturm counts on B^T B, to about 2^-60 relative."""
n = len(a)
a = [Fraction(x) for x in a]
b = [Fraction(x) for x in b]
diag = [a[k] ** 2 + (b[k - 1] ** 2 if k else 0) for k in range(n)]
off2 = [(a[k] * b[k]) ** 2 for k in range(n - 1)]
def below(x):
count, d = 0, diag[0] - x
count += d < 0
for k in range(1, n):
if d == 0:
d = Fraction(1, 10 ** 80)
d = diag[k] - x - off2[k - 1] / d
count += d < 0
return count
top = 2 * sum(diag) + 1
out = []
for i in range(n):
lo, hi = Fraction(0), Fraction(top)
while hi - lo > lo / 2 ** 60 and hi > Fraction(1, 10 ** 300):
if lo > 0 and hi > 4 * lo:
x = Fraction(math.sqrt(float(lo) * float(hi)))
elif lo == 0:
x = hi / 1024
else:
x = (lo + hi) / 2
if below(x) > i:
hi = x
else:
lo = x
out.append(math.sqrt(float((lo + hi) / 2)))
return sorted(out, reverse=True)
def assert_relative(got, want, c=8.0):
n = len(want)
for g, w in zip(got, want):
if w == 0:
assert g == 0
else:
assert abs(g - w) <= c * n * U * w, (g, w)
def test_page_example_graded():
a, b = [1.0, 1e-4, 1e-8, 1e-12], [1.0, 1e-4, 1e-8]
sigma, stats = singular_values(a, b)
assert stats["steps"] == 5 and stats["failures"] == 0
assert_relative(sigma, exact_singular_values(a, b), c=2.0)
assert abs(sigma[3] - 4.9999999859375e-13) <= 4 * U * 5e-13
q, e = [x * x for x in a], [x * x for x in b]
tau = newton_shift(q, e, 4)
assert abs(tau - 2.5e-25) <= 1e-6 * 2.5e-25
qh, eh, _ = dqds_step(q, e, 4, tau)
assert all(x > 0 for x in qh) and eh[2] < 1e-24 # e_3 fell from 1e-16 to 7.5e-25
def test_relative_accuracy_where_shifted_qr_is_absolute():
a, b = [1e-12, 1e-4, 1e-8, 1.0], [1e-6, 1e-6, 1e-2]
sigma, _ = singular_values(a, b)
exact = exact_singular_values(a, b)
assert_relative(sigma, exact)
assert abs(exact[3] - 7.070714e-13) < 1e-18
def test_known_spectra():
n = 12
sigma, _ = singular_values([1.0] * n, [1.0] * (n - 1))
assert_relative(sigma, [2 * math.cos(k * math.pi / (2 * n + 1)) for k in range(1, n + 1)])
n = 40 # Cholesky factor of tridiag(-1, 2, -1): singular values sqrt(2 - 2 cos(k pi / (n + 1)))
a = [math.sqrt((k + 1) / k) for k in range(1, n + 1)]
b = [-math.sqrt(k / (k + 1)) for k in range(1, n)]
sigma, _ = singular_values(a, b)
want = sorted((2 * math.sin(k * math.pi / (2 * (n + 1))) for k in range(1, n + 1)), reverse=True)
assert_relative(sigma, want, c=16.0)
def test_random_exact_invariants_and_rate():
rnd = random.Random(1)
n = 100
a = [rnd.uniform(0.1, 1) for _ in range(n)]
b = [rnd.uniform(0.1, 1) for _ in range(n - 1)]
sigma, stats = singular_values(a, b)
det = math.prod(a)
assert abs(math.prod(sigma) / det - 1) <= 8 * n * U
frob = float(sum(Fraction(x) ** 2 for x in a) + sum(Fraction(x) ** 2 for x in b))
assert abs(sum(s * s for s in sigma) - frob) <= 8 * n * U * frob
assert stats["steps"] <= 12 * n # observed about 9 transforms per singular value
a = [rnd.choice([-1, 1]) * 10 ** rnd.uniform(-6, 0) for _ in range(6)]
b = [rnd.choice([-1, 1]) * 10 ** rnd.uniform(-6, 0) for _ in range(5)]
assert_relative(singular_values(a, b)[0], exact_singular_values(a, b))
def test_trivial_and_zero_cases():
assert singular_values([-3.0], [])[0] == [3.0]
assert singular_values([0.0] * 3, [0.0] * 2)[0] == [0.0, 0.0, 0.0]
sigma, _ = singular_values([2.0, -5.0, 1.0], [0.0, 0.0]) # diagonal
assert sigma == [5.0, 2.0, 1.0]
for a, b in (([0.0, 1.0, 2.0], [1.0, 1.0]), ([1.0, 0.0, 2.0], [1.0, 1.0]), ([1.0, 2.0, 0.0], [1.0, 1.0])):
sigma, _ = singular_values(a, b)
assert sigma[-1] == 0.0
assert abs(sum(s * s for s in sigma) - (sum(x * x for x in a) + sum(x * x for x in b))) < 1e-14
a, b, c = 3.0, 4.0, 2.0 # 2 by 2: sigma1 sigma2 = |ac|, sigma1^2 + sigma2^2 = a^2 + b^2 + c^2
s1, s2 = singular_values([a, c], [b])[0]
assert abs(s1 * s2 - 6.0) < 1e-14 and abs(s1 * s1 + s2 * s2 - 29.0) < 1e-13
def test_repeated_and_close_values():
sigma, _ = singular_values([1.0] * 5, [0.0] * 4)
assert sigma == [1.0] * 5
sigma, _ = singular_values([1.0, 1.0], [1e-8]) # close pair 1 +- 5e-9
assert abs(sigma[0] - (1 + 5e-9)) < 4 * U and abs(sigma[1] - (1 - 5e-9)) < 4 * U
sigma, _ = singular_values([1.0, 1.0 + 1e-12, 1.0], [1e-300, 1e-300])
assert abs(sigma[0] - (1 + 1e-12)) < 4 * U and sigma[1] == 1.0 and sigma[2] == 1.0
What the next level needs
- the tests did not run in this build