QR (Gram–Schmidt, Householder, Givens) and least squares done right

Before this lesson: lesson 03 (inverse, determinant, dot product and orthogonal projection), lesson 04 (orthogonal matrices and change of basis) and lesson 08 (the SVD and the pseudoinverse). The condition number $\kappa(A)$ appears here and gets its full treatment in lesson 13. Time: about 90 minutes, plus the notebook.

---

Why this matters #

Least squares, $\min_x \|Ax - b\|_2$, is the most common computation in data analysis: linear regression, curve fitting, linear probes, every Gauss–Newton step, and many inverse problems in science. In exact arithmetic all standard methods give the same answer. In floating point they do not. The textbook formula $(A^\top A)^{-1}A^\top b$ squares the condition number: with $\kappa(A) = 10^7$ it keeps about 2–3 correct digits, while QR keeps about 9–10. QR is also a building block you will meet again: it orthonormalizes the sketch in randomized SVD (lesson 09), it drives the QR eigenvalue algorithm, and Arnoldi is Gram–Schmidt on a Krylov sequence (lesson 14).

1. Predict first #

2. The picture #

One column first. Fit $b \approx x\,a$ with one column $a \in \mathbb{R}^m$ and a scalar $x$. The best $x$ makes the error $b - xa$ orthogonal to $a$: $a^\top(b - \hat x a) = 0$, so $\hat x = a^\top b / a^\top a$. The fitted vector $\hat x a$ is the orthogonal projection of $b$ onto the line through $a$.

Many columns. $Ax = x_1a_1 + \dots + x_na_n$ ranges over the column space $C(A)$. The closest point to $b$ is the orthogonal projection $p = A\hat x$, and the residual $r = b - A\hat x$ is orthogonal to every column. In the language of lesson 02: $r \in N(A^\top)$.

Interactive figure: loading…

Orthonormal columns make this easy. If $A = Q$ with $Q^\top Q = I$, the coordinates are plain dot products, $\hat x = Q^\top b$: no system to solve. QR is the step that makes the columns orthonormal and records the change in an upper triangular matrix $R$.

Three ways to get there.

• Gram–Schmidt builds $q_1, q_2, \dots$ one at a time: take $a_k$, remove its components along $q_1, \dots, q_{k-1}$, normalize ("triangular orthogonalization").
• Householder multiplies $A$ on the left by reflections until it is upper triangular ("orthogonal triangularization").
• Givens does the same with 2-D rotations, each one zeroing a single entry.

3. Definitions and main results #

For $A \in \mathbb{R}^{m\times n}$, $m \ge n$, the reduced (thin) QR is $A = QR$ with $Q \in \mathbb{R}^{m\times n}$, $Q^\top Q = I_n$, and $R \in \mathbb{R}^{n\times n}$ upper triangular. The full QR is $A = \tilde Q \binom{R}{0}$ with $\tilde Q = [\,Q \;\; Q_2\,]$ orthogonal ($m \times m$). For full column rank, the columns of $Q_2$ are an orthonormal basis of $N(A^\top)$.

Existence and uniqueness of QR #

Gram–Schmidt by hand first. Take $a_1 = (1, 1, 0)$ and $a_2 = (1, 0, 1)$.

1. Normalize $a_1$: $r_{11} = \|a_1\| = \sqrt2$, so $q_1 = (1, 1, 0)/\sqrt2$.
2. The part of $a_2$ along $q_1$ has length $r_{12} = q_1^\top a_2 = (1 + 0 + 0)/\sqrt2 = 1/\sqrt2$.
3. Remove it: $v_2 = a_2 - r_{12}q_1$. Here $r_{12}q_1 = \tfrac{1}{\sqrt2}\cdot\tfrac{1}{\sqrt2}(1, 1, 0) = (\tfrac12, \tfrac12, 0)$, so $v_2 = (1, 0, 1) - (\tfrac12, \tfrac12, 0) = (\tfrac12, -\tfrac12, 1)$. Check: $q_1^\top v_2 = (\tfrac12 - \tfrac12 + 0)/\sqrt2 = 0$ ✓.
4. Normalize: $r_{22} = \|v_2\| = \sqrt{\tfrac14 + \tfrac14 + 1} = \sqrt{1.5}$, so $q_2 = v_2/\sqrt{1.5} = (1, -1, 2)/\sqrt6$.
5. Collect the coefficients in $R$. Column $k$ of $R$ holds the numbers that rebuild $a_k$: $a_1 = r_{11}q_1$ and $a_2 = r_{12}q_1 + r_{22}q_2$. So

