Prev
Unit 3 · Linear Systems: Direct Methods
Gaussian Elimination and LU Factorisation
60m100 XP
Next
🧮Linear Systems: Direct Methods

Gaussian Elimination and LU Factorisation

Row operations, triangular systems, A = LU, operation counts and many right-hand sides

Teddington, 1946
The Mathematics Division of Britain's National Physical Laboratory, a few months after the war. A system of 18 linear equations in 18 unknowns arrives to be solved. There is no electronic computer yet; the Division's tools are desk calculating machines, cranked by hand. Leslie Fox, Charles Goodwin and the young James Wilkinson — with Alan Turing, by some accounts, taking part — set to work with Gaussian elimination. It takes them about two weeks. Fox later recalled it as one of the experiences that convinced them to build an automatic computer. Two years later Turing published Rounding-off Errors in Matrix Processes. It recast elimination as the factorisation of a matrix into triangular pieces and introduced a term every engineer now uses: the condition number. Today a laptop solves those 18 equations in a few microseconds, and a structural model with a million unknowns overnight. The algorithm is the same one. What changed is that we now understand exactly what it costs, and exactly how accurate it is. (Cleve Moler tells the story.)
Linear systems are everywhere in engineering because linearity is the first approximation to everything. Nodal analysis of a circuit, the stiffness method for a truss, a finite-difference grid for heat flow, the Newton step of the last unit — all end in
Ax=b,A\mathbf x = \mathbf b,
with AA an n×nn \times n matrix of known coefficients, b\mathbf b the loads or sources, and x\mathbf x the unknown voltages, displacements or temperatures. Cramer's rule, taught in school, needs about (n+1)!(n+1)! operations: for n=20n = 20, longer than the age of the universe. Elimination needs about 23n3\tfrac23 n^3.
Gaussian elimination. Use row 1 to put zeros below the first diagonal entry: subtract mr1=ar1/a11m_{r1} = a_{r1}/a_{11} times row 1 from each row r>1r > 1. Then use the new row 2 to clear column 2 below the diagonal, and so on. After n−1n-1 stages the matrix is upper triangular, UU, and the right-hand side has been transformed the same way into c\mathbf c.
A triangular system is solved from the bottom up — back substitution:
xn=cnunn,xr=1urr(cr−∑j>rurj xj),r=n−1,…,1x_n = \frac{c_n}{u_{nn}}, \qquad x_r = \frac{1}{u_{rr}}\Big(c_r - \sum_{j>r} u_{rj}\,x_j\Big),\quad r = n-1, \dots, 1

Back substitution costs about n² operations — trivial next to the elimination that produced U.

💡Elimination is a factorisation: A = LU
Keep the multipliers mrkm_{rk} instead of throwing them away, and put them in a unit lower-triangular matrix LL (ones on the diagonal, mrkm_{rk} in position (r,k)(r, k)). Then, exactly, A=LU.A = LU. Elimination is the computation of LL and UU. Solving Ax=bA\mathbf x = \mathbf b then splits into two triangular solves: Ly=bL\mathbf y = \mathbf b (forward substitution) and Ux=yU\mathbf x = \mathbf y (back substitution). The polish astronomer Tadeusz Banachiewicz used triangular factorisations for hand computation in 1938; Turing's 1948 paper made the connection with elimination the standard way to think about it.
Derivation
Why the multipliers are exactly LL
Each stage of elimination is a multiplication on the left by a simple matrix. For a 3×33 \times 3 matrix, stage 1 subtracts m21m_{21} times row 1 from row 2 and m31m_{31} times row 1 from row 3:
E1=[100−m2110−m3101],E2=[1000100−m321],E2E1A=U.E_1 = \begin{bmatrix} 1 & 0 & 0 \\ -m_{21} & 1 & 0 \\ -m_{31} & 0 & 1 \end{bmatrix}, \qquad E_2 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -m_{32} & 1 \end{bmatrix}, \qquad E_2E_1A = U.
Undoing a stage means adding the rows back, so each inverse is the same matrix with the signs of the multipliers flipped. Multiply both sides by E1−1E2−1E_1^{-1}E_2^{-1}:
A=E1−1E2−1U,E1−1E2−1=[100m2110m3101][1000100m321]=[100m2110m31m321]=L.A = E_1^{-1}E_2^{-1}U, \qquad E_1^{-1}E_2^{-1} = \begin{bmatrix} 1 & 0 & 0 \\ m_{21} & 1 & 0 \\ m_{31} & 0 & 1 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & m_{32} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ m_{21} & 1 & 0 \\ m_{31} & m_{32} & 1 \end{bmatrix} = L.
In this order the product does no arithmetic at all: every multiplier simply lands in its own position. That is why storing the multipliers is computing LL, and it holds for any nn.
With A=LUA = LU, the system LUx=bLU\mathbf x = \mathbf b is solved in two triangular steps. First Ly=bL\mathbf y = \mathbf b from the top down — forward substitution, y1=b1y_1 = b_1 and yr=br−∑j<rlrjyjy_r = b_r - \sum_{j<r} l_{rj}y_j (no division, since LL has ones on its diagonal). Then Ux=yU\mathbf x = \mathbf y from the bottom up, by the back substitution formula above.
📐A resistor network
Four nodes, each joined to its neighbours and to ground by 1 Ω resistors, 1 A injected at node 1. Kirchhoff's current law gives Gv=iG\mathbf v = \mathbf i with conductance matrix G=[3−1−10−14−1−1−1−14−10−1−13].G = \begin{bmatrix} 3 & -1 & -1 & 0 \\ -1 & 4 & -1 & -1 \\ -1 & -1 & 4 & -1 \\ 0 & -1 & -1 & 3 \end{bmatrix}. The first stage uses multipliers −1/3,−1/3,0-1/3, -1/3, 0; the second −0.3636,−0.2727-0.3636, -0.2727; the last −0.6429-0.6429. The pivots — the diagonal of UU — are 3,3.6667,3.1818,2.14293, 3.6667, 3.1818, 2.1429. Back substitution gives node voltages 0.4667,0.2,0.2,0.13330.4667, 0.2, 0.2, 0.1333 V. The two middle voltages are equal, as the symmetry of the network demands.
We can only see a short distance ahead, but we can see plenty there that needs to be done.
– Alan Turing, Computing Machinery and Intelligence (1950)
Counting operations. At stage kk there are n−kn-k rows to update, each with about n−kn-k entries, at one multiply and one subtract each. Summing,
∑k=1n−12(n−k)2≈23n3versus2n2 for the two triangular solves\sum_{k=1}^{n-1} 2(n-k)^2 \approx \tfrac23 n^3 \qquad\text{versus}\qquad 2n^2 \text{ for the two triangular solves}

