Determinant by Kraut's method (LU decomposition) in Python
Computes a matrix determinant by factoring it into lower and upper triangular matrices and multiplying the diagonal of the factors.
- File
kraut_determinant.py, 58 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/kraut-determinant/. The authors are listed in prov/VFMNY1.jsonl.
from fractions import Fraction
from typing import Any, Sequence
def crout(a: Sequence[Sequence[Any]]) -> tuple[list[list[Any]], list[int], int] | None:
"""Crout (Kraut) factorization with row swaps over a field.
Returns (w, perm, s) such that, with P the rows of a in the order perm, P = L U, where L is
the lower triangle of w (with its diagonal) and U is the strict upper triangle of w with a unit
diagonal; s = +1 or -1 is the sign of the permutation. Returns None when a is singular.
Integers are converted to Fraction; other field elements are used as they are.
Cost: 2n^3/3 + O(n^2) field operations.
"""
w = [[Fraction(x) if isinstance(x, int) else x for x in row] for row in a]
n = len(w)
perm = list(range(n))
s = 1 # @line L1
for j in range(n): # @line L2
for i in range(j, n): # @line L3
acc = w[i][j]
for k in range(j):
acc = acc - w[i][k] * w[k][j]
w[i][j] = acc
p = next((i for i in range(j, n) if w[i][j] != 0), None) # @line L4
if p is None: # @line L5
return None
if p != j: # @line L6
w[j], w[p] = w[p], w[j]
perm[j], perm[p] = perm[p], perm[j]
s = -s
piv = w[j][j]
for t in range(j + 1, n): # @line L7
acc = w[j][t]
for k in range(j):
acc = acc - w[j][k] * w[k][t]
w[j][t] = acc / piv
return w, perm, s
def kraut_determinant(a: Sequence[Sequence[Any]]) -> Any:
"""Return det a for a square matrix over a field (ints become Fraction), by Crout's method.
Cost: 2n^3/3 + O(n^2) field operations.
"""
if not a:
return Fraction(1)
res = crout(a)
if res is None:
z = Fraction(a[0][0]) if isinstance(a[0][0], int) else a[0][0]
return z - z
w, _, s = res
d = w[0][0] if s == 1 else -w[0][0] # @line L8
for j in range(1, len(w)): # @line L9
d = d * w[j][j]
return d # @line L10
test_kraut_determinant.py 118 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/kraut-determinant/. The authors are listed in prov/VFMNY1.jsonl.
import random
from fractions import Fraction
from itertools import permutations
from kraut_determinant import crout, kraut_determinant
def det_leibniz(a):
n = len(a)
total = 0
for p in permutations(range(n)):
inv = sum(1 for i in range(n) for j in range(i + 1, n) if p[i] > p[j])
t = 1
for i in range(n):
t *= a[p[i]][i]
total += -t if inv % 2 else t
return total
class ModP:
"""An element of Z/pZ, enough of a field for the algorithm."""
def __init__(self, v: int, p: int) -> None:
self.v, self.p = v % p, p
def __sub__(self, o):
return ModP(self.v - o.v, self.p)
def __mul__(self, o):
return ModP(self.v * o.v, self.p)
def __truediv__(self, o):
return ModP(self.v * pow(o.v, self.p - 2, self.p), self.p)
def __neg__(self):
return ModP(-self.v, self.p)
def __eq__(self, o):
if isinstance(o, ModP):
return self.v == o.v
return self.v == o % self.p
def __ne__(self, o):
return not self == o
def check_factorization(a, res):
w, perm, s = res
n = len(a)
lo = [[w[i][j] if j <= i else 0 for j in range(n)] for i in range(n)]
up = [[1 if i == j else (w[i][j] if j > i else 0) for j in range(n)] for i in range(n)]
for i in range(n):
for j in range(n):
assert sum(lo[i][k] * up[k][j] for k in range(n)) == a[perm[i]][j]
inv = sum(1 for i in range(n) for j in range(i + 1, n) if perm[i] > perm[j])
assert s == (-1) ** inv
def test_examples() -> None:
a = [[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]]
assert kraut_determinant(a) == -12
w, perm, s = crout(a)
assert perm == [1, 0, 2, 3] and s == -1
assert [w[j][j] for j in range(4)] == [1, 1, -4, -3]
assert kraut_determinant([[2, 1, 0], [1, 3, 1], [0, 1, 4]]) == 18
def test_edges() -> None:
assert kraut_determinant([]) == 1
assert kraut_determinant([[5]]) == 5
assert kraut_determinant([[0]]) == 0
assert kraut_determinant([[0, 1], [1, 0]]) == -1
assert kraut_determinant([[1, 2], [2, 4]]) == 0
assert kraut_determinant([[0, 0, 0], [1, 2, 3], [4, 5, 6]]) == 0
assert crout([[0, 0], [0, 1]]) is None
def test_random_rational() -> None:
rng = random.Random(1941)
for _ in range(800):
n = rng.randint(1, 6)
a = [[rng.randint(-4, 4) for _ in range(n)] for _ in range(n)]
if rng.random() < 0.2:
a[rng.randrange(n)] = [0] * n
if rng.random() < 0.3:
a[0][0] = 0
if rng.random() < 0.2:
for i in range(n):
a[i][rng.randrange(n)] = 0
d = det_leibniz(a)
assert kraut_determinant(a) == d
res = crout(a)
assert (res is None) == (d == 0)
if res is not None:
check_factorization(a, res)
def test_random_mod_p() -> None:
rng = random.Random(23)
for p in (2, 3, 5, 7, 101):
for _ in range(150):
n = rng.randint(1, 5)
a = [[rng.randrange(p) for _ in range(n)] for _ in range(n)]
got = kraut_determinant([[ModP(x, p) for x in row] for row in a])
assert got.v == det_leibniz(a) % p
def test_fractions_and_large_entries() -> None:
rng = random.Random(8)
for _ in range(20):
n = 5
a = [[rng.randint(-10**30, 10**30) for _ in range(n)] for _ in range(n)]
assert kraut_determinant(a) == det_leibniz(a)
b = [[Fraction(rng.randint(-9, 9), rng.randint(1, 9)) for _ in range(4)] for _ in range(4)]
assert kraut_determinant(b) == det_leibniz(b)
What the next level needs
- the tests did not run in this build
Read the full page about Determinant by Kraut's method (LU decomposition)