Check one entry: row 3 of $Q$ times column 2 of $R$ is $0 \cdot \tfrac{1}{\sqrt2} + \tfrac{2}{\sqrt6}\sqrt{1.5} = \tfrac{2}{\sqrt6}\cdot\tfrac{\sqrt6}{2} = 1 = a_{32}$ ✓. $R$ is upper triangular because $a_1$ uses only $q_1$, and $a_2$ uses only $q_1$ and $q_2$.

Interactive figure: loading…

For the three columns of the figure, step 4 and 5 continue the same way: $r_{13} = q_1^\top a_3 = 1/\sqrt2$, $r_{23} = q_2^\top a_3 = 1/\sqrt6$, $v_3 = a_3 - r_{13}q_1 - r_{23}q_2 = (-\tfrac23, \tfrac23, \tfrac23)$, $r_{33} = 2/\sqrt3$ and $q_3 = (-1, 1, 1)/\sqrt3$.

Theorem (QR). If $A$ has independent columns $a_1, \dots, a_n$, then $A = QR$ with orthonormal columns in $Q$ and $R$ upper triangular with $r_{kk} > 0$, and this factorization is unique.

Proof of existence (Gram–Schmidt).

1. $a_1 \ne 0$ (independent columns). Set $r_{11} = \|a_1\|$ and $q_1 = a_1/r_{11}$.
2. Assume $q_1, \dots, q_{k-1}$ are orthonormal and $\operatorname{span}(q_1, \dots, q_{k-1}) = \operatorname{span}(a_1, \dots, a_{k-1})$.
3. Set $r_{jk} = q_j^\top a_k$ for $j < k$ and $v_k = a_k - \sum_{j<k} r_{jk}q_j$. Then $q_i^\top v_k = q_i^\top a_k - r_{ik} = 0$ for every $i < k$.
4. $v_k \ne 0$: otherwise $a_k \in \operatorname{span}(q_1, \dots, q_{k-1}) = \operatorname{span}(a_1, \dots, a_{k-1})$, which contradicts independence. Set $r_{kk} = \|v_k\| > 0$ and $q_k = v_k/r_{kk}$. The span property now holds for $k$.
5. Rearranged, $a_k = \sum_{j<k} r_{jk}q_j + r_{kk}q_k$: column $k$ of $A$ is $Q$ times column $k$ of $R$, and it uses only $q_1, \dots, q_k$, so $R$ is upper triangular. Hence $A = QR$. $\blacksquare$

Proof of uniqueness

1. Let $A = Q_1R_1 = Q_2R_2$, both of this form. Then $Q_1 = Q_2T$ with $T = R_2R_1^{-1}$, which is upper triangular with positive diagonal. Also $I = Q_1^\top Q_1 = T^\top Q_2^\top Q_2 T = T^\top T$, so $T$ is orthogonal.
2. $T^{-1} = T^\top$ is upper triangular (inverse of an upper triangular matrix) and lower triangular (transpose of one), so $T$ is diagonal with entries $\pm 1$. A positive diagonal forces $T = I$: $R_1 = R_2$ and $Q_1 = Q_2$. $\blacksquare$

Another route: $A^\top A = R^\top R$ is the Cholesky factorization of the positive definite $A^\top A$, and the Cholesky factor with positive diagonal is unique (lesson 07).

$r_{kk} = \|v_k\|$ is the distance from $a_k$ to the span of the earlier columns. A tiny $r_{kk}$ warns that $a_k$ is nearly a combination of earlier columns. Without independence, QR still exists (Householder computes it for any $A$), but some $r_{kk} = 0$ and uniqueness is lost.

Classical and modified Gram–Schmidt #

The two versions differ only in which vector each projection is computed from. Look at column 3:

• Classical Gram–Schmidt (CGS) computes all coefficients from the original column: $r_{13} = q_1^\top a_3$ and $r_{23} = q_2^\top a_3$. Then it subtracts both: $v_3 = a_3 - r_{13}q_1 - r_{23}q_2$.
• Modified Gram–Schmidt (MGS) subtracts one projection at a time and computes the next coefficient from what is left: $v \leftarrow a_3 - r_{13}q_1$ with $r_{13} = q_1^\top a_3$, then $r_{23} = q_2^\top v$ (the current, partly reduced vector) and $v_3 = v - r_{23}q_2$.

In exact arithmetic $q_2^\top q_1 = 0$, so $q_2^\top(a_3 - r_{13}q_1) = q_2^\top a_3$: the two methods give the same numbers. In floating point they do not:

