Numerical Linear Algebra: LU and QR Decompositions

⚡ Key Takeaways

LU and QR decompositions factor a matrix $A$ into structured pieces — $A = LU$ (lower × upper triangular, via Gaussian elimination) or $A = QR$ (orthogonal × upper triangular, via Gram–Schmidt or Householder reflections). LU solves square linear systems $Ax=b$ efficiently; QR solves least-squares and underpins eigenvalue algorithms. These are the two workhorse factorizations of numerical linear algebra, used everywhere from circuit simulation to statistics. Formalized in the 20th century (Turing, Householder, Givens) on foundations from Gauss.

💡 The Core Insight

Every elementary row operation is left-multiplication by a simple matrix. LU collects all the elimination steps into one lower-triangular matrix $L$ ; QR replaces "elimination" by "orthogonalization," collecting all the rotations/reflections into an orthogonal $Q$ . In both cases the right factor $R$ (or $U$ ) is upper triangular, so the original problem $Ax=b$ becomes two triangular systems, each solvable in $O(n^2)$ by back/forward substitution. Triangularity is the whole game: it is the structure that makes solving cheap.

DefinitionLU Decomposition

Let $A \in \mathbb{R}^{n\times n}$ . An LU decomposition is a factorization

1

\quad A = LU,

where $L$ is unit lower triangular ( $\ell_{ii}=1$ , $\ell_{ij}=0$ for $j>i$ ) and $U$ is upper triangular ( $u_{ij}=0$ for $i>j$ ). The entries below the diagonal of $L$ are the multipliers used in Gaussian elimination; $U$ is the row-echelon (upper) result. Derivation hint: count degrees of freedom — $A$ has $n^2$ entries; $L$ contributes $\binom{n}{2}$ free entries, $U$ contributes $\binom{n+1}{2}$ , summing to $n^2$ .

DefinitionQR Decomposition

Let $A \in \mathbb{R}^{m\times n}$ with $m \ge n$ . A QR decomposition is

2

\quad A = QR,

where $Q \in \mathbb{R}^{m\times m}$ is orthogonal ( $Q^{\mathsf T}Q = I$ ) and $R \in \mathbb{R}^{m\times n}$ is upper triangular. The reduced (thin) QR uses $\hat Q \in \mathbb{R}^{m\times n}$ with orthonormal columns and $\hat R \in \mathbb{R}^{n\times n}$ upper triangular, so $A = \hat Q \hat R$ . Derivation hint: apply Gram–Schmidt to the columns $a_1,\dots,a_n$ ; the orthonormal outputs form $\hat Q$ , and the projection coefficients fill $\hat R$ .

TheoremExistence and Uniqueness

(LU). If all leading principal minors of $A$ are nonzero, then $A$ has a unique LU decomposition with $L$ unit lower triangular. In general (any nonsingular $A$ ), a row permutation exists so that $PA = LU$ — the LU with partial pivoting.

(QR). Every $A \in \mathbb{R}^{m\times n}$ ( $m\ge n$ ) has a QR decomposition. If $A$ has full column rank, the reduced decomposition $A=\hat Q\hat R$ is unique once we require the diagonal entries of $\hat R$ to be positive.

ProofLU exists when leading minors are nonzero

Proceed by induction / direct elimination.

Step 1. The $(1,1)$ leading minor is $a_{11}\ne 0$ (it equals the first principal minor). Define multipliers $\ell_{i1}=a_{i1}/a_{11}$ for $i>1$ . Subtracting $\ell_{i1}$ times row 1 from row $i$ zeros the first column below the diagonal. This is left-multiplication by an elementary lower-triangular matrix $M_1$ with $-\ell_{i1}$ in its first column.

Step 2. The resulting matrix $A^{(2)}=M_1 A$ has the same trailing principal minors structure; its $(2,2)$ pivot is nonzero exactly because the second leading minor of $A$ is nonzero (the minor equals $a_{11}\cdot$ pivot $_2$ ). Repeat: form $M_2$ , then $M_3$ , ..., $M_{n-1}$ .

Step 3. After $n-1$ stages,

3

\quad M_{n-1}\cdots M_2 M_1 A = U,

Derivation hint: each $M_k$ clears column $k$ below the pivot; the product is upper-triangularization.

Step 4. Each $M_k$ is unit lower triangular, so its inverse is too, and

4

\quad L = M_1^{-1}M_2^{-1}\cdots M_{n-1}^{-1}

Derivation hint: inverse of a product reverses order; each $M_k^{-1}$ just flips the sign of its off-diagonal multipliers. is unit lower triangular. A clean miracle: $L$ 's subdiagonal entries are exactly the multipliers $\ell_{ik}$ (no interaction between stages), so $A = LU$ .

Step 5 (uniqueness). Suppose $L_1U_1=L_2U_2$ . Then $L_2^{-1}L_1 = U_2 U_1^{-1}$ . The left side is unit lower triangular, the right is upper triangular; a matrix that is both is the identity. Hence $L_1=L_2,\ U_1=U_2$ . $\blacksquare$

∎

ProofQR via Gram–Schmidt

