mathematics//linear algebra//matrix factorization//QR decomposition

The QR decomposition is the factorization of a matrix into a matrix \(Q\) with orthonormal columns and an upper triangular matrix \(R\), and it is the stable way to solve a least squares problem: the fit of a calibration, a regression or a parameter estimate, without the loss of digits that the textbook formula brings. For a tall matrix \(A\) with more rows than columns,

The QR decomposition is the factorization of a matrix into a matrix QQQ with orthonormal columns and an upper triangular matrix RRR, and it is the stable way to solve a least squares problem: the fit of a calibration, a regression or a parameter estimate, without the loss of digits that the textbook formula brings. For a tall matrix AAA with more rows than columns,

A=QR,x^=R−1QTb.A=QR,\qquad \hat x=R^{-1}Q^{\mathsf T}b .A=QR,x^=R−1QTb.

The second formula is never computed with an inverse: RRR is triangular, so Rx^=QTbR\hat x=Q^{\mathsf T}bRx^=QTb is solved by back substitution. Why it works is geometric. Multiplying by QTQ^{\mathsf T}QT is a rotation into coordinates aligned with the column space of AAA, and an orthogonal matrix has condition number 1, so that step neither amplifies errors nor loses digits. What is left is a triangular system whose conditioning is that of AAA itself.

QR works with κ(A)\kappa(A)κ(A); the normal equations work with κ(A)2\kappa(A)^2κ(A)2. Fitting a polynomial on poorly scaled data with κ(A)=108\kappa(A)=10^{8}κ(A)=108 loses about eight digits through QR and all sixteen through ATAA^{\mathsf T}AATA, which is the whole argument for never forming the normal matrix (condition number).

It is computed with Householder reflections, each one an orthogonal step that zeroes a column below the diagonal, at a cost of about 2mn22mn^22mn2 operations for an m×nm\times nm×n matrix: roughly twice the Cholesky route through ATAA^{\mathsf T}AATA, the price of the accuracy.

Among the least-squares solvers it sits in the middle. The Cholesky of the normal matrix is faster and fragile; the SVD is slower and also reports how close the problem is to rank loss, which QR does not, and NumPy's lstsq uses it. QR is the usual default in solvers that fit many well-posed problems, and column-pivoted QR adds a cheap estimate of the numerical rank.

QR can be updated as rows arrive, one orthogonal rotation per new measurement, which keeps a running fit numerically clean. Square-root information filters and some implementations of recursive least squares are built on that update, preferred on long missions where the plain recursion slowly loses the symmetry of its covariance.