Chapter 2 · Matrices

2.3 LU Factorisation

Record the elimination once as $A=LU$ — then every new right-hand side costs only two quick triangular solves, $L\mathbf{y}=\mathbf{b}$ and $U\mathbf{x}=\mathbf{y}$.

A structural engineer has a model of a bridge with a few thousand unknown joint displacements. The stiffness matrix $A$ depends only on the bridge — its beams, joints and materials. The load vector $\mathbf{b}$ changes all the time: rush-hour traffic, a lorry parked on one span, wind from the north, wind from the south, snow. Each load case is a new system $A\mathbf{x}=\mathbf{b}$ with the same matrix. Running Gaussian elimination from scratch for every load would repeat exactly the same row operations on $A$ hundreds of times.

The cure is to do the elimination on $A$ once and record it. The record takes the form of a product $A=LU$ of a lower triangular matrix $L$ (which stores the row operations) and an upper triangular matrix $U$ (the echelon form). After that, each new right-hand side is handled by two triangular systems, which are very cheap to solve. This LU factorisation is how computers actually solve linear systems — it is the engine inside the "solve" button of every engineering, statistics and graphics package.

Triangular Systems Are Easy

📘 Definition 2.9 — Triangular matrices

A square matrix $U$ is upper triangular if every entry below the main diagonal is zero ($u_{ij}=0$ whenever $i\gt j$). A square matrix $L$ is lower triangular if every entry above the main diagonal is zero ($\ell_{ij}=0$ whenever $i\lt j$). A lower triangular matrix with every diagonal entry equal to $1$ is called unit lower triangular. $$U=\begin{bmatrix}3&-1&2\\0&-2&1\\0&0&5\end{bmatrix}\qquad L=\begin{bmatrix}1&0&0\\2&1&0\\-3&-2&1\end{bmatrix}\ \text{(unit lower triangular)}$$

A triangular system needs no elimination at all. In $U\mathbf{x}=\mathbf{y}$ the last equation contains only $x_n$; solve it, substitute upwards, and each equation gives one new unknown. This is back substitution, which you met at the end of Gaussian elimination in Section 1.2. A lower triangular system $L\mathbf{y}=\mathbf{b}$ is solved the same way but from the top down — forward substitution. For example, $$\begin{bmatrix}1&0&0\\2&1&0\\-3&-2&1\end{bmatrix}\begin{bmatrix}y_1\\y_2\\y_3\end{bmatrix}=\begin{bmatrix}-1\\-7\\8\end{bmatrix}\quad\Longrightarrow\quad\begin{aligned}y_1&=-1,\\ y_2&=-7-2y_1=-5,\\ y_3&=8+3y_1+2y_2=-5.\end{aligned}$$ Because $L$ has $1$s on its diagonal there is not even any dividing to do.

Finding $L$ and $U$ by Elimination

📘 Definition 2.10 — LU factorisation

An LU factorisation of a square matrix $A$ is a factorisation $$A=LU$$ in which $L$ is unit lower triangular and $U$ is upper triangular. (The same idea works for an $m\times n$ matrix: then $L$ is $m\times m$ unit lower triangular and $U$ is an $m\times n$ echelon form.)

Where do $L$ and $U$ come from? Reduce $A$ to an echelon form using only the replacement operation "subtract a multiple of a higher row from a lower row", $R_i\to R_i-m_{ij}R_j$ with $j\lt i$. The number $m_{ij}$ is called the multiplier; it is $\dfrac{\text{entry to be eliminated}}{\text{pivot}}$. The echelon form we reach is $U$. Each operation is multiplication on the left by a unit lower triangular elementary matrix $E$ (Section 2.2), so $$E_k\cdots E_2E_1A=U\qquad\Longrightarrow\qquad A=\left(E_1^{-1}E_2^{-1}\cdots E_k^{-1}\right)U=LU .$$ The inverse of "subtract $m$ times row $j$ from row $i$" is "add $m$ times row $j$ to row $i$", and when we eliminate column by column the product of these inverses simply drops every multiplier into its own position. That is the remarkable fact behind the method:

📐 Theorem 2.12 — The multipliers form $L$

Let $A$ be an $n\times n$ matrix that can be reduced to an upper triangular matrix $U$ by row replacements $R_i\to R_i-m_{ij}R_j$ with $j\lt i$ only — no row interchanges and no scaling — carried out column by column from left to right. Then $A=LU$, where $L$ is the unit lower triangular matrix whose entry in position $(i,j)$ below the diagonal is the multiplier $m_{ij}$ (and $0$ where the entry was already zero and no operation was needed): $$L=\begin{bmatrix}1&0&0&\cdots&0\\ m_{21}&1&0&\cdots&0\\ m_{31}&m_{32}&1&\cdots&0\\ \vdots&\vdots&&\ddots&\vdots\\ m_{n1}&m_{n2}&\cdots&m_{n,n-1}&1\end{bmatrix}.$$

🧭 Problem-solving strategy — computing $A=LU$

Start with $L=I$ (diagonal of $1$s, unknown entries below).
Work down column $1$: for each row $i\gt 1$, compute $m_{i1}=a_{i1}/a_{11}$, do $R_i\to R_i-m_{i1}R_1$ and write $m_{i1}$ into position $(i,1)$ of $L$.
Move to column $2$ of the new matrix and repeat with pivot in position $(2,2)$, then column $3$, and so on. If a pivot is $0$, stop — a row interchange is needed (see the last part of this section).
The final matrix is $U$. Check by multiplying $LU$ — you must get $A$ back.

✏️ Example 1 — An LU factorisation of a $3\times3$ matrix

