Systems of Linear Equations

A system of m linear equations in n unknowns is written Ax = b, where A is the m × n coefficient matrix. The augmented matrix [A | b] carries the right-hand sides along.

Gaussian elimination

Use row operations (swap rows, multiply a row by a non-zero number, add a multiple of one row to another) to bring [A | b] to row echelon form, then solve from the bottom up (back substitution). Row operations never change the solution set.

Example:

x + y + z = 6 2x − y + 3z = 9 x + 2y − z = 2

from fractions import Fraction

def solve(A, b):
    n = len(A)
    M = [[Fraction(x) for x in row] + [Fraction(v)] for row, v in zip(A, b)]
    for c in range(n):
        pivot = next(r for r in range(c, n) if M[r][c] != 0)
        M[c], M[pivot] = M[pivot], M[c]
        for r in range(c + 1, n):
            f = M[r][c] / M[c][c]
            M[r] = [a - f * p for a, p in zip(M[r], M[c])]
    print("echelon form:")
    for row in M:
        print("  ", [str(x) for x in row])
    x = [Fraction(0)] * n
    for r in reversed(range(n)):                       # back substitution
        x[r] = (M[r][n] - sum(M[r][k] * x[k] for k in range(r + 1, n))) / M[r][r]
    return [str(v) for v in x]

print("solution:", solve([[1, 1, 1], [2, -1, 3], [1, 2, -1]], [6, 9, 2]))

echelon form:
   ['1', '1', '1', '6']
   ['0', '-3', '1', '-3']
   ['0', '0', '-5/3', '-5']
solution: ['1', '2', '3']

How many solutions? The rank test

Let r = rank(A), r′ = rank([A | b]), and n = number of unknowns:

For a square system, a unique solution exists exactly when det(A) ≠ 0 (and then x = A⁻¹b; Cramer's rule gives xᵢ = det(Aᵢ)/det(A), where Aᵢ has column i replaced by b).

Example with a parameter. For x + y + z = 1, x + 2y + 3z = 2, x + 4y + kz = 4, the determinant of the coefficients is k − 7.

• If k ≠ 7, there is exactly one solution.
• If k = 7, subtracting row 1 from rows 2 and 3 gives (0, 1, 2 | 1) and (0, 3, 6 | 3): the third is 3 times the second, so the system is consistent with rank 2 < 3, and it has infinitely many solutions.

Had the last right-hand side been anything other than 4, k = 7 would have made the system inconsistent.

Homogeneous systems

Ax = 0 always has the trivial solution x = 0. It has non-trivial solutions exactly when rank(A) < n, which for a square matrix means det(A) = 0. The solutions form a subspace (the null space) of dimension n − rank(A) (rank–nullity theorem).

LU decomposition

LU decomposition writes A = LU, with L lower triangular (1s on the diagonal in Doolittle's method) and U upper triangular. U is the echelon form from Gaussian elimination, and L holds the multipliers used. Solving Ax = b then takes two triangular solves, Ly = b (forward) and Ux = y (back), which is cheap to repeat for many right-hand sides b. (If a zero pivot appears, rows must be swapped: PA = LU.)

from fractions import Fraction

def lu(A):
    n = len(A)
    U = [[Fraction(x) for x in row] for row in A]
    L = [[Fraction(int(i == j)) for j in range(n)] for i in range(n)]
    for c in range(n):
        for r in range(c + 1, n):
            L[r][c] = U[r][c] / U[c][c]                  # store the multiplier
            U[r] = [a - L[r][c] * p for a, p in zip(U[r], U[c])]
    return L, U

L, U = lu([[2, 1, 1], [4, 3, 3], [8, 7, 9]])
print("L =", [[str(x) for x in row] for row in L])
print("U =", [[str(x) for x in row] for row in U])
print("det = product of U's diagonal =", U[0][0] * U[1][1] * U[2][2])

L = [['1', '0', '0'], ['2', '1', '0'], ['4', '3', '1']]
U = [['2', '1', '1'], ['0', '1', '1'], ['0', '0', '2']]
det = product of U's diagonal = 4

Practice questions