QR Decomposition

Factoring a matrix into an orthonormal part and an upper triangular part, the stable workhorse for least squares.

The factorization

QR decomposition writes a matrix A as A = QR, where Q has orthonormal columns and R is upper triangular. The columns of Q form an orthonormal basis for the column space of A, and R records how the original columns combine to produce them. It is the matrix form of the Gram-Schmidt process, though computed more stably.

How it is computed

Practical QR uses Householder reflections, each of which zeros out the entries below a diagonal position by reflecting a vector onto a coordinate axis. Applying a sequence of these reflections triangularizes A while keeping every operation orthogonal, so no error amplification occurs. Givens rotations are an alternative, well suited to sparse or streaming problems.

Solving least squares

For an overdetermined system Ax = b, the least-squares solution follows from QR without ever forming A^T A: since A = QR, the normal equations reduce to R x = Q^T b, a single triangular solve. This route is far more accurate than the normal equations because it avoids squaring the condition number of A.

The QR algorithm for eigenvalues

Repeatedly factoring a matrix as QR and re-multiplying as RQ produces a sequence that converges to a triangular matrix whose diagonal holds the eigenvalues. This QR algorithm, with shifts and deflation, is the standard method behind eigenvalue routines in numerical libraries.

import numpy as np
A = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
b = np.array([1.0, 2.0, 2.0])
Q, R = np.linalg.qr(A)
x = np.linalg.solve(R, Q.T @ b)   # least-squares solution
print(x)

Because orthogonal transformations never magnify error, QR is preferred wherever conditioning is a concern, including fits of noisy diagnostic data against physics models.