Find an LU factorisation of $A=\begin{bmatrix}3&-1&2\\6&-4&5\\-9&7&-3\end{bmatrix}$.

💡 Show solution▶

Column 1 The pivot is $3$. Multipliers $m_{21}=\tfrac63=2$ and $m_{31}=\tfrac{-9}{3}=-3$. Do $R_2\to R_2-2R_1$ and $R_3\to R_3-(-3)R_1=R_3+3R_1$: $$\begin{bmatrix}3&-1&2\\6&-4&5\\-9&7&-3\end{bmatrix}\longrightarrow\begin{bmatrix}3&-1&2\\0&-2&1\\0&4&3\end{bmatrix},\qquad L=\begin{bmatrix}1&0&0\\2&1&0\\-3&\ast&1\end{bmatrix}.$$

Column 2 The pivot is now $-2$. Multiplier $m_{32}=\tfrac{4}{-2}=-2$; do $R_3\to R_3-(-2)R_2=R_3+2R_2$: $$\begin{bmatrix}3&-1&2\\0&-2&1\\0&4&3\end{bmatrix}\longrightarrow\begin{bmatrix}3&-1&2\\0&-2&1\\0&0&5\end{bmatrix}=U,\qquad L=\begin{bmatrix}1&0&0\\2&1&0\\-3&-2&1\end{bmatrix}.$$

Answer $$A=LU=\begin{bmatrix}1&0&0\\2&1&0\\-3&-2&1\end{bmatrix}\begin{bmatrix}3&-1&2\\0&-2&1\\0&0&5\end{bmatrix}.$$

Check Row 3 of $LU$ is $-3(3,-1,2)-2(0,-2,1)+1(0,0,5)=(-9,\,3+4,\,-6-2+5)=(-9,7,-3)$, which is row 3 of $A$; rows 1 and 2 check the same way. ✓

⚠️ Watch the signs of the multipliers

1. $L$ stores $m$, not $-m$. In Example 1 we added $3R_1$ to $R_3$. Written in the standard form $R_3\to R_3-m_{31}R_1$, that is $m_{31}=-3$, and $-3$ is what goes into $L$. A quick test: $m_{ij}$ always equals (entry being eliminated) ÷ (pivot), signs included.

2. Use the current matrix. The multiplier $m_{32}$ is computed from the matrix after column 1 has been cleared ($4\div(-2)$), not from the original $A$ ($7\div(-4)$).

3. No scaling. Do not divide a row by its pivot to make a leading $1$ — that would change the diagonal of $L$. The pivots stay in $U$.

Each step performs one row replacement $R_i\to R_i-m_{ij}R_j$ on the right-hand matrix. The pivot is circled, the row being changed is highlighted, and the multiplier $m_{ij}$ drops into the same position of $L$. When the right-hand matrix is upper triangular it is $U$, and $LU=A$.

Matrix:

step—

row operation—

multiplier stored in $L$—

check $LU=A$—

Solving $A\mathbf{x}=\mathbf{b}$ with $L$ and $U$

Suppose $A=LU$. Then $A\mathbf{x}=\mathbf{b}$ reads $L(U\mathbf{x})=\mathbf{b}$. Give the bracket a name, $\mathbf{y}=U\mathbf{x}$, and the single system splits into two triangular ones.

📐 Theorem 2.13 — Solving with an LU factorisation

Let $A=LU$ with $L$ unit lower triangular and $U$ upper triangular with nonzero diagonal entries. Then for every $\mathbf{b}$ the system $A\mathbf{x}=\mathbf{b}$ has exactly one solution, found in two stages: $$\text{(1) solve } L\mathbf{y}=\mathbf{b}\ \text{ by forward substitution;}\qquad \text{(2) solve } U\mathbf{x}=\mathbf{y}\ \text{ by back substitution.}$$ Indeed $A\mathbf{x}=L(U\mathbf{x})=L\mathbf{y}=\mathbf{b}$.

✏️ Example 2 — Two triangular solves

Use the factorisation of Example 1 to solve $A\mathbf{x}=\mathbf{b}$ with $\mathbf{b}=\begin{bmatrix}-1\\-7\\8\end{bmatrix}$.

💡 Show solution▶

Forward $L\mathbf{y}=\mathbf{b}$: $$\begin{aligned}y_1&=-1\\ 2y_1+y_2&=-7\ \Rightarrow\ y_2=-7+2=-5\\ -3y_1-2y_2+y_3&=8\ \Rightarrow\ y_3=8-3-10=-5\end{aligned}\qquad \mathbf{y}=\begin{bmatrix}-1\\-5\\-5\end{bmatrix}.$$

Back $U\mathbf{x}=\mathbf{y}$: $$\begin{aligned}5x_3&=-5\ \Rightarrow\ x_3=-1\\ -2x_2+x_3&=-5\ \Rightarrow\ -2x_2=-4\ \Rightarrow\ x_2=2\\ 3x_1-x_2+2x_3&=-1\ \Rightarrow\ 3x_1=-1+2+2=3\ \Rightarrow\ x_1=1\end{aligned}\qquad \mathbf{x}=\begin{bmatrix}1\\2\\-1\end{bmatrix}.$$

Check $A\mathbf{x}=\begin{bmatrix}3-2-2\\6-8-5\\-9+14+3\end{bmatrix}=\begin{bmatrix}-1\\-7\\8\end{bmatrix}=\mathbf{b}$. ✓

✏️ Example 3 — A spring chain under several loads

