QR Decomposition

QR decomposition, also known as QR factorization, is a foundational procedure in numerical linear algebra and scientific computing that factors any real matrix A of dimension m × n (with m ≥ n) into the product A = QR, where Q is an orthogonal matrix satisfying Q^TQ = I and R is an upper triangular matrix. For complex-valued matrices, Q is a unitary matrix satisfying QQ = I, where Q denotes the conjugate transpose. Because orthogonal transformations preserve the Euclidean norm of vectors and leave the spectral condition number of a matrix unchanged, QR decomposition serves as the standard numerical backbone for solving overdetermined linear least-squares problems, constructing orthonormal bases of Krylov subspaces in methods such as the generalized minimal residual method, and computing matrix eigenvalues via the iterative QR algorithm.

Unlike LU decomposition, which relies on non-orthogonal row operations that can amplify rounding errors unless stabilized by pivoting, or Cholesky decomposition, which applies only to symmetric positive-definite matrices, QR factorization exists for every rectangular or square matrix. Computing Q via Householder reflections or Givens rotations achieves backward stability down to machine epsilon, avoiding the condition-number squaring of the normal equations A^TAx = A^Tb.

Let A ∈ ℝ^m×n be a rectangular matrix with m ≥ n. The full QR decomposition represents A as:

A = QR = [Q₁ Q₂] [R₁ ; 0]

Here, Q ∈ ℝ^m×m is a square orthogonal matrix partitioned into Q₁ ∈ ℝ^m×n, whose columns form an orthonormal basis for range(A) when A has full column rank, and Q₂ ∈ ℝ^m×(m−n), whose columns span the orthogonal complement null(A^T). The upper block R₁ ∈ ℝ^n×n is upper triangular, while the bottom block consists of m − n rows of zeros.

Dropping the zero rows of R and the columns of Q₂ yields the thin (or reduced / economy-size) QR decomposition A = Q₁R₁. If A has full column rank and the diagonal entries of R₁ are positive (r_ii > 0), the thin QR factorization is unique, and R₁ coincides with the upper triangular Cholesky factor of the Gram matrix A^TA = R₁^TR₁.

Three algorithmic families construct the QR factorization of a dense or sparse matrix, trading arithmetic operation counts against numerical orthogonality and parallelizability.

Gram–Schmidt Orthogonalization and Modified Gram–Schmidt

The classical Gram–Schmidt (CGS) process builds the columns q₁, …, q_n of Q₁ sequentially by subtracting from each column a_j of A its orthogonal projections onto the previously computed vectors q₁, …, q_j−1:

v_j = a_j − ∑_i=1^j−1 (q_i^Ta_j) q_i, r_ij = q_i^Ta_j, r_jj = ‖v_j‖₂, q_j = v_j / r_jj.

Although algebraically transparent and requiring 2mn² floating-point operations (flops), classical Gram–Schmidt suffers from substantial loss of orthogonality in finite-precision arithmetic whenever A is ill-conditioned: subtracting nearly parallel vectors triggers numerical cancellation, causing ‖Q₁^TQ₁ − I‖₂ to grow without bound. Modified Gram–Schmidt (MGS) reorganizes the inner loop as a sequence of rank-one projections against each freshly normalized vector q_i. MGS performs the same 2mn² flops as CGS but bounds the departure from orthogonality by O(ε_mach κ₂(A)), where κ₂(A) = σ_max(A)/σ_min(A) is the 2-norm condition number. When machine-precision orthogonality is required inside Krylov solvers such as the Lanczos algorithm or GMRES, MGS is paired with a second reorthogonalization pass ("twice is enough," as proved by Kahan and Parlett).

Householder Reflections

Introduced by Alston Scott Householder in 1958, Householder QR operates from the left by applying a sequence of n orthogonal elementary reflectors H₁, H₂, …, H_n to zero out all subdiagonal entries of A one column at a time. Each Householder matrix has the form:

H_k = I − 2 v_kv_k^T / (v_k^Tv_k),

