Numerical Linear Algebra: LU and QR Decompositions

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

A x = 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.

A matrix factorization rewrites a matrix as a product of simpler, structured matrices whose actions are easy to invert or analyze. The LU decomposition expresses a square matrix as $A = LU$ where $L$ is unit lower triangular and $U$ is upper triangular; it is Gaussian elimination recorded as matrix multiplication. The QR decomposition writes $A = QR$ where $Q$ has orthonormal columns and $R$ is upper triangular; it orthogonalizes the columns of $A$ . Both turn the hard problem of solving general linear systems into a sequence of trivial triangular solves or orthogonal multiplications. A surprising fact: although Gaussian elimination dates to antiquity, the realization that it is a factorization $A = LU$ — and that this view explains its numerical behavior — came only with Turing's 1948 analysis.

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 $A x = 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.

Formal Statement

LU Decomposition

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

A = LU,

where $L$ is unit lower triangular ( $ℓ_{ii} = 1$ , $ℓ_{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.

⟳count degrees of freedom —

A

has

n^{2}

entries;

L

contributes

(2 n)

free entries,

U

contributes

(2 n + 1)

, summing to

n^{2}

QR Decomposition

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

A = QR,

where $Q \in R^{m \times m}$ is orthogonal ( $Q^{T} Q = I$ ) and $R \in R^{m \times n}$ is upper triangular. The reduced (thin) QR uses $\hat{Q} \in R^{m \times n}$ with orthonormal columns and $\hat{R} \in R^{n \times n}$ upper triangular, so $A = \hat{Q} \hat{R}$ .

⟳apply Gram–Schmidt to the columns

a_{1}, \dots, a_{n}

; the orthonormal outputs form

\hat{Q}

, and the projection coefficients fill

\hat{R}

Existence 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

P A = LU

— the LU with partial pivoting. (QR). Every

A \in R^{m \times n}

(

m \geq 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.

LU exists when leading minors are nonzero

Proceed by induction / direct elimination.

Step 1. The

(1, 1)

leading minor is

a_{11} \neq = 0

(it equals the first principal minor). Define multipliers

ℓ_{i 1} = a_{i 1} / a_{11}

for

i > 1

. Subtracting

ℓ_{i 1}

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

- ℓ_{i 1}

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,

M_{n - 1} \dots M_{2} M_{1} A = U,

⟳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

L = M_{1}^{- 1} M_{2}^{- 1} \dots M_{n - 1}^{- 1}

⟳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 $ℓ_{ik}$ (no interaction between stages), so $A = LU$ .

Step 5 (uniqueness). Suppose

L_{1} U_{1} = L_{2} U_{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}

■

QR 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}

v_{k} = a_{k} - j = 1 \sum k - 1 (q_{j}^{T} a_{k}) q_{j}, r_{j k} = q_{j}^{T} a_{k} .

⟳

v_{k}

is the component of

a_{k}

orthogonal to

span (q_{1}, \dots, q_{k - 1})

Step 3. Normalize:

q_{k} = v_{k} /∥ v_{k} ∥

r_{k k} = ∥ v_{k} ∥ > 0

(nonzero by full rank). Step 4. Rearranging,

a_{k} = \sum_{j = 1}^{k} r_{j k} q_{j}

, which in matrix form is exactly

A = \hat{Q} \hat{R}

with

\hat{R}

upper triangular. Orthonormality gives

\hat{Q}^{T} \hat{Q} = I

. Step 5 (uniqueness). If

\hat{Q}_{1} \hat{R}_{1} = \hat{Q}_{2} \hat{R}_{2}

then

\hat{Q}_{2}^{T} \hat{Q}_{1} = \hat{R}_{2} \hat{R}_{1}^{- 1}

is both orthogonal and upper triangular with positive diagonal, hence the identity.

■

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.

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.

Building Intuition

Think of solving $A x = 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 (

∥ Q x ∥ = ∥ 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

∥ A x - b ∥

becomes minimizing

∥ R x - Q^{T} b ∥

, and a rigid motion of the residual changes nothing.

LU and QR as compositions of simpler maps

Worked Examples

LU decomposition and solving a system

Given. Solve

A x = b

with

A = 248137139, b = 51432 .

Step 1 — eliminate column 1. Pivot

2

. Multipliers

ℓ_{21} = 4/2 = 2

ℓ_{31} = 8/2 = 4

A^{(2)} = 200113115 .

Step 2 — eliminate column 2. Pivot

1

. Multiplier

ℓ_{32} = 3/1 = 3

U = 200110112, L = 124013001 .

Step 3 — forward solve $L y = 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 $U x = y$ .

2 x_{3} = 0 \Rightarrow x_{3} = 0

;

x_{2} + x_{3} = 4 \Rightarrow x_{2} = 4

;

2 x_{1} + x_{2} + x_{3} = 5 \Rightarrow x_{1} = \frac{1}{2}

x = (\frac{1}{2}, 4, 0)^{T} (dimensionless solution vector)

Check:

A x = (5, 14, 32)^{T} = b .

✓

QR decomposition by Gram–Schmidt

Given.

A = 110101 .

Step 1.

a_{1} = (1, 1, 0)^{T}

∥ a_{1} ∥ = 2

, so

q_{1} = \frac{1}{2} (1, 1, 0)^{T}

and

r_{11} = 2

. Step 2.

a_{2} = (1, 0, 1)^{T}

. Projection coefficient

r_{12} = q_{1}^{T} a_{2} = \frac{1}{2} (1) = \frac{1}{2}

v_{2} = a_{2} - r_{12} q_{1} = (1, 0, 1)^{T} - \frac{1}{2} (1, 1, 0)^{T} = (\frac{1}{2}, - \frac{1}{2}, 1)^{T} .

Step 3.

∥ v_{2} ∥ = \frac{1}{4} + \frac{1}{4} + 1 = 3/2

, so

r_{22} = 3/2

and

q_{2} = \frac{2}{3} (\frac{1}{2}, - \frac{1}{2}, 1)^{T} = \frac{1}{6} (1, - 1, 2)^{T} .

Result.

\hat{Q} = 1/ 2 1/ 2 0 1/ 6 - 1/ 6 2/ 6, \hat{R} = (20 1/ 2 3/2) .

A = \hat{Q} \hat{R}, \hat{Q}^{T} \hat{Q} = I_{2}

Check:

\hat{Q} \hat{R}

reproduces

A

, and the columns of

\hat{Q}

are orthonormal. ✓

Common Misconceptions

Pivoting is not optional

A nonzero determinant does not guarantee a stable (or even existing) LU without pivoting. The matrix $(0111)$ 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 $∣ ℓ_{ij} ∣ \leq 1$ and is the practical default.

Why classical Gram–Schmidt fails

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

a_{1} = (1, ε, 0)^{T}, a_{2} = (1, 0, ε)^{T}

with

ε = 1 0^{- 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.

Connections

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 = L L^{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.

Real-World Applications

Circuit simulation (SPICE). Each Newton iteration solves $A x = b$ for nodal voltages; sparse LU with pivoting is the standard solver.
Structural / finite-element analysis. Stiffness matrices are symmetric positive-definite — Cholesky (a specialized LU) solves $K u = f$ for displacements.
Statistics / least squares. Fitting $y \approx X β$ uses QR of the design matrix $X$ to solve the normal equations stably without forming $X^{T} X$ .
Computer graphics & robotics. QR (via Householder/Givens) orthonormalizes drifting rotation matrices and solves inverse-kinematics linearizations.
Eigenvalue computation. The shifted QR algorithm repeatedly factors and recombines $A = QR \to R Q$ , the basis of LAPACK's dgeev.
Geodesy and GPS. Overdetermined positioning systems are solved by QR least squares.
Recommender systems / data compression. Low-rank approximations begin with QR (column-pivoted QR) to reveal numerical rank before refining with SVD.

Historical Development

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.

Subtleties & Limits

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

The Core Insight

Formal Statement

LU Decomposition

QR Decomposition

Existence and Uniqueness

LU exists when leading minors are nonzero

QR via Gram–Schmidt

Classical Gram–Schmidt is numerically unstable

Alan Turing | 1912–1954

Building Intuition

LU and QR as compositions of simpler maps

Worked Examples

LU decomposition and solving a system

QR decomposition by Gram–Schmidt

Common Misconceptions

Pivoting is not optional

Why classical Gram–Schmidt fails

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also

Before this

Where this leads