• Why CGS fails. If $a_k$ is nearly in the span of the earlier columns, $v_k$ is a small difference of large, nearly equal vectors ($r_{kk} \ll \|a_k\|$). Rounding errors of size about $\varepsilon\|a_k\|$ ($\varepsilon \approx 2.2 \times 10^{-16}$, the machine epsilon of double precision, lesson 13) survive the subtraction and are divided by $r_{kk}$. So $q_k$ keeps a component of size about $\varepsilon\|a_k\|/r_{kk}$ along the earlier $q_j$, and the errors grow from column to column.
• Why MGS helps. MGS takes each coefficient from the vector that is actually being reduced. So each subtraction also removes the $q_j$-part of the rounding errors made in the earlier subtractions. CGS fixes all coefficients in advance and never corrects these errors.
• How large. The loss of orthogonality $\|Q^\top Q - I\|_2$ is about $\varepsilon\kappa(A)^2$ for CGS (order 1 from $\kappa \approx 10^8$ on), about $\varepsilon\kappa(A)$ for MGS, and about $\varepsilon$ for Householder, independent of $\kappa$. Orthogonalizing each column twice (CGS2) brings Gram–Schmidt back to order $\varepsilon$ when $\varepsilon\kappa(A) \ll 1$.

In the nearly dependent case of the figure, $a_1 = (1, \delta, 0)$, $a_2 = (1, 0, \delta)$, $a_3 = (1, 0, 0)$ with $\delta = 10^{-8}$ ($\kappa(A) \approx 3 \times 10^8$): CGS gives $\|Q^\top Q - I\|_2 \approx 0.71$ (the angle between $q_2$ and $q_3$ is $45°$), MGS gives about $10^{-8}$.

Householder reflections #

A Householder reflection is $H = I - 2\,vv^\top/(v^\top v)$ for a nonzero $v$. It reflects across the hyperplane orthogonal to $v$: $Hv = -v$, and $Hw = w$ for $w \perp v$. So $H$ is symmetric and orthogonal, and $H^2 = I$.

To map $x$ to a multiple of $e_1$, the target must be $\pm\|x\|e_1$ (reflections keep length), and $v = x - \text{target}$. Choose

with the convention $\operatorname{sign}(0) = +1$.

One step by hand. Let $x = (1, 2, 2)$ be the first column of a $3 \times 3$ matrix.

1. $\|x\| = \sqrt{1 + 4 + 4} = 3$ and $x_1 > 0$, so the target is $-3e_1$.
2. $v = x + 3e_1 = (4, 2, 2)$; here $v_1 = 1 + 3$ adds two positive numbers.
3. $v^\top v = 16 + 4 + 4 = 24$ and $v^\top x = 4 + 4 + 4 = 12$.
4. $Hx = x - 2\,\tfrac{v^\top x}{v^\top v}\,v = x - v = (-3, 0, 0)$ ✓: two entries are zeroed.
5. Apply the same $H$ to another column $y = (1, 0, 1)$, without building $H$: $v^\top y = 6$, so $Hy = y - 2 \cdot \tfrac{6}{24}\,v = y - \tfrac12 v = (-1, -1, 0)$. Check: $\|Hy\| = \sqrt2 = \|y\|$ ✓.

With the other sign, $v = x - 3e_1 = (-2, 2, 2)$ and $Hx = (3, 0, 0)$. That is also correct in exact arithmetic. The problem appears when $x$ is close to the $e_1$ axis: then $v_1 = x_1 - \|x\|$ subtracts two nearly equal numbers (cancellation). Example: $x = (1, 10^{-9}, 0)$. $\|x\|$ rounds to exactly 1, so the wrong sign gives $v_1 = 1 - 1 = 0$ (true value $-5 \times 10^{-19}$), and the computed $H$ maps $x$ to $(1, -10^{-9}, 0)$: nothing was zeroed. The right sign gives $v = (2, 10^{-9}, 0)$ and $Hx = (-1, 0, 0)$.

Interactive figure: loading…

Householder QR. For $k = 1, \dots, n$: let $x$ be column $k$ on and below the diagonal, build $v_k$, and update the remaining block $X$ by the rank-1 formula $X \leftarrow X - 2v_k(v_k^\top X)/(v_k^\top v_k)$. Never form the $m\times m$ matrix $H_k$. After $n$ steps $H_n\cdots H_1A = \binom{R}{0}$, so $\tilde Q = H_1\cdots H_n$. For least squares you do not need $Q$ itself: apply the same reflections to $b$. Cost: about $2mn^2 - \tfrac23 n^3$ flops.

