Chapter 16

Matrix Decompositions — LU

00 · Symbol Glossary

$L$L — lower triangular factor

A unit lower triangular matrix: ones on the main diagonal, zeros strictly above the diagonal. In the LU decomposition $A = LU$ , the entries below the diagonal of $L$ are the multipliers from Gaussian elimination — the scalars $c$ in row operations $R_i \to R_i - cR_j$ .

$U$U — upper triangular factor

An upper triangular matrix: zeros strictly below the main diagonal. In $A = LU$ , $U$ is the row echelon form of $A$ (before back-substitution scales pivots to 1). Solving $A\mathbf{x}=\mathbf{b}$ reduces to two easy triangular solves: $L\mathbf{y}=\mathbf{b}$ , then $U\mathbf{x}=\mathbf{y}$ .

$P$P — permutation matrix

A matrix obtained by permuting the rows of the identity $I_n$ . Each row and column has exactly one entry equal to 1, all others 0. When Gaussian elimination requires row swaps, the factorisation becomes $PA = LU$ — pivoting is encoded in $P$ , not in $L$ or $U$ .

$\ell_{ij}$ell sub i j — elimination multiplier

The multiplier stored at position $(i,j)$ of $L$ for $i > j$ . If elimination uses $R_i \to R_i - \ell_{ij} R_j$ , then $\ell_{ij}$ is the ratio of the entry being eliminated to the pivot above it: $\ell_{ij} = a_{ij}/a_{jj}$ at the moment of elimination.

01 · Why Factor a Matrix

Chapter 4 showed how to solve $A\mathbf{x}=\mathbf{b}$ once by Gaussian elimination. In quantitative work, the matrix $A$ is often fixed while the right-hand side $\mathbf{b}$ changes — stress scenarios, Monte Carlo paths, or successive time steps in a risk engine. Re-running full elimination for every new $\mathbf{b}$ wastes the work already invested in reducing $A$ .

The LU decomposition records elimination as a matrix product: $A = LU$ . The expensive step — reducing $A$ to triangular form — is done once. Each new $\mathbf{b}$ costs only two triangular solves, each $O(n^2)$ for an $n \times n$ system instead of $O(n^3)$ for full elimination.

Definition — LU Decomposition

A square matrix $A \in \mathbb{R}^{n \times n}$ has an LU decomposition if there exist matrices $L, U \in \mathbb{R}^{n \times n}$ such that:

A = LU

$L$ is unit lower triangular — $\ell_{ii} = 1$ for all $i$ , and $\ell_{ij} = 0$ for $i < j$ .

$U$ is upper triangular — $u_{ij} = 0$ for $i > j$ .

When no row interchanges are needed during elimination, $A$ admits a direct LU factorisation. When pivots require row swaps, the factorisation becomes $PA = LU$ with a permutation matrix $P$ .

✓ Example — Risk Engine with Many Right-Hand Sides

A $500 \times 500$ sensitivity matrix $A$ arises from a finite-difference approximation of a PDE used in counterparty exposure. For each of $10{,}000$ Monte Carlo scenarios, you must solve $A\mathbf{x}^{(k)} = \mathbf{b}^{(k)}$ .

Full Gaussian elimination per scenario: $10{,}000 \times O(500^3) \approx 1.25 \times 10^{12}$ flops.

LU once ( $O(500^3)$ ) plus $10{,}000$ triangular solves ( $2 \times 10{,}000 \times O(500^2) \approx 5 \times 10^9$ flops): dominated by the single factorisation. The LU approach is roughly $250\times$ faster — the difference between overnight and a week of compute.

❌ Failure — Re-Factoring When Only $\mathbf{b}$ Changes

Suppose $A = \begin{pmatrix}2&1\\4&3\end{pmatrix}$ is fixed and you need solutions for $\mathbf{b}_1 = \begin{pmatrix}3\\7\end{pmatrix}$ and $\mathbf{b}_2 = \begin{pmatrix}1\\5\end{pmatrix}$ .

Attempted computation: run Gaussian elimination from scratch on the augmented matrix $\left(\begin{array}{cc|c}2&1&3\\4&3&7\end{array}\right)$ , then again on $\left(\begin{array}{cc|c}2&1&1\\4&3&5\end{array}\right)$ .

