Prev
Unit 5 · Least Squares and Data Fitting
QR Factorisation
60m100 XP
Next
📐Least Squares and Data Fitting

QR Factorisation

Orthogonal matrices, Gram–Schmidt and its instability, Householder reflections, Givens rotations

Stanford, 1965
In the early 1960s, statistical software solves regression problems the way the textbooks say: form ATAA^{\mathsf T}A, solve the normal equations. On the computers of the day, with data that are often nearly collinear, the results are sometimes nonsense — and nothing warns the user. Gene Golub, a young numerical analyst, publishes a nine-page paper in Numerische Mathematik: Numerical methods for solving linear least squares problems. Its idea comes from Alston Householder's 1958 work on reflections. Transform AA by a sequence of mirror reflections — each one exactly orthogonal, so it cannot amplify any error — until it is upper triangular. Then the least-squares problem is a single back substitution, and the condition number is never squared. Within a decade, QR by Householder reflections had become the standard way numerical libraries solve least-squares problems. It is what MATLAB's backslash does with a rectangular matrix, and what LAPACK's QR-based least-squares drivers do.
A square matrix QQ is orthogonal if QTQ=IQ^{\mathsf T}Q = I: its columns are orthonormal. Orthogonal matrices are rotations and reflections, and they preserve length: ∥Qv∥2=∥v∥2\|Q\mathbf v\|_2 = \|\mathbf v\|_2. Two consequences make them the numerical analyst's favourite tool:
  • κ2(Q)=1\kappa_2(Q) = 1 — the best possible;
  • multiplying by QQ never magnifies errors already present.