Three masses sit in a horizontal chain on a smooth track between two fixed walls, joined by four identical springs of stiffness $k=1$ kN/cm. If forces $f_1,f_2,f_3$ (kN, positive to the right) act on the masses, the displacements $u_1,u_2,u_3$ (cm) satisfy $K\mathbf{u}=\mathbf{f}$ with $$K=\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{bmatrix}.$$ Factor $K=LU$ once, then find the displacements for the load cases (a) $\mathbf{f}=(1,0,0)$, (b) $\mathbf{f}=(0,1,0)$, (c) $\mathbf{f}=(1,1,1)$.

💡 Show solution▶

Factor $m_{21}=\tfrac{-1}{2}$: $R_2\to R_2+\tfrac12R_1$ gives row $\left(0,\tfrac32,-1\right)$; $m_{31}=0$ (nothing to do). Then $m_{32}=\dfrac{-1}{3/2}=-\tfrac23$: $R_3\to R_3+\tfrac23R_2$ gives row $\left(0,0,2-\tfrac23\right)=\left(0,0,\tfrac43\right)$. So $$L=\begin{bmatrix}1&0&0\\-\tfrac12&1&0\\0&-\tfrac23&1\end{bmatrix},\qquad U=\begin{bmatrix}2&-1&0\\0&\tfrac32&-1\\0&0&\tfrac43\end{bmatrix}.$$ Notice that $L$ and $U$ inherit the zeros of the "tridiagonal" $K$ — a big saving for long chains.

(a) Forward: $y_1=1$, $y_2=0+\tfrac12(1)=\tfrac12$, $y_3=0+\tfrac23\cdot\tfrac12=\tfrac13$. Back: $\tfrac43u_3=\tfrac13\Rightarrow u_3=\tfrac14$; $\tfrac32u_2-\tfrac14=\tfrac12\Rightarrow u_2=\tfrac12$; $2u_1-\tfrac12=1\Rightarrow u_1=\tfrac34$. So $\mathbf{u}=\left(\tfrac34,\tfrac12,\tfrac14\right)$ cm.

(b) Forward: $y_1=0$, $y_2=1$, $y_3=0+\tfrac23(1)=\tfrac23$. Back: $u_3=\tfrac23\div\tfrac43=\tfrac12$; $\tfrac32u_2-\tfrac12=1\Rightarrow u_2=1$; $2u_1-1=0\Rightarrow u_1=\tfrac12$. So $\mathbf{u}=\left(\tfrac12,1,\tfrac12\right)$ cm.

(c) Forward: $y_1=1$, $y_2=1+\tfrac12=\tfrac32$, $y_3=1+\tfrac23\cdot\tfrac32=2$. Back: $u_3=2\div\tfrac43=\tfrac32$; $\tfrac32u_2-\tfrac32=\tfrac32\Rightarrow u_2=2$; $2u_1-2=1\Rightarrow u_1=\tfrac32$. So $\mathbf{u}=\left(\tfrac32,2,\tfrac32\right)$ cm.

Check $K\left(\tfrac32,2,\tfrac32\right)=\left(3-2,\,-\tfrac32+4-\tfrac32,\,-2+3\right)=(1,1,1)$. ✓ Also, case (c) is not just (a)+(b): it also needs the response to $(0,0,1)$, which by symmetry is $\left(\tfrac14,\tfrac12,\tfrac34\right)$, and $\left(\tfrac34,\tfrac12,\tfrac14\right)+\left(\tfrac12,1,\tfrac12\right)+\left(\tfrac14,\tfrac12,\tfrac34\right)=\left(\tfrac32,2,\tfrac32\right)$ — linearity at work. The factorisation was done once; each load cost only two short substitutions.

The spring chain of Example 3 with $K=LU$ already factored. For the chosen load $\mathbf{f}$, watch forward substitution fill $\mathbf{y}$ from the top ($L\mathbf{y}=\mathbf{f}$), then back substitution fill $\mathbf{u}$ from the bottom ($U\mathbf{u}=\mathbf{y}$). As each displacement is found, the mass in the picture moves to its new position. While playing, the widget moves on to the next load case — the factorisation is never redone.

Load case $\mathbf{f}$ (kN):

stage—

equation used—

$\mathbf{u}$ (cm)—

arithmetic operations this load—

When Does LU Exist? Row Swaps and $PA=LU$

Elimination without row swaps breaks down as soon as a pivot is $0$. The simplest example is $$A=\begin{bmatrix}0&1\\1&1\end{bmatrix}.$$ Suppose $$A=\begin{bmatrix}1&0\\ \ell&1\end{bmatrix}\begin{bmatrix}u_{11}&u_{12}\\0&u_{22}\end{bmatrix}=\begin{bmatrix}u_{11}&u_{12}\\ \ell u_{11}&\ell u_{12}+u_{22}\end{bmatrix}.$$ Comparing $(1,1)$ entries, $u_{11}=0$, so the $(2,1)$ entry is $\ell\cdot0=0\ne1$. This invertible matrix has no LU factorisation. Swapping the two rows first fixes everything — and swaps are recorded by a permutation matrix.

📘 Definition 2.11 — Permutation matrix

A permutation matrix is a square matrix obtained by rearranging the rows of the identity matrix $I$. Multiplying $A$ on the left by a permutation matrix $P$ rearranges the rows of $A$ in the same way. For instance $P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}$ swaps rows 1 and 2. Every permutation matrix is invertible with $P^{-1}=P^{T}$.

📐 Theorem 2.14 — Existence of LU and of $PA=LU$

(a) Let $A$ be an invertible $n\times n$ matrix. Then $A$ has an LU factorisation (with $L$ unit lower triangular) if and only if elimination without row interchanges never meets a zero pivot; equivalently, if and only if every leading submatrix $A_k$ (the top-left $k\times k$ block, $k=1,\dots,n$) is invertible. When it exists, this factorisation is unique.

