Lecture 06: Factorizations: A = CR and A = LU

The idea in one sentence

A factorization writes $A$ as a product of simpler matrices, and the shape of the pieces answers a question the original block of numbers was keeping to itself.

Two questions a matrix will not answer directly

A matrix arrives as a rectangle of numbers. Two questions come up immediately, and neither one can be read off the rectangle.

The first: how much is actually in here? A table with four columns may only hold three columns worth of information, because one of them is a copy of the others in disguise. You cannot see that by looking.

The second: I need to solve $Ax = b$ a hundred times with a hundred different right-hand sides. Elimination worked the first time. Do I have to do all of it again?

Each question gets a factorization. $A = CR$ answers the first by pulling out the columns that carry real information. $A = LU$ answers the second by writing down what elimination did, so it never has to be done twice. Both are the same move: replace one hard matrix by a product of easy ones.

A = CR: the columns that matter, and the recipe

Here is a matrix with something to hide:

A = \begin{bmatrix} 2 & 1 & 3 & 0 \\ 3 & 4 & 7 & 0 \\ 1 & 2 & 3 & -1 \end{bmatrix}.

Walk left to right and ask each column one question: are you new, or can I already build you?

Column 1 is $(2, 3, 1)$ . Nothing came before it, so it is new. Column 2 is $(1, 4, 2)$ . To be a copy of column 1 it would have to be a single multiple of it, and $1/2$ , $4/3$ and $2/1$ disagree. New. Column 3 is $(3, 7, 3)$ , and $2 + 1 = 3$ , $3 + 4 = 7$ , $1 + 2 = 3$ . Column 3 is column 1 plus column 2. Old. Column 4 is $(0, 0, -1)$ . Building it from the first two would need $2a + b = 0$ together with $3a + 4b = 0$ , which force $a = b = 0$ , and then the last entry says $0 = -1$ . New.

Three new columns. Stand them side by side and call that matrix $C$ :

C = \begin{bmatrix} 2 & 1 & 0 \\ 3 & 4 & 0 \\ 1 & 2 & -1 \end{bmatrix}.

Now every column of $A$ is a combination of the columns of $C$ , and each combination is three numbers. Write those three numbers as a column and collect them:

R = \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}.

Read $R$ column by column. Column 1 of $R$ is $(1, 0, 0)$ : take one of the first ingredient and none of the others, which returns column 1 of $A$ . Column 3 is $(1, 1, 0)$ : one of the first, one of the second, which is the relation you found. Column 4 is $(0, 0, 1)$ : the third ingredient on its own.

$C$ holds the ingredients. $R$ holds the recipes. Their product is $A$ .

Fig. 1

Green marks $C$ , the columns that carry the information. Yellow marks the column that was a copy and the recipe column of $R$ that rebuilds it from the other two.

The number of ingredients has a name. It is the rank of $A$ , written $\rank(A)$ , and here it is 3. That is the point of this factorization. You never compute the rank from $A$ by a formula. You read it off as the width of $C$ .

Look at the shapes. $A$ is $m$ by $n$ . $C$ is $m$ by $r$ , where $r$ is the rank. $R$ is $r$ by $n$ . The product has to pass through a waist of width $r$ .

Fig. 2

The rank is the waist of the factorization. Green is $C$ , holding $r$ columns that span the column space; orange is $R$ , holding $r$ rows that span the row space. The same $r$ on both sides, which is why the two counts can never disagree.

That waist pays a dividend right away. Read the product the other way: row $i$ of $A$ is row $i$ of $C$ applied to $R$ , so every row of $A$ is a combination of the $r$ rows of $R$ . That means $r$ independent columns force at most $r$ independent rows. Run the same argument on $A$ transposed and the inequality comes back the other way. The number of independent rows and the number of independent columns have to be the same. Lecture 9 proves it properly.

python

import numpy as np

A = np.array([[2., 1, 3, 0], [3, 4, 7, 0], [1, 2, 3, -1]])
print(A[:, 2] - A[:, 0] - A[:, 1])  # [0. 0. 0.]  col 3 = col 1 + col 2

C = A[:, [0, 1, 3]]  # keep the columns that were new
R = np.array([[1., 0, 1, 0], [0, 1, 1, 0], [0, 0, 0, 1]])
print(C @ R - A)                 # all zeros: A = CR
print(C.shape, R.shape)          # (3, 3) (3, 4)
print(np.linalg.matrix_rank(A))  # 3  the width of C

python

import torch

A = torch.tensor([[2., 1, 3, 0], [3, 4, 7, 0], [1, 2, 3, -1]])
print(A[:, 2] - A[:, 0] - A[:, 1])  # tensor([0., 0., 0.])

C = A[:, [0, 1, 3]]  # keep the columns that were new
R = torch.tensor([[1., 0, 1, 0], [0, 1, 1, 0], [0, 0, 0, 1]])
print(C @ R - A)                    # all zeros: A = CR
print(C.shape, R.shape)             # torch.Size([3, 3]) torch.Size([3, 4])
print(torch.linalg.matrix_rank(A))  # tensor(3)  the width of C

python

import jax.numpy as jnp

A = jnp.array([[2., 1, 3, 0], [3, 4, 7, 0], [1, 2, 3, -1]])
print(A[:, 2] - A[:, 0] - A[:, 1])  # [0. 0. 0.]

C = A[:, [0, 1, 3]]  # keep the columns that were new
R = jnp.array([[1., 0, 1, 0], [0, 1, 1, 0], [0, 0, 0, 1]])
print(C @ R - A)                  # all zeros: A = CR
print(C.shape, R.shape)           # (3, 3) (3, 4)
print(jnp.linalg.matrix_rank(A))  # 3  the width of C

python

import tensorflow as tf

A = tf.constant([[2., 1, 3, 0], [3, 4, 7, 0], [1, 2, 3, -1]])
print(A[:, 2] - A[:, 0] - A[:, 1])  # [0. 0. 0.]

C = tf.gather(A, [0, 1, 3], axis=1)  # tf spells column pick "gather"
R = tf.constant([[1., 0, 1, 0], [0, 1, 1, 0], [0, 0, 0, 1]])
print(C @ R - A)                 # all zeros: A = CR
print(C.shape, R.shape)          # (3, 3) (3, 4)
print(tf.linalg.matrix_rank(A))  # 3  the width of C

This factorization pays for itself in a place you may already have used. When people fine-tune a large model they freeze the weight matrix $W$ and learn a thin $C$ and a wide $R$ beside it. The trick is called LoRA, and the thing they train is $CR$ , whose waist is the whole point.

Fig. 3

$W$ stays frozen and the learning happens in the two slivers beside it. The update $CR$ has the same $4096 \times 4096$ shape, so it can touch every entry of $W$ , yet its waist is 8: you train 65,536 numbers instead of 16,777,216. That is 0.39% of the matrix, and the rank is what sets the price.

A = LU: the receipt elimination leaves behind

Now the second question. Take the system from Lecture 4:

A = \begin{bmatrix} 2 & 4 & -2 \\ 4 & 9 & -3 \\ -2 & -3 & 7 \end{bmatrix}.

Elimination cleared the entries below the diagonal in three steps. Each step subtracted a multiple of one row from another, and that multiple is called a multiplier. Write $\ell_{ij}$ for the number that multiplied row $j$ when it was subtracted from row $i$ :

\ell_{21} = 2, \qquad \ell_{31} = -1, \qquad \ell_{32} = 1,

and the result is the triangular matrix

U = \begin{bmatrix} 2 & 4 & -2 \\ 0 & 1 & 1 \\ 0 & 0 & 4 \end{bmatrix}.

It is tempting to keep $U$ and throw the multipliers away. That wastes the work. Put the three multipliers into a matrix instead, each one in the slot it came from, with ones down the diagonal:

L = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 1 & 1 \end{bmatrix}.

Fig. 4

The receipt. Each elimination step spends one multiplier, and that multiplier is filed at the address of the entry it cleared. Nothing is mixed and nothing is rescaled on the way in.

Then $LU = A$ , exactly. Check the bottom row by hand. Row 3 of $LU$ is $-1$ times row 1 of $U$ , plus $1$ times row 2 of $U$ , plus $1$ times row 3 of $U$ :

-1\,(2, 4, -2) + 1\,(0, 1, 1) + 1\,(0, 0, 4) \;=\; (-2, -3, 7),

which is row 3 of $A$ . That computation is elimination read backwards. Steps 2 and 3 said

