LU factorisation

Elimination throws its multipliers away. Keep them, and the matrix splits into two triangles. Watch each one drop into place.

Keep the multipliers

On the page Elimination we turned

with the multipliers $\ell_{21}=3$, $\ell_{31}=0$ and $\ell_{32}=2$. Now put those multipliers in a lower triangle, below a diagonal of ones:

Then $\mathbf{A}=\mathbf{L}\mathbf{U}$. This is the LU factorisation: $\mathbf{L}$ is lower triangular, $\mathbf{U}$ is upper triangular, and their product is $\mathbf{A}$.

Why does it work? Elimination made each row of $\mathbf{U}$ by taking multiples of earlier rows away from a row of $\mathbf{A}$. Run that backwards, and each row of $\mathbf{A}$ is its row of $\mathbf{U}$ plus those same multiples put back. Row 2 of $\mathbf{A}$ is 3 times row 1 of $\mathbf{U}$ plus row 2 of $\mathbf{U}$:

Row 3 of $\mathbf{A}$ is 0 times row 1 plus 2 times row 2 plus row 3 of $\mathbf{U}$: $2\,(0,2,-2)+(0,0,5)=(0,4,1)$. Each row of $\mathbf{L}$ lists how much of each row of $\mathbf{U}$ to add, and that is what the product $\mathbf{L}\mathbf{U}$ does.

So LU is not a new method. It is elimination, written down.

Watch the signs. On the elimination page, the step matrix $\mathbf{E}_{21}$ held $-3$, because the step takes 3 of row 1 away. In $\mathbf{L}$ the 3 appears with a plus sign, because rebuilding puts it back. In matrices, $\mathbf{L}$ undoes the steps in reverse order, as the page The inverse matrix said an undo must:

Solving with two triangles

Now solve $\mathbf{A}\mathbf{x}=\mathbf{b}$, which is $\mathbf{L}(\mathbf{U}\mathbf{x})=\mathbf{b}$. Give the inside a name, $\mathbf{c}=\mathbf{U}\mathbf{x}$, and split the job in two:

Both are triangles, so both are quick. Take $\mathbf{b}=(2,12,2)$, the right-hand side from the elimination page. $\mathbf{L}\mathbf{c}=\mathbf{b}$ is solved from the top down, a forward substitution:

These are the right-hand sides that elimination produced, 2, 6 and $-10$. Then $\mathbf{U}\mathbf{x}=\mathbf{c}$ is solved from the bottom up, by back substitution: $5z=-10$ gives $z=-2$; $2y-2z=6$ gives $2y+4=6$, so $y=1$; and $x+2y+z=2$ gives $x+2-2=2$, so $x=2$.

Many right-hand sides

The work is in the factoring. Once $\mathbf{L}$ and $\mathbf{U}$ are known, a new $\mathbf{b}$ needs only the two sweeps.

Try the second one yourself before the picture finishes: $\mathbf{b}=(4,12,5)$. Forward: $c_1=4$, $c_2=12-3\cdot4=0$, $c_3=5-0-2\cdot0=5$. Back: $5z=5$, so $z=1$; $2y-2=0$, so $y=1$; $x+2+1=4$, so $x=1$. Check: $\mathbf{A}\,(1,1,1)=(4,12,5)$.

How much cheaper is it? Factoring an $n\times n$ matrix takes about $n^3/3$ multiply-and-subtract steps, as the elimination page counted. Each sweep through a triangle takes about $n^2/2$, so a solve takes about $n^2$. For $n=1000$, the factoring takes about $3.3\times10^8$ steps and each later solve about $10^6$: about 330 times less.

Row exchanges

Elimination sometimes needs a row exchange, when a pivot is 0. Take

The first pivot is 0. Swap the rows first, with the permutation matrix $\mathbf{P}=\begin{bmatrix}0&1\\1&0\end{bmatrix}$, which is the identity with its rows swapped. Then $\mathbf{P}\mathbf{A}=\begin{bmatrix}3&1\\0&2\end{bmatrix}$ is already upper triangular, so $\mathbf{U}=\mathbf{P}\mathbf{A}$ and $\mathbf{L}=\mathbf{I}$. In general, with the exchanges done first,

Library programs exchange rows more often than this page does. They move the largest available number into each pivot place, because that keeps rounding errors small. SciPy, for example, makes two exchanges on the $3\times3$ matrix above, even though none of its pivots is 0. Its $\mathbf{L}$ and $\mathbf{U}$ then differ from ours, but they solve the same system.

Worked example

1. Read $\mathbf{L}$ from the multipliers. Multipliers $3$, $0$, $2$ give $\ell_{21}=3$, $\ell_{31}=0$, $\ell_{32}=2$ in $\mathbf{L}$, below a diagonal of ones.
2. Rebuild $\mathbf{A}$. Row 2: $3\,(1,2,1)+(0,2,-2)=(3,8,1)$. Row 3: $2\,(0,2,-2)+(0,0,5)=(0,4,1)$.
3. First solve. $\mathbf{b}=(2,12,2)$ gives $\mathbf{c}=(2,6,-10)$ and $\mathbf{x}=(2,1,-2)$.
4. Second solve. $\mathbf{b}=(4,12,5)$ gives $\mathbf{c}=(4,0,5)$ and $\mathbf{x}=(1,1,1)$, with no new elimination.
5. Signs. In $\mathbf{E}_{21}$ the entry in row 2, column 1 is $-3$; in $\mathbf{L}$ it is $+3$.
6. A zero pivot. $\begin{bmatrix}0&2\\3&1\end{bmatrix}$ needs $\mathbf{P}=\begin{bmatrix}0&1\\1&0\end{bmatrix}$; then $\mathbf{U}=\begin{bmatrix}3&1\\0&2\end{bmatrix}$ and $\mathbf{L}=\mathbf{I}$.

Where you’ll meet this

A circuit simulator stepping through time often solves with the same matrix at every step, with a new right-hand side each time. It factors once and then only sweeps. Programs that compute how a structure bends under many different loads do the same: one factoring, one quick solve per load. In SciPy, a free Python library, the pair of functions lu_factor and lu_solve splits the work exactly this way.

Reference card