Why Householder QR is backward stable. Backward stable means: the computed $\hat R$ is the exact R-factor of a nearby matrix, $\tilde Q\hat R = A + \delta A$ with $\tilde Q$ exactly orthogonal and $\|\delta A\| \le c\,\varepsilon\|A\|$ ($c$ a modest constant depending on $m, n$).

1. Each step multiplies by an orthogonal matrix. Its rounding error is a perturbation of relative size $O(\varepsilon)$ of the current matrix.
2. Orthogonal maps keep 2-norms ($\|Hy\| = \|y\|$), so later steps carry earlier errors along without enlarging them. The errors only add up: no division by a small $r_{kk}$, no cancelling subtraction.
3. The computed $\hat Q$ is a product of reflections, so $\|\hat Q^\top\hat Q - I\| = O(\varepsilon)$ for every $\kappa(A)$. Orthogonality is built in, not produced by subtraction.

Gram–Schmidt also gives a small $\|A - \hat Q\hat R\|$, but for large $\kappa(A)$ its $\hat Q$ is far from every orthogonal matrix. So CGS (and MGS) are not backward stable as QR algorithms.

Interactive figure: loading…

Givens rotations #

A Givens rotation changes two coordinates and leaves the others alone. On the pair $(a, b)$:

Example: $(3, 4)$ gives $r = 5$, $c = 3/5$, $s = 4/5$, and $G(3, 4) = (9/5 + 16/5,\ -12/5 + 12/5) = (5, 0)$. Givens QR zeros the subdiagonal entries one by one.

Givens rotations in practice

• Compute $r$ with hypot(a, b) so that $a^2$ cannot overflow.
• Givens QR is as stable as Householder but costs about 1.5 times more on a dense matrix ($\approx 3mn^2 - n^3$ flops).
• It wins when only a few entries must be zeroed: Hessenberg matrices (one rotation per column in the QR eigenvalue algorithm, lesson 14), sparse matrices, and updating $R$ when a new data row arrives (one sweep of $n$ rotations instead of a new factorization).
• Rotations on disjoint pairs of rows commute, so many of them can run in parallel.

The normal equations and the QR solution #

Theorem (least squares). Let $A$ have full column rank. Then $\hat x$ minimizes $\|Ax - b\|_2$ if and only if $A^\top A\hat x = A^\top b$ (the normal equations), and $\hat x$ solves $R\hat x = Q^\top b$, where $A = QR$ is the reduced QR.

Proof.

1. Projection. The closest point $A\hat x \in C(A)$ to $b$ is the orthogonal projection, so the residual is orthogonal to every column: $A^\top(b - A\hat x) = 0$, i.e. $A^\top A\hat x = A^\top b$.
2. Global minimum. For any $x$: $b - Ax = (b - A\hat x) + A(\hat x - x)$. The first part is in $N(A^\top)$, the second in $C(A)$, so they are orthogonal. By Pythagoras $\|b - Ax\|^2 = \|b - A\hat x\|^2 + \|A(\hat x - x)\|^2$, which is smallest only at $x = \hat x$ ($A$ has full column rank, so $A(\hat x - x) = 0$ forces $x = \hat x$).
3. QR. $A^\top A = R^\top Q^\top QR = R^\top R$ and $A^\top b = R^\top Q^\top b$, so $R^\top R\hat x = R^\top Q^\top b$. $R^\top$ is invertible, hence $R\hat x = Q^\top b$: one back substitution.
4. Full-QR view, without normal equations. $\tilde Q^\top$ keeps norms, so $\|Ax - b\|^2 = \|Rx - Q^\top b\|^2 + \|Q_2^\top b\|^2$. The first term is zero at $R\hat x = Q^\top b$; the second is the minimal squared residual. $\blacksquare$

Interactive figure: loading…

Why QR and not the normal equations #

• Squared condition number. From the SVD $A = U\Sigma V^\top$: $A^\top A = V\Sigma^2V^\top$, so $\kappa(A^\top A) = \kappa(A)^2$ (proof in lesson 13).
• Error sizes. For a small residual, a backward-stable method has relative error in $\hat x$ of about $\varepsilon\kappa(A)$. Householder QR and SVD solvers are backward stable for least squares. The normal equations are not: their error is about $\varepsilon\kappa(A)^2$, even when the residual is zero. (For a large residual, the problem itself has sensitivity of order $\kappa + \kappa^2\tan\theta$, with $\theta$ the angle between $b$ and $C(A)$; QR reaches this, and no method does better.)
• Total failure. For $\kappa(A) \gtrsim 1/\sqrt{\varepsilon} \approx 10^8$, the computed $A^\top A$ can be exactly singular. Läuchli example: $A = \begin{pmatrix} 1 & 1 \\ \delta & 0 \\ 0 & \delta\end{pmatrix}$, $\delta = 10^{-8}$, has full rank, but $A^\top A = \begin{pmatrix} 1 + \delta^2 & 1 \\ 1 & 1 + \delta^2\end{pmatrix}$ and $1 + 10^{-16}$ rounds to $1$.
• Cost. Normal equations $\approx mn^2 + n^3/3$ flops, Householder QR $\approx 2mn^2 - \tfrac23 n^3$ (up to 2× more), SVD $\approx 2mn^2 + 11n^3$. For small $\kappa$ the normal equations are fine ($\kappa \le 10^3$ gives errors near $10^{-10}$).