Let $a_1,\dots,a_n$ be the columns of full-rank $A$ .

Step 1. Set $q_1 = a_1/\|a_1\|$ , with $r_{11}=\|a_1\|$ .

Step 2. For $k=2,\dots,n$ , subtract projections onto previous $q_j$ :

5

\quad v_k = a_k - \sum_{j=1}^{k-1} (q_j^{\mathsf T}a_k)\,q_j,\qquad r_{jk}=q_j^{\mathsf T}a_k.

Derivation hint: $v_k$ is the component of $a_k$ orthogonal to $\mathrm{span}(q_1,\dots,q_{k-1})$ .

Step 3. Normalize: $q_k = v_k/\|v_k\|$ , $r_{kk}=\|v_k\|>0$ (nonzero by full rank).

Step 4. Rearranging, $a_k = \sum_{j=1}^{k} r_{jk}q_j$ , which in matrix form is exactly $A=\hat Q\hat R$ with $\hat R$ upper triangular. Orthonormality gives $\hat Q^{\mathsf T}\hat Q=I$ .

Step 5 (uniqueness). If $\hat Q_1\hat R_1 = \hat Q_2\hat R_2$ then $\hat Q_2^{\mathsf T}\hat Q_1 = \hat R_2\hat R_1^{-1}$ is both orthogonal and upper triangular with positive diagonal, hence the identity. $\blacksquare$

∎

⚠️ Classical Gram–Schmidt is numerically unstable

The algorithm above (classical Gram–Schmidt, CGS) loses orthogonality catastrophically in finite precision when columns are nearly dependent. In practice use Modified Gram–Schmidt (MGS) — subtract projections one at a time, updating $a_k$ immediately — or, better, Householder QR. MGS and Householder are mathematically equivalent to CGS but vastly more accurate.

AT

Alan Turing

1912–1954

Turing's 1948 paper "Rounding-off errors in matrix processes" gave the first rigorous error analysis of Gaussian elimination and explicitly framed it as the matrix factorization $A=LU$ (he wrote it in terms of triangular factors). He introduced the condition number concept for matrices, quantifying how input errors amplify. This reconceptualization — elimination as factorization — is what allowed the subsequent stability theory of Wilkinson and others. Turing thus turned a centuries-old hand computation into an object of modern numerical mathematics.

🧠 Intuition

Think of solving $Ax=b$ as undoing a complicated transformation. LU says: $A$ stretches space in a way that decomposes into a "shear" ( $L$ ) followed by an "axis-aligned scaling with shear" ( $U$ ). To invert, undo $U$ then $L$ — each triangular, each cheap.

QR says: $A$ 's action equals "first rotate/reflect ( $Q$ ), then upper-triangular shear ( $R$ )." Because $Q$ is just a rigid motion ( $\|Qx\|=\|x\|$ ), it preserves lengths and angles — it never amplifies error. That is why QR is so stable, and why it shines for least-squares: minimizing $\|Ax-b\|$ becomes minimizing $\|Rx - Q^{\mathsf T}b\|$ , and a rigid motion of the residual changes nothing.

Figure

LU and QR as compositions of simpler maps

ExampleLU decomposition and solving a system

Given. Solve $Ax=b$ with

A=\begin{pmatrix}2&1&1\\4&3&3\\8&7&9\end{pmatrix},\qquad b=\begin{pmatrix}5\\14\\32\end{pmatrix}.

Step 1 — eliminate column 1. Pivot $2$ . Multipliers $\ell_{21}=4/2=2$ , $\ell_{31}=8/2=4$ .

A^{(2)}=\begin{pmatrix}2&1&1\\0&1&1\\0&3&5\end{pmatrix}.

Step 2 — eliminate column 2. Pivot $1$ . Multiplier $\ell_{32}=3/1=3$ .

U=\begin{pmatrix}2&1&1\\0&1&1\\0&0&2\end{pmatrix},\qquad L=\begin{pmatrix}1&0&0\\2&1&0\\4&3&1\end{pmatrix}.

Step 3 — forward solve $Ly=b$ . $y_1=5$ ; $\;2(5)+y_2=14\Rightarrow y_2=4$ ; $\;4(5)+3(4)+y_3=32\Rightarrow y_3=0$ .

Step 4 — back solve $Ux=y$ . $2x_3=0\Rightarrow x_3=0$ ; $\;x_2 + x_3 = 4\Rightarrow x_2=4$ ; $\;2x_1+x_2+x_3=5\Rightarrow x_1=\tfrac12$ .

\boxed{x = (\tfrac12,\;4,\;0)^{\mathsf T}}\quad\text{(dimensionless solution vector)}

Check: $A x = (5,14,32)^{\mathsf T} = b.$ ✓

ExampleQR decomposition by Gram–Schmidt

Given.

A=\begin{pmatrix}1&1\\1&0\\0&1\end{pmatrix}.

Step 1. $a_1=(1,1,0)^{\mathsf T}$ , $\|a_1\|=\sqrt2$ , so $q_1=\tfrac1{\sqrt2}(1,1,0)^{\mathsf T}$ and $r_{11}=\sqrt2$ .

