QR factorisation and numerical linear algebra
Practical computation depends on stable algorithms. The QR factorisation, built from orthogonal transformations, is the workhorse for least squares and eigenvalues.
The QR factorisation
Any m x n matrix A with independent columns can be written \(A = QR\), where Q is m x n with orthonormal columns spanning Col A, and R is n x n upper triangular with positive diagonal. Q comes from applying Gram-Schmidt to A's columns, and R records the coefficients: \(R = Q^TA\).
QR for least squares
With \(A = QR\), the normal equations \(A^TA\hat{\mathbf{x}} = A^T\mathbf{b}\) simplify to \(R\hat{\mathbf{x}} = Q^T\mathbf{b}\), a triangular system solved by back substitution. This avoids forming \(A^TA\), which squares the condition number, so QR gives more accurate least-squares solutions and is the standard approach in statistical software.
Householder reflections and Givens rotations
Rather than Gram-Schmidt, numerical libraries build QR from Householder reflections, introduced by Alston Householder in 1958, each reflecting a column to zero out entries below the diagonal, or Givens rotations, which zero one entry at a time and suit sparse matrices. Both use orthogonal matrices and are backward stable.
The QR algorithm for eigenvalues
The QR algorithm computes all eigenvalues: factor \(A_0 = Q_0R_0\), form \(A_1 = R_0Q_0\), and repeat. Each \(A_{k+1} = Q_k^TA_kQ_k\) is similar to A, so eigenvalues are preserved, and the iterates approach triangular form with the eigenvalues on the diagonal. With shifts and a preliminary reduction to Hessenberg form, it converges rapidly and is named among the top algorithms of the twentieth century.
Floating-point arithmetic
Computers store real numbers in floating point, typically IEEE 754 double precision with about 16 significant decimal digits, so nearly every operation rounds. A stable algorithm keeps the effect of rounding small relative to the problem's inherent sensitivity. Orthogonal transformations are favoured because they preserve lengths and never magnify errors.
Conditioning and accuracy
The condition number \(\kappa(A) = \|A\|\|A^{-1}\|\) bounds how much relative errors in data can be amplified in the solution of \(A\mathbf{x} = \mathbf{b}\). A rule of thumb: with condition number \(10^k\), expect to lose about k digits of accuracy. The Hilbert matrices, with entries \(\frac{1}{i + j - 1}\), are famously ill-conditioned: the 10 x 10 Hilbert matrix has a condition number around \(10^{13}\).
Software and scale
Libraries such as BLAS and LAPACK, developed from the 1970s onwards, implement these algorithms efficiently and underlie NumPy, MATLAB, R, Julia, and much scientific software. Iterative methods such as conjugate gradients (1952) and GMRES handle huge sparse systems, and GPUs accelerate dense linear algebra for machine learning, where matrix multiplication dominates computation.
Worked example
For \(A\) with columns \(\mathbf{a}_1 = (3, 4)\) and \(\mathbf{a}_2 = (1, 2)\): \(\mathbf{q}_1 = \frac{1}{5}(3, 4)\). The projection coefficient is \(r_{12} = \mathbf{q}_1 \cdot \mathbf{a}_2 = \frac{11}{5}\); then \(\mathbf{a}_2 - \frac{11}{5}\mathbf{q}_1 = (1, 2) - \frac{11}{25}(3, 4) = (-\frac{8}{25}, \frac{6}{25})\), with length \(\frac{10}{25} = \frac{2}{5}\), so \(\mathbf{q}_2 = (-0.8, 0.6)\) and \(R = \begin{bmatrix} 5 & 2.2 \\ 0 & 0.4 \end{bmatrix}\). Check: \(\det R = 2 = \det A\), since \(\det Q = 1\).
Key takeaways
A = QR with orthonormal-column Q and upper triangular R exists for matrices with independent columns; Gram-Schmidt gives Q and R = Q^T A. Least squares becomes R x-hat = Q^T b, avoiding the conditioning loss of A^T A. Libraries build QR from Householder reflections (1958) or Givens rotations, which are backward stable. The QR algorithm iterates A_{k+1} = R_k Q_k, preserving eigenvalues and converging to triangular form, especially with shifts and Hessenberg reduction. Floating point (IEEE 754 doubles, about 16 digits) rounds every step; orthogonal methods do not amplify errors. Condition number 10^k loses about k digits; Hilbert matrices are notoriously ill-conditioned. BLAS, LAPACK, iterative methods, and GPUs power modern computation.
Flashcard
A = Q R
Check your understanding
Why do statistical programs solve least squares with QR instead of the normal equations?