NumPy - QR Decomposition | Application Architect

Key Insights

QR decomposition factorizes any matrix A into an orthogonal matrix Q and upper triangular matrix R, enabling stable solutions for least squares problems and eigenvalue computations
NumPy’s np.linalg.qr() provides multiple modes including reduced QR for efficiency and complete QR for theoretical completeness, with pivoting available through SciPy for rank-deficient matrices
Understanding when to use QR over other decompositions (LU, SVD, Cholesky) significantly impacts numerical stability and performance in real-world applications like regression analysis and signal processing

Understanding QR Decomposition

QR decomposition breaks down an m×n matrix A into two components: Q (an orthogonal matrix) and R (an upper triangular matrix) such that A = QR. The orthogonal property of Q means Q^T Q = I, which preserves vector lengths and angles during transformations—a critical property for numerical stability.

import numpy as np

# Create a sample matrix
A = np.array([[12, -51, 4],
              [6, 167, -68],
              [-4, 24, -41]], dtype=float)

# Perform QR decomposition
Q, R = np.linalg.qr(A)

print("Original matrix A:")
print(A)
print("\nOrthogonal matrix Q:")
print(Q)
print("\nUpper triangular matrix R:")
print(R)
print("\nVerification A = QR:")
print(np.allclose(A, Q @ R))
print("\nVerification Q^T Q = I:")
print(np.allclose(Q.T @ Q, np.eye(3)))

The orthogonality of Q ensures that Q.T @ Q produces an identity matrix, while R contains all the information about the original matrix’s structure in a triangular form.

QR Modes: Complete vs Reduced

NumPy offers different QR modes depending on your needs. The reduced mode (default) produces economical decomposition where Q is m×n and R is n×n. The complete mode generates full-sized matrices where Q is m×m and R is m×n.

# Rectangular matrix (more rows than columns)
A = np.array([[1, 2],
              [3, 4],
              [5, 6],
              [7, 8]], dtype=float)

# Reduced QR (default, mode='reduced')
Q_reduced, R_reduced = np.linalg.qr(A, mode='reduced')
print(f"Reduced Q shape: {Q_reduced.shape}")  # (4, 2)
print(f"Reduced R shape: {R_reduced.shape}")  # (2, 2)

# Complete QR
Q_complete, R_complete = np.linalg.qr(A, mode='complete')
print(f"Complete Q shape: {Q_complete.shape}")  # (4, 4)
print(f"Complete R shape: {R_complete.shape}")  # (4, 2)

# Both reconstruct A correctly
print(f"Reduced reconstruction: {np.allclose(A, Q_reduced @ R_reduced)}")
print(f"Complete reconstruction: {np.allclose(A, Q_complete @ R_complete)}")

For most applications, reduced mode provides better memory efficiency and computational performance. Use complete mode when you need the full orthogonal basis for the space.

Solving Linear Systems with QR

QR decomposition excels at solving overdetermined systems (more equations than unknowns) that arise in least squares problems. Unlike normal equations which can be numerically unstable, QR-based solutions maintain accuracy.

# Overdetermined system: fit a line to noisy data
np.random.seed(42)
x = np.linspace(0, 10, 50)
y_true = 3 * x + 7
y_noisy = y_true + np.random.normal(0, 2, 50)

# Build design matrix for y = mx + b
A = np.column_stack([x, np.ones_like(x)])

# Solve using QR decomposition
Q, R = np.linalg.qr(A)
# Solve R * params = Q^T * y for params
params = np.linalg.solve(R, Q.T @ y_noisy)

print(f"Estimated slope: {params[0]:.4f} (true: 3.0)")
print(f"Estimated intercept: {params[1]:.4f} (true: 7.0)")

# Compare with numpy's lstsq (which uses SVD)
params_lstsq = np.linalg.lstsq(A, y_noisy, rcond=None)[0]
print(f"lstsq slope: {params_lstsq[0]:.4f}")
print(f"lstsq intercept: {params_lstsq[1]:.4f}")

The QR approach solves the normal equations A^T A x = A^T b by transforming them into R x = Q^T b, avoiding the explicit formation of A^T A which can amplify numerical errors.

Polynomial Fitting with QR

QR decomposition handles polynomial regression efficiently, especially for higher-degree polynomials where normal equations become ill-conditioned.

# Generate data with polynomial relationship
x = np.linspace(-3, 3, 100)
y_true = 2*x**3 - 5*x**2 + 3*x + 1
y_noisy = y_true + np.random.normal(0, 5, 100)

