Solving with PA = LU

Elementary Matrices and LU Factorization records elimination as $A = LU$. With Partial Pivoting the rows get swapped along the way; a permutation matrix $P$ (the identity with its rows reordered) records the swaps, and the factorization becomes $$PA = LU,$$ with $L$ unit lower triangular (its multipliers all at most 1 in size) and $U$ upper triangular.

Key idea

To solve $A\mathbf{x} = \mathbf{b}$ from $PA = LU$: reorder $\mathbf{b}$ to $P\mathbf{b}$, solve $L\mathbf{y} = P\mathbf{b}$ by forward substitution, then $U\mathbf{x} = \mathbf{y}$ by back substitution. Factoring costs $\frac{2}{3}n^{3}$ flops once; each new $\mathbf{b}$ costs only $2n^{2}$.

Why it works

Multiply $A\mathbf{x} = \mathbf{b}$ by $P$: $PA\mathbf{x} = P\mathbf{b}$, that is, $L(U\mathbf{x}) = P\mathbf{b}$. Call $U\mathbf{x} = \mathbf{y}$. Then $L\mathbf{y} = P\mathbf{b}$ is a lower triangular system, solved from the top down, and $U\mathbf{x} = \mathbf{y}$ is upper triangular, solved from the bottom up. Both are cheap because each unknown comes from one division once the earlier ones are known.

Worked example: Solving with the factors

For $A = \begin{bmatrix} 2 & 1 & -1 \\ -6 & 3 & 2 \\ 4 & -2 & 5 \end{bmatrix}$, pivoting gives $P$ = (swap rows 1 and 2), $$L = \begin{bmatrix} 1 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 \\ -\frac{2}{3} & 0 & 1 \end{bmatrix},\quad U = \begin{bmatrix} -6 & 3 & 2 \\ 0 & 2 & -\frac{1}{3} \\ 0 & 0 & \frac{19}{3} \end{bmatrix}.$$ Solve $A\mathbf{x} = (3, 2, 5)$.

$P\mathbf{b}$: swap the first two entries, $(2, 3, 5)$.
Forward: $y_1 = 2$; $y_2 = 3 - \left(-\frac{1}{3}\right)(2) = \frac{11}{3}$; $y_3 = 5 - \left(-\frac{2}{3}\right)(2) - 0 = \frac{19}{3}$.
Back: $x_3 = \frac{19/3}{19/3} = 1$; $x_2 = \frac{\frac{11}{3} + \frac{1}{3}(1)}{2} = 2$; $x_1 = \frac{2 - 3(2) - 2(1)}{-6} = 1$.
Check: $A(1, 2, 1) = (2 + 2 - 1,\ -6 + 6 + 2,\ 4 - 4 + 5) = (3, 2, 5)$.

Answer: $\mathbf{x} = (1, 2, 1)$

Worked example: A new right-hand side, no new elimination

Same factors; now solve $A\mathbf{x} = (3, -8, -1)$.

$P\mathbf{b} = (-8, 3, -1)$.
Forward: $y_1 = -8$; $y_2 = 3 + \frac{1}{3}(-8) = \frac{1}{3}$; $y_3 = -1 + \frac{2}{3}(-8) = -\frac{19}{3}$.
Back: $x_3 = -1$; $x_2 = \frac{\frac{1}{3} + \frac{1}{3}(-1)}{2} = 0$; $x_1 = \frac{-8 - 0 - 2(-1)}{-6} = 1$.
Only about $2n^{2}$ flops: the factorization was reused.

Answer: $\mathbf{x} = (1, 0, -1)$

Solving with PA = LU: reorder b, forward substitution, back substitution, and a check that Ax gives back b.

This split into "factor" and "solve" is how numerical libraries are organized, and it pays off whenever the same matrix meets many right-hand sides: Inverse and Shifted Power Methods solves with the same shifted matrix every step, implicit ODE methods with a fixed step solve with the same matrix every step, and the spline and boundary value systems later in the course are factored once and reused.

Reading $P$ is easy: row $i$ of $P$ has a single $1$, in column $j$, which means "entry $i$ of $P\mathbf{b}$ is entry $j$ of $\mathbf{b}$." Software usually stores $P$ as a list of row indices rather than as a matrix. In the example above, $P$ moves row 2 to the top, so $P\mathbf{b}$ starts with $b_2$.

Common mistake

Forgetting to apply the row swaps to $\mathbf{b}$ before forward substitution. $L$ and $U$ describe the swapped rows, so using $\mathbf{b}$ in its original order solves a different system. Check: multiply your answer by $A$; it must give back the original $\mathbf{b}$.

The same factors give $\det A = \pm u_{11}u_{22}\cdots u_{nn}$ (the sign set by the number of swaps), far cheaper than cofactor expansion.

Sources and further reading

Tobin A. Driscoll and Richard J. Braun, Fundamentals of Numerical Computation (University of Delaware, free online edition): chapter "Square linear systems", sections "LU factorization" and "Row pivoting"

MIT OpenCourseWare 18.335J, Introduction to Numerical Methods: lecture notes on LU factorization with pivoting

Solving with PA = LU

Practice

Builds on

Leads to

Sources and further reading