\text{row 3 of } U \;=\; \text{row 3 of } A \;-\; \ell_{31}\,(\text{row 1 of } U) \;-\; \ell_{32}\,(\text{row 2 of } U).

The rows being subtracted are rows of $U$ , because a pivot row is finished the moment it becomes a pivot row and never changes again. Move the two terms across the equals sign and row 3 of $A$ is a combination of the rows of $U$ with the multipliers as weights. That is exactly what $LU$ computes.

Fig. 5

The proof, drawn. Orange holds the finished rows of $U$ , yellow holds the weights, and green is what comes out: row 3 of $A$ , entry for entry.

Now the payoff on the second question. To solve $Ax = b$ , solve $LUx = b$ in two easy stages. Set $c = Ux$ and solve $Lc = b$ going down, then solve $Ux = c$ going up. Both matrices are triangular, so each stage is a walk through the entries with no searching. With $b = (2, 8, 10)$ this gives $c = (2, 4, 8)$ and then $x = (-1, 2, 2)$ . The expensive part, building $L$ and $U$ , happens once. Every right-hand side after that is cheap.

Fig. 6

Why the factorization pays. A triangular system solves itself: the first row of $Lc = b$ gives $c_1$ outright, and every row after it has one new unknown. The orange walk goes down, the green walk comes back up, and yellow is $c$ , the answer the first walk hands to the second.

Why the multipliers land unmixed

Something in the last section deserves suspicion. Matrix multiplication normally scrambles everything. Multiply two matrices carrying a 5 and a 4 and you expect a 20 to appear somewhere. Yet the multipliers walked into $L$ untouched. Here is the smallest example that shows why.

Let $E$ subtract 5 times row 1 from row 2, and let $F$ subtract 4 times row 2 from row 3:

E = \begin{bmatrix} 1 & 0 & 0 \\ -5 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \qquad F = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -4 & 1 \end{bmatrix}.

Elimination does $E$ first and $F$ second, so the combined step is the product $FE$ :

FE = \begin{bmatrix} 1 & 0 & 0 \\ -5 & 1 & 0 \\ 20 & -4 & 1 \end{bmatrix}.

There is the 20, sitting in the corner. It is not noise. $E$ replaced row 2 by row 2 minus 5 times row 1. Then $F$ subtracted 4 times that new row 2 from row 3, and the new row 2 had row 1 baked into it. So row 3 came out carrying $(-4)(-5) = 20$ times row 1. Going forward, row 3 feels row 1.

Now undo it. Undoing $E$ means adding 5 times row 1 back to row 2, and undoing $F$ means adding 4 times row 2 back to row 3. The last step taken has to be the first step undone, so $F$ is undone first and the product is $E^{-1}F^{-1}$ :

E^{-1}F^{-1} = \begin{bmatrix} 1 & 0 & 0 \\ 5 & 1 & 0 \\ 0 & 4 & 1 \end{bmatrix}.

Clean. The 20 is gone and the two multipliers sit in their own slots. In this order $F^{-1}$ acts first, and at that moment row 2 is still a finished row of $U$ with no row 1 in it. Then $E^{-1}$ repairs row 2 and never touches row 3. Row 3 never hears about row 1 at all.

Fig. 7

The blue arc is the leak. Going forward, row 1 reaches row 3 through the row 2 that changed underneath it, and the product records that with a 20, marked pink. Going backward in the reversed order that path does not exist, so the corner stays zero and the multipliers survive.

Order is the whole trick. Multiply the inverses in the other order, $F^{-1}E^{-1}$ , and the 20 comes straight back. Lecture 5 gave the rule: undoing a sequence reverses it, socks before shoes on the way in and shoes before socks on the way out. That is what makes $L$ readable.

Try it. Eliminate, and watch the receipt fill in

step 0 of 6

Press Next. Three steps clear the matrix into U, and each one files its multiplier in L. Three more steps replay L times U and rebuild A.

Blue is the row being changed, yellow is the multiplier just filed, green rings the pivots.

python

import numpy as np

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

