M5-T5 — QR Decomposition

Learning objectives

After this lesson, you should be able to:

• State the QR decomposition A = QR for a tall matrix and the shapes of Q and R.
• Derive that the least-squares solution satisfies R x̂ = Qᵀb and that the residual norm is the length of the leftover part of Qᵀb.
• Explain why QR is more accurate than the normal equations (orthogonal steps; κ(R) = κ(A)).
• Describe Householder QR as a sequence of reflections, and know to use a library routine.

Introduction and motivation

A robot fitting a trajectory, a calibration or a surface needs an answer it can trust to many digits even when the model is nearly redundant. The method behind nearly all production fitting code is QR. It is the natural continuation of the previous two modules: Gram-Schmidt showed that a matrix is an orthonormal part times a triangle, orthogonal matrices showed that orthonormal parts do no damage, and projections showed that least squares is a shadow. QR puts these together.

The intuition

Any tall matrix can be written A = QR, where Q has orthonormal columns and R is a small upper triangular matrix. Think of Q as a clean set of axes for the column space of A, and R as the coordinates of the original columns in those axes (this is what Gram-Schmidt computed).

Now the fit becomes easy. Rotating the data into the clean axes does not change any lengths; in those axes the part of b inside the column space is simply a list of numbers, Qᵀb, and what is left over is the residual. The unknowns solve a small triangular system R x = Qᵀb by back substitution.

Geometric and physical interpretation

Geometric view. Q rotates (and possibly reflects) the data space so that the column space of A lies along the first n coordinate axes. In the rotated frame the data vector has n components inside the column space (these determine the fit) and m − n components outside it (the residual). Rotations preserve lengths, so no error is amplified. Physical view. Imagine tilting a measurement table until the plane of possible predictions is flat on the table top. The fit is then read off from the position of the data over the table, and the residual is its height above it.

The mathematics

A = QR

the (economy) QR decomposition: Q is m × n with orthonormal columns, R is n × n upper triangular.

Qᵀb

the coordinates of the data in the clean axes: the first n entries of the rotated data.

κ(R)

the condition number of R, equal to κ(A).

The decomposition

A = Q R       QᵀQ = I (m × n Q),       R upper triangular (n × n),   rii ≠ 0 when the columns of A are independent

What it means. Q holds orthonormal axes for the column space; R holds the coordinates of the columns of A in those axes. R is triangular because column k of A needs only the first k axes.

The solution

minimise ‖A x − b‖²    ⇒    R x̂ = Qᵀ b          ‖r‖ = ‖ (I − QQᵀ) b ‖

What it means. Solve the small triangular system by back substitution. The residual is whatever part of b lies outside the span of the columns of Q.

Why it is accurate

κ(R) = κ(A)          (but κ(AᵀA) = κ(A)²)

What it means. Orthogonal factors do not change the singular values, so R has exactly the sensitivity of A. Solving with R faces κ(A); forming AᵀA faces κ(A)².

QR: A = QR, solve R x̂ = Qᵀb by back substitution. It avoids AᵀA, keeps the sensitivity at κ(A), and is the standard least-squares solver.

Step-by-step derivation: R x̂ = Qᵀb, and Householder reflections

From the cost to the triangular system. Let A = QR (economy size). From the normal equations AᵀAx̂ = Aᵀb:

1. Substitute A = QR: AᵀA = RᵀQᵀQR = RᵀR, since QᵀQ = I.
2. And Aᵀb = RᵀQᵀb. So the normal equations read RᵀRx̂ = RᵀQᵀb.
3. If the columns are independent, R is invertible, and so is Rᵀ; cancel it: R x̂ = Qᵀb.

We reached the same solution as the normal equations algebraically, but the computer never forms AᵀA.