(b) For every $n\times n$ matrix $A$ there is a permutation matrix $P$ such that $PA$ has an LU factorisation: $$PA=LU.$$ Then $A\mathbf{x}=\mathbf{b}$ is equivalent to $LU\mathbf{x}=P\mathbf{b}$, solved by forward and back substitution as before.

✏️ Example 4 — A zero pivot: finding $PA=LU$

Let $A=\begin{bmatrix}0&2&1\\1&1&1\\2&1&3\end{bmatrix}$. Find a permutation matrix $P$ and a factorisation $PA=LU$, then use it to solve $A\mathbf{x}=\begin{bmatrix}0\\2\\7\end{bmatrix}$.

💡 Show solution▶

Swap The $(1,1)$ entry is $0$, so swap rows 1 and 2 before eliminating: $$P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix},\qquad PA=\begin{bmatrix}1&1&1\\0&2&1\\2&1&3\end{bmatrix}.$$

Eliminate $m_{21}=0$, $m_{31}=2$: $R_3\to R_3-2R_1$ gives $(0,-1,1)$. Then $m_{32}=\tfrac{-1}{2}$: $R_3\to R_3+\tfrac12R_2$ gives $\left(0,0,\tfrac32\right)$. No further zero pivots, so $$L=\begin{bmatrix}1&0&0\\0&1&0\\2&-\tfrac12&1\end{bmatrix},\qquad U=\begin{bmatrix}1&1&1\\0&2&1\\0&0&\tfrac32\end{bmatrix},\qquad PA=LU.$$

Solve $P\mathbf{b}=(2,0,7)$. Forward: $y_1=2$, $y_2=0$, $y_3=7-2(2)+\tfrac12(0)=3$. Back: $\tfrac32x_3=3\Rightarrow x_3=2$; $2x_2+2=0\Rightarrow x_2=-1$; $x_1-1+2=2\Rightarrow x_1=1$. So $\mathbf{x}=(1,-1,2)$.

Check $A\mathbf{x}=(0-2+2,\ 1-1+2,\ 2-1+6)=(0,2,7)$. ✓

⚠️ Analysis — what software really does

1. Singular matrices can have an LU factorisation. $\begin{bmatrix}1&2\\2&4\end{bmatrix}=\begin{bmatrix}1&0\\2&1\end{bmatrix}\begin{bmatrix}1&2\\0&0\end{bmatrix}$ — the zero simply shows up on the diagonal of $U$. Theorem 2.13 then does not apply: $U\mathbf{x}=\mathbf{y}$ has no solution or infinitely many.

2. Tiny pivots are dangerous too. On a computer a pivot like $10^{-12}$ produces huge multipliers that magnify rounding errors. Library routines therefore use partial pivoting: in each column they swap up the entry of largest absolute value before eliminating. That is why a numerical package returns three matrices $P$, $L$, $U$ even when no swap was strictly necessary.

Counting the Work

Why go to the trouble of storing $L$? Count the arithmetic. A flop (floating-point operation) is one addition, subtraction, multiplication or division.

📐 Theorem 2.15 — Operation counts for an $n\times n$ system

Assume no row interchanges are needed.

Computing $A=LU$ by elimination takes $\tfrac23n^3-\tfrac12n^2-\tfrac16n$ flops — about $\tfrac23n^3$ for large $n$.
Forward substitution with unit lower triangular $L$ takes $n^2-n$ flops; back substitution with $U$ takes $n^2$ flops. Together: about $2n^2$.
Computing $A^{-1}$ takes about $2n^3$ flops, and then $A^{-1}\mathbf{b}$ costs another $2n^2$.

For $n=3$: factorising costs $18-4.5-0.5=13$ flops and each new right-hand side $6+9=15$ flops. For small $n$ the difference is minor — but $n^3$ grows much faster than $n^2$.

✏️ Example 5 — Fifty load cases on a large model

A finite-element model has $n=1000$ unknowns and must be solved for $50$ different load vectors. Estimate the number of flops if (a) elimination is repeated from scratch for each load, (b) $A=LU$ is computed once and reused.

💡 Show solution▶

(a) Each fresh elimination costs about $\tfrac23(1000)^3\approx6.67\times10^{8}$ flops (the extra $2n^2$ for the right-hand side is negligible). Fifty of them: $50\times6.67\times10^{8}\approx3.33\times10^{10}$ flops.

(b) One factorisation, $6.67\times10^{8}$ flops, plus $50$ pairs of triangular solves at $2(1000)^2=2\times10^{6}$ flops each, i.e. $10^{8}$: total $\approx7.67\times10^{8}$ flops.

Compare $\dfrac{3.33\times10^{10}}{7.67\times10^{8}}\approx43$. Reusing the factorisation is about $43$ times faster; on a machine doing $10^9$ flops per second that is $33$ seconds against under one second.

⚠️ Don't compute $A^{-1}$ just to solve a system

It is tempting to compute $A^{-1}$ once and then form $\mathbf{x}=A^{-1}\mathbf{b}$ for each load. But finding $A^{-1}$ costs about three times as much as finding $L$ and $U$, and multiplying by $A^{-1}$ costs the same $2n^2$ flops as the two triangular solves — so it is slower, and it usually adds more rounding error. For large sparse matrices (like the tridiagonal $K$ of Example 3) it is far worse: $L$ and $U$ keep the zeros, while $A^{-1}$ is typically completely full.

Exercises

Find an LU factorisation (no row interchanges) of each matrix. Check by multiplying $LU$.