where the reflection vector v_k = x + sign(x₁)‖x‖₂e₁ maps the active subcolumn x onto a multiple of the first unit coordinate vector e₁ without numerical cancellation. After n steps, H_n⋯H₁A = R, and Q = H₁⋯H_n is stored implicitly below the diagonal of R. Householder QR requires 2mn² − (2/3)n³ flops and achieves unconditional backward stability: ‖Q̂^TQ̂ − I‖₂ = O(ε_mach) independent of κ₂(A). Modern libraries such as LAPACK aggregate groups of b reflectors into a block WY form I − VTV^T to execute Level-3 BLAS matrix-matrix multiplications.

Givens Rotations

Whereas a Householder reflector zeros an entire subcolumn simultaneously, a Givens rotation G(i, k, θ) is a plane rotation matrix that differs from the identity only in rows and columns i and k, with cosine c = cos θ on the diagonal and sine s = ±sin θ on the off-diagonal. Choosing c = a / √(a² + b²) and s = −b / √(a² + b²) selectively zeroes a single entry b. For a dense m × n matrix, Givens QR requires 3mn² − n³ flops, but for upper Hessenberg matrices in Arnoldi iteration or tridiagonal matrices in the Lanczos algorithm, Givens rotations reduce the matrix to upper triangular form in only O(n²) or O(n) operations.

Linear Least Squares and Rank-Revealing QR

Given an overdetermined linear system Ax ≈ b with m > n, minimizing the Euclidean residual ‖Ax − b‖₂² via the normal equations A^TAx = A^Tb squares the condition number to κ₂(A^TA) = [κ₂(A)]². Using the orthogonal invariance of the 2-norm, ‖Ax − b‖₂² = ‖R₁x − Q₁^Tb‖₂² + ‖Q₂^Tb‖₂². The least-squares minimizer x̂ is therefore obtained by solving the n × n upper triangular system R₁x̂ = Q₁^Tb by back-substitution, retaining a condition sensitivity proportional to κ₂(A). When A is rank-deficient, QR with column pivoting (AP = QR, where P is a permutation matrix selected at step k to swap in the column with the largest remaining residual norm) exposes the numerical rank of A along the diagonal of R and underpins subset selection and principal component analysis preprocessing.

The QR Algorithm for Eigenvalue Problems

Developed independently by John G. F. Francis and Vera N. Kublanovskaya in 1961–1962, the QR algorithm computes all eigenvalues and eigenvectors of a square matrix A ∈ ℝ^n×n. After an initial orthogonal reduction A₀ = U^TAU to upper Hessenberg form, each iteration factors the shifted matrix A_k − μ_kI = Q_kR_k and forms the reverse product A_k+1 = R_kQ_k + μ_kI = Q_k^TA_kQ_k. With Wilkinson shifts μ_k, the subdiagonal entries decay at a quadratic or cubic rate, deflating A_k into real Schur form.

Why is QR decomposition preferred over normal equations for least squares?

Forming the Gram matrix A^TA squares the condition number of A: if κ₂(A) = 10⁸, then κ₂(A^TA) = 10¹⁶, exhausting IEEE 754 double-precision mantissa bits and causing Cholesky decomposition of A^TA to lose accuracy. Householder QR works directly on A without squaring κ₂(A).

Do we need to form the full matrix Q explicitly?

Rarely. In Householder QR, applying Q^Tb using the stored reflection vectors v₁, …, v_n takes only 2mn − n² flops, whereas multiplying out Q ∈ ℝ^m×m explicitly requires 4(m²n − mn² + n³/3) flops and m² memory words.

1. Householder, A. S. (1958). Unitary triangularization of a nonsymmetric matrix. Journal of the ACM, 5(4), 339–342. https://doi.org/10.1145/320941.320947
2. Francis, J. G. F. (1961). The QR transformation: A unitary analogue to the LR transformation—Part 1. The Computer Journal, 4(3), 265–271. https://doi.org/10.1093/comjnl/4.3.265
3. Björck, Å. (1967). Solving linear least squares problems by Gram–Schmidt orthogonalization. BIT Numerical Mathematics, 7(1), 1–21. https://doi.org/10.1007/BF01934122
4. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press. https://doi.org/10.56021/9781421407944