Step 2. $a_2=(1,0,1)^{\mathsf T}$ . Projection coefficient $r_{12}=q_1^{\mathsf T}a_2 = \tfrac1{\sqrt2}(1)=\tfrac1{\sqrt2}$ .

v_2 = a_2 - r_{12}q_1 = (1,0,1)^{\mathsf T} - \tfrac12(1,1,0)^{\mathsf T} = (\tfrac12,-\tfrac12,1)^{\mathsf T}.

Step 3. $\|v_2\| = \sqrt{\tfrac14+\tfrac14+1}=\sqrt{3/2}$ , so $r_{22}=\sqrt{3/2}$ and $q_2 = \sqrt{\tfrac{2}{3}}(\tfrac12,-\tfrac12,1)^{\mathsf T} = \tfrac1{\sqrt6}(1,-1,2)^{\mathsf T}.$

Result.

\hat Q=\begin{pmatrix}1/\sqrt2 & 1/\sqrt6\\[2pt]1/\sqrt2 & -1/\sqrt6\\[2pt]0 & 2/\sqrt6\end{pmatrix},\qquad \hat R=\begin{pmatrix}\sqrt2 & 1/\sqrt2\\[2pt]0 & \sqrt{3/2}\end{pmatrix}.

\boxed{A = \hat Q\hat R,\quad \hat Q^{\mathsf T}\hat Q = I_2}

Check: $\hat Q\hat R$ reproduces $A$ , and the columns of $\hat Q$ are orthonormal. ✓

⚠️ Pivoting is not optional

A nonzero determinant does not guarantee a stable (or even existing) LU without pivoting. The matrix $\begin{pmatrix}0&1\\1&1\end{pmatrix}$ has no LU factorization at all (its $(1,1)$ pivot is $0$ ), yet it is nonsingular. Even when a pivot is merely small, dividing by it amplifies rounding errors. Partial pivoting — swapping in the largest-magnitude entry of the active column as pivot — keeps all multipliers $|\ell_{ij}|\le 1$ and is the practical default.

CounterexampleWhy classical Gram–Schmidt fails

Take, in 3-digit arithmetic, the near-collinear columns

a_1=(1,\varepsilon,0)^{\mathsf T},\ a_2=(1,0,\varepsilon)^{\mathsf T}

with $\varepsilon = 10^{-3}$ . After normalizing $a_1$ and subtracting its projection from $a_2$ , the tiny orthogonal component is computed as the difference of two nearly-equal large numbers — catastrophic cancellation destroys it. The resulting $q_2$ can fail to be orthogonal to $q_1$ to working precision. Householder QR avoids this entirely because it never forms such differences; it works with norm-preserving reflections.

connection

LU is the matrix incarnation of Gaussian Elimination, and its determinant follows instantly: $\det A = \det L\cdot\det U = \prod u_{ii}$ . The symmetric positive-definite specialization is the Cholesky Decomposition $A=LL^{\mathsf T}$ , half the cost of LU. QR is built from Orthogonal Projection and the Gram–Schmidt Process, and it is the engine of the QR Algorithm for Eigenvalues and Eigenvectors. Both decompositions are quantitatively governed by the Condition Number, and QR's stability links to the Singular Value Decomposition, the most robust (but costliest) factorization. The least-squares use of QR is the computational core of Linear Regression.

🔬 Where this appears

⏳ History

- c. 179 CE — The Chinese Nine Chapters on the Mathematical Art describes elimination-style solution of linear systems.

- 1809–1810 — Gauss systematizes elimination while computing the orbit of Ceres (least squares), giving the method its modern name.

- 1888 — Wilhelm Jordan popularizes the Gauss–Jordan variant.

- 1907 — Erhard Schmidt formalizes the orthogonalization process (with roots in Gram, 1883).

- 1948 — Alan Turing publishes the first rounding-error analysis, framing elimination as the factorization

A=LU

and introducing matrix condition numbers.

- 1948 — Alston Householder and others develop systematic matrix-computation theory.

- 1958 — Householder introduces reflection-based QR, the stable standard.

- 1958–1961 — John Francis and Vera Kublanovskaya independently invent the QR algorithm for eigenvalues — one of the top algorithms of the 20th century.

- 1965 — James Wilkinson's The Algebraic Eigenvalue Problem consolidates backward-error analysis of LU and QR.

open-problem

Even partial pivoting can, in pathological cases, exhibit exponential element growth (the Wilkinson matrix), making LU unstable in principle — yet such matrices essentially never arise in practice, and why average-case growth is benign remains incompletely understood (the subject of "smoothed analysis," Sankar–Spielman–Teng, 2006). For sparse matrices, choosing a pivot order that minimizes fill-in while preserving stability is NP-hard, driving ongoing research in heuristic orderings. Communication cost, not arithmetic, now dominates on modern hardware, spurring "communication-avoiding" LU and QR variants (CALU, TSQR) whose stability–performance trade-offs are still being charted.

Numerical Linear Algebra: LU and QR Decompositions

Formal Statement

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also