E21 = np.array([[1., 0, 0], [-2, 1, 0], [0, 0, 1]])  # 2 x row1 off row2
E31 = np.array([[1., 0, 0], [0, 1, 0], [1, 0, 1]])   # -1 x row1 off row3
E32 = np.array([[1., 0, 0], [0, 1, 0], [0, -1, 1]])  # 1 x row2 off row3
U = E32 @ E31 @ E21 @ A
print(U)  # [[2. 4. -2.] [0. 1. 1.] [0. 0. 4.]]

L = np.array([[1., 0, 0], [2, 1, 0], [-1, 1, 1]])  # multipliers in place
print(L @ U - A)                                   # all zeros

b = np.array([2., 8, 10])
c = np.linalg.solve(L, b)  # [2. 4. 8.]   forward through L
x = np.linalg.solve(U, c)  # [-1. 2. 2.]  back through U
print(c, x)

python

import torch

A = torch.tensor([[2., 4, -2], [4, 9, -3], [-2, -3, 7]])

E21 = torch.tensor([[1., 0, 0], [-2, 1, 0], [0, 0, 1]])  # 2 x row1 off row2
E31 = torch.tensor([[1., 0, 0], [0, 1, 0], [1, 0, 1]])   # -1 x row1 off row3
E32 = torch.tensor([[1., 0, 0], [0, 1, 0], [0, -1, 1]])  # 1 x row2 off row3
U = E32 @ E31 @ E21 @ A
print(U)  # [[2. 4. -2.] [0. 1. 1.] [0. 0. 4.]]

L = torch.tensor([[1., 0, 0], [2, 1, 0], [-1, 1, 1]])  # multipliers in place
print(L @ U - A)                                       # all zeros

b = torch.tensor([2., 8, 10])
c = torch.linalg.solve(L, b)  # tensor([2., 4., 8.])   forward through L
x = torch.linalg.solve(U, c)  # tensor([-1., 2., 2.])  back through U
print(c, x)

python

import jax.numpy as jnp

A = jnp.array([[2., 4, -2], [4, 9, -3], [-2, -3, 7]])

E21 = jnp.array([[1., 0, 0], [-2, 1, 0], [0, 0, 1]])  # 2 x row1 off row2
E31 = jnp.array([[1., 0, 0], [0, 1, 0], [1, 0, 1]])   # -1 x row1 off row3
E32 = jnp.array([[1., 0, 0], [0, 1, 0], [0, -1, 1]])  # 1 x row2 off row3
U = E32 @ E31 @ E21 @ A
print(U)  # [[2. 4. -2.] [0. 1. 1.] [0. 0. 4.]]

L = jnp.array([[1., 0, 0], [2, 1, 0], [-1, 1, 1]])  # multipliers in place
print(L @ U - A)                                    # all zeros

b = jnp.array([2., 8, 10])
c = jnp.linalg.solve(L, b)  # [2. 4. 8.]   forward through L
x = jnp.linalg.solve(U, c)  # [-1. 2. 2.]  back through U
print(c, x)

python

import tensorflow as tf

A = tf.constant([[2., 4, -2], [4, 9, -3], [-2, -3, 7]])

E21 = tf.constant([[1., 0, 0], [-2, 1, 0], [0, 0, 1]])  # 2 x row1 off row2
E31 = tf.constant([[1., 0, 0], [0, 1, 0], [1, 0, 1]])   # -1 x row1 off row3
E32 = tf.constant([[1., 0, 0], [0, 1, 0], [0, -1, 1]])  # 1 x row2 off row3
U = E32 @ E31 @ E21 @ A
print(U)  # [[2. 4. -2.] [0. 1. 1.] [0. 0. 4.]]

L = tf.constant([[1., 0, 0], [2, 1, 0], [-1, 1, 1]])  # multipliers in place
print(L @ U - A)                                      # all zeros

b = tf.constant([[2.], [8.], [10.]])  # solve wants b as a column
c = tf.linalg.solve(L, b)             # [[2.] [4.] [8.]]    forward through L
x = tf.linalg.solve(U, c)             # [[-1.] [2.] [2.]]   back through U
print(c, x)

LDU: pulling the pivots out

$L$ and $U$ look like a mismatched pair. $L$ has ones on its diagonal, and $U$ has the pivots $2, 1, 4$ on its. That asymmetry is cosmetic, and one more factor removes it. Divide each row of $U$ by its own pivot and park the pivots in a diagonal matrix $D$ :