1 $\begin{bmatrix}2&5\\6&13\end{bmatrix}$
2 $\begin{bmatrix}4&-2\\-8&7\end{bmatrix}$
3 $\begin{bmatrix}1&2&3\\2&7&9\\-1&4&1\end{bmatrix}$
4 $\begin{bmatrix}2&-1&3\\-4&5&-4\\6&0&13\end{bmatrix}$
5 $\begin{bmatrix}3&1&0\\9&5&2\\0&4&7\end{bmatrix}$
6 $\begin{bmatrix}1&1&0&0\\2&3&1&0\\0&2&4&1\\0&0&3&5\end{bmatrix}$
7 The $2\times3$ matrix $\begin{bmatrix}2&1&4\\6&5&9\end{bmatrix}$ ($L$ is $2\times2$, $U$ is $2\times3$).
8 The general $2\times2$ matrix $\begin{bmatrix}a&b\\c&d\end{bmatrix}$ with $a\ne0$.

Check: multiply your factors back together

Part:

💡 Show solutions▶

1 $m_{21}=\tfrac62=3$; $R_2\to R_2-3R_1$ gives $(0,13-15)=(0,-2)$. $$L=\begin{bmatrix}1&0\\3&1\end{bmatrix},\quad U=\begin{bmatrix}2&5\\0&-2\end{bmatrix}.$$

2 $m_{21}=\tfrac{-8}{4}=-2$; $R_2\to R_2+2R_1$ gives $(0,7-4)=(0,3)$. $$L=\begin{bmatrix}1&0\\-2&1\end{bmatrix},\quad U=\begin{bmatrix}4&-2\\0&3\end{bmatrix}.$$

3 $m_{21}=2$, $m_{31}=-1$: rows become $(0,3,3)$ and $(0,6,4)$. Then $m_{32}=\tfrac63=2$: $R_3\to R_3-2R_2$ gives $(0,0,-2)$. $$L=\begin{bmatrix}1&0&0\\2&1&0\\-1&2&1\end{bmatrix},\quad U=\begin{bmatrix}1&2&3\\0&3&3\\0&0&-2\end{bmatrix}.$$

4 $m_{21}=-2$, $m_{31}=3$: rows become $(0,3,2)$ and $(0,3,4)$. Then $m_{32}=1$: $R_3\to R_3-R_2$ gives $(0,0,2)$. $$L=\begin{bmatrix}1&0&0\\-2&1&0\\3&1&1\end{bmatrix},\quad U=\begin{bmatrix}2&-1&3\\0&3&2\\0&0&2\end{bmatrix}.$$

5 $m_{21}=3$: row 2 becomes $(0,2,2)$; $m_{31}=0$. Then $m_{32}=\tfrac42=2$: $R_3\to R_3-2R_2$ gives $(0,0,3)$. $$L=\begin{bmatrix}1&0&0\\3&1&0\\0&2&1\end{bmatrix},\quad U=\begin{bmatrix}3&1&0\\0&2&2\\0&0&3\end{bmatrix}.$$

6 Only one entry below each pivot is nonzero. $m_{21}=2$: row 2 becomes $(0,1,1,0)$. $m_{32}=2$: row 3 becomes $(0,0,2,1)$. $m_{43}=\tfrac32$: row 4 becomes $\left(0,0,0,5-\tfrac32\right)=\left(0,0,0,\tfrac72\right)$. $$L=\begin{bmatrix}1&0&0&0\\2&1&0&0\\0&2&1&0\\0&0&\tfrac32&1\end{bmatrix},\quad U=\begin{bmatrix}1&1&0&0\\0&1&1&0\\0&0&2&1\\0&0&0&\tfrac72\end{bmatrix}.$$

7 $m_{21}=3$; $R_2\to R_2-3R_1$ gives $(0,2,-3)$. $$L=\begin{bmatrix}1&0\\3&1\end{bmatrix},\quad U=\begin{bmatrix}2&1&4\\0&2&-3\end{bmatrix}.$$

8 $m_{21}=\tfrac ca$; $R_2\to R_2-\tfrac caR_1$ gives $\left(0,d-\tfrac{bc}{a}\right)$. $$\begin{bmatrix}a&b\\c&d\end{bmatrix}=\begin{bmatrix}1&0\\ \tfrac ca&1\end{bmatrix}\begin{bmatrix}a&b\\0&d-\tfrac{bc}{a}\end{bmatrix}.$$ (The second pivot is $\tfrac{ad-bc}{a}$; it is $0$ exactly when the matrix is not invertible.)

Solving systems with an LU factorisation. Show the forward substitution ($L\mathbf{y}=\mathbf{b}$) and the back substitution ($U\mathbf{x}=\mathbf{y}$).

