Eric Campbell 6 min read Statistics

How to Perform QR Decomposition in Python

QR decomposition is a fundamental matrix factorization technique that decomposes any matrix A into the product of two matrices: Q (an orthogonal matrix) and R (an upper triangular matrix)....

Key Insights

  • QR decomposition factors any matrix A into an orthogonal matrix Q and upper triangular matrix R, providing superior numerical stability compared to methods like Gaussian elimination for solving linear systems.
  • NumPy’s qr() function offers ‘reduced’ and ‘complete’ modes—use reduced mode for overdetermined systems to save memory and computation time.
  • For solving least squares problems, QR decomposition outperforms the normal equations approach (A^T A x = A^T b) by avoiding the condition number squaring that leads to numerical instability.

Introduction to QR Decomposition

QR decomposition is a fundamental matrix factorization technique that decomposes any matrix A into the product of two matrices: Q (an orthogonal matrix) and R (an upper triangular matrix). Mathematically, this is expressed as A = QR.

The power of QR decomposition lies in its numerical stability. Unlike methods such as Gaussian elimination that can amplify rounding errors, QR decomposition maintains accuracy even when working with ill-conditioned matrices. This makes it the preferred choice for solving least squares problems, computing eigenvalues through the QR algorithm, and performing orthogonalization in machine learning applications like PCA and QR-based neural network training.

In practical terms, you’ll reach for QR decomposition when solving overdetermined linear systems (more equations than unknowns), computing matrix rank, or when you need a numerically stable alternative to the Gram-Schmidt process for orthogonalization.

Mathematical Foundation

An orthogonal matrix Q has a special property: its columns form an orthonormal set of vectors. This means Q^T Q = I (the identity matrix). Geometrically, orthogonal matrices represent rotations and reflections—they preserve lengths and angles.

An upper triangular matrix R has all zeros below its main diagonal. This structure makes it trivial to solve systems like Rx = b through back substitution, working from the bottom row upward.

When we decompose A = QR, we’re essentially finding an orthonormal basis (Q) and expressing our original matrix as a transformation in that basis (R). Let’s verify these properties with code:

import numpy as np

# Create a simple 3x3 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("Q matrix (orthogonal):")
print(Q)
print("\nR matrix (upper triangular):")
print(R)

# Verify Q is orthogonal: Q^T @ Q should equal I
print("\nQ^T @ Q (should be identity):")
print(Q.T @ Q)

# Verify decomposition: Q @ R should equal A
print("\nQ @ R (should equal original A):")
print(Q @ R)

# Check if R is upper triangular
print("\nIs R upper triangular?", np.allclose(R, np.triu(R)))

This code demonstrates that Q^T Q produces the identity matrix (within floating-point precision) and that multiplying Q and R reconstructs the original matrix A.

QR Decomposition with NumPy

NumPy provides numpy.linalg.qr() as the primary interface for QR decomposition. The function supports two modes that control the dimensions of the output matrices:

  • ‘reduced’ mode (default): For an m×n matrix where m ≥ n, returns Q as m×n and R as n×n
  • ‘complete’ mode: Returns full-sized Q (m×m) and R (m×n)

For most applications, especially overdetermined systems, reduced mode is preferred because it’s more memory-efficient and computationally faster.

import numpy as np

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

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

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

# Both should reconstruct A
print("\nReconstruction error (reduced):", 
      np.linalg.norm(A - Q_reduced @ R_reduced))
print("Reconstruction error (complete):", 
      np.linalg.norm(A - Q_complete @ R_complete))

# Verify orthogonality for reduced Q
print("\nQ_reduced^T @ Q_reduced:")
print(Q_reduced.T @ Q_reduced)

The reduced mode is particularly useful when you’re only interested in solving Ax = b, as the extra columns in complete mode don’t contribute to the solution.

QR Decomposition with SciPy

SciPy’s scipy.linalg.qr() offers additional functionality beyond NumPy, most notably column pivoting. Pivoting reorders columns to improve numerical stability and can reveal the rank of a matrix more reliably.

import numpy as np
from scipy import linalg

# Create a matrix with near-linear dependence
A = np.array([[1, 2, 2.001],
              [3, 4, 7.003],
              [5, 6, 11.005]], dtype=float)

# Standard QR
Q, R = linalg.qr(A)

# QR with column pivoting
Q_p, R_p, P = linalg.qr(A, pivoting=True)

print("Standard R matrix:")
print(R)
print("\nR with pivoting:")
print(R_p)
print("\nPivot indices:", P)