\begin{bmatrix} 2 & 4 & -2 \\ 0 & 1 & 1 \\ 0 & 0 & 4 \end{bmatrix} = \begin{bmatrix} 2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 4 \end{bmatrix} \begin{bmatrix} 1 & 2 & -1 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{bmatrix}.

So $A = LDU$ , with $L$ unit lower triangular, $D$ diagonal, and the new $U$ unit upper triangular. The two triangular factors are now built the same way and the pivots stand alone in the middle where you can read them. Their product is the determinant, which Lecture 14 collects.

This form is also the only one of its kind. Suppose $L_1 D_1 U_1 = L_2 D_2 U_2$ , with all four triangular factors carrying unit diagonals. Push the $L$ s to the left and the $U$ s to the right: $L_2^{-1} L_1 D_1 = D_2 U_2 U_1^{-1}$ . The left side is lower triangular, the right side is upper triangular, and a matrix that is both is diagonal. Unit diagonals then force $L_2^{-1}L_1 = I$ and $U_2U_1^{-1} = I$ . Uniqueness looks like bookkeeping. The next section cashes it in: a symmetric matrix has only one triangle to store.

Transposes, and what they are for

Flipping a matrix across its diagonal turns rows into columns. Written $A\T$ , the transpose obeys three rules:

(A + B)\T = A\T + B\T, \qquad (AB)\T = B\T A\T, \qquad (A^{-1})\T = (A\T)^{-1}.

The middle rule is the one worth deriving, because the reversal surprises people. Entry $(i,j)$ of $(AB)\T$ is entry $(j,i)$ of $AB$ , which is row $j$ of $A$ dotted with column $i$ of $B$ . Read that same dot product from the other side: it is row $i$ of $B\T$ dotted with column $j$ of $A\T$ , which is entry $(i,j)$ of $B\T A\T$ . The order flips because transposing turns the outside indices into the inside ones. The third rule follows from the second: transpose $AA^{-1} = I$ to get $(A^{-1})\T A\T = I$ .

Here is what those rules buy. Suppose $A$ is symmetric, so $A = A\T$ . Now transpose both sides of $A = LDU$ , using the product rule twice and $D\T = D$ because $D$ is diagonal:

A = A\T = (LDU)\T = U\T D L\T .

Look at $U\T D L\T$ . It is a unit lower triangular matrix, times a diagonal, times a unit upper triangular matrix, which is an $LDU$ factorization of $A$ . And $LDU$ is unique. So $U\T = L$ , or equivalently $U = L\T$ , and every symmetric matrix factors as

A = L D L\T .

Take the tridiagonal matrix $K$ from Lecture 5. Elimination gives multipliers $-1/2$ and $-2/3$ and pivots $2$ , $3/2$ , $4/3$ , and the factorization mirrors itself:

K = \begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{bmatrix} = L D L\T .

Fig. 8

A symmetric matrix folds along its diagonal. Cyan is one multiplier and its mirror image, yellow is the pivots. Software stores one triangle and the pivots, then gets the other triangle free.

That is worth real money. A symmetric matrix needs half the storage and half the elimination work, because computing $L$ hands you $L\T$ for nothing. Covariance matrices are symmetric, and so are the normal equations of Lecture 12. That shortcut runs underneath much of what machine learning code does.

python

import numpy as np

K = np.array([[2., -1, 0], [-1, 2, -1], [0, -1, 2]])
L = np.array([[1., 0, 0], [-0.5, 1, 0], [0, -2 / 3, 1]])
D = np.diag([2., 1.5, 4 / 3])

print(L @ D @ L.T - K)      # zeros: K = L D L transpose
print(np.prod(np.diag(D)))  # 4.0   product of the pivots
print(np.linalg.det(K))     # 4.0   and the determinant

python

import torch

K = torch.tensor([[2., -1, 0], [-1, 2, -1], [0, -1, 2]])
L = torch.tensor([[1., 0, 0], [-0.5, 1, 0], [0, -2 / 3, 1]])
D = torch.diag(torch.tensor([2., 1.5, 4 / 3]))

print(L @ D @ L.T - K)            # zeros: K = L D L transpose
print(torch.prod(torch.diag(D)))  # tensor(4.)   product of pivots
print(torch.linalg.det(K))        # tensor(4.)   and determinant