1 $\begin{bmatrix}2&5\\6&13\end{bmatrix}\mathbf{x}=\begin{bmatrix}1\\5\end{bmatrix}$ (use Ex. I.1).
2 $\begin{bmatrix}1&2&3\\2&7&9\\-1&4&1\end{bmatrix}\mathbf{x}=\begin{bmatrix}3\\6\\-5\end{bmatrix}$ (use Ex. I.3).
3 Solve $A\mathbf{x}=\mathbf{b}_1$ and $A\mathbf{x}=\mathbf{b}_2$ for the matrix of Ex. I.4, with $\mathbf{b}_1=\begin{bmatrix}-1\\0\\-7\end{bmatrix}$ and $\mathbf{b}_2=\begin{bmatrix}2\\1\\13\end{bmatrix}$.
4 $A=LU$ with $L=\begin{bmatrix}1&0&0\\-1&1&0\\2&-3&1\end{bmatrix}$ and $U=\begin{bmatrix}1&-2&1\\0&2&4\\0&0&-1\end{bmatrix}$. Solve $A\mathbf{x}=\begin{bmatrix}3\\-1\\4\end{bmatrix}$ without forming $A$, then form $A$ and check.
5 $\begin{bmatrix}3&1&0\\9&5&2\\0&4&7\end{bmatrix}\mathbf{x}=\begin{bmatrix}4\\16\\11\end{bmatrix}$ (use Ex. I.5).
6 Find the first column of $A^{-1}$ for the matrix of Ex. I.3 by solving $A\mathbf{x}=\mathbf{e}_1$.
7 Find $A^{-1}$ for $A=\begin{bmatrix}2&5\\6&13\end{bmatrix}$ by solving $A\mathbf{x}=\mathbf{e}_1$ and $A\mathbf{x}=\mathbf{e}_2$ with the factors of Ex. I.1.
8 Electrical circuit. The loop currents $i_1,i_2,i_3$ (amperes) in a three-loop circuit satisfy $R\mathbf{i}=\mathbf{v}$ with resistance matrix $R=\begin{bmatrix}4&-2&0\\-2&5&-1\\0&-1&3\end{bmatrix}$ (ohms) and source voltages $\mathbf{v}$ (volts). Factor $R=LU$ and find the currents for the two settings $\mathbf{v}=(8,3,1)$ and $\mathbf{v}=(2,1,5)$.

💡 Show solutions▶

1 Forward: $y_1=1$, $y_2=5-3(1)=2$. Back: $-2x_2=2\Rightarrow x_2=-1$; $2x_1+5(-1)=1\Rightarrow x_1=3$. $\mathbf{x}=(3,-1)$.

2 Forward: $y_1=3$, $y_2=6-2(3)=0$, $y_3=-5+1(3)-2(0)=-2$. Back: $-2x_3=-2\Rightarrow x_3=1$; $3x_2+3=0\Rightarrow x_2=-1$; $x_1-2+3=3\Rightarrow x_1=2$. $\mathbf{x}=(2,-1,1)$.

3 $\mathbf{b}_1$: $\mathbf{y}=(-1,\,0+2(-1),\,-7-3(-1)-(-2))=(-1,-2,-2)$; back: $x_3=-1$, $3x_2-2=-2\Rightarrow x_2=0$, $2x_1-0-3=-1\Rightarrow x_1=1$. $\mathbf{x}=(1,0,-1)$. $\mathbf{b}_2$: $\mathbf{y}=(2,\,1+4,\,13-6-5)=(2,5,2)$; back: $x_3=1$, $3x_2+2=5\Rightarrow x_2=1$, $2x_1-1+3=2\Rightarrow x_1=0$. $\mathbf{x}=(0,1,1)$. The same $L$ and $U$ served both.

4 Forward: $y_1=3$, $y_2=-1+3=2$, $y_3=4-2(3)+3(2)=4$. Back: $-x_3=4\Rightarrow x_3=-4$; $2x_2+4(-4)=2\Rightarrow x_2=9$; $x_1-18-4=3\Rightarrow x_1=25$. So $\mathbf{x}=(25,9,-4)$. Multiplying, $A=LU=\begin{bmatrix}1&-2&1\\-1&4&3\\2&-10&-11\end{bmatrix}$, and $A\mathbf{x}=(25-18-4,\,-25+36-12,\,50-90+44)=(3,-1,4)$. ✓

5 Forward: $y_1=4$, $y_2=16-12=4$, $y_3=11-0-2(4)=3$. Back: $3x_3=3\Rightarrow x_3=1$; $2x_2+2=4\Rightarrow x_2=1$; $3x_1+1=4\Rightarrow x_1=1$. $\mathbf{x}=(1,1,1)$.

6 Forward with $\mathbf{e}_1=(1,0,0)$: $y_1=1$, $y_2=-2$, $y_3=0+1-2(-2)=5$. Back: $-2x_3=5\Rightarrow x_3=-\tfrac52$; $3x_2+3\left(-\tfrac52\right)=-2\Rightarrow x_2=\tfrac{11}{6}$; $x_1+\tfrac{11}{3}-\tfrac{15}{2}=1\Rightarrow x_1=\tfrac{29}{6}$. The first column of $A^{-1}$ is $\left(\tfrac{29}{6},\tfrac{11}{6},-\tfrac52\right)$.

7 $\mathbf{e}_1$: $\mathbf{y}=(1,-3)$, $x_2=\tfrac32$, $2x_1+\tfrac{15}{2}=1\Rightarrow x_1=-\tfrac{13}{4}$. $\mathbf{e}_2$: $\mathbf{y}=(0,1)$, $x_2=-\tfrac12$, $2x_1-\tfrac52=0\Rightarrow x_1=\tfrac54$. $$A^{-1}=\begin{bmatrix}-\tfrac{13}{4}&\tfrac54\\ \tfrac32&-\tfrac12\end{bmatrix}.$$ Check: $\begin{bmatrix}2&5\\6&13\end{bmatrix}A^{-1}=\begin{bmatrix}-\tfrac{13}{2}+\tfrac{15}{2}&\tfrac52-\tfrac52\\-\tfrac{39}{2}+\tfrac{39}{2}&\tfrac{15}{2}-\tfrac{13}{2}\end{bmatrix}=I$. ✓

