LU Decomposition by Hand: Splitting a 3×3 Matrix into L·U

LU decomposition says that a square matrix A can be written as the product of two triangular matrices.

A   =   L   ·   U
3×3    3×3     3×3

Part	Shape	Meaning
`L`	lower triangular, 1s on the diagonal	records how many times each row was subtracted during elimination
`U`	upper triangular	the staircase matrix Gaussian elimination leaves behind

In one sentence: LU decomposition = Gaussian elimination, with the steps saved. U is the result of elimination, L is the log of what was done.

The problem

A = [ 2  1  1 ]
    [ 4  3  3 ]
    [ 8  7  9 ]

Find L and U such that A = LU.

The plan: eliminate, and keep the multipliers

Use row 1 to zero out the numbers below it in column 1, then use row 2 to zero out the number below it in column 2. Whenever you do “row i minus m times row j”, that m goes into position (i, j) of L.

m_ij = (the number to eliminate) / pivot      →  goes into L at (i, j)

Step 1: eliminate column 1

The pivot is a₁₁ = 2.

m₂₁ = 4 / 2 = 2     row2 − 2·row1 = [4 3 3] − [4 2 2] = [0 1 1]
m₃₁ = 8 / 2 = 4     row3 − 4·row1 = [8 7 9] − [8 4 4] = [0 3 5]

After that:

[ 2  1  1 ]
[ 0  1  1 ]
[ 0  3  5 ]

Step 2: eliminate column 2

The pivot is the 1 in row 2, column 2.

m₃₂ = 3 / 1 = 3     row3 − 3·row2 = [0 3 5] − [0 3 3] = [0 0 2]

That leaves U:

U = [ 2  1  1 ]
    [ 0  1  1 ]
    [ 0  0  2 ]

Step 3: put the multipliers into L

1s on the diagonal, and the three multipliers below it:

L = [ 1  0  0 ]     ← m₂₁ = 2 at (2,1)
    [ 2  1  0 ]     ← m₃₁ = 4 at (3,1), m₃₂ = 3 at (3,2)
    [ 4  3  1 ]

Check by multiplying back:

row 1 of LU = 1·[2 1 1]                          = [2 1 1]   ✓
row 2 of LU = 2·[2 1 1] + 1·[0 1 1]              = [4 3 3]   ✓
row 3 of LU = 4·[2 1 1] + 3·[0 1 1] + 1·[0 0 2]  = [8 7 9]   ✓

Why can the multipliers go straight into L? Row 3 of LU = 4·(row 1 of U) + 3·(row 2 of U) + 1·(row 3 of U). That is the elimination run backwards: elimination took 4 times row 1 and 3 times row 2 out of row 3, and now they are added back, which gives the original row 3. Each row of L says “this original row is made of these rows of U”.

Solving Ax = b with LU

Since A = LU, the system Ax = b becomes L(Ux) = b. Let y = Ux and do two steps:

Solve Ly = b. L is lower triangular, so substitute from the top down. This is forward substitution.
Solve Ux = y. U is upper triangular, so substitute from the bottom up. This is back substitution.

Take b = (4, 10, 24)ᵀ:

Forward, Ly = b:
  y₁ = 4
  2·y₁ + y₂ = 10          →  y₂ = 10 − 8 = 2
  4·y₁ + 3·y₂ + y₃ = 24   →  y₃ = 24 − 16 − 6 = 2

Back, Ux = y:
  2·x₃ = 2                →  x₃ = 1
  x₂ + x₃ = 2             →  x₂ = 1
  2·x₁ + x₂ + x₃ = 4      →  x₁ = 1

x = (1, 1, 1)ᵀ. Check: A·(1, 1, 1) = (2+1+1, 4+3+3, 8+7+9) = (4, 10, 24). Correct.

A free bonus: the determinant

det(A) = det(L) · det(U) = 1 · (2·1·2) = 4. The diagonal of L is all 1s, so det(A) is just the product of the diagonal of U.

When the pivot is 0: swap rows

A = [ 0  2 ]
    [ 1  3 ]

The first pivot is 0, so m₂₁ = 1/0 can’t be computed. The fix is to swap the two rows first, then eliminate:

PA = [ 0  1 ] [ 0  2 ]  =  [ 1  3 ]  =  [ 1  0 ] [ 1  3 ]  =  L · U
     [ 1  0 ] [ 1  3 ]     [ 0  2 ]     [ 0  1 ] [ 0  2 ]

P is a permutation matrix that records which rows were swapped. The decomposition is now written PA = LU.

A pivot that is tiny but not zero is also a problem, because dividing by it magnifies rounding error. So software by default picks the entry with the largest absolute value in the current column as the pivot (partial pivoting). In practice, the decomposition you actually get is almost always PA = LU.

Where it goes wrong