python

import jax.numpy as jnp

K = jnp.array([[2., -1, 0], [-1, 2, -1], [0, -1, 2]])
L = jnp.array([[1., 0, 0], [-0.5, 1, 0], [0, -2 / 3, 1]])
D = jnp.diag(jnp.array([2., 1.5, 4 / 3]))

print(L @ D @ L.T - K)        # zeros: K = L D L transpose
print(jnp.prod(jnp.diag(D)))  # 4.0   product of the pivots
print(jnp.linalg.det(K))      # 4.0   and the determinant

python

import tensorflow as tf

K = tf.constant([[2., -1, 0], [-1, 2, -1], [0, -1, 2]])
L = tf.constant([[1., 0, 0], [-0.5, 1, 0], [0, -2 / 3, 1]])
D = tf.linalg.diag(tf.constant([2., 1.5, 4 / 3]))

print(L @ D @ tf.transpose(L) - K)             # zeros: K = L D L transpose
print(tf.reduce_prod(tf.linalg.diag_part(D)))  # 4.0   prod is reduce_prod
print(tf.linalg.det(K))                        # 4.0   and the determinant

When a pivot is zero: PA = LU

Everything above assumed elimination never got stuck. It gets stuck. Try this matrix:

A = \begin{bmatrix} 1 & 2 & 1 \\ 2 & 4 & 3 \\ 1 & 1 & 2 \end{bmatrix}.

Clear the first column. Row 2 minus 2 times row 1 gives $(0, 0, 1)$ . Row 3 minus row 1 gives $(0, -1, 1)$ . Now look at the second pivot position, the middle of the matrix. It holds a zero, and you cannot divide by zero. Elimination has nothing to work with.

Nothing is wrong with $A$ . Its rank is 3 and it is invertible. The rows arrived in an unlucky order. Exchange rows 2 and 3 and the problem disappears.

Exchanging rows is itself a matrix. A permutation matrix is the identity with its rows reordered, and multiplying by it reorders the rows of whatever it hits. To swap rows 2 and 3, take the identity and swap rows 2 and 3:

P = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}, \qquad PA = \begin{bmatrix} 1 & 2 & 1 \\ 1 & 1 & 2 \\ 2 & 4 & 3 \end{bmatrix}.

Permutation matrices have a pleasant property that falls out of the transpose rules. The columns of $P$ are the standard basis vectors in some order, so entry $(i,j)$ of $P\T P$ is column $i$ of $P$ dotted with column $j$ of $P$ , which is 1 when $i = j$ and 0 otherwise. So $P\T P = I$ , and

P^{-1} = P\T .

To undo a shuffle, transpose it. Now eliminate $PA$ . Row 2 minus 1 times row 1 gives $(0, -1, 1)$ , and row 3 minus 2 times row 1 gives $(0, 0, 1)$ , which already has its zero in the right place, so the third multiplier is 0. The multipliers are $\ell_{21} = 1$ , $\ell_{31} = 2$ , $\ell_{32} = 0$ :

L = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 2 & 0 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 1 & 2 & 1 \\ 0 & -1 & 1 \\ 0 & 0 & 1 \end{bmatrix}, \qquad PA = LU .

Fig. 9

The zero in the middle comes from the order the rows arrived in. The matrix itself is fine. Watch the orange and cyan rows trade places: after the swap, the same elimination runs to the end and lands three yellow pivots.

So the honest general statement is $PA = LU$ . For every square matrix there is a row order that makes elimination work, and $P$ records it.

python

import numpy as np

A = np.array([[1., 2, 1], [2, 4, 3], [1, 1, 2]])
E = np.array([[1., 0, 0], [-2, 1, 0], [-1, 0, 1]])  # clear column 1
print(E @ A)                 # [[1. 2. 1.] [0. 0. 1.] [0. -1. 1.]]  pivot is 0

P = np.array([[1., 0, 0], [0, 0, 1], [0, 1, 0]])  # swap rows 2 and 3
print(P.T @ P)                            # identity: P inverse is P transpose

L = np.array([[1., 0, 0], [1, 1, 0], [2, 0, 1]])
U = np.array([[1., 2, 1], [0, -1, 1], [0, 0, 1]])
print(L @ U - P @ A)  # all zeros: PA = LU

