Gauss-Jordan elimination in Python
Gaussian elimination that continues to reduced row echelon form, giving the inverse or solution directly.
- File
gauss_jordan_elimination.py, 68 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/gauss-jordan-elimination/. The authors are listed in prov/NK577G.jsonl.
from fractions import Fraction
from typing import Sequence
def gauss_jordan_elimination(a: Sequence[Sequence[int | Fraction]], n: int | None = None) -> tuple[list[list[Fraction]], int, list[int]]:
"""Return (e, r, pivots): the reduced row echelon form e of the m by n matrix a over Q,
its rank r and its pivot columns (0-based). The input is not changed. Pass n when a has no rows.
Cost: O(m n min(m, n)) field operations.
"""
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
c = e[r][j] # @line L6
row_r = e[r]
for t in range(j, n):
row_r[t] /= c
for i in range(m): # @line L7
if i == r:
continue
f = e[i][j] # @line L8
if f != 0:
row_i = e[i]
for t in range(j, n): # @line L9
row_i[t] -= f * row_r[t]
pivots.append(j) # @line L10
r += 1
return e, r, pivots # @line L11
def inverse(a: Sequence[Sequence[int | Fraction]]) -> list[list[Fraction]] | None:
"""Return the inverse of the n by n matrix a over Q, or None when a is singular.
Runs Gauss-Jordan on [a | I]; a is invertible exactly when the pivots are the columns 0..n-1.
Cost: O(n^3) field operations.
"""
n = len(a)
aug = [list(a[i]) + [int(i == j) for j in range(n)] for i in range(n)]
e, r, pivots = gauss_jordan_elimination(aug, 2 * n)
if pivots != list(range(n)):
return None
return [row[n:] for row in e]
def solve(a: Sequence[Sequence[int | Fraction]], b: Sequence[int | Fraction]) -> list[Fraction] | None:
"""Return one solution x of a x = b over Q (free variables set to 0), or None when there is none.
Runs Gauss-Jordan on [a | b]; the system is inconsistent exactly when the last column is a pivot.
Cost: O(m n min(m, n)) field operations for an m by n matrix a.
"""
m = len(a)
n = len(a[0]) if m else 0
e, r, pivots = gauss_jordan_elimination([list(a[i]) + [b[i]] for i in range(m)], n + 1)
if pivots and pivots[-1] == n:
return None
x = [Fraction(0)] * n
for i, j in enumerate(pivots):
x[j] = e[i][n]
return x
test_gauss_jordan_elimination.py 135 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/gauss-jordan-elimination/. The authors are listed in prov/NK577G.jsonl.
import random
from fractions import Fraction
from itertools import combinations, permutations
from gauss_jordan_elimination import gauss_jordan_elimination as gj
from gauss_jordan_elimination import inverse, 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):
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_rref(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] and e[i][pivots[i]] == 1
else:
assert not nz
for k, j in enumerate(pivots):
assert all(e[i][j] == (1 if i == k else 0) for i in range(m))
def matmul(a, b):
return [[sum(a[i][k] * b[k][j] for k in range(len(b))) for j in range(len(b[0]))] for i in range(len(a))]
def test_example() -> None:
a = [[0, 2, 1], [1, 1, 1], [2, 3, 1]]
inv = inverse(a)
assert inv == [[Fraction(-2, 3), Fraction(1, 3), Fraction(1, 3)],
[Fraction(1, 3), Fraction(-2, 3), Fraction(1, 3)],
[Fraction(1, 3), Fraction(4, 3), Fraction(-2, 3)]]
e, r, p = gj([[0, 2, 4, 2], [1, 1, 1, 0], [2, 4, 6, 2]])
assert e == [[1, 0, -1, -1], [0, 1, 2, 1], [0, 0, 0, 0]] and r == 2 and p == [0, 1]
def test_edges() -> None:
assert gj([]) == ([], 0, [])
assert gj([[], []]) == ([[], []], 0, [])
assert gj([[4]]) == ([[1]], 1, [0])
assert gj([[0, 0, 3]]) == ([[0, 0, 1]], 1, [2])
assert inverse([]) == []
assert inverse([[2]]) == [[Fraction(1, 2)]]
assert inverse([[0]]) is None
assert inverse([[1, 2], [2, 4]]) is None
assert inverse([[0, 1], [1, 0]]) == [[0, 1], [1, 0]]
def test_random_rref() -> None:
rng = random.Random(3307)
for _ in range(600):
m, n = rng.randint(0, 4), rng.randint(0, 4)
a = [[rng.randint(-3, 3) for _ in range(n)] for _ in range(m)]
if m >= 2 and rng.random() < 0.3:
a[0] = [x + y for x, y in zip(a[1], a[-1])]
e, r, p = gj(a, n)
check_rref(e, r, p, m, n)
assert r == rank_brute(a, n)
assert rank_brute(a + e, n) == r # same row space
# uniqueness: a row permutation of a has the same reduced form
b = a[:]
rng.shuffle(b)
assert gj(b, n)[0] == e
def test_random_inverse() -> None:
rng = random.Random(611)
for _ in range(600):
n = rng.randint(1, 5)
a = [[rng.randint(-5, 5) for _ in range(n)] for _ in range(n)]
if rng.random() < 0.25:
a[rng.randrange(n)] = [0] * n
if rng.random() < 0.25:
for i in range(n):
a[i][i] = 0
inv = inverse(a)
if det_leibniz(a) == 0:
assert inv is None
else:
ident = [[int(i == j) for j in range(n)] for i in range(n)]
assert matmul(a, inv) == ident and matmul(inv, a) == ident
def test_large_entries() -> None:
rng = random.Random(31)
for _ in range(10):
n = 5
a = [[rng.randint(-10**30, 10**30) for _ in range(n)] for _ in range(n)]
inv = inverse(a)
ident = [[int(i == j) for j in range(n)] for i in range(n)]
assert matmul(a, inv) == ident
def test_random_solve() -> None:
rng = random.Random(4242)
for _ in range(600):
m, n = rng.randint(1, 4), rng.randint(1, 4)
a = [[rng.randint(-3, 3) for _ in range(n)] for _ in range(m)]
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(a, n) == rank_brute([a[i] + [b[i]] for i in range(m)], n + 1)
if consistent:
assert x is not None
assert all(sum(a[i][j] * x[j] for j in range(n)) == b[i] for i in range(m))
else:
assert x is None
assert solve([[0, 2, 1], [1, -1, 0], [2, 0, 3]], [3, -1, -1]) == [1, 2, -1]
assert solve([[1, 1], [1, 1]], [1, 2]) is None
What the next level needs
- the tests did not run in this build