Incomplete LU factorization in Python
LU with dropped fill-in, used as a sparse preconditioner.
- File
incomplete_lu_factorization.py, 150 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/incomplete-lu-factorization/. The authors are listed in prov/YZPFA2.jsonl.
"""Incomplete LU factorization ILU(p) by levels of fill, in CSR form.
ILU(0) keeps exactly the sparsity pattern of A; ILU(p) keeps fill entries
of level at most p (Saad, Section 10.3). The row-wise IKJ variant is used.
Works with floats or fractions.Fraction. Standard library only.
"""
import heapq
class CSR:
"""Compressed sparse row matrix with sorted column indices per row."""
def __init__(self, n, ptr, ind, data):
self.n, self.ptr, self.ind, self.data = n, ptr, ind, data
@classmethod
def from_dense(cls, a):
ptr, ind, data = [0], [], []
for row in a:
for j, v in enumerate(row):
if v != 0:
ind.append(j)
data.append(v)
ptr.append(len(ind))
return cls(len(a), ptr, ind, data)
@classmethod
def from_rows(cls, n, rows):
"""rows[i] is a dict column -> value."""
ptr, ind, data = [0], [], []
for i in range(n):
for j in sorted(rows[i]):
ind.append(j)
data.append(rows[i][j])
ptr.append(len(ind))
return cls(n, ptr, ind, data)
def row(self, i):
return {self.ind[p]: self.data[p] for p in range(self.ptr[i], self.ptr[i + 1])}
def matvec(self, x):
out = []
for i in range(self.n):
acc = 0 * x[0]
for p in range(self.ptr[i], self.ptr[i + 1]):
acc += self.data[p] * x[self.ind[p]]
out.append(acc)
return out
def to_dense(self):
z = 0 * self.data[0] if self.data else 0
out = [[z] * self.n for _ in range(self.n)]
for i in range(self.n):
for p in range(self.ptr[i], self.ptr[i + 1]):
out[i][self.ind[p]] = self.data[p]
return out
def nnz(self):
return len(self.ind)
class ZeroPivot(ArithmeticError):
pass
class ILU:
"""L (unit lower, strict part stored) and U (upper with diagonal) in CSR."""
def __init__(self, lower, upper):
self.L, self.U = lower, upper
def solve(self, r):
"""Return (L U)^{-1} r by a forward and a backward sweep."""
n = self.L.n
y = list(r)
lp, li, ld = self.L.ptr, self.L.ind, self.L.data
for i in range(n):
acc = y[i]
for p in range(lp[i], lp[i + 1]):
acc -= ld[p] * y[li[p]]
y[i] = acc
up, ui, ud = self.U.ptr, self.U.ind, self.U.data
for i in range(n - 1, -1, -1):
acc = y[i]
diag = None
for p in range(up[i], up[i + 1]):
j = ui[p]
if j == i:
diag = ud[p]
else:
acc -= ud[p] * y[j]
y[i] = acc / diag
return y
def nnz(self):
return self.L.nnz() + self.U.nnz()
def ilu(a, p=0):
"""ILU(p) of the CSR matrix a. Raises ZeroPivot on a zero pivot u_kk."""
n = a.n
INF = float("inf")
urows = [] # dict column -> value for row k of U (columns >= k)
ulev = [] # level of each entry of U
lrows = []
for i in range(n): # @line L1
w = a.row(i)
lev = {j: 0 for j in w}
heap = [j for j in w if j < i]
heapq.heapify(heap)
done = set()
while heap: # @line L2
k = heapq.heappop(heap)
if k in done:
continue
done.add(k)
if lev[k] > p:
continue # dropped entry: not used as a multiplier
ukk = urows[k][k]
if ukk == 0:
raise ZeroPivot("zero pivot in row %d" % k)
w[k] = w[k] / ukk # @line L3
lik, levk = w[k], lev[k]
for j, ukj in urows[k].items(): # @line L4
if j <= k:
continue
newlev = levk + ulev[k][j] + 1
if j in w:
w[j] -= lik * ukj # @line L5
if newlev < lev[j]:
lev[j] = newlev
elif p > 0: # @line L6
# fill (level >= 1, so never kept by ILU(0)); dropped in L7
# if its final level exceeds p
w[j] = -lik * ukj
lev[j] = newlev
if j < i:
heapq.heappush(heap, j)
lrow = {j: v for j, v in w.items() if j < i and lev[j] <= p} # @line L7
urow = {j: v for j, v in w.items() if j >= i and lev[j] <= p}
if i not in urow or urow[i] == 0:
raise ZeroPivot("zero pivot in row %d" % i)
lrows.append(lrow)
urows.append(urow)
ulev.append({j: lev[j] for j in urow})
return ILU(CSR.from_rows(n, lrows), CSR.from_rows(n, urows))
test_incomplete_lu_factorization.py 169 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/incomplete-lu-factorization/. The authors are listed in prov/YZPFA2.jsonl.
import math
import random
from fractions import Fraction as F
import pytest
from incomplete_lu_factorization import CSR, ZeroPivot, ilu
def product(fac):
n = fac.L.n
lo, up = fac.L.to_dense(), fac.U.to_dense()
for i in range(n):
lo[i][i] = 1
return [[sum(lo[i][k] * up[k][j] for k in range(n)) for j in range(n)] for i in range(n)]
def exact_lu(a):
n = len(a)
u = [row[:] for row in a]
lo = [[F(int(i == j)) for j in range(n)] for i in range(n)]
for k in range(n):
for i in range(k + 1, n):
m = u[i][k] / u[k][k]
lo[i][k] = m
u[i] = [s - m * t for s, t in zip(u[i], u[k])]
return lo, u
def random_sparse(rng, n, density, dom=True):
a = [[F(0)] * n for _ in range(n)]
for i in range(n):
for j in range(n):
if i != j and rng.random() < density:
a[i][j] = F(rng.randint(-4, 4))
a[i][i] = F(sum(abs(t) for t in a[i]) + rng.randint(1, 3)) if dom else F(rng.randint(1, 5))
return a
def test_page_example():
a = [[F(4), F(-1), F(-1), F(0)],
[F(-1), F(4), F(0), F(-1)],
[F(-1), F(0), F(4), F(-1)],
[F(0), F(-1), F(-1), F(4)]]
fac = ilu(CSR.from_dense(a), 0)
lu = product(fac)
r = [[lu[i][j] - a[i][j] for j in range(4)] for i in range(4)]
# the only nonzero of R = LU - A is the dropped fill at (2,3) and (3,2), 1-based
assert r[1][2] == F(1, 4) and r[2][1] == F(1, 4)
assert sum(1 for row in r for t in row if t != 0) == 2
assert fac.U.to_dense()[3][3] == F(52, 15)
def test_one_by_one_and_diagonal():
fac = ilu(CSR.from_dense([[F(3)]]))
assert fac.solve([F(6)]) == [F(2)]
d = [[F(i + 1) if i == j else F(0) for j in range(4)] for i in range(4)]
fac = ilu(CSR.from_dense(d))
assert fac.solve([F(1)] * 4) == [F(1), F(1, 2), F(1, 3), F(1, 4)]
def test_random_pattern_property_500():
# (LU)_ij = a_ij for every (i, j) in the pattern of A (ILU(0)).
rng = random.Random(1234)
for trial in range(500):
n = rng.randint(1, 7)
a = random_sparse(rng, n, 0.35)
fac = ilu(CSR.from_dense(a), 0)
lu = product(fac)
for i in range(n):
for j in range(n):
if a[i][j] != 0:
assert lu[i][j] == a[i][j]
# pattern of L + U is the pattern of A
assert fac.nnz() == sum(1 for row in a for t in row if t != 0)
def test_no_fill_patterns_are_exact():
# tridiagonal and arrowhead-down matrices have no fill: ILU(0) = LU
rng = random.Random(5)
for trial in range(50):
n = rng.randint(2, 8)
a = [[F(0)] * n for _ in range(n)]
for i in range(n):
a[i][i] = F(rng.randint(5, 9))
if i:
a[i][i - 1] = F(rng.randint(-2, 2))
a[i - 1][i] = F(rng.randint(-2, 2))
assert product(ilu(CSR.from_dense(a))) == a
def test_large_level_is_exact_lu():
rng = random.Random(77)
for trial in range(60):
n = rng.randint(2, 7)
a = random_sparse(rng, n, 0.3)
fac = ilu(CSR.from_dense(a), n)
lo, up = exact_lu(a)
assert fac.U.to_dense() == up
assert product(fac) == a
def test_levels_monotone_fill():
n = 8
a = [[F(0)] * (n * n) for _ in range(n * n)]
for i in range(n):
for j in range(n):
k = i * n + j
a[k][k] = F(4)
for di, dj in ((1, 0), (-1, 0), (0, 1), (0, -1)):
if 0 <= i + di < n and 0 <= j + dj < n:
a[k][(i + di) * n + j + dj] = F(-1)
c = CSR.from_dense([[float(t) for t in row] for row in a])
counts = [ilu(c, p).nnz() for p in range(4)]
assert counts[0] == c.nnz()
assert counts == sorted(counts) and counts[3] > counts[0]
def test_zero_pivot_reported():
with pytest.raises(ZeroPivot):
ilu(CSR.from_dense([[F(0), F(1)], [F(1), F(0)]]))
with pytest.raises(ZeroPivot): # singular
ilu(CSR.from_dense([[F(1), F(1)], [F(1), F(1)]]))
def test_m_matrix_richardson_converges():
# For an M-matrix ILU(0) is a regular splitting: rho(I - M^{-1} A) < 1.
n = 12
rows = []
for i in range(n * n):
rows.append([0.0] * (n * n))
for i in range(n):
for j in range(n):
k = i * n + j
rows[k][k] = 4.0
for di, dj, c in ((1, 0, -1.3), (-1, 0, -0.7), (0, 1, -1.0), (0, -1, -1.0)):
if 0 <= i + di < n and 0 <= j + dj < n:
rows[k][(i + di) * n + j + dj] = c
a = CSR.from_dense(rows)
fac = ilu(a, 0)
rng = random.Random(8)
b = [rng.uniform(-1, 1) for _ in range(n * n)]
x = [0.0] * (n * n)
norm = lambda v: math.sqrt(sum(t * t for t in v))
res = []
for it in range(400):
r = [bi - t for bi, t in zip(b, a.matvec(x))]
res.append(norm(r))
if res[-1] < 1e-10 * norm(b):
break
x = [xi + di for xi, di in zip(x, fac.solve(r))]
assert res[-1] < 1e-10 * norm(b)
# compare Jacobi splitting: ILU(0) needs far fewer sweeps
xj = [0.0] * (n * n)
for it in range(len(res)):
r = [bi - t for bi, t in zip(b, a.matvec(xj))]
xj = [xi + ri / 4.0 for xi, ri in zip(xj, r)]
rj = norm([bi - t for bi, t in zip(b, a.matvec(xj))])
assert rj > 1e3 * res[-1]
def test_hilbert_ilu_is_full_lu():
n = 6
h = [[F(1, i + j + 1) for j in range(n)] for i in range(n)]
fac = ilu(CSR.from_dense(h), 0) # dense pattern: no dropping
assert product(fac) == h
What the next level needs
- the tests did not run in this build