QR decomposition

The formula

The core identity, for an m×n matrix with m ≥ n:

A = Q · R

Q has orthonormal columns (Qᵀ·Q = I), R is upper triangular. Reading: the columns of Q are a clean perpendicular reference frame for the space your data columns actually span; R records how to rebuild each original column from that frame.

Two sizes show up in library docs:

full QR:  Q is m×m,  R is m×n with zeros below row n
thin QR:  Q is m×n,  R is n×n            (what you want for least squares)

Reading: same factorization, with or without the extra columns of Q that point out of your data's span. For a matrix with independent columns, the thin factorization is unique once you require positive diagonal entries in R ^[1].

The reason engineers meet it — the overdetermined system A·x ≈ b, more measurements than unknowns:

minimize ‖A·x − b‖   →   solve  R·x = Qᵀ·b   by back substitution

Reading: multiplying by Qᵀ rotates the problem without stretching anything, so the fit collapses to one triangular solve. Compare what your stats textbook taught:

normal equations:  (Aᵀ·A)·x = Aᵀ·b     with  cond(Aᵀ·A) = cond(A)²

Reading: forming Aᵀ·A squares the condition number before you've solved anything. QR never forms it. That single line is the entire reason this page exists.

Where you meet it

Every calibration fit you've ever run. MATLAB's backslash on a tall system, LAPACK's least-squares drivers, the fit button in your DAQ software — the workhorse underneath is QR (or its close cousin the SVD), not the normal equations ^[1]. When you fit a fifth-order polynomial through thermocouple calibration points and the answer comes back sane, a Householder QR quietly earned its keep.
Force-balance and multi-channel strain calibration on the stand. A six-component balance calibration applies dozens of known load combinations and records every bridge output, cross-talk included. That's an overdetermined system — more loading runs than sensitivity coefficients — and the calibration matrix falls out of one QR solve. The redundancy is the point: extra runs average down the noise instead of being discarded.
Navigation software, every measurement update. Square-root Kalman filters carry a triangular factor of the covariance instead of the covariance itself, and each update is a fresh QR re-triangularization. Working with the factor means the effective precision doubles — the same trick as least squares: never square what you can keep as a square root.
At the review board, when the chart says x = (AᵀA)⁻¹Aᵀb. As notation, fine. As the implementation, it's a finding. Ask what the condition number of the design matrix is; if the answer is 10⁵ or worse, the normal equations are running with 10¹⁰ and the coefficient error bars on the next chart are decoration.

How it works

QR is Gram–Schmidt with the receipts kept. Walk the columns of A left to right; from each column subtract its projections onto the clean directions already built, normalize what's left, and file the subtracted amounts into R. Small worked case:

A = | 3  1 |      Q = | 0.6  −0.8 |      R = | 5  3 |
    | 4  3 |          | 0.8   0.6 |          | 0  1 |

Column 1 of A has length 5 → r₁₁ = 5, q₁ = (0.6, 0.8). Column 2 leans on q₁ by amount q₁·a₂ = 3 → r₁₂ = 3; the perpendicular remainder has length 1 → q₂ = (−0.8, 0.6). Multiply Q·R and A comes back exactly, and det(A) = 5·1 = 5 falls out of the diagonal of R for free. One caution when checking against software: signs are only pinned down by convention — NumPy returns this example with the first column of Q and first row of R negated. Same factorization, different handshake; compare Q·R, not the factors.

Why it beats the normal equations. Orthogonal matrices have condition number 1 — multiplying by Qᵀ adds no damage. So QR solves the fit at sensitivity cond(A); the normal equations solve it at cond(A)². In double precision (~16 digits) that's the difference between losing 8 digits and losing all 16 on a cond(A) = 10⁸ problem. It isn't hypothetical: take A = |1 1; δ 0; 0 δ| with δ = 10⁻⁸ and data consistent with the answer (1, 1). Form Aᵀ·A in double precision and the 10⁻¹⁶ terms round away, leaving |1 1; 1 1| — exactly singular; the solver throws. Thin QR on the same data returns (1.0, 1.0). Verified: cond(A) ≈ 1.4×10⁸, comfortably solvable, but not after squaring.

Not all QRs are equal. Three routes to the same factorization:

