QR Decomposition

Fitting a regression by literally computing $(X^TX)^{-1}X^Ty$ works on paper but is numerically fragile: squaring $X$ into $X^TX$ squares any near-collinearity problems in the data before you even start inverting. QR decomposition sidesteps this entirely — it factors the design matrix directly, without ever forming $X^TX$, and is the actual method most serious regression libraries use under the hood.

An analogy: a rotate-then-stretch recipe

Any matrix's effect on space can be described as first rotating (and possibly reflecting) everything rigidly, without any distortion, and then stretching or shearing along the new axes. QR decomposition writes exactly this recipe down explicitly: $Q$ is the pure rotation/reflection part — it preserves lengths and angles perfectly, like a rigid camera turn — and $R$ is the triangular "stretch" part that does all the actual distorting. Because $Q$ never distorts anything, it never amplifies numerical errors, which is precisely why algorithms built on QR tend to be far more numerically trustworthy than ones that square the matrix first.

The mathematics, one symbol at a time

QR decomposition factors a matrix $A$ (say, $m \times n$ with $m \geq n$, a tall design matrix) as

$$A = QR,$$

where $Q$ is an $m \times n$ matrix with orthonormal columns (each column has length 1, and every pair of columns is orthogonal, so $Q^TQ = I$) and $R$ is an $n \times n$ upper-triangular matrix. In words: $Q$'s columns form a "clean," non-overlapping coordinate system for the same space $A$'s columns span, built via Gram-Schmidt (see Gram-Schmidt Orthogonalization), and $R$ records how the original columns of $A$ combine those clean directions to reconstruct themselves. For solving least squares, $Ax \approx b$ becomes $QRx \approx b$; multiplying both sides by $Q^T$ (which, crucially, exactly undoes $Q$ since $Q^TQ=I$, with no numerical error introduced) gives $Rx = Q^Tb$ — a small, easy, upper-triangular system, solved directly by back-substitution, with none of the error-amplifying squaring that $(X^TX)^{-1}X^Ty$ requires.

Worked example 1: QR by hand via Gram-Schmidt

Take $A$'s columns $v_1=(3,1)$, $v_2=(2,2)$ — the same pair from the Gram-Schmidt example. Orthogonalizing gave $u_1=(3,1)$, $u_2=(-0.4,1.2)$. Normalize: $\|u_1\|=\sqrt{10}\approx3.162$, so $q_1=(0.949,0.316)$; $\|u_2\|=\sqrt{0.16+1.44}=\sqrt{1.6}\approx1.265$, so $q_2=(-0.316,0.949)$. These form $Q$. Then $R = Q^TA$: $r_{11}=\langle q_1,v_1\rangle = \|u_1\| = 3.162$, $r_{12}=\langle q_1,v_2\rangle=0.949(2)+0.316(2)\approx2.530$, $r_{22}=\|u_2\|=1.265$ (and $r_{21}=0$ by construction). So $R\approx\begin{pmatrix}3.162 & 2.530\\0&1.265\end{pmatrix}$ — upper-triangular, exactly as required.

Worked example 2: solving a small regression via QR

Reusing the earlier one-variable regression setup, $A=\begin{pmatrix}1\\2\\3\end{pmatrix}$, $b=\begin{pmatrix}2\\3\\5\end{pmatrix}$: $\|A\|=\sqrt{1+4+9}=\sqrt{14}\approx3.742$, so $q_1 = A/3.742 \approx (0.267,0.535,0.802)$, and $R = [\,\|A\|\,] = [3.742]$ (a $1\times1$ matrix here, since there's only one column). Then $Q^Tb = 0.267(2)+0.535(3)+0.802(5) = 0.535+1.604+4.008\approx6.147$. Solving $R\beta = Q^Tb$: $3.742\beta = 6.147 \Rightarrow \beta\approx1.643$ — matching the earlier pseudoinverse answer of $23/14\approx1.643$ exactly, but reached without ever squaring $A$ into $A^TA$.

What this means in practice

QR decomposition is the default method inside numpy.linalg.lstsq, R's lm(), and most production regression code precisely because it avoids ever forming $X^TX$, keeping numerical error manageable even with moderately collinear factors. It's also the standard way to compute an orthonormal basis for a subspace, and shows up as a building block inside eigenvalue algorithms and the singular value decomposition (see Singular Value Decomposition).