LU Decomposition

The last lesson promised that the row-elimination from lesson four could be bottled and reused. Here is the small annoyance that makes that worth doing: Gaussian elimination does a pile of work to triangulate a matrix, and then we throw every bit of it away. LU decomposition is elimination that keeps the receipt. You split A into two triangular pieces, A = L · U. The U is the upper-triangular echelon form — the staircase shape with nothing but zeros below the diagonal — that elimination already hands you; the L is a lower-triangular matrix (with 1s on its diagonal) that records the multipliers you used on the way.

The factorisation

The construction adds nothing to the elimination you already run. When you clear an entry with R2 → R2 − m·R1, the number m is a multiplier; U is what A becomes once all such operations are done, and L is the identity with each multiplier m slotted into the exact position it cleared. To then solve Ax = b, rewrite it as L(Ux) = b and split into two triangular solves:

• Forward-substitution on Ly = b (top row down) recovers y.
• Back-substitution on Ux = y (bottom row up) recovers x.

Each substitution costs about O(n²), while the factorisation costs O(n³) and is done once. So for k different right-hand sides the bill is one O(n³) factorisation plus k cheap O(n²) solves — far better than re-eliminating k times. That reuse is the whole reason LU exists, and why a solver like scipy.linalg.solve factors once under the hood.

A worked example

Factor A = [[2, 1], [4, 3]]:

A = [ 2  1 ]     clear (2,1): multiplier m = 4 / 2 = 2,  R2 → R2 − 2·R1
    [ 4  3 ]

U = [ 2  1 ]     L = [ 1  0 ]     (slot m = 2 into position (2,1); keep 1s on L's diagonal)
    [ 0  1 ]         [ 2  1 ]

verify L·U:  [ 1 0 ][ 2 1 ]   [ 2          1        ]   [ 2  1 ]
             [ 2 1 ][ 0 1 ] = [ 2·2+1·0   2·1+1·1   ] = [ 4  3 ] = A  ✓

So L = [[1, 0], [2, 1]] and U = [[2, 1], [0, 1]] — the single multiplier 4/2 = 2 lives in L, the pivots stay in U.

One sign catches almost everyone, so name it plainly: the row operation was a subtraction, R2 → R2 − 2·R1, yet the entry stored in L is +2, not −2. The reason is what L is for. U is A with the subtraction already carried out, so L has to be the matrix that puts it back. Read the product row by row and you can watch it happen: row 2 of L·U is 2·(row 1 of U) + 1·(row 2 of U) = 2·[2, 1] + [0, 1] = [4, 3], which is row 2 of the original A. L re-adds exactly what elimination took away. So store the multiplier with the sign you used to divide (4/2 = 2), never the sign you used to subtract.

A question to carry forward

LU broke a matrix into two triangular factors to tame it. Sometimes, though, it is the matrix itself we want to break — chop a big one into a grid of smaller rectangular blocks and treat each block as a single entry. Here is the thread onward: can you add and multiply matrices block by block as if the blocks were ordinary numbers, and when does that shortcut hold?

In one breath

• LU: A = L·U — elimination that keeps the receipt; U = the upper-triangular echelon form, L = lower-triangular multipliers with a unit diagonal.
• Build it free: the multiplier m from R2 → R2 − m·R1 slots into L at the position it cleared; U is the eliminated A.
• The win is reuse: factor once at O(n³), then each b is a forward solve Ly=b + back solve Ux=y, both O(n²). Great for many right-hand sides.
• Don’t put pivots in L (they live in U); L’s diagonal is all 1s.
• Zero pivot ⇒ plain LU fails; pivot (swap rows) and factor PA = LU — also used for stability.

Practice