NLA MT25, LU factorisation
Flashcards
@Describe how the LU factorisation of $A \in \mathbb R^{n \times n}$ can be interpreted as writing $A$ as a sum of $n$ rank-$1$ outer products.
If $A = LU$ with $L = [L _ 1 \mid \cdots \mid L _ n]$ (columns) and $U = [U _ 1 \; ; \; \cdots \; ; \; U _ n]$ (rows), then matrix multiplication gives
\[A = LU = \sum^n _ {i = 1} L _ i U _ i\]Each $L _ i U _ i$ is a rank-$1$ outer product (column vector times row vector). The triangular structure of $L$ and $U$ means we have
\[ L _ i = \begin{pmatrix} 0 _ {i-1} \\ * \\ \vdots \\ * \end{pmatrix}, \qquad U _ i = \begin{pmatrix} 0 _ {i-1}^\top & * & \cdots & * \end{pmatrix} \]so $L _ i U _ i$ has zeros in its first $i-1$ rows and first $i-1$ columns, with a non-zero $(n-i+1) \times (n-i+1)$ block in the bottom-right. For $n = 5$, this looks like
\[A = \underbrace{\begin{pmatrix} * \\ * \\ * \\ * \\ * \end{pmatrix} \begin{pmatrix} * & * & * & * & * \end{pmatrix}} _ {L _ 1 U _ 1} + \underbrace{\begin{pmatrix} 0 \\ * \\ * \\ * \\ * \end{pmatrix} \begin{pmatrix} 0 & * & * & * & * \end{pmatrix}} _ {L _ 2 U _ 2} + \cdots + \underbrace{\begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ * \end{pmatrix} \begin{pmatrix} 0 & 0 & 0 & 0 & * \end{pmatrix}} _ {L _ 5 U _ 5}. \]@Describe the pivoted Gaussian elimination @algorithm for computing the LU factorisation $PA = LU$ of $A \in \mathbb R^{n \times n}$, in terms of the rank-1 decomposition view of LU (∆rank-1-view-of-lu).
Initialise $P^{(0)} = I _ n$, $A^{(1)} = A$. At step $k = 1, 2, \ldots, n$:
Pick the pivot. Look at the first column of the current $(n - k + 1) \times (n - k + 1)$ sub-block $A^{(k)}$. Choose the row index $r \in \{k, k+1, \ldots, n\}$ for which $ \vert a _ {r1}^{(k)} \vert $ is maximal (if this maximum is $0$, then the algorithm fails. However, this can’t happen when $A$ is nonsingular, see ∆pivoted-lu-decomposition-existence).
Swap. If $r \ne k$, swap rows $k$ and $r$ of the entire current matrix. This is equivalent to left-multiplying by a transposition permutation matrix. Record this in the running $P^{(k)}$.
Extract a rank-1 outer product. Write the active block (after the swap) as
\[A^{(k)} = \begin{pmatrix} \alpha & u^\top \\ v & B \end{pmatrix}\]where $\alpha \in \mathbb R$ is the pivot, $u \in \mathbb R^{n - k}$ is the rest of the pivot row, $v \in \mathbb R^{n - k}$ is the rest of the pivot column, and $B \in \mathbb R^{(n-k) \times (n-k)}$ is the remainder. Then
\[\begin{pmatrix} \alpha & u^\top \\ v & B \end{pmatrix} = \underbrace{\begin{pmatrix} 1 \\ v/\alpha \end{pmatrix} \begin{pmatrix} \alpha & u^\top \end{pmatrix}} _ {\text{rank-}1} + \begin{pmatrix} 0 & 0 \\ 0 & B - v u^\top / \alpha \end{pmatrix}\]So we may set
\[L _ k = \begin{pmatrix} 0 _ {k-1} \\ 1 \\ v/\alpha \end{pmatrix} \in \mathbb R^n, \qquad U _ k = \begin{pmatrix} 0 _ {k-1}^\top & \alpha & u^\top \end{pmatrix} \in \mathbb R^n,\](where we pad with zeros so $L _ k U _ k$ acts on the right rows and columns of the full $n \times n$ matrix). This gives the partial rank-1 decomposition at step $k$ as
\[P^{(k)} A = \sum^k _ {i = 1} L _ i U _ i + \begin{pmatrix} 0 & 0 \\ 0 & A^{(k+1)} \end{pmatrix}.\]where next active block $A^{(k+1)}$ is the Schur complement
\[A^{(k+1)} = B - \frac{vu^\top}{\alpha} \in \mathbb R^{(n-k) \times (n-k)}.\]Iterate. Increment $k$ and repeat for $A^{(k+1)}$.
After $n$ steps, the accumulated rank-$1$ outer products give $L = [L _ 1 \mid \cdots \mid L _ n]$ (unit lower triangular), and $U = [U _ 1 \; ; \; \cdots \; ; \; U _ n]$ (upper triangular), with running permutation $P = P^{(n)}$ satisfying
\[PA = LU.\]@Prove that any nonsingular matrix $A$ has a pivoted LU decomposition, i.e. that there exists a permutation matrix $P$ such that
\[PA = LU\]
We know that pivoted Gaussian elimination will produce an $LU$ decomposition when it succeeds (∆pivoted-lu-decomposition-algorithm). Hence we may argue that pivoted Gaussian elimination cannot fail when $A$ is nonsingular.
We argue by contradiction. Suppose that the algorithm fails at step $k$. By definition of the algorithm, this means it has run successfully through steps $1, \ldots, k-1$, and hence
\[P^{(k-1)} A = \sum^{k-1} _ {i = 1} L _ i U _ i + \begin{pmatrix} 0 & 0 \\ 0 & A^{(k)} \end{pmatrix}\]where $A^{(k)} \in \mathbb R^{(n-k+1) \times (n-k+1)}$ is the active sub-block at the start of step $k$.
Failure at step $k$ means that there is no non-zero entry in the first column of $A^{(k)}$. Hence the first column of $A^{(k)}$ is a zero vector, and $A^{(k)}$ is rank deficient with $\text{rank}(A^{(k)}) \le n - k$. Each $L _ i U _ i$ is rank $1$, so $\sum^{k-1} _ {i=1} L _ i U _ i$ has rank at most $k - 1$. Combining via sub-additivity of the rank, this implies
\[\text{rank}(P^{(k-1)}A) \le (k-1) + (n-k) = n-1,\]but since permutation matrices are invertible, this implies $A$ is rank-deficient, and hence singular. This is a contradiction.
Therefore the algorithm completes at step $n$, producing $PA = LU$.
Symmetric matrices: the $LDL^\top$ form
@Prove that if $A \in \mathbb R^{n \times n}$ is symmetric and admits an LU factorisation $A = LU$, then it also admits an $LDL^\top$ factorisation: $A = M D M^\top$ for some unit lower triangular $M$ and diagonal $D$.
Setup: $A = LU$ with $L$ unit lower triangular, $U$ upper triangular; $A = A^\top$. Since $L$ is unit lower triangular, $\det L = 1$, so $L^{-1}$ (and $L^{-\top}$) exists.
Transpose $LU = A$ and use symmetry:
\[\begin{aligned} LU &= A = A^\top = U^\top L^\top &&(\star 1) \\ \Rightarrow U L^{-\top} &= L^{-1} U^\top &&(\star 2) \end{aligned}\]Now the opposite-triangularity trick:
- LHS $U L^{-\top}$ is a product of two upper triangular matrices, hence upper triangular.
- RHS $L^{-1} U^\top$ is a product of two lower triangular matrices, hence lower triangular.
A matrix that is both upper and lower triangular must be diagonal; call this common value $D$:
\[U L^{-\top} = L^{-1} U^\top = D. \qquad (\star 3)\]Rearranging $UL^{-\top} = D$ gives $U = D L^\top$, so
\[A = LU = L D L^\top = M D M^\top \qquad \text{with } M = L. \qquad \blacksquare\]- $(\star 1)$ Apply $\,\cdot^\top\,$ to $A = LU$: transposing reverses order and swaps upper/lower triangularity.
- $(\star 2)$ Pre-multiply by $L^{-1}$ and post-multiply by $L^{-\top}$.
- $(\star 3)$ Both sides equal yet one is upper triangular and the other lower triangular ⇒ both diagonal. This is the load-bearing trick of the proof.
Why useful: $LDL^\top$ exposes the symmetric structure in the factorisation itself ($M$ on both sides; only $D$ stored separately), so storage and the per-step cost halve compared to plain LU. It also handles indefinite symmetric $A$, with the signs of the $d _ i$ encoding the inertia of $A$. For $A \succ 0$ all $d _ i > 0$ and setting $R = D^{1/2} M^\top$ refines $LDL^\top$ into the Cholesky factorisation $A = R^\top R$ (∆cholesky-existence-conditions); for $A \preceq 0$ all $d _ i \le 0$ and setting $N = M (-D)^{1/2}$ gives $A = -NN^\top$ (NLA 2019 Q1©).
Stability analysis
@State a (nonexaminable) result about the backward stability of computing the LU factorisation of a matrix.
where $u$ is the unit roundoff and $c _ n$ is a polynomial in the dimension $n$.
The relevant ratio is $\rho _ n := \|U\| / \|A\|$, the growth factor. Empirically $\rho _ n = O(\sqrt n)$ for partial pivoting, but in the worst case it can be exponentially large.
We have the result that the LU factorisation of the matrix has the following backwards stability properties:
\[\frac{ \vert \vert \hat L \hat U - A \vert \vert }{ \vert \vert L \vert \vert \, \vert \vert U \vert \vert } = \epsilon\]
Under what conditions will $\hat x$ be a backward stable solution to $LUx = b$?
Worked example
Compute the pivoted LU factorisation $PA = LU$ of $A = \begin{bmatrix} 1 & 2 & 2 \\ 4 & 4 & 2 \\ 4 & 6 & 4 \end{bmatrix}$ by hand, using the rank-1 outer-product algorithm (∆pivoted-lu-decomposition-algorithm).
Step 1: column 1 is $(1, 4, 4)^\top$; the largest-magnitude entry is $4$ (row 2), so swap rows 1↔2. Pivot $\alpha = 4$; multipliers $\tfrac14$ (row 2), $1$ (row 3). Eliminating leaves the active block $\begin{bmatrix} 1 & 1.5 \\ 2 & 2 \end{bmatrix}$ (rows 2–3, cols 2–3).
Step 2: its first column is $(1, 2)^\top$; the larger magnitude is in the lower row, so swap the two working rows — this also swaps the stored column-1 multipliers of $L$ (they become $1$ and $\tfrac14$). Pivot $2$, multiplier $\tfrac12$; one more elimination finishes.
\[P A = \begin{bmatrix} 4 & 4 & 2 \\ 4 & 6 & 4 \\ 1 & 2 & 2 \end{bmatrix} = \underbrace{\begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 0.25 & 0.5 & 1 \end{bmatrix}} _ {L}\underbrace{\begin{bmatrix} 4 & 4 & 2 \\ 0 & 2 & 2 \\ 0 & 0 & 0.5 \end{bmatrix}} _ {U},\]with $P$ reordering $A$’s rows to $(2, 3, 1)$. $L$ is unit lower triangular with $ \vert L _ {ij} \vert \le 1$ thanks to partial pivoting. Check: row 3 of $LU = 0.25(4,4,2) + 0.5(0,2,2) + (0,0,0.5) = (1, 2, 2)$ = row 3 of $PA$ $\checkmark$.
Bite-sized
LU factorisation of $A \in \mathbb R^{n \times n}$ costs $\tfrac{2}{3} n^3$ flops. Once $A = LU$ (or $PA = LU$) is computed, solving $Ax = b$ takes only $O(n^2)$ additional flops via two triangular substitutions.
Give a $2 \times 2$ matrix that has no LU factorisation without pivoting, and explain why.
$A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$. The first LU step requires dividing the first column by the pivot $A _ {1,1} = 0$, so the elimination fails immediately. Partial pivoting solves this by swapping rows: $PA = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = I \cdot I = LU$.
The LU growth factor $\rho _ n$ quantifies element growth during elimination: the lecture takes $\rho _ n := \ \vert U\ \vert / \ \vert A\ \vert $ (equivalently, the standard textbook form is $\rho _ n = \max _ {i,j,k} \vert a^{(k)} _ {ij} \vert \big/ \max _ {i,j} \vert a _ {ij} \vert $, the largest entry occurring at any elimination step over the largest original entry). It is the factor governing LU backward error, $\ \vert \hat L \hat U - A\ \vert \lesssim c _ n\, \rho _ n\, u\, \ \vert A\ \vert $ (since $\ \vert L\ \vert = O(1)$ with partial pivoting), so LU + partial pivoting is backward stable exactly when $\rho _ n$ stays modest. Empirically $\rho _ n$ is $O(\sqrt n)$ with partial pivoting, but the worst case (Wilkinson’s matrix) reaches $2^{n-1}$.
Despite QR giving guaranteed backward stability for $Ax = b$, why is LU + partial pivoting still the default solver?
LU+PP is roughly $2\times$ faster than QR ($\tfrac{2}{3} n^3$ vs $\tfrac{4}{3} n^3$). The pathological matrices that cause exponential growth in $\|U\| / \|A\|$ are essentially never encountered in real-world dense problems, so LU+PP is reliably backward stable in practice. Resolving the gap between worst-case (exponential growth) and practice is one of the biggest open problems in NLA.
Once an LU factorisation $A = LU$ is in hand, solving $A x _ i = b _ i$ for $k$ different right-hand sides $b _ 1, \ldots, b _ k$ takes only $O(k n^2)$ additional flops — far cheaper than $k$ independent solves at $O(k n^3)$.
Strategy for proving any nonsingular matrix has a pivoted LU decomposition (∆pivoted-lu-decomposition-existence).
Show the algorithm cannot fail when $A$ is nonsingular, by contradiction.
- Suppose pivoted Gaussian elimination fails at step $k$. The first column of the active block $A^{(k)}$ must be all zero — so $A^{(k)}$ has rank $\le n - k$.
- Use the partial-decomposition $P^{(k-1)}A = \sum _ {i < k} L _ i U _ i + \mathrm{diag}(0, A^{(k)})$. Each $L _ i U _ i$ is rank 1, so the sum has rank $\le k - 1$.
- Add ranks: $\mathrm{rank}(P^{(k-1)} A) \le (k-1) + (n-k) = n - 1$. Since permutations are invertible, $A$ is rank-deficient.
- Contradicts the assumption that $A$ is nonsingular.