Answers to section 1. (1) Normal equations: $\varepsilon\kappa^2 \approx 10^{-16}\cdot 10^{14} = 10^{-2}$, so about 2–3 digits (a test with a random $A$ gave $10^{-3}$). QR: $\varepsilon\kappa \approx 10^{-9}$, about 9–10 digits (the test gave $7 \times 10^{-11}$). (2) No. Classical Gram–Schmidt's $Q$ has lost orthogonality completely ($\|Q^\top Q - I\|$ of order 1); Householder's is at $10^{-15}$. A small $\|A - QR\|$ says nothing about orthogonality.

The minimum-norm solution #

If $A$ is $m \times n$ with $m < n$ and full row rank, $Ax = b$ has infinitely many solutions $x_p + N(A)$.

Theorem. The solution of smallest $\|x\|_2$ is the unique solution in the row space: $x^+ = A^\top(AA^\top)^{-1}b = A^+b$. We call it the minimum-norm solution.

1. Try $x = A^\top y$: then $AA^\top y = b$, and $AA^\top$ is invertible (full row rank). So $x^+$ solves $Ax = b$ and lies in $C(A^\top)$.
2. Any other solution is $x^+ + z$ with $z \in N(A) \perp C(A^\top)$, so $\|x^+ + z\|^2 = \|x^+\|^2 + \|z\|^2 \ge \|x^+\|^2$, with equality only for $z = 0$. $\blacksquare$

Simplest case: the plane $x_1 + 2x_2 + 2x_3 = 9$. The closest point to the origin lies along the normal $a = (1, 2, 2)$: $x^+ = 9a/(a^\top a) = (1, 2, 2)$, norm 3. The solution $(9, 0, 0)$ has norm 9. To avoid forming $AA^\top$ (again $\kappa^2$), use $A^\top = QR$: then $A = R^\top Q^\top$; set $x = Qy$ and solve $R^\top y = b$ by forward substitution.

Rank-deficient $A$ (any shape). The least-squares solutions are $\hat x + N(A)$; the minimum-norm one is $A^+b = \sum_{\sigma_i > 0}(u_i^\top b/\sigma_i)\,v_i$ (the pseudoinverse, lesson 08). Plain QR fails here: back substitution divides by $r_{kk} \approx 0$. Use instead

• a truncated SVD (singular values below a tolerance count as 0), as np.linalg.lstsq does with rcond;
• QR with column pivoting, $AP = QR$: at each step the remaining column of largest norm goes first, so $|r_{11}| \ge |r_{22}| \ge \dots$, and a drop to about $\varepsilon\|A\|$ usually reveals the numerical rank;
• ridge, $(A^\top A + \lambda I)x = A^\top b$, which tends to $A^+b$ as $\lambda \to 0$.

4. Worked example #

This is the worked example case of the first figure.

Gram–Schmidt.

• $r_{11} = \|a_1\| = \sqrt{4 + 1 + 4} = 3$, so $q_1 = (2, 1, 2)/3$.
• $r_{12} = q_1^\top a_2 = (2 + 1 + 0)/3 = 1$, and $v_2 = (1, 1, 0) - (2, 1, 2)/3 = (1, 2, -2)/3$.
• $r_{22} = \|v_2\| = \sqrt{1 + 4 + 4}/3 = 1$, so $q_2 = (1, 2, -2)/3$. Check: $q_1^\top q_2 = (2 + 2 - 4)/9 = 0$.

Least squares by QR. $Q^\top b = \big((2 + 2 + 6)/3,\; (1 + 4 - 6)/3\big) = (10/3, -1/3)$. Back substitution, from the last row up: row 2 says $\hat x_2 = -1/3$; row 1 says $3\hat x_1 + \hat x_2 = 10/3$, so $3\hat x_1 = 10/3 + 1/3 = 11/3$ and $\hat x_1 = 11/9$.