Where it breaks: the coefficient matrix operations — eliminating the $(2,1)$ entry, updating rows 2 — are identical in both runs. Only the augmented column differs. Re-doing elimination on $A$ duplicates $O(n^3)$ work that LU stores permanently in $L$ and $U$ .

Consequence: for $m$ right-hand sides, naive elimination costs $O(mn^3)$ . LU plus $m$ solves costs $O(n^3 + mn^2)$ . When $m \gg n$ , this gap is decisive.

02 · LU from Gaussian Elimination

Every multiplier written down during elimination becomes an entry of $L$ . The upper triangular result of elimination is $U$ .

Definition — Elimination Multiplier

When Gaussian elimination eliminates the entry at position $(i,j)$ using pivot row $j$ (with $j < i$ ), the multiplier is:

\ell_{ij} = \frac{a_{ij}^{(\text{before})}}{a_{jj}^{(\text{pivot})}}

$a_{ij}^{(\text{before})}$ — the entry in row $i$ , column $j$ immediately before the elimination step.

$a_{jj}^{(\text{pivot})}$ — the pivot in row $j$ , column $j$ , which must be nonzero.

The row operation is $R_i \to R_i - \ell_{ij} R_j$ . The multiplier $\ell_{ij}$ is stored at position $(i,j)$ of $L$ .

Step-by-step — LU decomposition of $A = \begin{pmatrix}1&2&3\\4&5&6\\7&8&10\end{pmatrix}$

Identify the matrix: $A \in \mathbb{R}^{3 \times 3}$ . No zero pivot appears in column 1 — row swaps are not needed. Initialise $L = I_3$ (ones on the diagonal, zeros elsewhere) and work on a copy of $A$ to build $U$ .

Eliminate below pivot $(1,1) = 1$ : target column 1, rows 2 and 3.

$\ell_{21} = a_{21}/a_{11} = 4/1 = 4$ . Operation: $R_2 \to R_2 - 4R_1$ . Row 2 becomes $(4-4,\ 5-8,\ 6-12) = (0,\ -3,\ -6)$ .

$\ell_{31} = a_{31}/a_{11} = 7/1 = 7$ . Operation: $R_3 \to R_3 - 7R_1$ . Row 3 becomes $(7-7,\ 8-14,\ 10-21) = (0,\ -6,\ -11)$ .

Store $\ell_{21} = 4$ at $L_{21}$ , $\ell_{31} = 7$ at $L_{31}$ .

Eliminate below pivot $(2,2) = -3$ : target column 2, row 3 only.

$\ell_{32} = a_{32}/a_{22} = (-6)/(-3) = 2$ . Operation: $R_3 \to R_3 - 2R_2$ . Row 3 becomes $(0,\ -6-2(-3),\ -11-2(-6)) = (0,\ 0,\ 1)$ .

Store $\ell_{32} = 2$ at $L_{32}$ .

Read off $U$ from the reduced matrix:

U = \begin{pmatrix}1 & 2 & 3 \\ 0 & -3 & -6 \\ 0 & 0 & 1\end{pmatrix}

The three pivots sit on the diagonal:

1

-3

1

. No free variables —

A

is invertible.

Assemble $L$ from stored multipliers:

L = \begin{pmatrix}1 & 0 & 0 \\ 4 & 1 & 0 \\ 7 & 2 & 1\end{pmatrix}

Ones on the diagonal by construction. Below-diagonal entries are exactly the

\ell_{ij}

computed above.

Verify $A = LU$ entry by entry:

$(LU)_{11} = 1 \cdot 1 = 1 = a_{11}$ ✓.

$(LU)_{12} = 1 \cdot 2 = 2 = a_{12}$ ✓.

$(LU)_{21} = 4 \cdot 1 + 1 \cdot 0 = 4 = a_{21}$ ✓.

$(LU)_{22} = 4 \cdot 2 + 1 \cdot (-3) = 8 - 3 = 5 = a_{22}$ ✓.

$(LU)_{33} = 7 \cdot 3 + 2 \cdot (-6) + 1 \cdot 1 = 21 - 12 + 1 = 10 = a_{33}$ ✓.

All nine entries match. The factorisation is correct.

What $L$ Encodes

$L$ is not mysterious — it is a bookkeeping matrix for the elimination multipliers. Row $i$ of $L$ (below the diagonal) records every scalar by which row $i$ was reduced during elimination. $U$ is the echelon form. Together they reconstruct $A$ exactly.

