Gaussian elimination with Python

The lesson on linear systems presented Gaussian elimination as the procedure np.linalg.solve runs internally. This one is about what that procedure leaves behind. Elimination does not merely produce a solution: it produces a factorisation of the matrix, and that factorisation is reusable.

What elimination produces

Applying Gaussian elimination with partial pivoting to $A$ amounts to constructing three matrices:

PA = LU

where $P$ is a permutation recording the row interchanges, $L$ is lower triangular with ones on the diagonal and holds the multipliers used, and $U$ is upper triangular and is the visible result of the elimination.

In the figure, $U$ takes shape through elimination while $L$ stores each multiplier $\ell_{rk}$ in the position that has just been zeroed:

100010001

21-1-3-12-212

P A = L U

Step 1 / 7

Nothing is discarded. Every row operation is recorded: the interchanges in $P$ , the multipliers in $L$ , and the result in $U$ . That is the difference between running the elimination and factorising.

from scipy.linalg import lu

P, L, U = lu(A)
np.allclose(P @ L @ U, A)     # True

SciPy's convention differs from the usual one. It returns $P$ such that $A = PLU$ , whereas the textbooks, and the statement above, write $PA = LU$ . Both describe the same factorisation: since $P$ is a permutation, $P^{-1} = P^\top$ , so $A = PLU$ is equivalent to $P^\top A = LU$ . SciPy's $P$ is the transpose of the $P$ in the statement.

An immediate by-product: the determinant is read off the diagonal of $U$ .

\det A = (-1)^{s} \prod_{i=1}^{n} u_{ii}

where $s$ is the number of row interchanges. For the matrix in the figure, two interchanges and a diagonal of $(-3,\; 5/3,\; 1/5)$ give $\det A = -1$ . This is how np.linalg.det obtains its result: not by cofactor expansion, which would cost $O(n!)$ , but by factorising.

Factor once, solve many times

The factorisation costs roughly $\tfrac{2}{3}n^3$ operations. Solving with it, by contrast, is two triangular substitutions — forward with $L$ , backward with $U$ — costing $2n^2$ together.

That asymmetry is the reason to keep the factors. When the same $A$ is solved against several right-hand sides, a common situation in practice, repeating the elimination throws away all the expensive work:

from scipy.linalg import lu_factor, lu_solve

lu_piv = lu_factor(A)          # once: O(n³)
x1 = lu_solve(lu_piv, b1)      # each one: O(n²)
x2 = lu_solve(lu_piv, b2)

solve k times

1.7·10⁹

factor once

93.3·10⁶

speed-up ×18.0

Flop counts: 2n³/3 to factor, 2n² per triangular substitution.

n500k20

matrix 500×500 · right-hand sides k = 20

The advantage grows with $n$ and with the number of systems. For large $n$ the cost approaches that of a single factorisation, regardless of how many right-hand sides have to be processed.

If all the right-hand sides are known in advance, np.linalg.solve accepts a matrix as its second argument and factorises only once:

B = np.column_stack([b1, b2, b3])
X = np.linalg.solve(A, B)      # one factorisation, three solutions

When the matrix is symmetric: Cholesky

The lesson on the inverse and transpose established that $A^\top A$ is symmetric and positive semidefinite. For matrices with that structure a cheaper factorisation exists:

A = LL^\top

with $L$ lower triangular with positive diagonal. This is the Cholesky factorisation, and it requires $A$ to be symmetric and positive definite. It costs $\tfrac{1}{3}n^3$ operations, half of LU, because it exploits the symmetry rather than ignoring it.

L = np.linalg.cholesky(K)          # fails if K is not positive definite
alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))

The failure is informative rather than inconvenient: if Cholesky does not go through, the matrix is not positive definite, which in a Gaussian process usually signals an ill-conditioned covariance matrix. Standard practice is to add a regularisation term to the diagonal, $K + \sigma^2 I$ , which shifts every eigenvalue and restores positive definiteness.

When the system is rectangular: QR

For least squares, the previous lesson noted that forming $A^\top A$ squares the condition number. The alternative is to factorise $A$ directly:

A = QR

with $Q$ having orthonormal columns and $R$ upper triangular. Substituting into the normal equations, the problem reduces to $R\vec{w} = Q^\top\vec{b}$ , a triangular system, without ever constructing $A^\top A$ .

Q, R = np.linalg.qr(A)
w = np.linalg.solve(R, Q.T @ b)

This is, in essence, what np.linalg.lstsq does internally.

Factorisation	Requires	Cost	Used for
LU	square, non-singular	$\tfrac{2}{3}n^3$	general systems, determinant
Cholesky	symmetric positive definite	$\tfrac{1}{3}n^3$	covariances, Gaussian processes
QR	full column rank	$2mn^2$	least squares
SVD	none	$\sim 20n^3$	rank, null space, pseudoinverse

The reduced row echelon form

Elimination can be carried past $U$ to the reduced row echelon form, with pivots equal to one and zeros above them as well. NumPy does not provide it, and the omission is deliberate: for solving a system it adds nothing that $LU$ does not already give, and in floating point the decision of which entry counts as an exact zero is ambiguous.

Where it does belong is in exact algebra, with sympy:

import sympy as sp

M = sp.Matrix([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]])
M.rref()          # exact rational arithmetic, no rounding

The distinction is one of domain, not of quality. sympy works over exact rationals and suits determining structure — rank, a basis for the null space, dependencies among rows — in small matrices. numpy and scipy work in floating point over optimised implementations and are the tool for numerical computation. Using the second to reason about exact structure, or the first to solve a system of any size, inverts both.

Why it is not hand-coded

The routines shown delegate to LAPACK, a library with decades of debugging whose implementations account for the memory hierarchy: they operate in blocks to exploit the cache, which a direct transcription of the pseudocode does not.

Writing one's own elimination is an instructive exercise and a poor production decision. It will be one to two orders of magnitude slower and, more likely than not, less stable: the pivoting, the criterion for detecting singularity and the handling of edge cases are where the difficulty concentrates, and they are precisely the parts a first attempt leaves out.

Exercise. Factorise $A$ with scipy.linalg.lu and check that the product of the diagonal of U, signed by the number of interchanges, reproduces np.linalg.det(A). Then determine how many right-hand sides are needed for factorising once to be cheaper than calling solve repeatedly, with $n = 1000$ .