Gaussian elimination in Python
Reduces a matrix to echelon form by row operations to solve systems and read off rank, nullspace and determinant.
- File
gaussian_elimination.py, 52 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/gaussian-elimination/. The authors are listed in prov/82ZZ65.jsonl.
from fractions import Fraction
from typing import Sequence
def gaussian_elimination(a: Sequence[Sequence[int | Fraction]], n: int | None = None) -> tuple[list[list[Fraction]], int, list[int]]:
"""Bring an m by n matrix over Q to row echelon form by row operations.
Returns (e, r, pivots): e is an echelon form of a (a new matrix, a is not changed), r is the rank,
and pivots lists the r pivot columns (0-based, increasing). Pass n when a has no rows.
Cost: O(m n min(m, n)) field operations; with Fraction the entries stay quotients of minors of a.
"""
e = [[Fraction(x) for x in row] for row in a]
m = len(e)
if n is None:
n = len(e[0]) if m else 0
r, pivots = 0, [] # @line L1
for j in range(n): # @line L2
p = next((i for i in range(r, m) if e[i][j] != 0), None) # @line L3
if p is None: # @line L4
continue
e[r], e[p] = e[p], e[r] # @line L5
for i in range(r + 1, m): # @line L6
f = e[i][j] / e[r][j] # @line L7
if f != 0:
row_i, row_r = e[i], e[r]
for t in range(j, n): # @line L8
row_i[t] -= f * row_r[t]
pivots.append(j) # @line L9
r += 1 # @line L10
return e, r, pivots # @line L11
def solve(a: Sequence[Sequence[int | Fraction]], b: Sequence[int | Fraction]) -> list[Fraction] | None:
"""Return one solution x of a x = b over Q, or None when there is none. Free variables are set to 0.
Cost: the elimination of the m by (n + 1) augmented matrix, then O(n r) for back substitution.
"""
m = len(a)
n = len(a[0]) if m else 0
aug = [list(a[i]) + [b[i]] for i in range(m)]
e, r, pivots = gaussian_elimination(aug, n + 1)
if r > 0 and pivots[-1] == n:
return None # a row 0 = nonzero
x = [Fraction(0)] * n
for k in range(r - 1, -1, -1):
j = pivots[k]
s = e[k][n] - sum((e[k][t] * x[t] for t in range(j + 1, n)), Fraction(0))
x[j] = s / e[k][j]
return x
test_gaussian_elimination.py 123 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/gaussian-elimination/. The authors are listed in prov/82ZZ65.jsonl.
import random
from fractions import Fraction
from itertools import combinations, permutations
from gaussian_elimination import gaussian_elimination, solve
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
def rank_brute(a, n):
"""The largest t such that some t by t minor is nonzero."""
m = len(a)
for t in range(min(m, n), 0, -1):
for rows in combinations(range(m), t):
for cols in combinations(range(n), t):
if det_leibniz([[a[i][j] for j in cols] for i in rows]) != 0:
return t
return 0
def check_echelon(e, r, pivots, m, n):
assert len(pivots) == r and pivots == sorted(set(pivots))
for i in range(m):
nz = [j for j in range(n) if e[i][j] != 0]
if i < r:
assert nz and nz[0] == pivots[i]
else:
assert not nz
def random_matrix(rng, m, n, lo=-3, hi=3):
a = [[rng.randint(lo, hi) for _ in range(n)] for _ in range(m)]
if m and n and rng.random() < 0.3: # force rank deficiency: copy a combination of rows
i, k = rng.randrange(m), rng.randrange(m)
c = rng.randint(-2, 2)
a[i] = [c * x for x in a[k]]
if n and rng.random() < 0.2: # a zero column
j = rng.randrange(n)
for row in a:
row[j] = 0
return a
def test_example() -> None:
a = [[0, 2, 1], [1, -1, 0], [2, 0, 3]]
e, r, pivots = gaussian_elimination([row + [b] for row, b in zip(a, [3, -1, -1])])
assert e == [[1, -1, 0, -1], [0, 2, 1, 3], [0, 0, 2, -2]]
assert r == 3 and pivots == [0, 1, 2]
assert solve(a, [3, -1, -1]) == [1, 2, -1]
def test_edges() -> None:
assert gaussian_elimination([]) == ([], 0, [])
assert gaussian_elimination([], 3) == ([], 0, [])
assert gaussian_elimination([[], []]) == ([[], []], 0, [])
assert gaussian_elimination([[5]]) == ([[5]], 1, [0])
assert gaussian_elimination([[0]]) == ([[0]], 0, [])
assert gaussian_elimination([[0, 0], [0, 0]])[1] == 0
assert solve([], []) == []
assert solve([[0, 0]], [1]) is None
assert solve([[0, 0]], [0]) == [0, 0]
# zero first pivot that needs a swap
e, r, p = gaussian_elimination([[0, 1], [1, 0]])
assert e == [[1, 0], [0, 1]] and r == 2 and p == [0, 1]
# input is not changed
a = [[0, 1], [1, 0]]
gaussian_elimination(a)
assert a == [[0, 1], [1, 0]]
def test_large_entries() -> None:
rng = random.Random(77)
for _ in range(20):
n = 6
a = [[rng.randint(-10**30, 10**30) for _ in range(n)] for _ in range(n)]
x0 = [Fraction(rng.randint(-10**20, 10**20), rng.randint(1, 10**10)) for _ in range(n)]
b = [sum(a[i][j] * x0[j] for j in range(n)) for i in range(n)]
x = solve(a, b)
assert x == x0 # the random matrix is nonsingular with high probability; then x is unique
def test_random_rank_and_shape() -> None:
rng = random.Random(4021)
for _ in range(600):
m, n = rng.randint(0, 4), rng.randint(0, 4)
a = random_matrix(rng, m, n)
e, r, pivots = gaussian_elimination(a, n)
check_echelon(e, r, pivots, m, n)
assert r == rank_brute(a, n)
# e has the same row space as a: stacking does not raise the rank
assert rank_brute(a + [[int(x) if x.denominator == 1 else x for x in row] for row in e], n) == r
def test_random_solve() -> None:
rng = random.Random(918)
for _ in range(600):
m, n = rng.randint(1, 4), rng.randint(1, 4)
a = random_matrix(rng, m, n)
if rng.random() < 0.5:
x0 = [rng.randint(-3, 3) for _ in range(n)]
b = [sum(a[i][j] * x0[j] for j in range(n)) for i in range(m)]
else:
b = [rng.randint(-3, 3) for _ in range(m)]
x = solve(a, b)
consistent = rank_brute([row + [bi] for row, bi in zip(a, b)], n + 1) == rank_brute(a, n)
if x is None:
assert not consistent
else:
assert consistent
assert all(sum(a[i][j] * x[j] for j in range(n)) == b[i] for i in range(m))
What the next level needs
- the tests did not run in this build