03 · Solving Systems via Forward and Back Substitution

Once $A = LU$ is known, solving $A\mathbf{x} = \mathbf{b}$ splits into two triangular systems.

Definition — Forward and Back Substitution

Given $A = LU$ and $\mathbf{b} \in \mathbb{R}^n$ , solve $A\mathbf{x} = \mathbf{b}$ in two stages:

Stage 1 — Forward substitution: solve $L\mathbf{y} = \mathbf{b}$ for $\mathbf{y}$ . Since $L$ is lower triangular with unit diagonal, solve top to bottom: $y_1 = b_1$ , then $y_i = b_i - \sum_{j=1}^{i-1}\ell_{ij} y_j$ .

Stage 2 — Back substitution: solve $U\mathbf{x} = \mathbf{y}$ for $\mathbf{x}$ . Since $U$ is upper triangular, solve bottom to top: $x_n = y_n/u_{nn}$ , then $x_i = (y_i - \sum_{j=i+1}^{n} u_{ij} x_j)/u_{ii}$ .

Step-by-step — Solving $A\mathbf{x}=\mathbf{b}$ for $A$ above, $\mathbf{b}=\begin{pmatrix}6\\15\\28\end{pmatrix}$

Forward substitution — solve $L\mathbf{y}=\mathbf{b}$ :

$y_1 = b_1 = 6$ .

$y_2 = b_2 - \ell_{21} y_1 = 15 - 4(6) = 15 - 24 = -9$ .

$y_3 = b_3 - \ell_{31} y_1 - \ell_{32} y_2 = 28 - 7(6) - 2(-9) = 28 - 42 + 18 = 4$ .

\mathbf{y} = \begin{pmatrix}6\\-9\\4\end{pmatrix}

Back substitution — solve $U\mathbf{x}=\mathbf{y}$ :

$x_3 = y_3/u_{33} = 4/1 = 4$ .

$x_2 = (y_2 - u_{23} x_3)/u_{22} = (-9 - (-6)(4))/(-3) = (-9+24)/(-3) = 15/(-3) = -5$ .

$x_1 = (y_1 - u_{12} x_2 - u_{13} x_3)/u_{11} = (6 - 2(-5) - 3(4))/1 = (6+10-12)/1 = 4$ .

\mathbf{x} = \begin{pmatrix}4\\-5\\4\end{pmatrix}

Verify in the original system: $A\mathbf{x}$ :

Row 1: $1(4)+2(-5)+3(4) = 4-10+12 = 6 = b_1$ ✓.

Row 2: $4(4)+5(-5)+6(4) = 16-25+24 = 15 = b_2$ ✓.

Row 3: $7(4)+8(-5)+10(4) = 28-40+40 = 28 = b_3$ ✓.

❌ Failure — Solving $U\mathbf{x}=\mathbf{y}$ Top-Down

With $U = \begin{pmatrix}1&2&3\\0&-3&-6\\0&0&1\end{pmatrix}$ and $\mathbf{y}=\begin{pmatrix}6\\-9\\4\end{pmatrix}$ , attempting $x_1 = y_1/u_{11} = 6$ first is wrong.

Why it breaks: $u_{12}=2$ and $u_{13}=3$ mean $x_2$ and $x_3$ contribute to row 1. Until $x_2$ and $x_3$ are known, $x_1$ cannot be determined. Upper triangular systems must be solved from the bottom row upward.

Consequence: solving top-down gives a vector that satisfies row 1 by accident but violates rows 2 and 3. Back substitution is mandatory.

04 · Pivoting and the $PA = LU$ Factorisation

When a pivot is zero (or dangerously small), elimination must swap rows. The permutation is factored separately.

Definition — $PA = LU$ Factorisation

If Gaussian elimination on $A$ requires row interchanges, there exists a permutation matrix $P$ such that:

PA = LU

$P$ records the sequence of row swaps. Equivalently, $A = P^{-1}LU$ . Since $P$ is a permutation matrix, $P^{-1} = P^T$ — inversion is transposition.

To solve $A\mathbf{x}=\mathbf{b}$ : compute $PA = LU$ , then solve $L\mathbf{y} = P\mathbf{b}$ (forward), then $U\mathbf{x}=\mathbf{y}$ (back). The permutation $P$ reorders $\mathbf{b}$ to match the pivoted rows of $A$ .

