QR Decomposition Explained: Gram-Schmidt, Orthogonal Factors, and Numerical Stability

QR decomposition is one of the most important and widely used matrix factorizations in numerical linear algebra. It decomposes any matrix A into the product A = QR, where Q is an orthogonal matrix (its columns form an orthonormal set) and R is an upper triangular matrix. This factorization is the computational backbone of the QR algorithm for eigenvalues, the standard method for solving least squares problems, and a key ingredient in many signal processing and data analysis algorithms.

The elegance of QR decomposition lies in its construction through the Gram-Schmidt orthogonalization process. Just as you can orthogonalize a set of vectors to form a basis, you can decompose any matrix into an orthonormal factor and a triangular factor that records the "coefficients" of the original columns expressed in terms of the orthonormal basis. This connection between matrix factorization and basis construction gives QR decomposition its geometric meaning and computational power.

1. What Is QR Decomposition?

QR decomposition expresses a matrix A ∈ ℝ^m×n (with m ≥ n) as the product A = QR, where Q ∈ ℝ^m×n has orthonormal columns (Q^TQ = I_n) and R ∈ ℝ^n×n is upper triangular with positive diagonal entries. The columns of Q span the same subspace as the columns of A, and R encodes how the original columns relate to the orthonormal basis.

The Q Factor (Orthonormal Columns)

The matrix Q has orthonormal columns, meaning each column has unit length and any two distinct columns are orthogonal (their dot product is zero). Geometrically, Q^TQ = I means that applying Q preserves norms and angles — it is an isometry restricted to the column space. The columns of Q form an orthonormal basis for the column space of A (also called the range of A).

The R Factor (Upper Triangular)

The matrix R is upper triangular, meaning all entries below the main diagonal are zero. The diagonal entries R_ii are positive and represent the "lengths" of the projections of the original columns of A onto the orthonormal directions. The off-diagonal entries R_ij (for i < j) represent the components of column j of A that lie along the first i orthonormal directions. In essence, R records the Gram-Schmidt coefficients.

2. The Gram-Schmidt Process

The classical Gram-Schmidt process constructs the QR decomposition by orthogonalizing the columns of A one at a time. For each column a_j, you subtract its projections onto the previously computed orthonormal vectors q₁, q₂, ..., q_j−1, and then normalize the result. The projection coefficients become the entries of R.

Step-by-Step Algorithm

1. Initialize: Set q₁ = a₁ / ‖a₁‖ and R₁₁ = ‖a₁‖.
2. For j = 2 to n: Compute the projection of a_j onto each previous q_i: R_ij = q_i^T a_j for i = 1, ..., j−1.
3. Subtract projections: v = a_j − Σ_i=1^j−1 R_ij q_i.
4. Normalize: R_jj = ‖v‖ and q_j = v / R_jj.

Worked Example: 3×2 Matrix

Let A = [[1, 1], [1, 0], [0, 1]]. The columns are a₁ = [1, 1, 0]^T and a₂ = [1, 0, 1]^T.

Step 1: q₁ = a₁ / ‖a₁‖ = [1, 1, 0]^T / √2 ≈ [0.707, 0.707, 0]^T. R₁₁ = √2 ≈ 1.414.

Step 2: R₁₂ = q₁^T a₂ = 0.707 × 1 + 0.707 × 0 + 0 × 1 = 0.707.

Step 3: v₂ = a₂ − R₁₂ q₁ = [1, 0, 1] − 0.707 × [0.707, 0.707, 0] = [0.5, −0.5, 1].

Step 4: R₂₂ = ‖v₂‖ = √(0.25 + 0.25 + 1) = √1.5 ≈ 1.225. q₂ = v₂ / R₂₂ ≈ [0.408, −0.408, 0.816]^T.

The final result is Q = [[0.707, 0.408], [0.707, −0.408], [0, 0.816]] and R = [[1.414, 0.707], [0, 1.225]]. You can verify: QR = [[0.707, 0.408], [0.707, −0.408], [0, 0.816]] × [[1.414, 0.707], [0, 1.225]] ≈ [[1, 1], [1, 0], [0, 1]] = A.

3. Computing QR: Full Example with a 3×3 Matrix

For a full 3×3 matrix A, the QR decomposition proceeds identically — you orthogonalize three columns. Consider A = [[1, 1, 0], [1, 0, 1], [0, 1, 1]]. After applying Gram-Schmidt:

q₁ = [1/√2, 1/√2, 0]^T, q₂ = [1/√6, −1/√6, 2/√6]^T, q₃ = [1/√3, −1/√3, −1/√3]^T.

