LU-Factorization & Input–Output Models
Splitting a matrix into triangular pieces for fast solving — and a Nobel-winning application to economics
§1Why Factor a Matrix at All?
Imagine you run a factory. Every week the same machine $A$ processes a different batch of orders $\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3, \dots$ — same $A$, new $\mathbf{b}$ each time. Solving $A\mathbf{x} = \mathbf{b}$ from scratch every week, by full Gaussian elimination, repeats the same expensive work over and over. There has to be a smarter way.
§2Triangular Matrices — the Stars of the Show
Everything rests on triangular matrices, so let us define them carefully and see exactly which entries matter. First, the diagonal.
The main diagonal of a matrix runs from the top-left corner downward to the right: the entries $a_{11}, a_{22}, a_{33}, \dots$ where the row index equals the column index.
A matrix is lower triangular if every entry above the main diagonal is zero. All the nonzero action sits on or below the diagonal — the lower-left wedge.
A matrix is upper triangular if every entry below the main diagonal is zero. The nonzero entries fill the upper-right wedge, on or above the diagonal.
Square examples — the easy case
$$\begin{bmatrix} 2 & 0 & 0 \\ 5 & 1 & 0 \\ -3 & 4 & 7 \end{bmatrix}$$
Zeros fill the whole region above the diagonal.
$$\begin{bmatrix} 1 & 2 & -1 \\ 0 & 3 & 5 \\ 0 & 0 & 4 \end{bmatrix}$$
Zeros fill the whole region below the diagonal.
Non-square examples — same rule
Triangular-ness is about the diagonal, and the definition works for rectangular matrices too: "lower" means zeros above the diagonal, "upper" means zeros below it, wherever that diagonal runs.
$$\begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 1 & 4 & 0 \\ 0 & 0 & 2 & 5 \end{bmatrix}$$
Every entry below the diagonal (positions $(2,1),(3,1),(3,2)$) is zero. This is exactly the shape of a row-echelon matrix.
$$\begin{bmatrix} 2 & 0 & 0 \\ 1 & 3 & 0 \\ 4 & -1 & 5 \\ 0 & 2 & 6 \end{bmatrix}$$
Every entry above the diagonal is zero. Square matrices are easiest, but the idea extends cleanly.
§3Upper Triangular ⟹ Back Substitution
When the coefficient matrix is upper triangular, the system practically solves itself. The bottom equation has only one unknown; solve it, then climb upward, substituting as you go.
Solve $U\mathbf{x} = \mathbf{b}$ where
$$\begin{bmatrix} 1 & 2 & -1 \\ 0 & 1 & 3 \\ 0 & 0 & 2 \end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 3 \\ 4 \\ 6 \end{bmatrix}.$$
The system, written out, is $\begin{cases} x_1 + 2x_2 - x_3 = 3 \\ x_2 + 3x_3 = 4 \\ 2x_3 = 6 \end{cases}$. Start at the bottom and work up:
Row 3: $2x_3 = 6 \Rightarrow x_3 = 3$.
Row 2: $x_2 + 3(3) = 4 \Rightarrow x_2 = 4 - 9 = -5$.
Row 1: $x_1 + 2(-5) - 3 = 3 \Rightarrow x_1 = 3 + 10 + 3 = 16$.
So $\mathbf{x} = (16, -5, 3)$. No elimination needed — just substitution from the bottom up. This is back substitution.
§4Lower Triangular ⟹ Forward Substitution
A lower triangular system is just as easy, but the other way round. The top equation has only one unknown; solve it, then descend.
Solve $L\mathbf{y} = \mathbf{b}$ where
$$\begin{bmatrix} 2 & 0 & 0 \\ 1 & 3 & 0 \\ -1 & 2 & 1 \end{bmatrix}\begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 4 \\ 7 \\ 1 \end{bmatrix}.$$
The system is $\begin{cases} 2y_1 = 4 \\ y_1 + 3y_2 = 7 \\ -y_1 + 2y_2 + y_3 = 1 \end{cases}$. Start at the top and work down:
Row 1: $2y_1 = 4 \Rightarrow y_1 = 2$.
Row 2: $2 + 3y_2 = 7 \Rightarrow y_2 = \tfrac{5}{3}$.
Row 3: $-2 + 2(\tfrac{5}{3}) + y_3 = 1 \Rightarrow y_3 = 1 + 2 - \tfrac{10}{3} = -\tfrac{1}{3}$.
So $\mathbf{y} = (2, \tfrac{5}{3}, -\tfrac{1}{3})$. Top-down substitution — this is forward substitution.
§5The Two-Stage Method — Solving Ax = b via LU
Now we combine the two. Suppose $A = LU$. We want to solve $A\mathbf{x} = \mathbf{b}$, i.e. $LU\mathbf{x} = \mathbf{b}$. The trick is to give the middle piece $U\mathbf{x}$ a name.
To solve $A\mathbf{x} = \mathbf{b}$ when $A = LU$:
$\textbf{Stage 1.}$ Solve $L\mathbf{y} = \mathbf{b}$ for $\mathbf{y}$ by forward substitution ($L$ is lower triangular).
$\textbf{Stage 2.}$ Solve $U\mathbf{x} = \mathbf{y}$ for $\mathbf{x}$ by back substitution ($U$ is upper triangular).
Then $\mathbf{x}$ solves the original system.
The resulting $\mathbf{x}$ is a solution because
$$A\mathbf{x} = (LU)\mathbf{x} = L(U\mathbf{x}) = L\mathbf{y} = \mathbf{b}. \;\checkmark$$
Moreover, every solution arises this way: given any solution $\mathbf{x}$ of $A\mathbf{x}=\mathbf{b}$, set $\mathbf{y} = U\mathbf{x}$; then $L\mathbf{y} = LU\mathbf{x} = A\mathbf{x} = \mathbf{b}$, so this $\mathbf{y}$ is exactly the one Stage 1 produces. Nothing is lost. And because each stage is pure substitution, the method adapts beautifully to a computer.
Solve $A\mathbf{x} = \mathbf{b}$ where $A = LU$ with $L = \begin{bmatrix}1&0&0\\2&1&0\\-1&3&1\end{bmatrix}$, $U = \begin{bmatrix}2&1&-1\\0&1&2\\0&0&3\end{bmatrix}$, and $\mathbf{b} = \begin{bmatrix}1\\4\\6\end{bmatrix}$.
§6Two Facts About Triangular Matrices (Lemma 2.7.1)
Let $A$ and $B$ be matrices.
$\textbf{1.}$ If $A$ and $B$ are both lower (upper) triangular, so is their product $AB$.
$\textbf{2.}$ If $A$ is $n\times n$ and lower (upper) triangular, then $A$ is invertible if and only if every main-diagonal entry is nonzero. In that case $A^{-1}$ is also lower (upper) triangular.
Part 1: when you multiply two lower-triangular matrices, every entry above the diagonal of the product is a sum of terms each containing a zero factor — so it stays zero. The triangular shape is preserved.
Part 2: the diagonal of a triangular matrix is its set of pivots. All diagonal entries nonzero means a pivot in every row — full rank — hence invertible. Inverting by the $[A\mid I]$ method never disturbs the triangular shape, so $A^{-1}$ comes out triangular too.
§7Finding L and U — Where They Come From
Here is the beautiful part. We already know how to carry $A$ to a row-echelon (upper triangular) matrix $U$ using elementary matrices. That very process hands us $L$ for free.
Reduce $A$ to row-echelon form $U$ with elementary matrices:
$$E_k E_{k-1} \cdots E_2 E_1 A = U.$$
Then $A = LU$ where
$$L = (E_k \cdots E_1)^{-1} = E_1^{-1} E_2^{-1} \cdots E_k^{-1}.$$
If we use no row interchanges (and never add a row to a row above it), every $E_i$ is lower triangular — so by Lemma 2.7.1, $L$ is lower triangular and invertible. This is the LU-factorization.
If $A$ can be lower reduced to a row-echelon matrix $U$ (that is, reduced using no row interchanges), then $A = LU$ where $L$ is lower triangular and invertible and $U$ is upper triangular and row-echelon. Such a factorization $A = LU$ is called an LU-factorization of $A$.
Let $A$ be $m\times n$ of rank $r$, lower-reducible to row-echelon $U$. Then $A = LU$, where $L$ is built as follows:
$\textbf{1.}$ If $A = 0$, take $L = I_m$ and $U = 0$.
$\textbf{2.}$ If $A \neq 0$, write $A_1 = A$ and let $\mathbf{c}_1$ be its leading column. Use $\mathbf{c}_1$ to create the first leading $1$ and zeros below it (by lower reduction). Let $A_2$ be the matrix of rows $2$ to $m$ of the result.
$\textbf{3.}$ If $A_2 \neq 0$, let $\mathbf{c}_2$ be its leading column and repeat Step 2 on $A_2$ to create $A_3$.
$\textbf{4.}$ Continue until $U$ is reached (all rows below the last leading $1$ are zero). This takes $r$ steps.
$\textbf{5.}$ Build $L$ by placing $\mathbf{c}_1, \mathbf{c}_2, \dots, \mathbf{c}_r$ at the bottom of the first $r$ columns of $I_m$.
§8A Full 4×4 Walkthrough
Let us run the book's algorithm completely on a $4\times4$ matrix, watching the leading column at each stage — that column is exactly what gets stored into $L$.
$$A = \begin{bmatrix} 2 & 4 & -2 & 6 \\ 1 & 5 & 5 & 0 \\ 3 & 5 & -4 & 12 \\ -1 & 0 & 9 & 5 \end{bmatrix}.$$
The first column $\mathbf{c}_1 = (2,1,3,-1)$ is the leading column. Store it — it becomes column 1 of $L$. Then scale row 1 by $\tfrac12$ and clear below:
$$A \to \begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 3 & 6 & -3 \\ 0 & -1 & -1 & 3 \\ 0 & 2 & 8 & 8 \end{bmatrix}$$
Look at rows 2–4. The leading column is $\mathbf{c}_2 = (3,-1,2)$ (from rows 2,3,4 of column 2). Store it as column 2 of $L$ (placed at the bottom). Scale and clear below:
$$\to \begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 1 & 2 & -1 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 4 & 10 \end{bmatrix}$$
Rows 3–4 now. Leading column $\mathbf{c}_3 = (1,4)$. Store it as column 3 of $L$. Clear below (row 4 $-\,4\times$ row 3):
$$\to \begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 1 & 2 & -1 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 0 & 2 \end{bmatrix}$$
The last submatrix is just the single entry $2$, so $\mathbf{c}_4 = (2)$. Store it. Scale row 4 to a leading $1$:
$$U = \begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 1 & 2 & -1 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 0 & 1 \end{bmatrix}$$
Now assemble $L$: drop the four stored leading columns into the bottom of columns 1–4 of the identity. Column 1 gets $(2,1,3,-1)$, column 2 gets $(3,-1,2)$ (starting at row 2), column 3 gets $(1,4)$ (starting at row 3), column 4 gets $(2)$ (row 4):
$$L = \begin{bmatrix} 2 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 \\ 3 & -1 & 1 & 0 \\ -1 & 2 & 4 & 2 \end{bmatrix}, \qquad U = \begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 1 & 2 & -1 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 0 & 1 \end{bmatrix}.$$
You can multiply $LU$ and confirm it reproduces $A$ exactly. Notice the pattern: each stored leading column slots straight into $L$, and $U$ is just the row-echelon form. The factorization fell out of ordinary elimination.
§9Exercises — Nicholson §2.7
(a) $A = \begin{bmatrix} 2 & 6 & -2 & 0 & 2 \\ 3 & 9 & -3 & 3 & 1 \\ -1 & -3 & 1 & -3 & 1 \end{bmatrix}$.
(b) $A = \begin{bmatrix} 2 & 4 & 2 \\ 1 & -1 & 3 \\ -1 & 7 & -7 \end{bmatrix}$.
$A = \begin{bmatrix} 2 & 0 & 0 \\ 0 & -1 & 0 \\ 1 & 1 & 3 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 0 & 1 \end{bmatrix} = LU$, $\;\mathbf{b} = \begin{bmatrix} 1 \\ -1 \\ 2 \end{bmatrix}$. Solve by finding $\mathbf{y}$ with $L\mathbf{y} = \mathbf{b}$, then $\mathbf{x}$ with $U\mathbf{x} = \mathbf{y}$.
Show that $\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} = LU$ is impossible, where $L$ is lower triangular and $U$ is upper triangular.
§10Section 2.8 — An Application to Input–Output Economics
Leontief's "input–output" tables eventually catalogued hundreds of sectors of the U.S. economy, and the computations were among the first serious industrial uses of early computers. A whole branch of economics grew from organising production as a matrix.
Let $E$ be the input–output matrix: the entry $e_{ij}$ is the fraction of industry $j$'s output consumed by industry $i$. A price vector $\mathbf{p}$ is an equilibrium (each industry breaks even) exactly when
$$E\mathbf{p} = \mathbf{p}, \qquad\text{equivalently}\qquad (I - E)\mathbf{p} = \mathbf{0}.$$
We seek a nonzero, nonnegative solution $\mathbf{p}$. Since this is homogeneous, equilibrium prices are determined only up to a positive scale factor — only the ratios of prices matter, which makes economic sense (currency units are arbitrary).
Find the equilibrium price structure for the input–output matrix $E = \begin{bmatrix} 0.1 & 0.2 & 0.3 \\ 0.6 & 0.2 & 0.3 \\ 0.3 & 0.6 & 0.4 \end{bmatrix}$.
(Each column sums to $1$ — every industry's entire output is distributed somewhere, as it should be.)
Find the equilibrium price structures for three industries whose input–output matrix is $E = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}$, and discuss why the answer is unusual.
§11Exercises — Nicholson §2.8
(b) $E = \begin{bmatrix} 0.5 & 0 & 0.5 \\ 0.1 & 0.9 & 0.2 \\ 0.4 & 0.1 & 0.3 \end{bmatrix}$.
Industries $A, B, C$ are such that all output of $A$ is used by $B$, all output of $B$ is used by $C$, and all output of $C$ is used by $A$. Find the possible equilibrium price structures.
We have split matrices into triangular pieces and watched economies balance. Next, the determinant takes centre stage.
The determinant is the single number that decides invertibility, measures how a matrix scales volume, and — pleasingly — is trivial to read off a triangular matrix (just multiply the diagonal). Everything you learned about triangular shapes today pays off immediately in the next chapter.