QR Decomposition

What is QR Decomposition?

Every matrix A with linearly independent columns can be factored into a product A = QR, where Q has orthonormal columns and R is upper triangular with positive diagonal entries. This factorization exists for any m × n matrix with m ≥ n and rank n — that is, for any matrix whose columns are linearly independent. Unlike LU decomposition, which requires a square matrix, QR works for tall rectangular matrices too, making it the natural tool for overdetermined least-squares problems.

The orthogonal factor Q captures the geometry of A's column space: its columns form an orthonormal basis for that space. The upper triangular factor R captures the "scaling and shearing" needed to go from the standard orthonormal basis to A's actual columns. Together they encode everything about A in a form that is exceptionally well-conditioned for numerical computation.

QR decomposition appears in three major applications: solving least-squares problems, computing eigenvalues via the QR algorithm, and orthogonalizing a set of vectors. In all three, the orthogonality of Q is the key — it means Q⁻¹ = Qᵀ, so matrix-vector products with Q or Qᵀ are perfectly conditioned and can never amplify errors.

The Q and R Factors

For an m × n matrix A with m ≥ n and linearly independent columns, the QR factorization writes A = QR where:

• Q is an m × n matrix with orthonormal columns: QᵀQ = I_n. The columns of Q are an orthonormal basis for the column space of A.
• R is an n × n upper triangular matrix with positive diagonal entries. R encodes how to reconstruct A's columns from Q's columns.

There is also a "full" QR factorization where Q is extended to a full m × m orthogonal matrix by appending m − n additional orthonormal columns (spanning the left null space of A), and R is extended to an m × n matrix by padding with rows of zeros. Both forms are useful, but the "thin" or "economy" QR (with Q being m × n) is the one most commonly used in practice.

Computing QR: Gram-Schmidt Process

The conceptually simplest way to compute QR is the Gram-Schmidt process, which directly constructs the orthonormal columns of Q one at a time. Recall from Module 5 that Gram-Schmidt takes a set of linearly independent vectors and produces an orthonormal basis for their span. Applied to the columns a₁, a₂, …, aₙ of A, it produces the columns q₁, q₂, …, qₙ of Q.

The process is iterative. At step k, we take aₖ, subtract off its projections onto all previously constructed basis vectors q₁, …, q_{k−1}, and normalize the remainder:

Reading off the entries of R from the Gram-Schmidt steps: the diagonal entry rₖₖ = ‖ẽₖ‖ is the norm of the residual at step k, and the off-diagonal entries rᵢₖ = qᵢᵀaₖ are the projection coefficients. The column a₁, a₂, …, aₙ can be recovered from q₁, q₂, …, qₙ via the upper triangular matrix R, confirming that A = QR.

Classical vs. Modified Gram-Schmidt

The classical Gram-Schmidt process as described above is mathematically correct but numerically unstable when run in floating-point arithmetic: rounding errors in the projection steps can cause the computed vectors to lose orthogonality, especially when columns are nearly parallel. The modified Gram-Schmidt reorders the projections to subtract each projection as soon as it is computed, rather than all at once. It produces the same result in exact arithmetic but is much more numerically stable, because errors introduced in one projection do not propagate into subsequent ones.

For most practical work, the modified Gram-Schmidt process offers a good balance of simplicity and numerical quality. But when very high accuracy is required — as in eigenvalue computation — Householder reflections provide even stronger guarantees.

Computing QR: Householder Reflections

The Householder reflection approach is the method used in production numerical libraries (LAPACK, NumPy, MATLAB). Instead of building Q column by column as in Gram-Schmidt, it builds Q as a product of elementary orthogonal matrices applied from the left to A.

A Householder reflector H is an orthogonal, symmetric matrix of the form H = I − 2vvᵀ / (vᵀv), where v is called the Householder vector. Geometrically, H reflects any vector across the hyperplane perpendicular to v. By choosing v appropriately at each step, we can zero out all entries below the diagonal in a given column of A.

