Week 8: Least Squares and QR

Least Squares Fitting#

Week 7 showed that when $b$ does not lie in the column space $C(A)$, the equation $Ax = b$ has no exact solution. The closest you can get is the orthogonal projection of $b$ onto $C(A)$: find $\hat{x}$ so that

minimizes the residual length $\|b - Ax\|$ over all $x$. That is the least-squares problem

The geometry is the same as last week: the residual $e = b - A\hat{x}$ must be orthogonal to every column of $A$, which forces the normal equations

When $A$ is $m \times n$ with $m > n$ and full column rank $n$, the Gram matrix $A^\top A$ is $n \times n$, symmetric, and invertible, so a unique least-squares solution exists:

The fitted values are $p = A\hat{x}$; the residual $e = b - p$ lives in the left null space $N(A^\top)$.

When this shows up. Overdetermined systems are the everyday case: more measurements than parameters. Fitting a line through noisy sensor readings, estimating robot pose from many landmarks, regressing a policy on a batch of trajectories — all are least squares.

Worked example: fit a line. Suppose we want $y = c + d t$ for three data points

Stack the model as $A x = b$ with $x = (d, c)^\top$ (slope first, then intercept):

The first column is the $t$-values; the second is the constant column of ones. No exact line hits all three points (check: the slopes between consecutive pairs differ). Form the normal equations:

Solve $A^\top A\, \hat{x} = A^\top b$:

Subtract the second row from the first after scaling: from $3\hat{d} + 3\hat{c} = 5$ we get $\hat{d} + \hat{c} = 5/3$. From $5\hat{d} + 3\hat{c} = 6$, eliminate $\hat{c}$ to find $\hat{d} = 1/2$, then $\hat{c} = 5/3 - 1/2 = 7/6$. So

The fitted values are $A\hat{x} = (7/6,\; 5/3,\; 13/6)^\top$, and the residual $b - A\hat{x} = (-1/6,\; 1/3,\; -1/6)^\top$ is orthogonal to both columns of $A$ (dot products zero). Among all lines, this one minimizes the sum of squared vertical errors.

Connection to Week 7. The normal equations are exactly the subspace-projection equations $A^\top(b - A\hat{x}) = 0$. Least squares is projection of $b$ onto $C(A)$, rephrased as “best fit in the 2-norm.”

Gram-Schmidt Orthogonalization#

Projections become cleanest when the basis is orthonormal: vectors $q_1, \ldots, q_n$ with

(unit length, pairwise orthogonal). The projection of $b$ onto $\operatorname{span}\{q_1, \ldots, q_n\}$ collapses to a sum of line projections:

with no matrix inverse. Gram-Schmidt turns any independent set $a_1, \ldots, a_n$ into such an orthonormal set.

Algorithm. Start with independent vectors $a_1, \ldots, a_n \in \mathbb{R}^m$.

1. $v_1 = a_1$, then $q_1 = v_1 / \|v_1\|$.
2. For $j = 2, \ldots, n$:
   • Subtract projections onto the previous $q_i$'s: $$v_j = a_j - \sum_{i=1}^{j-1} (q_i^\top a_j)\, q_i.$$
   • Normalize: $q_j = v_j / \|v_j\|$.

Each $v_j$ is the component of $a_j$ orthogonal to $\operatorname{span}\{q_1, \ldots, q_{j-1}\} = \operatorname{span}\{a_1, \ldots, a_{j-1}\}$. Independence of the $a_i$'s guarantees $v_j \neq 0$, so normalization is always defined.

Worked example: two vectors in $\mathbb{R}^3$. Take

First vector: $\|a_1\| = \sqrt{2}$, so

Second vector: the coefficient of $q_1$ is

Then

Norm: $\|v_2\| = \sqrt{(1/2)^2 + (-1/2)^2 + 1^2} = \sqrt{3/2} = \sqrt{6}/2$. So

Check: $q_1^\top q_2 = (1/\sqrt{12})(1 - 1 + 0) = 0$, and both have unit length. The pair $\{q_1, q_2\}$ is an orthonormal basis for $\operatorname{span}\{a_1, a_2\}$.