Putting m in the wrong place. m₃₂ is “row 3 minus a multiple of row 2”, so it goes at (3,2), not (2,3).
Computing m from the original matrix. The 3 in the numerator of step 2 is a₃₂ after step 1, not the 7 in the original.
Forgetting that the diagonal of L is 1. The diagonal values have already been absorbed into U.
Mixing up the direction of substitution. L goes top-down, U goes bottom-up.

Checking with NumPy

import numpy as np
from scipy.linalg import lu

A = np.array([[2, 1, 1], [4, 3, 3], [8, 7, 9]], dtype=float)
P, L, U = lu(A)          # scipy's convention: A = P @ L @ U

P    # [[0. 1. 0.]
     #  [0. 0. 1.]
     #  [1. 0. 0.]]
L    # [[1.     0.     0.]
     #  [0.25   1.     0.]
     #  [0.5    0.6667 1.]]
U    # [[ 8.    7.    9.   ]
     #  [ 0.   -0.75 -1.25 ]
     #  [ 0.    0.   -0.6667]]

The result is different from the hand calculation, and that’s expected: scipy uses partial pivoting, the largest number in column 1 is 8, so it moves row 3 to the top first. The decomposition isn’t unique, but both multiply back to A, and the determinant agrees: 8·(−0.75)·(−2/3) = 4, times the sign of P (+1, since the three rows are cyclically shifted, an even permutation).

Two small traps:

scipy’s convention is A = P @ L @ U, not PA = LU. The two differ by a transpose: P becomes Pᵀ.
If you want L and U without any row swaps, you have to write it yourself. It’s a few lines: for each column j and each row i below it, compute m = U[i,j] / U[j,j], store it in L[i,j], then do U[i] -= m * U[j].

What it’s for

A = LU is not the goal by itself. Its value is this: do the expensive part, elimination, once, and reuse it.

1. The same A, many b’s

The direct way to solve Ax = b is to eliminate from scratch every time, about n³/3 multiplications each. With LU, elimination happens once, and each new b costs one forward and one back substitution, about n² in total.

n = 1000, 100 different b’s:

Approach	Work (order of magnitude, multiplications)
Re-eliminate every time	100 × 3.3×10⁸ ≈ 3.3×10¹⁰
One LU + 100 substitutions	3.3×10⁸ + 100 × 10⁶ ≈ 4.3×10⁸

That’s roughly 80× less. The situation is common: in a time-stepping loop the matrix stays the same and only the right-hand side changes, and iterative algorithms often solve the same system every round.

2. Inverting a matrix is solving n systems

Column j of A⁻¹ is the solution of Ax = eⱼ (column j of the identity). So inverting = one LU + n substitutions.

The more important lesson runs the other way: if all you want is A⁻¹b, don’t compute the inverse. Solve with LU. It’s faster and numerically more stable. np.linalg.solve(A, b) does LU internally, not invert-then-multiply.

3. The determinant is nearly free

det(A) is the product of the diagonal of U, times +1 or −1 depending on how many row swaps happened. Cofactor expansion is O(n!); LU is O(n³).

4. It’s a building block for other algorithms

The normal equations in linear regression, (XᵀX)β = Xᵀy, are a linear system. XᵀX is symmetric positive definite, so you can use Cholesky, the symmetric version of LU, which is about twice as fast.
Newton’s method solves J·Δ = −f at every step, with J the Jacobian. If the same J is reused for several steps (variants such as the chord method), it is factored only once.
Gaussian processes and Kalman filters both have a “solve a system with the covariance matrix” step, and underneath it’s Cholesky or LU.
Implicit time-stepping for differential equations solves a system of the same shape at every step.

In one sentence: whenever you see “solve Ax = b” and A shows up more than once, think LU.

Practice

By hand. For A = [[1,2,1],[3,8,7],[2,7,9]], find L and U, multiply back to verify, and compute det(A).
Solve. Use the L and U from problem 1 to solve Ax = (1, 2, 3)ᵀ: forward, then back, then plug into A to check.
Code. Write my_lu(A) in plain NumPy that returns (L, U) without row swaps, and test it on this article’s example with np.allclose(L @ U, A). Then try the matrix from the “When the pivot is 0” section and see what happens.

Answers (try first, then look)

m₂₁ = 3, m₃₁ = 2; after column 1, rows 2 and 3 are [0,2,4] and [0,3,7]; m₃₂ = 3/2. L = [[1,0,0],[3,1,0],[2,1.5,1]], U = [[1,2,1],[0,2,4],[0,0,1]], det = 1·2·1 = 2.
Forward gives y = (1, −1, 2.5), back gives x = (9.5, −5.5, 2.5). Check: A·x = (1, 2, 3).
The zero-pivot matrix divides by zero at U[j,j] == 0 and produces inf or nan.