Check with the normal equations. $A^\top A = \begin{pmatrix} 9 & 3 \\ 3 & 2\end{pmatrix} = R^\top R$ and $A^\top b = (10, 3)$. Solving $9\hat x_1 + 3\hat x_2 = 10$, $3\hat x_1 + 2\hat x_2 = 3$ gives the same $\hat x = (11/9, -1/3)$.

Residual. $r = b - A\hat x = (-10, 10, 5)/9$, and $A^\top r = 0$ ✓. $\|r\| = \sqrt{225}/9 = 5/3$. The full QR gives this without $\hat x$: $q_3 = (2, -2, -1)/3$ completes the orthonormal basis, and $|q_3^\top b| = |2 - 4 - 3|/3 = 5/3$ ✓.

Householder on the same $A$. $x = a_1 = (2, 1, 2)$, $\|x\| = 3$, $x_1 > 0$, so $v_1 = (5, 1, 2)$, $v_1^\top v_1 = 30$, $H_1a_1 = a_1 - 2v_1(15/30) = a_1 - v_1 = (-3, 0, 0)$ and $H_1a_2 = (1, 1, 0) - \tfrac{2\cdot 6}{30}(5, 1, 2) = (-1, 3/5, -4/5)$. Step 2 on $x = (3/5, -4/5)$, $\|x\| = 1$: $v_2 = (8/5, -4/5)$ and $H_2x = (-1, 0)$. So Householder gives $\begin{pmatrix} -3 & -1 \\ 0 & -1\end{pmatrix} = -R$: the same factorization up to the signs $D = -I$, as uniqueness predicts.

5. In code #

NumPy checks the hand computation. The notebook for this lesson has these cells and the exercises.

import numpy as np

A = np.array([[2., 1.], [1., 1.], [2., 0.]])
b = np.array([1., 2., 3.])

Q, R = np.linalg.qr(A)                  # Householder (LAPACK geqrf); the signs of diag(R) are free
print(R.round(12))                      # [[-3. -1.] [ 0. -1.]]: our R with D = -I
d = np.sign(np.diag(R))
print((d[:, None] * R).round(12) + 0.0) # [[3. 1.] [0. 1.]]: the unique R with positive diagonal

x, res, rank, sv = np.linalg.lstsq(A, b, rcond=None)   # SVD-based LAPACK solver (gelsd)
print(x, res)                           # [ 1.2222 -0.3333] [2.7778]: (11/9, -1/3), ‖r‖² = 25/9
print(A.T @ (b - A @ x))                # ~1e-15: the residual is orthogonal to C(A)

Orthogonality that is built in (Householder) against orthogonality computed through $A^\top A$ ("Cholesky QR": $R$ from the Cholesky factor of $A^\top A$, then $Q = AR^{-1}$):

import numpy as np
from scipy.linalg import solve_triangular

rng = np.random.default_rng(0)
def with_condition(kappa, m=80, n=20):  # random A with singular values from 1 down to 1/kappa
    U, _ = np.linalg.qr(rng.standard_normal((m, n)))
    V, _ = np.linalg.qr(rng.standard_normal((n, n)))
    return (U * np.logspace(0, -np.log10(kappa), n)) @ V.T

for kappa in [1e2, 1e4, 1e6]:
    A = with_condition(kappa)
    I = np.eye(A.shape[1])
    Qh, _ = np.linalg.qr(A)                         # Householder
    L = np.linalg.cholesky(A.T @ A)                 # Cholesky QR: R = Lᵀ, Q = A R⁻¹
    Qc = solve_triangular(L, A.T, lower=True).T
    print(f"kappa {kappa:.0e}  cond(AᵀA) {np.linalg.cond(A.T @ A):.1e}  "
          f"Householder {np.abs(Qh.T @ Qh - I).max():.0e}  Cholesky-QR {np.abs(Qc.T @ Qc - I).max():.0e}")
# cond(AᵀA) = 1e4, 1e8, 1e12 = kappa²; Householder stays near 1e-15; Cholesky-QR grows like
# eps·kappa² (2e-13, 1e-9, 2e-5)

The Läuchli matrix, a rank-deficient matrix with pivoted QR, and the minimum-norm solution:

import numpy as np
from scipy.linalg import qr

d = 1e-8
A = np.array([[1., 1.], [d, 0.], [0., d]])
print(A.T @ A)                                      # [[1. 1.] [1. 1.]]: 1 + d² rounded to 1
print(np.linalg.matrix_rank(A), np.linalg.matrix_rank(A.T @ A))   # 2 1
print(np.linalg.lstsq(A, np.array([2., d, d]), rcond=None)[0])     # [1. 1.]: fine on A itself