# Pivoting reveals near-rank deficiency through small diagonal values
print("\nDiagonal of R_p:", np.diag(R_p))

The pivoting option returns an additional permutation array P, where A[:, P] = Q @ R. This reordering places columns with larger norms first, making it easier to detect rank deficiency by examining the diagonal elements of R.

For most standard applications, NumPy’s implementation is sufficient. Use SciPy when you need pivoting for rank detection or when working with matrices that might be rank-deficient.

Practical Applications

The most common application of QR decomposition is solving overdetermined linear least squares problems. Given Ax = b where A is m×n (m > n), we want to find x that minimizes ||Ax - b||.

import numpy as np

# Create an overdetermined system
np.random.seed(42)
A = np.random.randn(100, 5)  # 100 equations, 5 unknowns
x_true = np.array([1, 2, 3, 4, 5])
b = A @ x_true + 0.1 * np.random.randn(100)  # Add noise

# Solve using QR decomposition
Q, R = np.linalg.qr(A, mode='reduced')
# Transform: Ax = b -> QRx = b -> Rx = Q^T b
x_qr = np.linalg.solve(R, Q.T @ b)

print("True x:", x_true)
print("Estimated x:", x_qr)
print("Residual norm:", np.linalg.norm(A @ x_qr - b))

# Compare with normal equations (less stable)
x_normal = np.linalg.solve(A.T @ A, A.T @ b)
print("Normal equations x:", x_normal)

QR decomposition also provides an efficient way to compute matrix rank:

import numpy as np

def matrix_rank_qr(A, tol=1e-10):
    """Compute matrix rank using QR decomposition"""
    Q, R = np.linalg.qr(A)
    # Count diagonal elements above tolerance
    return np.sum(np.abs(np.diag(R)) > tol)

# Test with rank-deficient matrix
A = np.array([[1, 2, 3],
              [2, 4, 6],
              [3, 5, 7]], dtype=float)

print("Matrix rank:", matrix_rank_qr(A))
print("NumPy rank (for verification):", np.linalg.matrix_rank(A))

Performance Comparison and Best Practices

Understanding when to use QR versus other decompositions is crucial for writing efficient numerical code. Here’s a practical comparison:

import numpy as np
import time

def benchmark_qr(n, m):
    """Benchmark QR decomposition for different sizes"""
    A = np.random.randn(n, m)
    
    # NumPy QR
    start = time.time()
    Q, R = np.linalg.qr(A, mode='reduced')
    numpy_time = time.time() - start
    
    # SciPy QR
    from scipy import linalg
    start = time.time()
    Q, R = linalg.qr(A, mode='economic')
    scipy_time = time.time() - start
    
    return numpy_time, scipy_time

# Test different sizes
sizes = [(100, 50), (500, 100), (1000, 200)]
print("Matrix Size | NumPy (s) | SciPy (s)")
print("-" * 40)
for n, m in sizes:
    np_t, sp_t = benchmark_qr(n, m)
    print(f"{n}x{m:4d}    | {np_t:.5f}   | {sp_t:.5f}")

Best practices:

  1. Use reduced mode for overdetermined systems to save memory
  2. Prefer QR over normal equations for least squares—it’s more stable and doesn’t square the condition number
  3. Use SVD instead of QR when you need the complete rank-revealing decomposition or when dealing with severely ill-conditioned matrices
  4. Use LU decomposition for square systems where numerical stability isn’t a primary concern—it’s faster
  5. Enable pivoting when rank detection is important or when working with potentially rank-deficient matrices

For large matrices, consider that QR decomposition has O(mn²) complexity for an m×n matrix. If you’re repeatedly solving systems with the same matrix, factor once and reuse Q and R.

Conclusion

QR decomposition is an essential tool in numerical linear algebra, offering superior stability for solving overdetermined systems and least squares problems. NumPy’s qr() function handles most use cases efficiently, while SciPy’s version provides pivoting for rank detection.

Use reduced mode for memory efficiency, leverage the decomposition for stable least squares solving, and remember that while QR is more expensive than LU decomposition, its numerical stability often justifies the cost. For production code dealing with real-world data, QR decomposition should be your default choice for solving rectangular systems.

For deeper exploration, investigate the QR algorithm for eigenvalue computation, Householder reflections for implementing QR from scratch, and the relationship between QR decomposition and the Gram-Schmidt process.

Liked this? There's more.

Every week: one practical technique, explained simply, with code you can use immediately.