✓ Example — Pivot Required

$A = \begin{pmatrix}0&1\\2&3\end{pmatrix}$ . Pivot in position $(1,1)$ is $0$ — elimination stalls.

Swap rows: $P = \begin{pmatrix}0&1\\1&0\end{pmatrix}$ , so $PA = \begin{pmatrix}2&3\\0&1\end{pmatrix}$ .

Now $\ell_{21} = 0/2 = 0$ , giving $L = I_2$ and $U = \begin{pmatrix}2&3\\0&1\end{pmatrix}$ .

$PA = LU$ ✓. Without the swap, no LU factorisation of $A$ exists without pivoting.

❌ Failure — Ignoring a Zero Pivot

For $A = \begin{pmatrix}0&2\\1&3\end{pmatrix}$ , attempting $\ell_{21} = 1/0$ is undefined.

Why it breaks: division by zero in the multiplier formula. The matrix may still be invertible — here $\det(A) = 0\cdot3 - 2\cdot1 = -2 \neq 0$ — but elimination must swap rows first so a nonzero pivot sits in position $(1,1)$ .

Consequence: LU without pivoting fails on matrices that are perfectly well-conditioned. Production solvers always use $PA=LU$ .

05 · Complexity and Quant Application

Definition — LU Computational Cost

For $A \in \mathbb{R}^{n \times n}$ :

LU factorisation: $O(n^3)$ flops — same order as one Gaussian elimination.

Each triangular solve: $O(n^2)$ flops.

$m$ right-hand sides: $O(n^3 + mn^2)$ total — factorise once, solve $m$ times.

Determinant via LU: $\det(A) = \det(L)\det(U) = 1 \cdot \prod_{i=1}^n u_{ii}$ (times $\det(P) = \pm 1$ if pivoting was used).

✓ Example — Greeks via Multiple Linear Systems

A delta-gamma approximation for a book of $n = 200$ risk factors requires solving $A\mathbf{x} = \mathbf{b}^{(k)}$ for each of $m = 50$ bump scenarios, where $A$ is the Hessian of mark-to-market value and each $\mathbf{b}^{(k)}$ is a gradient under a different stress.

LU factorisation of $A$ : one $O(200^3) = 8 \times 10^6$ flop pass.

50 solves: $50 \times 2 \times O(200^2) = 50 \times 80{,}000 = 4 \times 10^6$ flops.

Total $\approx 1.2 \times 10^7$ flops. Naive elimination: $50 \times 8 \times 10^6 = 4 \times 10^8$ flops — roughly $33\times$ more work. In a nightly risk batch processing thousands of instruments, LU is the standard backend for repeated solves.

06 · Practice Exercises

EXERCISE 16.1

Run Gaussian elimination without row swaps. Store each multiplier $\ell_{ij} = a_{ij}/a_{jj}$ in $L_{ij}$ . The final reduced matrix is $U$ .

$A = \begin{pmatrix}2&1&0\\4&-2&2\\-2&2&1\end{pmatrix}$ .

$\ell_{21} = 4/2 = 2$ . $R_2 \to R_2 - 2R_1$ : row 2 becomes $(0, -4, 2)$ .

$\ell_{31} = -2/2 = -1$ . $R_3 \to R_3 + R_1$ : row 3 becomes $(0, 3, 1)$ .

$\ell_{32} = 3/(-4) = -3/4$ . $R_3 \to R_3 + \frac{3}{4}R_2$ : row 3 becomes $(0, 0, 1 + \frac{3}{4}(2)) = (0, 0, 2.5)$ .

L = \begin{pmatrix}1&0&0\\2&1&0\\-1&-\frac{3}{4}&1\end{pmatrix}, \quad U = \begin{pmatrix}2&1&0\\0&-4&2\\0&0&2.5\end{pmatrix}

Verify $(LU)_{22} = 2(1) + 1(-4) = -2$ ✓. $(LU)_{33} = -1(0) - \frac{3}{4}(2) + 2.5 = -1.5 + 2.5 = 1$ ✓.

Find the LU decomposition of $A = \begin{pmatrix}2&1&0\\4&-2&2\\-2&2&1\end{pmatrix}$ . Show all three multipliers and verify $A = LU$ for at least three entries.

