BrevFeed
Sign in

6. Elementary matrices and LU decomposition

📖 Reading · 11 min
💡 Every code box below is live — edit it and hit Run.

Every lesson in this module has treated row operations as actions performed on a matrix. This closing lesson reveals they're secretly multiplications by specific matrices — which both explains several claims taken on faith earlier (why row operations affect the determinant the way §16.4 said, why Gauss-Jordan produces A^{-1}) and builds the factorization, LU decomposition, that real solvers use instead of raw elimination.

Elementary matrices

An elementary matrix is what you get by applying a single row operation (§16.1) to the identity matrix I. There are three kinds, matching the three operations:

  • Row swap: I with two rows swapped.
  • Row scale: I with one diagonal entry replaced by c\neq0.
  • Row addition: I with a single off-diagonal entry c added.

The claim: left-multiplying A by an elementary matrix E performs that row operation on A. For example, in 2\times2,

E=\begin{pmatrix}1&0\\-3&1\end{pmatrix}\quad\Longrightarrow\quad EA=\begin{pmatrix}1&0\\-3&1\end{pmatrix}\begin{pmatrix}a&b\\c&d\end{pmatrix}=\begin{pmatrix}a&b\\c-3a&d-3b\end{pmatrix}

— row 2 minus 3\timesrow 1, exactly the row operation E encodes, applied to every column of A at once (consistent with §16.2's column-by-column view of multiplication). Every elementary matrix is invertible (its own inverse undoes the same operation: swap undoes itself, scale by c undoes with 1/c, add c\timesrow j undoes by subtracting it back out) — which is why row operations never change a system's solution set (§16.1) or a matrix's invertibility (§16.3): each step just left-multiplies by something invertible, and invertible matrices don't create or destroy solutions.

Gaussian elimination, restated: reducing A to row-echelon form U is exactly

E_k\cdots E_2E_1A=U

for some sequence of elementary matrices — and §16.3's Gauss-Jordan inverse computation is precisely this: apply the same E_k\cdots E_1 to I, and since E_k\cdots E_1A=I when A row-reduces fully, that product E_k\cdots E_1 is A^{-1}.

LU decomposition

Collect every elementary matrix used during forward elimination only (no row swaps, and only the "add a multiple of one row to another" step that clears below the pivots) into L=(E_k\cdots E_1)^{-1}. Since each E_i was lower-triangular in this no-swap case, so is its product, and so is L — and it turns out to have 1s on its diagonal with the negatives of the multipliers used sitting directly below, at exactly the positions that were zeroed out. This gives the LU decomposition:

A=LU

where U is the upper-triangular row-echelon form from elimination, and L is unit lower-triangular (diagonal all 1s). Every matrix that can be row-reduced without row swaps has one (a matrix needing a swap has a closely related PA=LU, with P a permutation matrix recording the swap — not developed further here).

Why bother, when Gaussian elimination already works? Because solving A\vec x=\vec b becomes two cheap triangular solves instead of one elimination:

A\vec x=\vec b\ \Longrightarrow\ LU\vec x=\vec b\ \Longrightarrow\ \begin{cases}L\vec y=\vec b&\text{(forward substitution)}\\U\vec x=\vec y&\text{(back substitution)}\end{cases}

and critically, L and U don't depend on \vec b at all — compute them once, then solve for as many different right-hand sides as needed almost for free, which is exactly the situation software faces (e.g. re-solving a system with many different loads or boundary conditions). This is also, in practice, how \det A gets computed for anything larger than 3\times3: \det A=\det L\cdot\det U=1\cdot(\text{product of }U\text{'s diagonal}), since a triangular matrix's determinant is just its diagonal product — far cheaper than cofactor expansion.

Doing it in Python

import numpy as np
import scipy.linalg as la

A = np.array([[2., 1., 1.], [4., 3., 3.], [8., 7., 9.]])

P, L, U = la.lu(A)
print("L =")
for row in L.tolist():
    print([round(v, 4) for v in row])
print("U =")
for row in U.tolist():
    print([round(v, 4) for v in row])

reconstructed = P @ L @ U
print("\nP @ L @ U == A:", bool(np.allclose(reconstructed, A)))

L =
[1.0, 0.0, 0.0]
[0.25, 1.0, 0.0]
[0.5, 0.6667, 1.0]
U =
[8.0, 7.0, 9.0]
[0.0, -0.75, -1.25]
[0.0, 0.0, -0.6667]

P @ L @ U == A: True

(scipy's lu pivots for numerical stability, so L/U here don't match a strict top-to-bottom, no-swap elimination exactly — but PLU always reconstructs A, and the two triangular-solve idea is identical.)

Solving with a precomputed factorization, several right-hand sides at almost no extra cost:

import numpy as np
import scipy.linalg as la

A = np.array([[2., 1., 1.], [4., 3., 3.], [8., 7., 9.]])
lu, piv = la.lu_factor(A)  # compute once

for b in ([4., 10., 24.], [1., 2., 3.], [0., 0., 6.]):
    x = la.lu_solve((lu, piv), np.array(b))
    print(f"b={b} -> x={[round(v, 4) or 0.0 for v in x.tolist()]}")

