Mathematical foundations of Householder QR Factorization.
Householder QR factorization—introduced in 1958 by the American mathematician Alston Scott Householder in his seminal paper 'Unitary Triangularization of a Nonsymmetric Matrix'—is one of the most magnificent algorithmic monuments in numerical mathematics. In linear algebra, orthogonal transformations are the gold standard of numerical computation: because orthogonal matrices satisfy QᵀQ = I, multiplying a vector by Q preserves its Euclidean length (‖Qx‖₂ = ‖x‖₂) and inner products, meaning orthogonal transformations never amplify round-off errors. For decades, mathematicians relied on the Gram-Schmidt orthogonalization process to convert matrix columns into an orthonormal basis. However, in finite-precision floating-point arithmetic, Classical Gram-Schmidt (CGS) suffers from a catastrophic pathology: accumulated round-off errors in inner products cause the computed basis vectors to rapidly lose orthogonality, degrading to linear dependence. Alston Householder solved this by abandoning vector-by-vector projections in favor of elementary hyperplane reflections. A Householder reflector H = I − 2(vvᵀ)/(vᵀv) is a symmetric orthogonal matrix that reflects an entire vector across a hyperplane. By choosing v = x + sign(x₁)‖x‖₂ e₁, Householder demonstrated that an entire column below the diagonal can be annihilated in a single stroke, projecting the column directly onto the coordinate axis without subtractive cancellation. By applying n − 1 successive reflections, matrix A is reduced to upper-triangular form: H_{n-1} ... H₂ H₁ A = R, yielding the orthogonal factorization A = QR with Q = H₁ H₂ ... H_{n-1}. Because every individual reflector is perfectly orthogonal, the orthogonality of Q (‖QᵀQ − I‖_F ≤ O(ε_mach)) is unconditionally guaranteed down to machine precision. Although Householder QR requires (4/3)n³ operations—twice the cost of Gaussian elimination—its absolute stability makes it the undisputed method of choice for solving overdetermined least-squares problems, computing the QR algorithm for matrix eigenvalues, and stabilizing ill-conditioned systems. In this computational laboratory, students evaluate Householder QR factorization across dense matrices from n = 10 to n = 60, verify in-place Householder vector storage, measure orthogonality preservation, and extract the asymptotic scaling constant C = 4/3 ≈ 1.333 FLOPs/n³ via ordinary least-squares regression.
Viva voce questions & analytical theory
What is the geometric and algebraic derivation of a Householder reflector?
Introduced in 1958 by American mathematician Alston S. Householder, a Householder reflection (or elementary reflector) is an orthogonal matrix of the form H = I − 2(vvᵀ)/(vᵀv), where v is a non-zero vector orthogonal to a reflecting hyperplane. Geometrically, H reflects any vector x across the hyperplane: Hx = x − 2(vᵀx / vᵀv)v. To zero out all entries of a column vector x below the first component, one reflects x onto the coordinate axis vector e₁ = [1, 0, ..., 0]ᵀ: Hx = ±‖x‖₂ e₁. The vector v is chosen as v = x + sign(x₁)‖x‖₂ e₁, where matching the sign of x₁ prevents catastrophic subtractive cancellation.
Why is Householder QR vastly superior to Classical Gram-Schmidt (CGS) in numerical stability?
In Classical Gram-Schmidt (CGS), orthogonalizing column a_k against previously computed vectors q_1, ..., q_{k-1} subtracts projected components. In finite-precision arithmetic, small rounding errors in the inner products fail to subtract the components completely. As k increases, the computed vectors lose orthogonality catastrophically: ‖QᵀQ − I‖ can grow to O(1), meaning the columns of Q become linearly dependent! Modified Gram-Schmidt (MGS) improves stability, but still suffers from orthogonality loss on ill-conditioned matrices. Householder QR, by contrast, applies exact orthogonal reflections: because the product of orthogonal matrices is unconditionally orthogonal, ‖QᵀQ − I‖ ≤ O(ε_mach) is maintained to machine precision.
Why does Householder QR require (4/3)n³ FLOPs—twice the cost of Gaussian elimination?
At step k, reflecting an (n − k) vector requires: (1) Computing the inner product vᵀA (costing 2(n−k)² FLOPs); and (2) Subtracting the outer product βv(vᵀA) from A (costing 2(n−k)² FLOPs). Total per step is 4(n − k)² FLOPs—exactly twice the 2(n − k)² FLOPs of Gaussian elimination! Summing over all steps: ∑_{k=1}^{n-1} 4(n − k)² = 4 · (n³/3) = (4/3)n³ FLOPs. The extra factor of 2 is the price paid for unconditional numerical stability and perfect orthogonality.
How are Householder reflectors stored compactly in-place without forming n×n matrices?
A Householder matrix H_k is an n×n dense matrix, but it is completely determined by the single vector v_k of length n − k. Normalizing v_k so that its first component is 1, the remaining n − k − 1 components of v_k can be stored directly in the zeroed-out entries of column k below the diagonal of A! The diagonal and upper triangle of A store the factor R. Thus, the entire QR factorization is stored in-place in the original memory of A, requiring only an auxiliary vector of length n to store the scalar multipliers β_k.
How does Householder QR compare to Givens rotations?
Householder reflectors annihilate an entire column below the diagonal in one step, costing (4/3)n³ FLOPs. Givens rotations annihilate one single element at a time using 2×2 planar rotations, costing (4/3)n³ · (3/2) = 2n³ FLOPs (50% more expensive for dense matrices). However, Givens rotations are far superior for banded matrices, sparse matrices, and updating/downdating QR factorizations when rows or columns are added or deleted.
How does QR decomposition solve overdetermined linear least-squares problems Ax ≈ b?
Given an m×n matrix A with m > n, factoring A = QR partitions Q into [Q₁ Q₂] and R into [R₁; 0]. The least-squares problem min ‖Ax − b‖₂² = min ‖Qᵀ(Ax − b)‖₂² = min (‖R₁x − Q₁ᵀb‖₂² + ‖Q₂ᵀb‖₂²). The norm is minimized when the first term is zero: R₁x = Q₁ᵀb. Because R₁ is upper-triangular, x is solved trivially via back-substitution in O(n²) operations, with the residual norm given directly by ‖r‖₂ = ‖Q₂ᵀb‖₂.
How is Householder QR blocked in high-performance libraries like LAPACK (DGEQRF)?
The compact WY representation (Schreiber and Van Loan, 1989) expresses the product of r Householder reflectors as Q = I − W Yᵀ, where W and Y are n×r matrices. This transforms Level-2 BLAS rank-1 updates into Level-3 BLAS matrix-matrix multiplications (GEMM), enabling full saturation of CPU cache hierarchies and GPU memory bandwidth.
Algorithm precautions and stability rules
- Always include the sign function in v = x + sign(x₁)‖x‖₂ e₁ to eliminate subtractive cancellation.
- Never form the explicit n×n matrix H = I − 2vvᵀ/vᵀv; always apply it as a rank-1 update A ← A − βv(vᵀA).
- Store Householder vectors v_k in the vacated lower triangle of A to achieve zero-memory in-place factorization.
- For least-squares problems, accumulate Qᵀb simultaneously during factorization to avoid computing explicit Q.
- If column rank deficiency is suspected, use column pivoting (Businger-Golub QR with pivoting).
- Maintain 64-bit precision to preserve orthogonality down to ‖QᵀQ − I‖_F ≈ 10⁻¹⁶.
Theoretical limitations
- Requires (4/3)n³ operations—twice as expensive as Gaussian elimination ((2/3)n³).
- Computing the explicit orthogonal matrix Q requires an additional (4/3)n³ operations if needed.
- More complex memory access patterns than LU decomposition.