# Create Vandermonde matrix for degree 3 polynomial
degree = 3
A = np.vander(x, degree + 1, increasing=True)

# QR-based solution
Q, R = np.linalg.qr(A)
coeffs = np.linalg.solve(R, Q.T @ y_noisy)

print("Estimated coefficients:", coeffs)
print("True coefficients: [1, 3, -5, 2]")

# Evaluate fit quality
y_fit = A @ coeffs
residual = np.linalg.norm(y_noisy - y_fit)
print(f"Residual norm: {residual:.4f}")

Computing Matrix Rank and Null Space

QR with column pivoting (available through SciPy) reveals matrix rank and helps identify linear dependencies. While NumPy doesn’t provide pivoting directly, understanding the R matrix’s diagonal elements indicates rank.

from scipy.linalg import qr

# Rank-deficient matrix
A = np.array([[1, 2, 3],
              [2, 4, 6],
              [1, 1, 2],
              [3, 5, 8]], dtype=float)

# QR with column pivoting
Q, R, P = qr(A, pivoting=True)

# Determine rank from R diagonal
tolerance = 1e-10
rank = np.sum(np.abs(np.diag(R)) > tolerance)
print(f"Matrix rank: {rank}")
print(f"R diagonal: {np.diag(R)}")

# The permutation matrix P reorders columns
print(f"Permutation: {P}")
print(f"Verification: {np.allclose(A[:, P], Q @ R)}")

Gram-Schmidt Orthogonalization Implementation

QR decomposition is essentially the Gram-Schmidt process in matrix form. Understanding this connection helps debug numerical issues and implement custom variants.

def gram_schmidt_qr(A):
    """Manual QR using classical Gram-Schmidt"""
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))
    
    for j in range(n):
        # Start with original column
        v = A[:, j].copy()
        
        # Subtract projections onto previous columns
        for i in range(j):
            R[i, j] = Q[:, i] @ A[:, j]
            v -= R[i, j] * Q[:, i]
        
        # Normalize
        R[j, j] = np.linalg.norm(v)
        Q[:, j] = v / R[j, j]
    
    return Q, R

# Test implementation
A = np.random.randn(5, 3)
Q_gs, R_gs = gram_schmidt_qr(A)
Q_np, R_np = np.linalg.qr(A)

print("Reconstruction error (Gram-Schmidt):", 
      np.linalg.norm(A - Q_gs @ R_gs))
print("Orthogonality error (Gram-Schmidt):", 
      np.linalg.norm(Q_gs.T @ Q_gs - np.eye(3)))

Note that classical Gram-Schmidt can suffer from numerical instability. NumPy uses Householder reflections, which provide better numerical properties.

Performance Considerations

QR decomposition has O(mn²) complexity for an m×n matrix. For large matrices, choose the appropriate mode and consider whether you need the full decomposition.

import time

sizes = [100, 500, 1000, 2000]
for n in sizes:
    A = np.random.randn(n, n)
    
    # Time QR decomposition
    start = time.time()
    Q, R = np.linalg.qr(A)
    qr_time = time.time() - start
    
    # Time LU decomposition for comparison
    from scipy.linalg import lu
    start = time.time()
    P, L, U = lu(A)
    lu_time = time.time() - start
    
    print(f"n={n}: QR={qr_time:.4f}s, LU={lu_time:.4f}s")

When to Use QR Over Alternatives

Use QR decomposition when:

Solving least squares problems where numerical stability matters
Working with tall matrices (m » n) in overdetermined systems
Computing orthonormal bases for subspaces
Implementing iterative eigenvalue algorithms (QR algorithm)

Avoid QR when:

Solving square systems with unique solutions (LU is faster)
Working with symmetric positive definite matrices (Cholesky is faster)
Needing singular values or dealing with severe ill-conditioning (use SVD)

# Comparison: solving Ax = b with different methods
n = 1000
A = np.random.randn(n, n)
b = np.random.randn(n)

# QR method
Q, R = np.linalg.qr(A)
x_qr = np.linalg.solve(R, Q.T @ b)

# Direct solve (uses LU)
x_direct = np.linalg.solve(A, b)

# Both should give same result
print(f"Solution difference: {np.linalg.norm(x_qr - x_direct)}")
print(f"Residual (QR): {np.linalg.norm(A @ x_qr - b)}")
print(f"Residual (LU): {np.linalg.norm(A @ x_direct - b)}")

QR decomposition provides a robust middle ground between the speed of LU decomposition and the comprehensive analysis of SVD, making it the go-to choice for stable solutions to overdetermined linear systems.