b=[4.0, 10.0, 24.0] -> x=[1.0, 1.0, 1.0]
b=[1.0, 2.0, 3.0] -> x=[0.5, 0.5, -0.5]
b=[0.0, 0.0, 6.0] -> x=[0.0, -3.0, 3.0]

Worked example

Find the LU decomposition (no pivoting needed) of A=\begin{pmatrix}1&2\\3&8\end{pmatrix}, then use it to solve A\vec x=(5,19).

Eliminate: R_2\leftarrow R_2-3R_1 turns A into U=\begin{pmatrix}1&2\\0&2\end{pmatrix}, using multiplier 3. L has 1s on the diagonal and the multiplier (not its negative — L records what must be added back to reverse the elimination) below:

L=\begin{pmatrix}1&0\\3&1\end{pmatrix},\qquad U=\begin{pmatrix}1&2\\0&2\end{pmatrix}

Check: LU=\begin{pmatrix}1&0\\3&1\end{pmatrix}\begin{pmatrix}1&2\\0&2\end{pmatrix}=\begin{pmatrix}1&2\\3&6+2\end{pmatrix}=\begin{pmatrix}1&2\\3&8\end{pmatrix}=A ✓.

Solve L\vec y=(5,19) (forward substitution): y_1=5; 3(5)+y_2=19\Rightarrow y_2=4.

Solve U\vec x=(5,4) (back substitution): 2x_2=4\Rightarrow x_2=2; x_1+2(2)=5\Rightarrow x_1=1.

\boxed{\vec x=(1,2)}

Sanity check. A\vec x=(1(1)+2(2),\,3(1)+8(2))=(5,19) ✓ — matching the target exactly, and confirming both triangular solves were done correctly without ever re-deriving U from scratch.

Your turn

1. Write the elementary matrix E that performs "R_2\leftarrow R_2+2R_1" on a 2\times2 matrix, and compute E\begin{pmatrix}1&3\\-2&0\end{pmatrix} to confirm.

2. For A=\begin{pmatrix}2&0\\4&3\end{pmatrix}, find L and U (one elimination step; identify the multiplier directly from A's entries).

3. True or false: if A=LU, then \det A equals the product of U's diagonal entries only, regardless of L.

Solutions

1. E=\begin{pmatrix}1&0\\2&1\end{pmatrix}. E\begin{pmatrix}1&3\\-2&0\end{pmatrix}=\begin{pmatrix}1&3\\2(1)+(-2)&2(3)+0\end{pmatrix}=\begin{pmatrix}1&3\\0&6\end{pmatrix} — row 2 is now the original row 2 plus 2\timesrow 1, exactly as claimed: (-2,0)+2(1,3)=(-2+2,\,0+6)=(0,6) ✓.

2. Multiplier to clear the 4: R_2\leftarrow R_2-2R_1 gives U=\begin{pmatrix}2&0\\0&3\end{pmatrix}, with L=\begin{pmatrix}1&0\\2&1\end{pmatrix} (the multiplier 2 recorded, not negated). Check: LU=\begin{pmatrix}2&0\\4&0+3\end{pmatrix}=\begin{pmatrix}2&0\\4&3\end{pmatrix}=A ✓.

3. True, since L is unit lower-triangular. \det A=\det L\cdot\det U, and a triangular matrix's determinant is its diagonal product — L's diagonal is all 1s by construction, so \det L=1 always, leaving \det A=\det U regardless of what L's off-diagonal multipliers happen to be.

Check yourself in code

Compute the LU decomposition of A=\begin{pmatrix}6&2\\3&4\end{pmatrix} (no pivoting needed — column 1's largest entry, 6, is already on top), and use it to compute \det A as the product of U's diagonal.

Print exactly this:

L =
[1.0, 0.0]
[0.5, 1.0]
U =
[6.0, 2.0]
[0.0, 3.0]
det via LU: 18.0
det via numpy: 18.0

import numpy as np
import scipy.linalg as la

A = np.array([[6., 2.], [3., 4.]])
P, L, U = la.lu(A)

print("L =")
for row in L.tolist():
    print(row)
# print U the same way, then det via the product of U's diagonal,
# then np.linalg.det(A) for comparison

import numpy as np
import scipy.linalg as la

A = np.array([[6., 2.], [3., 4.]])
P, L, U = la.lu(A)

print("L =")
for row in L.tolist():
    print(row)
print("U =")
for row in U.tolist():
    print(row)

det_lu = U[0, 0] * U[1, 1]
print("det via LU:", round(float(det_lu), 4))
print("det via numpy:", round(float(np.linalg.det(A)), 4))

Every row operation is left-multiplication by an invertible elementary matrix, which is why row operations preserve a system's solutions and why Gauss-Jordan produces A^{-1}. Collecting the operations from forward elimination gives A=LU: a unit lower-triangular L and an upper-triangular U, turning any solve into two cheap triangular passes and reusable across many right-hand sides.

This closes Module 16. Every operation here — elimination, matrix arithmetic, inverses, determinants — was performed on lists of numbers in \mathbb{R}^n. Module 17 steps back and asks what makes those operations work at all, by stripping away the numbers and defining a vector space from the properties alone.

← Previous Next →