The R matrix is computed as R = Q^TA, which is the most direct way to obtain R once Q is known:

Note that R_ij = 0 for i > j (upper triangular), and all diagonal entries are positive — this is the standard QR decomposition with positive diagonal R. The factorization is unique when A has full column rank.

4. Properties and Applications

Solving Least Squares Problems

The primary application of QR decomposition is solving the least squares problem: minimize ‖Ax − b‖₂ for an overdetermined system (m > n). Using A = QR, the normal equations become:

The key advantage over the normal equations (x = (A^TA)⁻¹A^Tb) is numerical stability. Computing A^TA squares the condition number of the problem, amplifying rounding errors. QR decomposition avoids this squaring, making it the preferred method for ill-conditioned systems.

The QR Algorithm for Eigenvalues

The QR algorithm computes all eigenvalues of a matrix by repeatedly applying QR decomposition and reversing the factors. Starting with A₀ = A, compute A_k+1 = R_k Q_k (reversed order from A_k = Q_k R_k). The matrices A_k converge to an upper triangular (or Schur) form whose diagonal entries are the eigenvalues. With shifts and deflation, this algorithm is the standard method for computing eigenvalues in production software.

Solving Linear Systems

For a square invertible system Ax = b, QR decomposition gives x = R⁻¹Q^Tb. Back substitution on the triangular system Rx = Q^Tb solves for x in O(n²) time. While LU decomposition is generally faster for square systems, QR is preferred when the matrix is ill-conditioned or when you also need the condition number (which comes free from the R factor).

• Numerical stability: QR decomposition does not square the condition number, making it safer than normal equations.
• Minimum-norm solutions: For underdetermined systems, QR decomposition of A^T gives the minimum-norm solution.
• Rank determination: The diagonal entries of R reveal the numerical rank of A — small R_ii values indicate near-dependence.
• Updating solutions: QR decomposition can be efficiently updated when a row or column is added (Givens rotations), making it useful for recursive least squares.

5. Common Mistakes

• Using normal equations instead of QR: Computing (A^TA)⁻¹A^Tb squares the condition number. Always prefer QR decomposition for ill-conditioned problems.
• Confusing Q with Q^T: When solving Ax = b via QR, you compute Q^Tb (not Qb). Q^T rotates the right-hand side into the coordinate system defined by Q.
• Assuming R has positive diagonal: The standard QR decomposition requires positive diagonal entries in R. If you use Householder reflections, some R_ii may be negative — multiply the corresponding column of Q and row of R by −1 to fix this.
• Using Gram-Schmidt for large matrices: Classical Gram-Schmidt is numerically unstable for large or ill-conditioned matrices. Modified Gram-Schmidt or Householder reflections should be used instead.
• Forgetting QR is for m ≥ n: QR decomposition is defined for matrices with at least as many rows as columns. For m < n, decompose A^T instead.

6. Frequently Asked Questions

Is QR decomposition unique?

Yes, when A has full column rank and we require positive diagonal entries in R, the QR decomposition is unique. If A is rank-deficient, R may have zero diagonal entries and Q is not uniquely determined.

What is the difference between Gram-Schmidt, Householder, and Givens QR?

Gram-Schmidt orthogonalizes columns sequentially. Householder uses reflection matrices to zero out subdiagonal entries, working column by column from left to right. Givens rotations zero out individual entries using 2×2 rotation matrices. Householder is the most numerically stable and efficient for dense matrices. Givens rotations are preferred for sparse or banded matrices because they can target specific entries.

How does QR decomposition relate to LU decomposition?

Both are matrix factorizations used for solving linear systems. LU factors A = LU (lower × upper triangular) and is faster for square systems. QR factors A = QR (orthogonal × upper triangular) and is more stable for ill-conditioned and overdetermined systems. QR is preferred for least squares; LU is preferred for square systems with multiple right-hand sides.

Can QR decomposition be used for non-square matrices?

Yes. For an m × n matrix with m > n (more rows than columns), the "thin" QR gives Q as m × n and R as n × n. For m < n, decompose A^T to get a "thin" QR of the transpose. The full (economy) QR uses the complete orthonormal basis but is less efficient.

What is the computational complexity of QR decomposition?

For an m × n matrix with m ≥ n: Gram-Schmidt is O(mn²), Householder is O(2mn² − 2n³/3), and Givens rotations are O(3mn²). Householder is the standard for dense matrices. For sparse matrices, the complexity depends on the sparsity pattern and can be significantly lower.