B = np.array([[1., 0., 1.], [0., 1., 1.], [1., 1., 2.], [1., 0., 1.]])   # col 3 = col 1 + col 2
Q, R, perm = qr(B, pivoting=True)                   # column-pivoted QR: B[:, perm] = Q R
print(np.abs(np.diag(R)))                           # [2.6 0.85 ~2e-16]: numerical rank 2
x = np.linalg.lstsq(B, np.array([1., 2., 3., 4.]), rcond=None)[0]   # minimum-norm LS solution
print(x @ np.array([1., 1., -1.]))                  # ~0: x is orthogonal to N(B) = span(1, 1, -1)

6. Where it shows up #

• Regression software. R's lm uses Householder QR; NumPy's lstsq and scikit-learn's LinearRegression use the SVD-based LAPACK driver gelsd. None of them forms $(X^\top X)^{-1}$.
• Random orthogonal matrices. torch.nn.init.orthogonal_ takes the QR of a Gaussian matrix and multiplies the columns of $Q$ by $\operatorname{sign}(\operatorname{diag} R)$. Without this sign fix, $Q$ is not uniformly (Haar) distributed.
• ML for science. Gauss–Newton and Levenberg–Marquardt steps in inverse problems and data assimilation are linear least-squares problems. Square-root Kalman filters update a triangular factor with orthogonal transformations, so the covariance stays positive semidefinite. Streaming sensor data can update $R$ with Givens rotations, one sweep per new row.
• Overparameterized models. Gradient descent on $\|Ax - b\|^2$ from $x_0 = 0$ moves only along gradients $A^\top(Ax - b) \in C(A^\top)$, so with a small enough step it converges to the minimum-norm solution $A^+b$: the simplest example of the implicit bias of gradient descent.

7. Common mistakes #

• Computing inv(A.T @ A) @ A.T @ b. It squares $\kappa$ and adds an explicit inverse. Use QR or lstsq.
• The wrong Householder sign. $v = x - \operatorname{sign}(x_1)\|x\|e_1$ cancels when $x$ is near the $e_1$ axis.
• Building the $m \times m$ matrix $H$. Apply $H$ as a rank-1 update instead.
• Testing a QR code only with $\|A - QR\|$. Gram–Schmidt passes this test even when $Q$ is far from orthogonal. Check $\|Q^\top Q - I\|$ too, and fix the signs with $D = \operatorname{diag}(\operatorname{sign} r_{kk})$ before comparing with a library QR.
• Computing every CGS coefficient from the original column and calling it MGS. MGS takes $r_{jk}$ from the partly reduced vector; this one change is the difference between $\varepsilon\kappa^2$ and $\varepsilon\kappa$.
• Plain QR on a rank-deficient matrix. Use pivoted QR, a truncated SVD or ridge.

8. Check yourself #

1. Give the Householder vector $v$ that maps $x = (3, 4)$ to a multiple of $e_1$, and $Hx$. Give the Givens $(c, s)$ for the same vector.
   Answer

   $x_1 > 0$, so $v = (3 + 5, 4) = (8, 4)$. $v^\top x = 40$, $v^\top v = 80$, so $Hx = x - v = (-5, 0)$. Givens: $r = 5$, $c = 3/5$, $s = 4/5$, and $G x = (5, 0)$.
2. A regression has $\kappa(A) = 10^5$ and a small residual. About how many correct digits do the normal equations and Householder QR give in double precision?
   Answer

   Normal equations: $\varepsilon\kappa^2 \approx 2 \times 10^{-6}$, about 5–6 digits. QR: $\varepsilon\kappa \approx 2 \times 10^{-11}$, about 10–11 digits.
3. Find the minimum-norm solution of $x_1 + x_3 = 1$, $x_2 + x_3 = 2$, and check it against $N(A)$.
   Answer

   $AA^\top = \begin{pmatrix} 2 & 1 \\ 1 & 2\end{pmatrix}$, $y = (AA^\top)^{-1}(1, 2) = (0, 1)$, $x^+ = A^\top y = (0, 1, 1)$. $N(A) = \operatorname{span}(1, 1, -1)$ and $(0, 1, 1)\cdot(1, 1, -1) = 0$. All solutions: $(0, 1, 1) + t(1, 1, -1)$, with $\|x\|^2 = 2 + 3t^2$, smallest at $t = 0$.
4. CGS and MGS give the same numbers in exact arithmetic. How do they compare with each other and with Householder in floating point, for $\kappa(A) = 10^8$?
   Answer

   Loss of orthogonality: CGS $\approx \varepsilon\kappa^2$, i.e. order 1 (orthogonality gone); MGS $\approx \varepsilon\kappa \approx 10^{-8}$; Householder $\approx 10^{-16}$, independent of $\kappa$. MGS is better than CGS but not as good as Householder.