So transforming a problem by orthogonal matrices does not change its sensitivity — unlike forming ATAA^{\mathsf T}A.
💡The QR factorisation
Any m×nm \times n matrix AA (m≥nm \ge n) can be written A=QR=[Q1Q2][R10]=Q1R1,A = QR = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix}\begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1R_1, with QQ orthogonal (m×mm \times m) and R1R_1 upper triangular (n×nn \times n). Because QQ preserves length, ∥b−Ax∥2=∥QTb−Rx∥2=∥Q1Tb−R1x∥2+∥Q2Tb∥2.\|\mathbf b - A\mathbf x\|^2 = \|Q^{\mathsf T}\mathbf b - R\mathbf x\|^2 = \|Q_1^{\mathsf T}\mathbf b - R_1\mathbf x\|^2 + \|Q_2^{\mathsf T}\mathbf b\|^2. The second term does not depend on x\mathbf x. The first is zero when R1x=Q1TbR_1\mathbf x = Q_1^{\mathsf T}\mathbf b — one triangular solve. And κ(R1)=κ(A)\kappa(R_1) = \kappa(A), not κ(A)2\kappa(A)^2.
Gram–Schmidt is the obvious way to compute QQ: take the columns of AA in turn and subtract from each its projections onto the orthonormal vectors already found, then normalise. It is taught in every linear algebra course — and in floating point it can fail badly. When columns are nearly parallel, the subtraction removes almost everything and the small remainder is dominated by rounding error; the computed QQ is no longer orthogonal. Modified Gram–Schmidt reorders the same operations (subtract each new direction from all later columns immediately) and loses orthogonality only in proportion to κ\kappa, not κ2\kappa^2.
Derivation
Gram–Schmidt produces QQ and RR together
Take the columns a1,a2,…\mathbf a_1, \mathbf a_2, \dots of AA in turn. Normalise the first: r11=∥a1∥r_{11} = \|\mathbf a_1\|, q1=a1/r11\mathbf q_1 = \mathbf a_1/r_{11}. For column jj, remove its components along the directions already found, then normalise what is left:
rij=qiTaj (i<j),vj=aj−∑i<jrij qi,rjj=∥vj∥,qj=vj/rjj.r_{ij} = \mathbf q_i^{\mathsf T}\mathbf a_j \ (i < j), \qquad \mathbf v_j = \mathbf a_j - \sum_{i<j} r_{ij}\,\mathbf q_i, \qquad r_{jj} = \|\mathbf v_j\|, \qquad \mathbf q_j = \mathbf v_j/r_{jj}.
Rearranged, each column is a combination of the q\mathbf qs found so far:
aj=r1jq1+r2jq2+⋯+rjjqj.\mathbf a_j = r_{1j}\mathbf q_1 + r_{2j}\mathbf q_2 + \dots + r_{jj}\mathbf q_j.
Column jj uses only q1,…,qj\mathbf q_1, \dots, \mathbf q_j, so the coefficients form an upper-triangular matrix: A=QRA = QR.
Where it goes wrong. If aj\mathbf a_j is nearly a combination of earlier columns, vj\mathbf v_j is a tiny difference of nearly equal vectors, ∥vj∥≪∥aj∥\|\mathbf v_j\| \ll \|\mathbf a_j\|. The rounding errors in that subtraction are of size u∥aj∥u\|\mathbf a_j\|, so relative to vj\mathbf v_j they are u∥aj∥/∥vj∥u\|\mathbf a_j\|/\|\mathbf v_j\| — catastrophic cancellation (lesson 1-2) — and they point in arbitrary directions, including along the earlier q\mathbf qs. Normalising divides by the small ∥vj∥\|\mathbf v_j\|, so those errors survive in qj\mathbf q_j at that relative size — close to 1 when the columns are nearly dependent — and qj\mathbf q_j is no longer orthogonal to its predecessors.
💡Householder reflections
For a vector x\mathbf x, the reflection H=I−2vvTvTv,v=x+sign⁡(x1)∥x∥e1,H = I - 2\dfrac{\mathbf v\mathbf v^{\mathsf T}}{\mathbf v^{\mathsf T}\mathbf v}, \qquad \mathbf v = \mathbf x + \operatorname{sign}(x_1)\|\mathbf x\|\mathbf e_1, maps x\mathbf x to −sign⁡(x1)∥x∥e1-\operatorname{sign}(x_1)\|\mathbf x\|\mathbf e_1: it zeroes every entry below the first in one step. HH is symmetric and orthogonal; it is a mirror. Apply one to column 1 of AA, the next to column 2 of what remains, and so on: after nn reflections AA is upper triangular, and QQ is the product of the mirrors (never formed explicitly — the vectors v\mathbf v are stored instead). The sign choice avoids cancellation in forming v\mathbf v. Cost: 2mn2−23n32mn^2 - \tfrac23 n^3 operations, about twice the normal equations — the price of stability.
Derivation
The mirror is orthogonal, and sends x\mathbf x onto the axis
Orthogonal. HH is clearly symmetric. Square it, writing P=vvT/vTvP = \mathbf v\mathbf v^{\mathsf T}/\mathbf v^{\mathsf T}\mathbf v; then P2=v(vTv)vT/(vTv)2=PP^2 = \mathbf v(\mathbf v^{\mathsf T}\mathbf v)\mathbf v^{\mathsf T}/(\mathbf v^{\mathsf T}\mathbf v)^2 = P, so
HTH=H2=(I−2P)2=I−4P+4P2=I.H^{\mathsf T}H = H^2 = (I - 2P)^2 = I - 4P + 4P^2 = I.
Geometrically, PP projects onto the line of v\mathbf v and I−2PI - 2P reflects through the plane perpendicular to v\mathbf v; reflecting twice changes nothing.
Onto the axis. Let α=−sign⁡(x1)∥x∥\alpha = -\operatorname{sign}(x_1)\|\mathbf x\|, so v=x−αe1\mathbf v = \mathbf x - \alpha\mathbf e_1. Using α2=∥x∥2\alpha^2 = \|\mathbf x\|^2:
vTv=∥x∥2−2αx1+α2=2(∥x∥2−αx1),vTx=∥x∥2−αx1.\mathbf v^{\mathsf T}\mathbf v = \|\mathbf x\|^2 - 2\alpha x_1 + \alpha^2 = 2\big(\|\mathbf x\|^2 - \alpha x_1\big), \qquad \mathbf v^{\mathsf T}\mathbf x = \|\mathbf x\|^2 - \alpha x_1.
The second is exactly half the first, so
Hx=x−2v vTxvTv=x−v=α e1.H\mathbf x = \mathbf x - 2\mathbf v\,\frac{\mathbf v^{\mathsf T}\mathbf x}{\mathbf v^{\mathsf T}\mathbf v} = \mathbf x - \mathbf v = \alpha\,\mathbf e_1.
The sign. With α\alpha of the opposite sign to x1x_1, the quantity ∥x∥2−αx1=∥x∥2+∣x1∣ ∥x∥\|\mathbf x\|^2 - \alpha x_1 = \|\mathbf x\|^2 + |x_1|\,\|\mathbf x\| is a sum of two positive numbers. The other choice would subtract them — and when x\mathbf x already points almost along e1\mathbf e_1 they are nearly equal: cancellation.
📐The polynomial fit, done right
The degree-9 fit of lesson 5-1, solved by a hand-written Householder QR: relative coefficient error 1.0×10−101.0 \times 10^{-10}, against 2.2×10−42.2 \times 10^{-4} from the normal equations — six more correct digits from the same data, the same arithmetic and the same answer in exact arithmetic. The orthogonality of the computed QQ tells the same story. For 12 columns 1,t,…,t111, t, \dots, t^{11} (κ=1.2×108\kappa = 1.2 \times 10^8), ∥QTQ−I∥\|Q^{\mathsf T}Q - I\| is 2.42.4 for classical Gram–Schmidt — not orthogonal at all — 5×10−95 \times 10^{-9} for modified Gram–Schmidt and 1.4×10−151.4 \times 10^{-15} for Householder.
Symmetry is a vast subject, significant in art and nature. Mathematics lies at its root, and it would be hard to find a better one on which to demonstrate the working of the mathematical intellect.
– Hermann Weyl, Symmetry (1952)
💡Givens rotations
A Givens rotation acts on two rows only, rotating them so that one chosen entry becomes zero: [cs−sc][ab]=[a2+b20],c=aa2+b2, s=ba2+b2.\begin{bmatrix} c & s \\ -s & c \end{bmatrix}\begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sqrt{a^2+b^2} \\ 0 \end{bmatrix}, \qquad c = \frac{a}{\sqrt{a^2+b^2}},\ s = \frac{b}{\sqrt{a^2+b^2}}. Where Householder clears a whole column at once, Givens clears one entry at a time — which is exactly right for sparse matrices and for updating: when a new measurement row arrives, a few rotations fold it into RR without revisiting the old data. A sensor calibration refined over 200 incoming readings this way, storing only the 2×22 \times 2 factor RR and two numbers, matches the batch solution to 10−1210^{-12}.
💡Where you will use this
GNSS Navigation and Attitude Control: square-root information filters propagate RR by Householder and Givens transformations instead of covariance matrices, because they stay accurate in single precision on flight computers. System identification: recursive least squares with Givens updates tracks changing parameters online. Structures: QR inside eigenvalue solvers for vibration (Unit 6). Data analysis: every regression routine in a serious statistics package.
📚Further reading
  • Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lectures 7, 8, 10 and 11 — QR, Gram–Schmidt, Householder triangularisation and least squares.
  • Gene Golub, 'Numerical methods for solving linear least squares problems', Numerische Mathematik 7, 206–216 (1965) — the paper of the cold open.
  • Åke Björck, Numerical Methods for Least Squares Problems (SIAM, 1996), chapters 2–3 — QR methods, modified Gram–Schmidt and updating a factorisation as data arrive.
↓But there is a problem
QR solves full-rank least-squares problems stably. But what if the columns of AA are nearly dependent — two sensors that measure almost the same thing, a model with one parameter too many? Then R1R_1 has a tiny diagonal entry, the solution is enormous and meaningless, and QR does not tell you which combination of parameters the data cannot determine. In October 2006 Netflix offered one million dollars to anyone who could predict its customers' film ratings 10% better than its own system. The data were a matrix of 480,000 customers by 17,770 films — almost entirely empty. The ideas that won were built on one factorisation.