The Householder approach has two key advantages over Gram-Schmidt. First, it is numerically stable: each Householder reflector is an exact orthogonal transformation, so no orthogonality is lost due to rounding. Second, it is efficient: working with vvᵀ implicitly (without forming the full n×n matrix H) makes each step cost O(mn) rather than O(mn²). In practice, LAPACK stores the Householder vectors and applies them lazily, making the overall factorization cost about 2mn² − 2n³/3 flops.

Computing QR: Givens Rotations

A third approach builds Q as a product of Givens rotations — elementary orthogonal matrices that rotate two coordinates at a time. A Givens rotation G(i, j, θ) affects only rows i and j, rotating them by angle θ so that a target entry becomes zero. This is the surgical scalpel version: while Householder reflections zero out an entire column below the diagonal in one step, Givens rotations zero out one entry at a time.

Givens rotations shine when A is sparse: if most entries of A are already zero, it is wasteful to apply a full Householder reflector (which modifies the entire column). Givens rotations touch only the two rows involved in each zero-creation, preserving sparsity much better. They are the method of choice for banded matrices and for updating an existing QR factorization when one row of A changes.

Application: Least-Squares via QR

The primary motivation for QR decomposition in applications is solving overdetermined least-squares problems: given an m × n matrix A with m > n and a vector b ∈ ℝᵐ, find x ∈ ℝⁿ that minimizes ‖Ax − b‖².

The standard approach — the normal equations AᵀAx = Aᵀb — computes AᵀA, which doubles the condition number of the problem (κ(AᵀA) = κ(A)²) and can be catastrophically inaccurate when A is ill-conditioned. The QR approach avoids forming AᵀA entirely.

The QR approach costs 2mn² − 2n³/3 flops to factorize (compared to mn² for forming AᵀA and n³/3 for its Cholesky factorization), but the reduction in condition number from κ(A)² to κ(A) is usually worth the extra cost. When A is nearly rank-deficient or the data is noisy, the QR approach can give correct answers where the normal equations produce garbage.

Application: The QR Algorithm for Eigenvalues

QR decomposition is also the engine of the most widely used algorithm for computing all eigenvalues of a matrix: the QR algorithm (not to be confused with QR decomposition itself). The basic iteration is remarkably simple:

1. Start with A₀ = A.
2. At each step k: compute the QR factorization Aₖ = QₖRₖ.
3. Form Aₖ₊₁ = RₖQₖ (reverse the product).
4. Repeat until Aₖ converges to an upper triangular (or nearly triangular) matrix.

The QR algorithm with shifts is the method of choice for dense eigenvalue problems and is implemented in LAPACK's dgeev (general matrices) and dsyev (symmetric matrices). For symmetric matrices, the sequence converges to a diagonal matrix, directly revealing all eigenvalues. The QR algorithm requires O(n³) work per iteration but typically converges in O(1) iterations per eigenvalue after deflation, giving an overall cost of O(n³).

Numerical Stability: Why QR Beats the Normal Equations

The condition number of a linear system measures how much the solution can change relative to perturbations in the data. For the least-squares problem min ‖Ax − b‖², the relevant condition number is κ(A) — the ratio of the largest to smallest singular value of A. Solving via the normal equations introduces AᵀA, whose condition number is κ(A)². This squaring can be catastrophic: if κ(A) = 10⁶ (a mild degree of ill-conditioning for real data), then κ(AᵀA) = 10¹², and about 12 digits of accuracy are lost before you even start solving.

QR decomposition avoids this squaring. Because Q is orthogonal, the system Rx = Qᵀb has condition number κ(R) = κ(A), not κ(A)². The transformation Qᵀb cannot make anything worse — orthogonal transformations preserve all lengths and angles exactly. This is the fundamental reason why QR decomposition is the preferred method for least-squares computation in any numerically sensitive application.