8 $m_{21}=-\tfrac12$: row 2 becomes $(0,4,-1)$; $m_{31}=0$; $m_{32}=-\tfrac14$: row 3 becomes $\left(0,0,3-\tfrac14\right)=\left(0,0,\tfrac{11}{4}\right)$. So $L=\begin{bmatrix}1&0&0\\-\tfrac12&1&0\\0&-\tfrac14&1\end{bmatrix}$, $U=\begin{bmatrix}4&-2&0\\0&4&-1\\0&0&\tfrac{11}{4}\end{bmatrix}$. For $\mathbf{v}=(8,3,1)$: $\mathbf{y}=\left(8,\,7,\,\tfrac{11}{4}\right)$, so $i_3=1$, $4i_2-1=7\Rightarrow i_2=2$, $4i_1-4=8\Rightarrow i_1=3$: currents $(3,2,1)$ A. For $\mathbf{v}=(2,1,5)$: $\mathbf{y}=\left(2,\,2,\,\tfrac{11}{2}\right)$, so $i_3=2$, $4i_2-2=2\Rightarrow i_2=1$, $4i_1-2=2\Rightarrow i_1=1$: currents $(1,1,2)$ A.

III

Existence, row swaps and permutation matrices.

1 Show that $\begin{bmatrix}0&1\\1&1\end{bmatrix}$ has no LU factorisation, and find a permutation matrix $P$ with $PA=LU$.
2 Does the singular matrix $\begin{bmatrix}1&2\\2&4\end{bmatrix}$ have an LU factorisation? If so, find it.
3 Find $P$, $L$ and $U$ with $PA=LU$ for $A=\begin{bmatrix}0&3&1\\2&4&6\\1&5&2\end{bmatrix}$.
4 For which values of $k$ does $A=\begin{bmatrix}1&2&0\\2&k&1\\0&1&3\end{bmatrix}$ have an LU factorisation (no row interchanges)?
5 Write down the $3\times3$ permutation matrix $P$ that swaps rows 1 and 3. Compute $P^2$ and explain why $P^{-1}=P=P^{T}$.
6 Suppose $A=LU$ is invertible. Explain why no diagonal entry of $U$ can be $0$.

💡 Show solutions▶

1 As shown before Definition 2.11: $LU$ has $(1,1)$ entry $u_{11}$, so $u_{11}=0$; then its $(2,1)$ entry is $\ell u_{11}=0\ne1$. Contradiction. With $P=\begin{bmatrix}0&1\\1&0\end{bmatrix}$, $PA=\begin{bmatrix}1&1\\0&1\end{bmatrix}$, which is already upper triangular: $L=I$, $U=PA$.

2 Yes. $m_{21}=2$ and $R_2\to R_2-2R_1$ gives $(0,0)$: $\begin{bmatrix}1&2\\2&4\end{bmatrix}=\begin{bmatrix}1&0\\2&1\end{bmatrix}\begin{bmatrix}1&2\\0&0\end{bmatrix}$. A zero pivot in the last column causes no trouble for the factorisation; it just shows the matrix is singular.

3 Swap rows 1 and 2: $P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}$, $PA=\begin{bmatrix}2&4&6\\0&3&1\\1&5&2\end{bmatrix}$. Then $m_{21}=0$, $m_{31}=\tfrac12$ (row 3 becomes $(0,3,-1)$), $m_{32}=1$ (row 3 becomes $(0,0,-2)$): $$L=\begin{bmatrix}1&0&0\\0&1&0\\ \tfrac12&1&1\end{bmatrix},\quad U=\begin{bmatrix}2&4&6\\0&3&1\\0&0&-2\end{bmatrix}.$$

4 $R_2\to R_2-2R_1$ gives $(0,k-4,1)$. If $k\ne4$ the second pivot is nonzero, $m_{32}=\tfrac{1}{k-4}$, and elimination finishes, so an LU factorisation exists. If $k=4$, the second pivot is $0$ while the entry below it is $1$: comparing entries of $A=LU$ gives $u_{11}=1$, $u_{12}=2$, $\ell_{21}=2$, $\ell_{31}=0$, $u_{22}=k-4=0$, and then the $(3,2)$ entry requires $\ell_{31}u_{12}+\ell_{32}u_{22}=0=1$ — impossible. Answer: all $k\ne4$. (When $k=\tfrac{13}{3}$ the last pivot $3-\tfrac{1}{k-4}$ is $0$: $A$ is singular but still has an LU factorisation.)

5 $P=\begin{bmatrix}0&0&1\\0&1&0\\1&0&0\end{bmatrix}$. Swapping rows 1 and 3 twice restores the original order, so $P^2=I$, hence $P^{-1}=P$. $P$ is symmetric, so $P^{T}=P$ as well.

6 $L$ is invertible (forward substitution always works because its diagonal entries are $1$), so $U=L^{-1}A$ is a product of invertible matrices and is invertible. An upper triangular matrix with a $0$ on the diagonal has fewer than $n$ pivots, so it is not invertible. Hence every $u_{ii}\ne0$.

Counting the work (use the approximate counts $\tfrac23n^3$ for LU and $2n^2$ for a pair of triangular solves).

1 For $n=500$, estimate the flops for computing $A=LU$ and for one pair of triangular solves. How many solves cost as much as one factorisation?
2 A system with $n=2000$ must be solved for $200$ right-hand sides. Compare the cost of repeating elimination each time with the cost of factorising once.
3 A computer performs $10^{9}$ flops per second. How long does the LU factorisation of a $10\,000\times10\,000$ matrix take, and how long does each further right-hand side take?
4 If $n$ is doubled, by what factor does the cost of the factorisation grow? And the cost of a pair of triangular solves?
5 A classmate suggests computing $A^{-1}$ once and then forming $A^{-1}\mathbf{b}$ for every load. Compare the costs and give two reasons to prefer $LU$.
6 Long chains. For a tridiagonal matrix (nonzero entries only on the diagonal and directly above and below it), such as the spring matrix of Example 3 with $n$ masses, explain why $L$ and $U$ are also mostly zeros and why the whole solve then costs a number of flops proportional to $n$ rather than $n^3$.