A cleaner derivation of the residual. Extend Q to a full orthogonal matrix [Q Q⊥] (square). Rotation does not change length, so ‖Ax − b‖² = ‖[Q Q⊥]ᵀ(Ax − b)‖². Now QᵀAx = Rx and Q⊥ᵀAx = 0 (A's columns lie in the span of Q), so

‖Ax − b‖² = ‖R x − Qᵀb‖² + ‖Q⊥ᵀ b‖²

What it means. The first piece can be made exactly zero by choosing x, which solves Rx = Qᵀb. The second piece, the size of the part of b outside the column space, does not depend on x at all: it is the unavoidable residual.

So the minimum is attained at Rx̂ = Qᵀb, and the minimum cost is ‖Q⊥ᵀb‖² = ‖r‖².

Why not Gram-Schmidt in practice. Classical Gram-Schmidt can lose orthogonality completely (last module). Householder QR builds Q as a product of reflections H = I − 2uuᵀ: each one is exactly orthogonal, chosen to zero the entries of one column below the diagonal. Because every step is an orthogonal transformation, errors cannot grow, and the computed Q is orthonormal to machine precision (about 3·10⁻¹⁶ on the Hilbert example against 1 for classical Gram-Schmidt).

A worked example: the line fit again, through QR

Use A = [[1, 1], [2, 1], [3, 1]] and b = (1, 2, 2).

1. Factor. Gram-Schmidt on the two columns gives q₁ = (1, 2, 3)/√14 = (0.267, 0.535, 0.802), then the remainder of column 2 after removing its shadow gives q₂ = (0.873, 0.218, −0.436), and R = [[3.742, 1.604], [0, 0.655]].
2. Rotate the data. Qᵀb = (q₁·b, q₂·b) = (0.267 + 1.069 + 1.604, 0.873 + 0.436 − 0.873) = (2.940, 0.436).
3. Back substitute R x = (2.940, 0.436): from the second row, 0.655 c = 0.436 gives c = 0.667. From the first, 3.742 m + 1.604·0.667 = 2.940 gives m = 0.5.
4. Residual size. ‖b‖² − ‖Qᵀb‖² = 9 − (8.644 + 0.190) = 0.167, matching ‖r‖² = 0.1667 from the normal equations.
5. Conditioning. The singular values of A are the same as those of R, so κ(R) = κ(A) = 6.79, while the normal equations deal with 46.1.
6. The hard case. For the degree-9 polynomial fit, this solve returns a relative error near 7·10⁻¹¹, versus 3·10⁻⁴ for the normal equations.

Robotics application

Default solver. numpy.linalg.lstsq, MATLAB's backslash on tall matrices and most calibration toolboxes use QR (or the SVD) internally. When you call a library to fit a plane, a polynomial or a calibration, this is what is happening.

Plane and line fitting with data scaling. Fitting a plane z = ax + by + c to LiDAR points far from the origin gives a design matrix whose columns are nearly parallel; QR (after centring the data) handles it, where the normal equations may not. The SLAM lessons on LiDAR point clouds and point-cloud preprocessing use such fits.

Square-root filtering. Kalman filters in square-root form carry a triangular factor (the Cholesky or QR factor of the covariance) instead of the covariance itself, for exactly the same conditioning reason: square roots halve the dynamic range of the numbers. This is standard in flight code.

Rank detection. A very small diagonal entry rkk shows the k-th column is nearly dependent on the earlier ones, the same early warning as a small Gram-Schmidt remainder.

Engineering implications and limitations

Practical rules

• Call a library QR (Householder), not your own Gram-Schmidt. The results are orthonormal to machine precision.
• Use the economy factorisation. For a tall matrix store Q as m × n, not m × m. A full Q for a million points would be a trillion entries.
• Check the diagonal of R. Entries near zero relative to the largest mean dependent columns. Use column-pivoted QR or the SVD when the rank is in doubt.
• The cost is about 2mn² − 2n³/3 operations, roughly twice that of forming and solving the normal equations, a small price for accuracy; for extremely large sparse problems sparse QR or Cholesky on regularised normal equations may be used knowingly.
• Reuse the factorisation for several right-hand sides, as with LU: factor A once, apply Qᵀ to each b and back-substitute.

Limits

QR assumes independent columns. For rank-deficient problems it needs pivoting, and for the most delicate cases (an explicit rank decision, a minimum-norm answer) the SVD of Module 6 is the right tool. QR improves the algorithm; it cannot repair an ill-conditioned problem: if κ(A) is 10¹⁵ the answer is uncertain by nature.

Common misconceptions

Engineering takeaways

1. A tall matrix factors as A = QR with orthonormal columns in Q and an upper triangular R.
2. Least squares reduces to R x̂ = Qᵀb, solved by back substitution; the residual norm is the size of what lies outside the span of Q.
3. Orthogonal steps do not amplify error: κ(R) = κ(A), whereas the normal equations face κ(A)².
4. Householder QR (reflections) is the stable way to compute it; classical Gram-Schmidt is not.
5. Use library QR, centre your data, and keep the SVD in mind for rank questions.

References and further reading

Primary papers and authoritative documentation.

• MIT 18.06SC Linear Algebra (OCW): lectures on orthogonal bases and Gram-Schmidt
• Boyd and Vandenberghe, Introduction to Applied Linear Algebra (VMLS): the chapter on QR factorisation and least squares; the site did not load when we last checked
• Mathematics for Machine Learning: the chapter on matrix decompositions
• Michigan ROB 101: Computational Linear Algebra: QR and least squares in a computational course

Knowledge check

Three conceptual questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

So far every measurement counted equally. Real sensors differ in quality. Weighted least squares lets a precise measurement count for more than a sloppy one, and connects the weights to covariance.

Related content

Nothing below is required to finish this lesson.

• Before
  M5-T4: Normal Equations

  Derive and solve the normal equations, and state when they are valid and when they lose accuracy.

• Before
  M4-T4: The Gram-Schmidt Process

  Turn independent vectors into an orthonormal basis step by step, and recognise its numerical weakness.

• Lab
  M5-L1: LiDAR Plane Fitting

  Fit the ground plane to a noisy LiDAR scan, find out why a naive fit fails on a wall, and fix it with a stable method.