Week 3: Elimination and LU Factorization

Solving Systems by Elimination#

A linear system $Ax = b$ asks for a vector $x$ whose linear combination of the columns of $A$ equals the right-hand side $b$. Gaussian elimination is the algorithm that turns this into an equivalent triangular system you can solve by back-substitution. The idea is simple: use each pivot to clear every entry below it, one column at a time, until only an upper triangular matrix $U$ remains.

Worked $3 \times 3$ example. Solve $Ax = b$ with

Write the augmented matrix $[A \mid b]$:

Column 1. The first pivot is $a_{11} = 2$. Multipliers for the rows below are

Replace row 2 by $\text{row}_2 - 2 \cdot \text{row}_1$ and row 3 by $\text{row}_3 - 4 \cdot \text{row}_1$:

Column 2. The second pivot is now $1$. The multiplier for row 3 is $\ell_{32} = 3/1 = 3$. Replace row 3 by $\text{row}_3 - 3 \cdot \text{row}_2$:

The coefficient matrix is now upper triangular:

Back-substitution. Solve $Ux = c$ from the bottom up:

So

Check: $A x = \begin{pmatrix} 2-1 \\ 4-3 \\ 8-7 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} = b$. The pivots of $U$ are $2$, $1$, and $2$ — the successive diagonal entries used to divide when eliminating below. If any pivot is zero, this process stalls; we return to that case in the last section.

Elimination did not change the solution set: each row operation is reversible (add back the multiple you subtracted), so $Ax = b$ and $Ux = c$ have exactly the same $x$.

Elimination as Matrix Multiplication#

Every elimination step is left-multiplication by an elementary matrix. Week 2 treated matrix multiplication as composition of linear maps; here each elementary matrix encodes one row operation, and the product of those matrices encodes the whole elimination.

To subtract $\ell$ times row $j$ from row $i$, multiply on the left by the matrix $E_{ij}$ that looks like the identity except for a single nonzero off-diagonal entry $-\ell$ in position $(i,j)$. For our $3 \times 3$ example:

In general, when the current pivot is $a_{jj}$ (after previous steps) and the entry below it is $a_{ij}$, the multiplier is

and $E_{ij}$ has $-\ell$ in position $(i,j)$. Applying all three steps to $A$ yields

Write $E = E_{32} E_{31} E_{21}$ for the product of all elementary matrices. Then the full elimination is the single matrix equation

The same $E$ acts on $b$: $E b = c$, so $Ux = c$ is equivalent to $Ax = b$. Because each $E_{ij}$ is lower triangular with $1$'s on the diagonal, so is $E$; each is invertible (put $+\ell$ where $-\ell$ was), so $E$ is invertible.

Small $2 \times 2$ check. For $A = \begin{pmatrix} 2 & 1 \\ 4 & 3 \end{pmatrix}$, the only multiplier is $\ell = a_{21}/a_{11} = 4/2 = 2$. Then

That multiplier is the quantity the exercise below asks for.

The LU Factorization#

From $EA = U$ we can recover $A$ by inverting $E$:

Define $L = E^{-1}$. Then

This is the LULower-Upper (matrix factorization) factorization of $A$ (when no row exchanges are needed). Because $E$ is unit lower triangular, so is $L = E^{-1}$: it has $1$'s on the diagonal and the elimination multipliers below the diagonal (with positive signs). For the $3 \times 3$ example,

The subdiagonal entries of $L$ are exactly $\ell_{21} = 2$, $\ell_{31} = 4$, $\ell_{32} = 3$ — the multipliers we computed while eliminating. You never need to form $E$ and invert it by hand: record the multipliers as you go and place them under the diagonal of $L$; the upper triangular result is $U$.

Verify $LU = A$. The $(3,3)$ entry of $LU$ is $4\cdot 1 + 3\cdot 1 + 1\cdot 2 = 9$, matching $A$; the rest of the product recovers $A$ column by column the same way.

Why factor once? Solving $Ax = b$ for many different right-hand sides reuses the same $L$ and $U$:

1. Forward substitution: solve $Ly = b$ for $y$ (cheap, because $L$ is triangular with unit diagonal).
2. Back-substitution: solve $Ux = y$ for $x$.

Factoring $A = LU$ costs $O(n^3)$ for an $n \times n$ matrix; each subsequent solve costs only $O(n^2)$. In robotics, simulation, and least-squares pipelines where $A$ is fixed and $b$ changes often, this amortization is the reason $LU$ is a workhorse factorization. Week 4 will also use the pivots on the diagonal of $U$: when no row exchanges occur, $\det(A)$ equals the product of those pivots.

Row Exchanges and PA=LU#

Plain elimination fails when a pivot position becomes zero (or is so small that dividing by it is numerically unsafe). You cannot form $\ell = a_{ij}/a_{jj}$ if $a_{jj} = 0$. The fix is a row exchange: swap the current row with a lower row that has a nonzero entry in the pivot column, then continue.

A row exchange is left-multiplication by a permutation matrix $P$ — the identity with two rows swapped (or a product of such swaps). For example, swapping rows 1 and 2 of a $2 \times 2$ system uses

After the necessary permutations, elimination proceeds without division by zero and produces $L$ and $U$ for the permuted matrix. The general factorization is

where $P$ is a permutation matrix, $L$ is unit lower triangular, and $U$ is upper triangular. Equivalently, $A = P^{-1} LU = P^\top LU$ (since $P^{-1} = P^\top$ for a permutation matrix).

When do you need $P$? Whenever a pivot is zero (or tiny) while a nonzero candidate sits below it. If an entire pivot column from the diagonal downward is zero, the matrix is singular: no sequence of row exchanges produces a full set of nonzero pivots, and $Ax = b$ either has no solution or infinitely many. Week 5 will connect that failure mode to rank and linear dependence.

Practical note. Production libraries often use partial pivoting (swap in the largest available entry below the diagonal) and return $P$, $L$, $U$ with $PA = LU$. That can reorder rows even when every pivot is already nonzero, so the printed $L$ and $U$ may differ from a hand calculation without row swaps — both factorizations are valid. The browser lab below builds $L$ and $U$ with plain elimination (no pivoting) in NumPy so the factors match the worked example above; for this $A$, $P = I$ and $PA = LU$ still holds.

Knowledge check#

Browser lab#

Compute an $LU$ factorization with NumPy-only elimination (same $A$ as the worked example) and verify $PA = LU$ with $P = I$.

import numpy as np

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

n = A.shape[0]
U = A.copy()
L = np.eye(n)
P = np.eye(n)  # no row exchanges needed for this A

for k in range(n - 1):
    for i in range(k + 1, n):
        L[i, k] = U[i, k] / U[k, k]
        U[i, k:] -= L[i, k] * U[k, k:]

print("L =\n", L)
print("U =\n", U)
print("P (identity) =\n", P)
print("L @ U =\n", L @ U)
print("PA ≈ LU:", np.allclose(P @ A, L @ U))
print("LU ≈ A:", np.allclose(L @ U, A))