💡 Show solutions▶

1 LU: $\tfrac23(500)^3\approx8.33\times10^{7}$ flops. Solves: $2(500)^2=5\times10^{5}$ flops. Ratio $\approx167$: one factorisation costs as much as about $167$ pairs of triangular solves (in general the ratio is $\tfrac{n}{3}$).

2 Repeating: $200\times\tfrac23(2000)^3\approx200\times5.33\times10^{9}\approx1.07\times10^{12}$ flops. Once: $5.33\times10^{9}+200\times2(2000)^2=5.33\times10^{9}+1.6\times10^{9}\approx6.93\times10^{9}$ flops — about $154$ times less work.

3 LU: $\tfrac23\times10^{12}$ flops $\div10^{9}$ flops/s $\approx667$ s, about $11$ minutes. Each new right-hand side: $2\times10^{8}$ flops, i.e. $0.2$ s.

4 $\tfrac23(2n)^3=8\cdot\tfrac23n^3$: the factorisation becomes $8$ times as expensive. $2(2n)^2=4\cdot2n^2$: the triangular solves become $4$ times as expensive.

5 Computing $A^{-1}$ costs about $2n^3$ flops — three times the $\tfrac23n^3$ of $LU$ — and each product $A^{-1}\mathbf{b}$ costs $2n^2$, exactly the same as forward plus back substitution. So the inverse never saves work. In addition, (i) it typically accumulates more rounding error, and (ii) for sparse matrices $L$ and $U$ keep the zeros while $A^{-1}$ is usually full, wasting memory and time.

6 Below each pivot there is only one nonzero entry, so each column needs one multiplier, and the row operation changes only the one row beneath, in just the entry to the right of the eliminated one (about $3$ flops per column). No new nonzeros are created: $L$ has only its diagonal and the entries just below it; $U$ has only its diagonal and the entries just above it. Forward and back substitution then use about $2$ and $3$ flops per row. Everything is a fixed number of flops per row — proportional to $n$. A chain with a million masses is solved in a few million flops instead of about $10^{18}$.

🔬 Research–Teaching NexusResearch-tutored · optional, about 30 minutes

Research question. When engineers model a long structure — a rail track on sleepers, a pipeline, a heated metal rod, a chain of floors in a tall building — the matrix is tridiagonal. What pattern do the pivots and multipliers of its LU factorisation follow, and how fast can such a system really be solved?

Why it matters. Finite-difference models of heat flow, beam bending and vibration in civil, mechanical and biomedical engineering all lead to the matrix $K_n$ with $2$ on the diagonal and $-1$ beside it. Solvers for these "banded" systems run inside weather models, bridge-monitoring software and heart-tissue simulations, often millions of times.

What to do.

Factor $K_4$ and $K_5$ by hand, extending Example 3. List the pivots $u_{11},u_{22},\dots$ and the multipliers. Guess the general pattern for $K_n$ and test it on $K_6$ with Interactive 1's method.
Use your pattern to predict the product of the pivots of $K_n$ (this number is $\det K_n$, Chapter 3).
In a spreadsheet, solve $K_{10}\mathbf{u}=\mathbf{f}$ for the uniform load $\mathbf{f}=(1,1,\dots,1)$ by forward and back substitution (one column for $\mathbf{y}$, one for $\mathbf{u}$). Plot $u_i$ against $i$ — what shape do you see?
Optional: in Python (scipy.linalg.lu_factor / lu_solve) or Octave, time solving with $n=1000,2000,4000$ for one and for $100$ right-hand sides; compare with the predictions $\tfrac23n^3$ and $2n^2$.

Tools: Pen and paper, spreadsheet; optionally Python or Octave

Hand in: A one-page note: the pivot and multiplier patterns with a short justification, the displacement plot, and one sentence on why banded solvers make very large engineering models possible.

Going further. Look up the "Thomas algorithm" used for tridiagonal systems in computational fluid dynamics, and compare its operation count with your answer to Exercise IV.6.

Instructor note: the pivots of $K_n$ are $2,\tfrac32,\tfrac43,\dots,\tfrac{n+1}{n}$ (pivot $k$ is $\tfrac{k+1}{k}$), the multipliers are $-\tfrac12,-\tfrac23,\dots,-\tfrac{n-1}{n}$, and the product of the pivots telescopes to $n+1$. For the uniform load the displacements are $u_i=\tfrac{i(n+1-i)}{2}$ — a parabola (for $n=10$: $5,9,12,14,15,15,14,12,9,5$).

🧪 Self-check quiz

0/5 correct

Q1. Elimination on a $3\times3$ matrix uses, in order, $R_2\to R_2-4R_1$, $R_3\to R_3+2R_1$, $R_3\to R_3-5R_2$ and no other operations. Then $L=$

Q2. Given $A=LU$, the first step in solving $A\mathbf{x}=\mathbf{b}$ is to solve

Q3. Which matrix has no LU factorisation without row interchanges?

Q4. For large $n$ the LU factorisation costs about $\tfrac23n^3$ flops. If the number of unknowns grows from $n$ to $10n$, the factorisation cost is multiplied by about

Q5. If $PA=LU$ with a permutation matrix $P$, then $A\mathbf{x}=\mathbf{b}$ is equivalent to

← 2.2 Inverses and Elementary Matrices 2.4 Markov Chains and Input–Output Models →