Block Lanczos algorithm in Python
Block Lanczos iteration for finding nullspace vectors of large sparse matrices over GF(2) and other finite fields.
- File
block_lanczos_algorithm.py, 194 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/block-lanczos-algorithm/. The authors are listed in prov/ZXGZZS.jsonl.
"""Montgomery's block Lanczos algorithm over GF(2).
A block of s vectors of length n is a list of n integers; bit c of entry r is the
entry in row r and column c. An s x s matrix is a list of s integers (its rows).
The sparse matrix B is a list of rows, each a list of the column indices of its ones.
"""
import random
from typing import Sequence
SparseRows = Sequence[Sequence[int]]
def b_times(b: SparseRows, v: Sequence[int]) -> list[int]:
"""B V for a sparse m x n matrix B and an n x s block V."""
out = []
for row in b:
acc = 0
for j in row:
acc ^= v[j]
out.append(acc)
return out
def bt_times(b: SparseRows, u: Sequence[int], n: int) -> list[int]:
"""B^T U for a sparse m x n matrix B and an m x s block U."""
out = [0] * n
for k, row in enumerate(b):
x = u[k]
if x:
for j in row:
out[j] ^= x
return out
def inner(v: Sequence[int], w: Sequence[int], s: int) -> list[int]:
"""V^T W, an s x s matrix, for two blocks with the same number of rows."""
res = [0] * s
for x, y in zip(v, w):
while x:
low = x & -x
res[low.bit_length() - 1] ^= y
x ^= low
return res
def times_small(v: Sequence[int], m: Sequence[int]) -> list[int]:
"""V M for a block V and an s x s matrix M (also M1 M2 for two s x s matrices)."""
out = []
for x in v:
acc = 0
while x:
low = x & -x
acc ^= m[low.bit_length() - 1]
x ^= low
out.append(acc)
return out
def add(m1: Sequence[int], m2: Sequence[int]) -> list[int]:
return [x ^ y for x, y in zip(m1, m2)]
def select(t: Sequence[int], prev_mask: int, s: int) -> tuple[list[int], int] | None:
"""Montgomery's choice of S_i and W_i^inv = S (S^T T S)^{-1} S^T from T = V^T A V.
The columns not in the previous selection are taken first. Returns (winv, mask),
or None if the step breaks down.
"""
full = (1 << s) - 1
order = [c for c in range(s) if not prev_mask >> c & 1] + [c for c in range(s) if prev_mask >> c & 1]
m = [t[c] | (1 << (s + c)) for c in range(s)] # row c of [T | I], the identity in bits s..2s-1
mask = 0
for idx, c in enumerate(order):
k = next((order[q] for q in range(idx, s) if m[order[q]] >> c & 1), None)
if k is not None:
m[c], m[k] = m[k], m[c]
mask |= 1 << c
for r in range(s):
if r != c and m[r] >> c & 1:
m[r] ^= m[c]
else:
k = next((order[q] for q in range(idx, s) if m[order[q]] >> (s + c) & 1), None)
if k is None:
return None
m[c], m[k] = m[k], m[c]
for r in range(s):
if r != c and m[r] >> (s + c) & 1:
m[r] ^= m[c]
m[c] = 0
if (mask | prev_mask) != full:
return None
return [(m[c] >> s) & full for c in range(s)], mask
def block_lanczos_kernel(b: SparseRows, n: int, s: int = 16, seed: int = 0,
tries: int = 3, trace: list | None = None) -> list[list[int]]:
"""Independent nonzero vectors x in GF(2)^n with B x = 0, for a sparse m x n matrix B.
Runs block Lanczos on A = B^T B with block width s. Every returned vector is checked;
the list may be empty if all tries break down or find nothing. If trace is a list,
each iteration appends (V_i, mask of S_i) to it, for testing.
Cost: about n / (s - 0.76) iterations, each with one product by B and one by B^T on
an n x s block, and O(n s) word operations.
"""
rng = random.Random(seed)
full = (1 << s) - 1
for _ in range(tries):
if trace is not None:
trace.clear()
found = _one_run(b, n, s, rng, full, trace)
if found:
return found
return []
def _one_run(b: SparseRows, n: int, s: int, rng: random.Random, full: int,
trace: list | None) -> list[list[int]]:
if n == 0:
return []
y = [rng.getrandbits(s) for _ in range(n)] # @line L1
v = bt_times(b, b_times(b, y), n) # V_0 = A Y
v0 = v
x = [0] * n
zero_block = [0] * n
zero_small = [0] * s
v1 = v2 = zero_block # V_{i-1}, V_{i-2}
winv1 = winv2 = zero_small
t1 = t2_1 = zero_small # T and T2 of step i-1
mask1 = full
for _ in range(n + 10): # @line L2
av = bt_times(b, b_times(b, v), n) # @line L3
t = inner(v, av, s)
t2 = inner(av, av, s)
if not any(t): # @line L4
break
sel = select(t, mask1, s) # @line L5
if sel is None:
break # breakdown: finish with the current block, the final step checks everything
winv, mask = sel
if trace is not None:
trace.append((v, mask))
x = add(x, times_small(v, times_small(winv, inner(v, v0, s)))) # @line L6
ident = [1 << c for c in range(s)]
t2m = [r & mask for r in t2]
d = add(ident, times_small(winv, add(t2m, t))) # @line L7
e = times_small(winv1, [r & mask for r in t])
f = times_small(times_small(winv2, add(ident, times_small(t1, winv1))),
[r & mask for r in add([r & mask1 for r in t2_1], t1)])
new = add(add([r & mask for r in av], times_small(v, d)), add(times_small(v1, e), times_small(v2, f))) # @line L8
v2, v1, v = v1, v, new
winv2, winv1 = winv1, winv
t1, t2_1, mask1 = t, t2, mask
# @line L9
z = [(xi ^ yi) | (vi << s) for xi, yi, vi in zip(x, y, v)] # columns: X + Y, then V_m
bz = b_times(b, z)
cols = 2 * s
zc = [_column(z, c) for c in range(cols)]
bzc = [_column(bz, c) for c in range(cols)]
basis: dict[int, tuple[int, int]] = {}
kernel = []
for c in range(cols): # @line L10
img, pre = bzc[c], zc[c]
while img:
top = img.bit_length() - 1
if top not in basis:
break
bi, bp = basis[top]
img, pre = img ^ bi, pre ^ bp
if img:
basis[img.bit_length() - 1] = (img, pre)
elif pre:
kernel.append(pre)
out: dict[int, int] = {}
for vec in kernel:
while vec:
top = vec.bit_length() - 1
if top not in out:
out[top] = vec
break
vec ^= out[top]
return [[vec >> r & 1 for r in range(n)] for vec in out.values()] # @line L11
def _column(block: Sequence[int], c: int) -> int:
"""Column c of a block as an integer with bit r for row r."""
col = 0
for r, x in enumerate(block):
if x >> c & 1:
col |= 1 << r
return col
test_block_lanczos_algorithm.py 207 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/block-lanczos-algorithm/. The authors are listed in prov/ZXGZZS.jsonl.
import random
from block_lanczos_algorithm import (add, b_times, block_lanczos_kernel, bt_times, inner,
select, times_small)
# Dense GF(2) references: a matrix is a list of rows, each a list of 0/1.
def dense(b, n):
return [[1 if j in set(row) else 0 for j in range(n)] for row in b]
def matmul(p, q):
if not p:
return []
k = len(q)
cols = len(q[0]) if q else 0
return [[sum(p[i][t] & q[t][j] for t in range(k)) & 1 for j in range(cols)] for i in range(len(p))]
def transpose(p, cols):
return [[p[i][j] for i in range(len(p))] for j in range(cols)]
def to_dense(block, s):
return [[x >> c & 1 for c in range(s)] for x in block]
def from_dense(rows):
return [sum(bit << c for c, bit in enumerate(r)) for r in rows]
def rank(rows):
"""Rank over GF(2) of a list of 0/1 vectors, by elimination on bit masks."""
basis = {}
for r in rows:
vec = sum(bit << c for c, bit in enumerate(r))
while vec:
top = vec.bit_length() - 1
if top not in basis:
basis[top] = vec
break
vec ^= basis[top]
return len(basis)
def nullity(b, n):
return n - rank(dense(b, n))
def random_sparse(rng, m, n, per_row):
return [sorted(rng.sample(range(n), min(per_row, n))) for _ in range(m)]
def check_output(b, n, vecs):
"""Every vector is nonzero with B x = 0, and the vectors are independent."""
d = dense(b, n)
for x in vecs:
assert len(x) == n and any(x)
assert all(sum(r[j] & x[j] for j in range(n)) % 2 == 0 for r in d)
assert rank(vecs) == len(vecs)
assert len(vecs) <= nullity(b, n)
def test_helpers_against_dense():
rng = random.Random(1)
for _ in range(60):
m, n, s = rng.randint(0, 9), rng.randint(1, 9), rng.randint(1, 6)
b = random_sparse(rng, m, n, rng.randint(0, n))
bd = dense(b, n)
v = [rng.getrandbits(s) for _ in range(n)]
u = [rng.getrandbits(s) for _ in range(m)]
w = [rng.getrandbits(s) for _ in range(n)]
sm = [rng.getrandbits(s) for _ in range(s)]
vd, ud, wd, smd = to_dense(v, s), to_dense(u, s), to_dense(w, s), to_dense(sm, s)
if m:
assert to_dense(b_times(b, v), s) == matmul(bd, vd)
assert to_dense(bt_times(b, u, n), s) == matmul(transpose(bd, n), ud)
assert to_dense(inner(v, w, s), s) == matmul(transpose(vd, s), wd)
assert to_dense(times_small(v, sm), s) == matmul(vd, smd)
assert add(v, w) == [a ^ c for a, c in zip(v, w)]
def test_select_defining_conditions():
rng = random.Random(2)
seen_ok = 0
for _ in range(400):
s = rng.randint(1, 6)
full = (1 << s) - 1
g = [[rng.getrandbits(1) for _ in range(s)] for _ in range(rng.randint(1, 8))]
t = from_dense(matmul(transpose(g, s), g)) # symmetric, like V^T A V
prev = rng.getrandbits(s)
res = select(t, prev, s)
if res is None:
continue
seen_ok += 1
winv, mask = res
assert mask | prev == full
sel = [c for c in range(s) if mask >> c & 1]
td, wd = to_dense(t, s), to_dense(winv, s)
# winv vanishes outside the selected rows and columns
assert all(wd[i][j] == 0 for i in range(s) for j in range(s) if i not in sel or j not in sel)
# on the selected part, winv is the inverse of S^T T S
sub = [[td[i][j] for j in sel] for i in sel]
wsub = [[wd[i][j] for j in sel] for i in sel]
assert matmul(sub, wsub) == [[int(i == j) for j in range(len(sel))] for i in range(len(sel))]
# S^T T S has full rank, and adding any other column of T makes it singular
r_t = rank(td)
assert len(sel) == r_t
assert seen_ok > 100
def test_worked_example():
b = [[0, 1], [1, 2], [0, 2], [3, 4], [4, 5]]
trace = []
vecs = block_lanczos_kernel(b, 6, s=2, seed=0, tries=1, trace=trace)
assert sorted(vecs) == [[1, 1, 1, 0, 0, 0], [1, 1, 1, 1, 1, 1]]
assert len(trace) == 1
v0, mask = trace[0]
assert to_dense(v0, 2) == [[0, 1], [0, 0], [0, 1], [0, 1], [1, 1], [1, 0]] and mask == 3
check_output(b, 6, vecs)
def test_edge_cases():
assert block_lanczos_kernel([], 0) == []
assert block_lanczos_kernel([[], []], 0) == []
# no rows: every vector is in the kernel
vecs = block_lanczos_kernel([], 5, s=8)
check_output([], 5, vecs)
assert len(vecs) == 5
# zero matrix
vecs = block_lanczos_kernel([[], [], []], 4, s=8)
check_output([[], [], []], 4, vecs)
assert len(vecs) == 4
# 1 x 1
assert block_lanczos_kernel([[0]], 1, s=4) == []
vecs = block_lanczos_kernel([[]], 1, s=4)
assert vecs == [[1]]
# identity: nullity 0, nothing to find
for n in (1, 2, 7, 30):
assert block_lanczos_kernel([[i] for i in range(n)], n, s=8) == []
# zero columns: e_j for an unused column j is in the kernel
b = [[0, 1], [1, 3], [0, 3]]
vecs = block_lanczos_kernel(b, 5, s=8)
check_output(b, 5, vecs)
assert len(vecs) == nullity(b, 5) == 3
# more columns than rows
rng = random.Random(3)
b = random_sparse(rng, 10, 25, 4)
vecs = block_lanczos_kernel(b, 25, s=16)
check_output(b, 25, vecs)
assert len(vecs) >= 1
# repeated rows
b = [[0, 2, 5]] * 6 + [[1, 3]]
vecs = block_lanczos_kernel(b, 6, s=8)
check_output(b, 6, vecs)
assert len(vecs) == 4
def test_orthogonality_of_w_blocks():
rng = random.Random(4)
for _ in range(5):
n, m = 60, 58
b = random_sparse(rng, m, n, 4)
trace = []
block_lanczos_kernel(b, n, s=8, seed=rng.randrange(10**6), tries=1, trace=trace)
assert len(trace) >= 3
ws = [[x & mask for x in v] for v, mask in trace] # W_i spanned by V_i S_i
aws = [bt_times(b, b_times(b, w), n) for w in ws]
for i in range(len(ws)):
for j in range(i):
assert not any(inner(ws[j], aws[i], 8))
def test_random_correctness_and_success():
rng = random.Random(5)
with_kernel = hit = found_total = null_total = 0
for case in range(500):
n = rng.randint(1, 24)
m = rng.randint(0, n + 3)
b = random_sparse(rng, m, n, rng.randint(1, 4))
s = rng.choice((4, 8, 16))
vecs = block_lanczos_kernel(b, n, s=s, seed=case, tries=3)
check_output(b, n, vecs)
k = nullity(b, n)
if k:
with_kernel += 1
hit += bool(vecs)
found_total += len(vecs)
null_total += min(k, s) # one run yields at most about s vectors
else:
assert vecs == []
assert with_kernel > 200
assert hit >= 0.95 * with_kernel
assert found_total >= 0.85 * null_total
def test_large_sparse():
rng = random.Random(6)
n, m = 300, 290
b = random_sparse(rng, m, n, 5)
vecs = block_lanczos_kernel(b, n, s=32, seed=1)
check_output(b, n, vecs)
assert len(vecs) >= 1
What the next level needs
- the tests did not run in this build