EXERCISE 16.2

After factorising, solve $L\mathbf{y}=\mathbf{b}$ top to bottom, then $U\mathbf{x}=\mathbf{y}$ bottom to top.

From Exercise 16.1: $L = \begin{pmatrix}1&0&0\\2&1&0\\-1&-\frac{3}{4}&1\end{pmatrix}$ , $U = \begin{pmatrix}2&1&0\\0&-4&2\\0&0&2.5\end{pmatrix}$ , $\mathbf{b}=\begin{pmatrix}1\\0\\3\end{pmatrix}$ .

Forward: $y_1 = 1$ . $y_2 = 0 - 2(1) = -2$ . $y_3 = 3 - (-1)(1) - (-\frac{3}{4})(-2) = 3 + 1 - 1.5 = 2.5$ .

Back: $x_3 = 2.5/2.5 = 1$ . $x_2 = (-2 - 2(1))/(-4) = -4/(-4) = 1$ . $x_1 = (1 - 1(1))/2 = 0/2 = 0$ .

$\mathbf{x} = \begin{pmatrix}0\\1\\1\end{pmatrix}$ . Check: $2(0)+1(1)+0(1)=1=b_1$ ✓; $4(0)+(-2)(1)+2(1)=0=b_2$ ✓.

Using $A = \begin{pmatrix}2&1&0\\4&-2&2\\-2&2&1\end{pmatrix}$ with LU factors from elimination, solve $A\mathbf{x}=\mathbf{b}$ where $\mathbf{b}=\begin{pmatrix}1\\0\\3\end{pmatrix}$ . Show forward and back substitution steps.

EXERCISE 16.3

Row swap first: $P$ exchanges rows 1 and 2. Factorise $PA$ , not $A$ . Then $\det(A) = \det(P)\det(L)\det(U) = (-1)\prod u_{ii}$ .

$A = \begin{pmatrix}0&1\\2&3\end{pmatrix}$ . $P = \begin{pmatrix}0&1\\1&0\end{pmatrix}$ , $PA = \begin{pmatrix}2&3\\0&1\end{pmatrix}$ .

$\ell_{21} = 0/2 = 0$ . $L = I_2$ , $U = \begin{pmatrix}2&3\\0&1\end{pmatrix}$ .

$\det(U) = 2 \cdot 1 = 2$ . $\det(P) = -1$ (one row swap). $\det(A) = (-1)(2) = -2$ .

Direct check: $\det(A) = 0\cdot3 - 1\cdot2 = -2$ ✓.

Find $PA=LU$ for $A = \begin{pmatrix}0&1\\2&3\end{pmatrix}$ . State $P$ , $L$ , $U$ , and compute $\det(A)$ from the factorisation.

EXERCISE 16.4

LU exists without pivoting only if every leading principal minor is nonzero. Equivalently, every upper-left $k \times k$ submatrix must have nonzero determinant.

$A = \begin{pmatrix}1&2\\2&4\end{pmatrix}$ . Leading $1 \times 1$ minor: $\det(1) = 1 \neq 0$ ✓. Leading $2 \times 2$ minor: $\det(A) = 1\cdot4 - 2\cdot2 = 0$ .

Zero determinant means rows are linearly dependent: row 2 = $2 \times$ row 1. Elimination gives $U = \begin{pmatrix}1&2\\0&0\end{pmatrix}$ — zero pivot on the diagonal. No LU factorisation with nonzero diagonal of $U$ exists.

For $\mathbf{b}=\begin{pmatrix}3\\6\end{pmatrix}$ : system is consistent ( $b_2 = 2b_1$ ) with infinitely many solutions. For $\mathbf{b}=\begin{pmatrix}3\\7\end{pmatrix}$ : inconsistent ( $7 \neq 2\cdot3$ ). LU cannot be used to solve either case reliably because $U$ is singular.

Explain why $A = \begin{pmatrix}1&2\\2&4\end{pmatrix}$ has no usable LU factorisation without pivoting. What happens when you try to solve $A\mathbf{x}=\mathbf{b}$ for $\mathbf{b}=\begin{pmatrix}3\\6\end{pmatrix}$ versus $\mathbf{b}=\begin{pmatrix}3\\7\end{pmatrix}$ ?