Classical Gram–Schmidt — the textbook loop. Numerically treacherous: on a 50×12 polynomial design matrix with cond ≈ 10⁸, the computed Q satisfies ‖QᵀQ − I‖ ≈ 2.3. That is not a typo; the "orthonormal" columns have completely lost orthogonality ^[1]^[2].
Modified Gram–Schmidt — same arithmetic count, projections subtracted one at a time against the running remainder. Same matrix: ‖QᵀQ − I‖ ≈ 5×10⁻⁹. Better, and good enough for many jobs ^[2].
Householder reflections — build Q from mirror operations instead of projections. Same matrix: ‖QᵀQ − I‖ ≈ 10⁻¹⁵. Backward stable, and what LAPACK runs by default ^[1]. Givens rotations do the same job one matrix entry at a time — the tool of choice for sparse problems and for cheaply updating a factorization when the next data row arrives mid-test ^[1].

The mistakes people make:

Forming AᵀA in a script because the formula in the textbook says so. The cost savings are trivial — Householder QR runs about twice the arithmetic of the normal equations — and the accuracy cost is your condition number, squared.
Trusting QR to rescue a bad experiment. QR avoids adding damage; it can't remove the cond(A) your test design baked in. Nearly collinear predictors — say, load and temperature that always moved together during the run — stay nearly collinear. Check conditioning; better, fix the test matrix.
Rank deficiency. If a column of A is a combination of the others, plain QR hits a zero on R's diagonal. Column-pivoted QR (A·P = Q·R) reorders columns to expose the rank and tells you which channel is redundant ^[1].

History

The algorithm is far older than the name. Laplace was manipulating observation equations with what amounts to sequential orthogonalization in a 1820 supplement to his treatise on probability — least-squares astronomy, the original noisy test data ^[2]^[3]. Jørgen Pedersen Gram, a Danish actuary, published the process in 1883 in a paper on series expansion by least squares; Erhard Schmidt gave the clean inner-product formulation in 1907 while building the theory of integral equations, and noted himself that the procedure was in essence Gram's ^[3]^[4]. The double name stuck anyway.

The modern, machine-era story is compressed into a few years around 1958–1965. Alston Householder at Oak Ridge showed in a four-page 1958 note in the Journal of the ACM that a matrix could be triangularized by unitary reflections — more accurate than elimination-style transforms and cheaper than the rotation schemes then in use ^[5]^[6]. Wallace Givens, whose career ran through Argonne, published the plane-rotation route to triangular form the same year ^[7]^[8]. Gene Golub closed the loop in 1965 in Numerische Mathematik, showing Householder's factorization was the right production tool for least squares — the paper that made "don't form the normal equations" official doctrine ^[3]^[9].

The strangest chapter belongs to the QR algorithm — the eigenvalue iteration built on repeated QR factorization, the method inside every eig() call and every modal-frequency extraction. John Francis, a programmer at the National Research Development Corporation in London with no PhD, submitted his first paper on it in October 1959, at age 25; it appeared in The Computer Journal in 1961–62 ^[10]^[11]. Vera Kublanovskaya at the Steklov Institute in Leningrad submitted the same idea independently in July 1960 ^[10]^[12]. Francis left the field by 1962 and spent the next four and a half decades entirely unaware that his method had become a pillar of scientific computing — it was named one of the ten most important algorithms of the twentieth century in 2000 ^[11]^[13]. Gene Golub and Frank Uhlig finally tracked him down in 2007, retired on the English south coast ^[10]^[11]. He got a standing ovation at the algorithm's 50th-anniversary meeting in 2009 and an honorary doctorate from Sussex in 2015 — roughly 54 years after the fact ^[11].

Related tools

/tools/strain-gauge-bridge — one channel of the multi-run calibration whose full sensitivity matrix a QR least-squares solve extracts
/tools/loop-4-20ma — the two-point exactly-determined scaling; add calibration points beyond coefficients and QR is the safe way to use them
/tools/snr-enob — the acquisition noise floor that decides how many digits your fit had to lose before QR vs. normal equations mattered
/tools/vibration-natural-freq — the single-mode estimate; the multi-mode version comes from an eigensolver running Francis's QR iteration

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →