Solving linear systems

Solve A·x = b by LUP decomposition and two substitutions: L·y = P·b forward, then U·x = y backward; in exact fractions.

With the factorization P·A = L·U, the equation A·x = b becomes two triangular ones:

Read the full explanation ↓

Time: O(n³), then O(n²) for each new b
Space: O(n²)
CLRS, section 28.1

Input:

A 3×3 · b 3×1PDF export

A = [[1, 2, 0], [3, 4, 4], [5, 6, 3]] and b = (3, 7, 8); the solution is x = (−7/5, 11/5, 3/5).

from fractions import Fraction
def solve(A, b):
n = len(A)
if len(b) != n:
return None
a = [[Fraction(x) for x in row] for row in A]
pi = list(range(n))
# 1. P·A = L·U (as on the LUP decomposition page).
for k in range(n):
p = max(range(k, n), key=lambda i: abs(a[i][k]))
if a[p][k] == 0:
return None
pi[k], pi[p] = pi[p], pi[k]
a[k], a[p] = a[p], a[k]
for i in range(k + 1, n):
a[i][k] = a[i][k] / a[k][k]
for j in range(k + 1, n):
a[i][j] = a[i][j] - a[i][k] * a[k][j]
lu = a
# 2. Forward substitution: L·y = P·b.
y = [Fraction(0)] * n
for i in range(n):
y[i] = b[pi[i]][0] - sum(lu[i][j] * y[j] for j in range(i))
# 3. Back substitution: U·x = y.
x = [Fraction(0)] * n
for i in range(n - 1, -1, -1):
x[i] = (y[i] - sum(lu[i][j] * x[j] for j in range(i + 1, n))) / lu[i][i]
return [str(v) for v in x]

Invariant: In back substitution, x[i+1] to x[n−1] are already known, so row i's equation has only one unknown.

Start

pi—

y—

x—

The table is still empty.

Local variables

Orange = changed in this step

`A` list

[1, 2, 0][3, 4, 4][5, 6, 3]

`b` list

[3][7][8]

Call stack

solve()line 4

Step 1 of 61

Explanation

The idea

With the factorization P·A = L·U, the equation A·x = b becomes two triangular ones:

L·y = P·b: since L is lower triangular, y[0] is known right away, then y[1], and so on (forward substitution).
U·x = y: since U is upper triangular, start from the last row: x[n−1], then upward (back substitution).

Cost

The factorization costs O(n³), but each substitution only O(n²). For several bs with the same A, factor once. (Computing A⁻¹ and multiplying by b is slower and numerically worse.)

What to watch in the trace

The lu table is L and U in one matrix; the y and x cards fill top down and bottom up, respectively.