python

import torch

A = torch.tensor([[1., 2, 1], [2, 4, 3], [1, 1, 2]])
E = torch.tensor([[1., 0, 0], [-2, 1, 0], [-1, 0, 1]])  # clear column 1
print(E @ A)                 # [[1. 2. 1.] [0. 0. 1.] [0. -1. 1.]]  pivot is 0

P = torch.tensor([[1., 0, 0], [0, 0, 1], [0, 1, 0]])  # swap rows 2 and 3
print(P.T @ P)                            # identity: P inverse is P transpose

L = torch.tensor([[1., 0, 0], [1, 1, 0], [2, 0, 1]])
U = torch.tensor([[1., 2, 1], [0, -1, 1], [0, 0, 1]])
print(L @ U - P @ A)  # all zeros: PA = LU

python

import jax.numpy as jnp

A = jnp.array([[1., 2, 1], [2, 4, 3], [1, 1, 2]])
E = jnp.array([[1., 0, 0], [-2, 1, 0], [-1, 0, 1]])  # clear column 1
print(E @ A)                 # [[1. 2. 1.] [0. 0. 1.] [0. -1. 1.]]  pivot is 0

P = jnp.array([[1., 0, 0], [0, 0, 1], [0, 1, 0]])  # swap rows 2 and 3
print(P.T @ P)                            # identity: P inverse is P transpose

L = jnp.array([[1., 0, 0], [1, 1, 0], [2, 0, 1]])
U = jnp.array([[1., 2, 1], [0, -1, 1], [0, 0, 1]])
print(L @ U - P @ A)  # all zeros: PA = LU

python

import tensorflow as tf

A = tf.constant([[1., 2, 1], [2, 4, 3], [1, 1, 2]])
E = tf.constant([[1., 0, 0], [-2, 1, 0], [-1, 0, 1]])  # clear column 1
print(E @ A)                 # [[1. 2. 1.] [0. 0. 1.] [0. -1. 1.]]  pivot is 0

P = tf.constant([[1., 0, 0], [0, 0, 1], [0, 1, 0]])  # swap rows 2 and 3
print(tf.transpose(P) @ P)                # identity: P inverse is P transpose

L = tf.constant([[1., 0, 0], [1, 1, 0], [2, 0, 1]])
U = tf.constant([[1., 2, 1], [0, -1, 1], [0, 0, 1]])
print(L @ U - P @ A)  # all zeros: PA = LU

Where this is going

Two factorizations, two answers. $A = LU$ freezes elimination so it never has to run twice, and $PA = LU$ patches the one case where it stalls. $A = CR$ found the columns that carry information and named their count the rank.

But $A = CR$ leaned on a word that has not been defined. Column 3 was called a copy because it was a combination of two others, and column 4 was called new because it was not. Behind that test sits a bigger object: the set of every vector you can build out of the columns of $A$ . It has a shape and a dimension, and it is the reason the rank means anything at all. Naming it takes the next chapter.

Check yourself

Factor $B = \begin{bmatrix} 1 & 3 & 4 \\ 2 & 1 & 3 \\ 1 & 2 & 3 \end{bmatrix}$ as $CR$ . Which column is the copy? What is $\rank(B)$ , and what shapes do $C$ and $R$ have?
Eliminate $\begin{bmatrix} 2 & 1 \\ 6 & 8 \end{bmatrix}$ by hand. Write down $U$ , the single multiplier, and $L$ . Multiply $LU$ and check you get back what you started with.
Let $E$ subtract 3 times row 1 from row 2, and let $F$ subtract 2 times row 2 from row 3. Compute $FE$ and $E^{-1}F^{-1}$ . Where does the cross term appear, and what number is it?
Factor $\begin{bmatrix} 4 & 2 \\ 2 & 3 \end{bmatrix}$ as $LDL\T$ . Check your answer by multiplying the three factors back together.
Which of these needs a row exchange before elimination can start: $\begin{bmatrix} 0 & 2 \\ 3 & 1 \end{bmatrix}$ or $\begin{bmatrix} 3 & 1 \\ 0 & 2 \end{bmatrix}$ ? Write the $P$ that fixes the one that does.
Redraw Figure 7 from memory, both panels and every label. The drawing is the test.