Failure mode. If the input vectors are linearly dependent, some $v_j$ is zero and you cannot normalize — Gram-Schmidt is ill-defined (or you drop that vector and reduce rank). Orthogonal or already unit-length independent vectors are fine; the algorithm simply skips the “subtract” step when coefficients are zero.

The QR Factorization#

Pack the independent columns of $A$ as $a_1, \ldots, a_n$ and the Gram-Schmidt outputs as $q_1, \ldots, q_n$. By construction, each $a_j$ is a linear combination of $q_1, \ldots, q_j$ only (earlier $q_i$'s appear with coefficient $q_i^\top a_j$; later ones do not appear). In matrix form that is the QROrthogonal-Upper (matrix factorization) factorization

where

• $Q$ is $m \times n$ with orthonormal columns ($Q^\top Q = I_n$),
• $R$ is $n \times n$ upper triangular.

Explicitly, the entries of $R$ are the Gram-Schmidt coefficients:

Column $j$ of $A = QR$ reads $a_j = r_{1j} q_1 + \cdots + r_{jj} q_j$, which is exactly the reverse of the Gram-Schmidt subtractions.

Example continued. With $a_1, a_2$ from above,

Check $R_{11} = \|a_1\| = \sqrt{2}$, $R_{12} = q_1^\top a_2 = 1/\sqrt{2}$, $R_{22} = \|v_2\| = \sqrt{6}/2$, and $R_{21} = 0$. Multiplying $QR$ recovers $A = [a_1\; a_2]$.

Full vs thin QR. When $A$ is $m \times n$ with $m \geq n$ and full column rank, the form above is the thin (or reduced) QR: $Q$ is $m \times n$, $R$ is $n \times n$. A full QR extends $Q$ to an $m \times m$ orthogonal matrix by adding an orthonormal basis for the left null space; the extra rows of $R$ are zero. For least squares the thin form is enough.

Why care. $QR$ is a stable way to encode “orthonormal basis for the column space” plus “change-of-basis coefficients.” Week 13 reuses orthonormal columns in the SVD construction; Gram-Schmidt is the constructive prototype for that idea.

Least Squares via QR#

Substitute $A = QR$ into the normal equations $A^\top A\, \hat{x} = A^\top b$:

Since $Q^\top Q = I$ and $R$ is invertible (upper triangular with positive diagonal entries from full column rank), multiply both sides by $(R^\top)^{-1}$:

That is a triangular system. Solve it by back-substitution — no need to form $A^\top A$ or invert it.

Why this is better numerically. Forming $A^\top A$ squares the condition number: $\kappa(A^\top A) = \kappa(A)^2$. If $A$ is already moderately ill-conditioned, $A^\top A$ can lose all accurate digits before you solve. Working with $R\hat{x} = Q^\top b$ keeps the condition number of order $\kappa(A)$, not $\kappa(A)^2$. Library routines (numpy.linalg.lstsq, np.linalg.qr followed by a triangular solve) use QR or related orthogonal factorizations for this reason.

Algorithm summary.

1. Compute thin QR: $A = QR$.
2. Form the projected right-hand side $c = Q^\top b$.
3. Solve $R\hat{x} = c$ by back-substitution.
4. Fitted values: $p = A\hat{x} = Q(Q^\top b)$; residual $e = b - p$ is orthogonal to $C(A)$.

For the line-fit design matrix from the first section, the same $\hat{x}$ appears whether you solve normal equations or use QR; the QR path is the one you should prefer on real data.

Knowledge check#

Browser lab#

Fit a noisy line with least squares, solved via QR.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
x = np.linspace(0, 10, 20)
y_true = 2.0 * x + 1.0
y = y_true + rng.normal(scale=2.0, size=x.shape)

A = np.column_stack([x, np.ones_like(x)])
Q, R = np.linalg.qr(A)
x_hat = np.linalg.solve(R, Q.T @ y)

print("slope, intercept =", x_hat)

fig, ax = plt.subplots(figsize=(6, 4))
ax.scatter(x, y, color="tab:blue", label="data")
ax.plot(x, A @ x_hat, color="tab:orange", label="least-squares fit")
ax.legend()
ax.set_title("Least squares line fit via QR")
plt.show()