5. Run Gram–Schmidt on $a_1 = (3, 4)$, $a_2 = (1, 0)$. Give $Q$ and $R$.
   Answer

   $r_{11} = 5$, $q_1 = (3, 4)/5$. $r_{12} = q_1^\top a_2 = 3/5$, $v_2 = (1, 0) - \tfrac35 \cdot \tfrac15(3, 4) = (16, -12)/25$, $r_{22} = \|v_2\| = 20/25 = 4/5$, $q_2 = (4, -3)/5$. So $Q = \tfrac15\begin{pmatrix} 3 & 4 \\ 4 & -3 \end{pmatrix}$ and $R = \begin{pmatrix} 5 & 3/5 \\ 0 & 4/5 \end{pmatrix}$. Check: $r_{11}r_{22} = 4 = |\det A|$.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. State the QR theorem for a matrix with independent columns.

   Show answer

2. Write one Gram–Schmidt step for column $k$.

   Show answer

3. Gram–Schmidt on $a_1 = (1, 1, 0)$, $a_2 = (1, 0, 1)$: what are $q_2$ and $r_{22}$?

   Show answer

4. Which vector $v$ defines the Householder reflection that maps $x$ to a multiple of $e_1$, and why this sign?

   Show answer

5. Give the Givens rotation that maps $(a, b)$ to $(r, 0)$.

   Show answer

6. Why do the normal equations $A^\top A\hat x = A^\top b$ hold at the least-squares solution?

   Show answer

7. How do you solve least squares with a QR factorization?

   Show answer

8. What is $\kappa(A^\top A)$, and what does it cost the normal equations?

   Show answer

9. What is the minimum-norm solution of $Ax = b$ for a wide $A$ with full row rank?

   Show answer

10. Why is Householder QR backward stable while Gram–Schmidt is not?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you fit linear regression in practice? Why not $(X^\top X)^{-1}X^\top y$?" — QR or an SVD solver (lstsq). The closed form squares $\kappa$ and fails near $\kappa \approx 10^8$. For small $\kappa$, the normal equations (or ridge with Cholesky) are fine and cheaper.
• "Gram–Schmidt or Householder?" — Gram–Schmidt builds orthogonality by subtraction (CGS $\approx\varepsilon\kappa^2$, MGS $\approx\varepsilon\kappa$). Householder reaches triangular form with reflections: backward stable, orthogonal to $O(\varepsilon)$. Gram–Schmidt (with reorthogonalization) stays useful when vectors arrive one at a time, as in Krylov methods.
• "What does np.linalg.lstsq return for a rank-deficient $X$?" — the minimum-norm least-squares solution; singular values below rcond·$\sigma_{\max}$ count as zero.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Write classical and modified Gram–Schmidt.
2. Write Householder QR with the stable sign choice, applying each reflection as a rank-1 update.
3. Write a Givens rotation with hypot.
4. Solve least squares by QR ($R\hat x = Q^\top b$) and compare its accuracy with the normal equations as $\kappa$ grows.
5. Compute the minimum-norm solution of an underdetermined system through the QR of $A^\top$.

The notebook then runs the stability shootout of the fourth figure with your own Gram–Schmidt and Householder codes, and shows the Läuchli matrix.

11. Watch and read #

16. Projection Matrices and Least SquaresGilbert Strang · MIT 18.06 Linear Algebra, lecture 16 · 48 min

17. Orthogonal Matrices and Gram-SchmidtGilbert Strang · MIT 18.06 Linear Algebra, lecture 17 · 49 min

9. Four Ways to Solve Least Squares ProblemsGilbert Strang · MIT 18.065 (2018) · 50 min

MATH426: Householder QRToby Driscoll · 15 min

• Driscoll and Braun, Fundamentals of Numerical Computation, free online, chapter 3: §3.2 The normal equations, §3.3 The QR factorization and §3.4 Computing QR factorizations.
• Trefethen and Bau, Numerical Linear Algebra (SIAM 1997), Lectures 7–8 (QR, Gram–Schmidt), 10 (Householder), 11 (least squares), 16 and 19 (stability of Householder QR and of least-squares algorithms).
• Gilbert Strang, MIT 18.065 (2018), lectures 9–11 (four ways to solve least squares, difficulties with $Ax = b$, minimum norm).
• Cleve Moler, "Compare Gram-Schmidt and Householder orthogonalization algorithms", Cleve's Corner (MathWorks blog, 2016): the same experiment as the fourth figure, in MATLAB.