Doubling n multiplies the factorisation cost by 8 and the solve cost by 4. For n = 1000: about 7 × 10⁸ operations to factorise, 2 × 10⁶ to solve — a fraction of a second, then microseconds.

Derivation
Adding up the operations
At stage kk the rows below the pivot are k+1,…,nk + 1, \dots, n — that is n−kn - k rows — and each needs its entries in columns k+1,…,nk + 1, \dots, n updated: another n−kn - k numbers, one multiplication and one subtraction each. Substitute j=n−kj = n - k and use the formula for a sum of squares:
∑k=1n−12(n−k)2=2∑j=1n−1j2=2 (n−1) n (2n−1)6≈2n33.\sum_{k=1}^{n-1} 2(n - k)^2 = 2\sum_{j=1}^{n-1} j^2 = 2\,\frac{(n-1)\,n\,(2n-1)}{6} \approx \frac{2n^3}{3}.
(The multipliers themselves cost one division per row, only about n2/2n^2/2 in total.) A triangular solve is cheaper: row rr of back substitution needs n−rn - r multiplications and additions, and ∑r2(n−r)≈n2\sum_r 2(n - r) \approx n^2. Two triangular solves: about 2n22n^2.
⚠Never compute the inverse
x=A−1b\mathbf x = A^{-1}\mathbf b is notation, not an algorithm. Forming A−1A^{-1} costs about 2n32n^3 operations — three times an LU — and multiplying by it is less accurate than two triangular solves. For a sparse matrix it is far worse: A−1A^{-1} is generally completely full even when AA is almost empty. In NumPy use np.linalg.solve, or scipy.linalg.lu_factor and lu_solve when the matrix is reused.
🕰Older than Gauss
Chapter eight of The Nine Chapters on the Mathematical Art, a Chinese text parts of which date to about 150 BC, solves systems of up to five equations by exactly this sequence of row operations on a counting board. Gauss used the method in the early 1800s to fit the orbit of the asteroid Pallas by least squares, and his name stuck to it in the twentieth century.
💡Where you will use this
Structures and Geotechnics: the stiffness method assembles Ku=fK\mathbf u = \mathbf f and factorises KK once for every load case. Electrical circuits and Power Systems: nodal analysis and every Newton step of a load-flow solution. Attitude Control and GNSS Navigation: Kalman filters solve small dense systems thousands of times per second, where the 23n3\tfrac23 n^3 count decides what fits on the flight computer.
📚Further reading
  • Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lecture 20, 'Gaussian Elimination' — elimination as a product of lower-triangular matrices, exactly as derived here.
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 2, 'Systems of Linear Equations'.
  • Alan Turing, 'Rounding-off errors in matrix processes', Quarterly Journal of Mechanics and Applied Mathematics 1, 287–308 (1948) — the paper that named the LU decomposition and the condition number.
↓But there is a problem
Elimination divides by the pivot ukku_{kk} at every stage. Nothing in the algorithm checks how big it is. Try [10−17111]x=[12]\begin{bmatrix} 10^{-17} & 1 \\ 1 & 1 \end{bmatrix}\mathbf x = \begin{bmatrix} 1 \\ 2 \end{bmatrix}. The answer is almost exactly x1=x2=1x_1 = x_2 = 1. Elimination as written returns x1=0x_1 = 0. In 1943 the statistician Harold Hotelling estimated that rounding errors in elimination could grow so fast that large systems would be hopeless. He was wrong — but the reason he was wrong took twenty years to understand.