EXERCISE 16.5

Count flops: one factorisation is $\frac{2}{3}n^3$ ; each triangular solve is $n^2$ . Compare $m \cdot \frac{2}{3}n^3$ to $\frac{2}{3}n^3 + 2mn^2$ .

Let $n = 1000$ , $m = 500$ .

Naive: $m \cdot O(n^3) \approx 500 \times \frac{2}{3}(10^9) \approx 3.33 \times 10^{11}$ flops.

LU: $O(n^3) + 2m \cdot O(n^2) \approx \frac{2}{3}(10^9) + 2(500)(10^6) \approx 6.67 \times 10^8 + 10^9 \approx 1.67 \times 10^9$ flops.

Ratio: $3.33 \times 10^{11} / 1.67 \times 10^9 \approx 200$ . LU is roughly 200 times cheaper for these parameters.

Break-even: LU wins when $m > 1$ (one factorisation plus $m$ solves beats $m$ factorisations). The advantage grows linearly with $m$ .

A risk system has $n = 1000$ factors and $m = 500$ scenario right-hand sides. Argue quantitatively why storing an LU factorisation beats re-running Gaussian elimination for each scenario. State the break-even value of $m$ .

EXERCISE 16.6

Factorise $A$ once. For each $\mathbf{b}^{(k)}$ , forward then back substitute. Compare total flop count to $3m$ full eliminations.

$A = \begin{pmatrix}3&1\\2&4\end{pmatrix}$ . $\ell_{21} = 2/3$ . $U = \begin{pmatrix}3&1\\0&10/3\end{pmatrix}$ , $L = \begin{pmatrix}1&0\\2/3&1\end{pmatrix}$ .

$\mathbf{b}^{(1)} = \begin{pmatrix}4\\6\end{pmatrix}$ : $y_1=4$ , $y_2 = 6 - \frac{2}{3}(4) = 6 - 8/3 = 10/3$ . $x_2 = (10/3)/(10/3) = 1$ . $x_1 = (4-1)/3 = 1$ . $\mathbf{x}^{(1)} = \begin{pmatrix}1\\1\end{pmatrix}$ .

$\mathbf{b}^{(2)} = \begin{pmatrix}5\\14\end{pmatrix}$ : $y_1=5$ , $y_2 = 14 - \frac{2}{3}(5) = 14 - 10/3 = 32/3$ . $x_2 = (32/3)/(10/3) = 32/10 = 16/5$ . $x_1 = (5 - 16/5)/3 = (9/5)/3 = 3/5$ . $\mathbf{x}^{(2)} = \begin{pmatrix}3/5\\16/5\end{pmatrix}$ .

Factorisation cost: $O(8)$ flops for $n=2$ . Each solve: $O(4)$ flops. Total: $8 + 2(4) = 16$ flops. Two full eliminations: $2 \times 8 = 16$ flops — break-even at $m=2$ for tiny $n$ ; for $n=500$ , factorisation dominates once and solves are cheap.

A $2 \times 2$ sensitivity matrix $A = \begin{pmatrix}3&1\\2&4\end{pmatrix}$ is fixed. Solve $A\mathbf{x}=\mathbf{b}^{(k)}$ for $\mathbf{b}^{(1)}=\begin{pmatrix}4\\6\end{pmatrix}$ and $\mathbf{b}^{(2)}=\begin{pmatrix}5\\14\end{pmatrix}$ using one LU factorisation. Show both solutions and comment on the cost advantage for large $n$ .

07 · Summary

Term	Definition
LU decomposition	$A = LU$ with $L$ unit lower triangular, $U$ upper triangular
Multiplier $\ell_{ij}$	Scalar $a_{ij}/a_{jj}$ stored in $L$ during elimination
Forward substitution	Solve $L\mathbf{y}=\mathbf{b}$ top to bottom
Back substitution	Solve $U\mathbf{x}=\mathbf{y}$ bottom to top
$PA = LU$	Pivoting version; $P$ records row swaps
Complexity	Factorise $O(n^3)$ once; each solve $O(n^2)$
Determinant	$\det(A) = \det(P)\prod_i u_{ii}$

Next: Matrix Decompositions — QR — orthogonal-triangular factorisation via Gram-Schmidt, and why it delivers numerically stable least squares where normal equations fail.