30. Least squares approximation and orthogonal projections
§16.1 solved A\vec x=\vec b when a solution exists. Most systems that arise from real data don't have one — more equations than unknowns, with noisy or inconsistent measurements (§17.5 already flagged this as typical for tall matrices). This lesson uses §20.4's projection to find the best possible answer anyway, and it's the direct linear-algebra engine behind the statistics course's regression chapter.
The problem
A\vec x=\vec b has no solution exactly when \vec b\notin\operatorname{Col}(A) (§17.4). Rather than giving up, ask for the \vec x that makes A\vec x as close to \vec b as possible — minimizing \|\vec b-A\vec x\|.
The normal equations
The closest point in \operatorname{Col}(A) to \vec b is \operatorname{proj}_{\operatorname{Col}(A)}\vec b (§20.4's projection, applied to the specific subspace \operatorname{Col}(A)). The key geometric fact: the residual \vec r=\vec b-A\hat{\vec x} at the best-fit \hat{\vec x} must be orthogonal to all of \operatorname{Col}(A) — if it weren't, the component of \vec r still pointing into \operatorname{Col}(A) could be chased down further, shrinking \|\vec r\| more, contradicting that \hat{\vec x} was already optimal.
\vec r\perp\operatorname{Col}(A) means \vec r\in\operatorname{Col}(A)^\perp, and by §20.2's identity \operatorname{Col}(A)^\perp=\operatorname{Row}(A^T)^\perp=\operatorname{Null}(A^T) (§20.2 stated this for a matrix's own rows; here it's applied to A^T, whose rows are A's columns). So A^T\vec r=\vec0:
A^T(\vec b-A\hat{\vec x})=\vec0\ \Longrightarrow\ \boxed{A^TA\hat{\vec x}=A^T\vec b}
the normal equations — solvable by ordinary elimination (§16.1), now applied to the new system (A^TA)\hat{\vec x}=A^T\vec b instead of the original inconsistent one.
When is the solution unique?
A^TA is invertible exactly when A's columns are independent (a fact usable here without full proof: A^TA\vec x=\vec0\iff \|A\vec x\|^2=\vec x^TA^TA\vec x=0\iff A\vec x=\vec0 by positive definiteness of the dot product, so A^TA and A share a null space — trivial exactly when A's columns are independent, §17.1). When it is,
\hat{\vec x}=(A^TA)^{-1}A^T\vec b
— a single closed-form solution, and exactly the formula behind §20.5's namesake application, ordinary least-squares regression (statistics course §6.0's simple linear regression is the case where A's columns are a constant column and one predictor column; multiple regression is the identical formula with more columns).
Doing it in Python
import numpy as np
# An inconsistent system: fit y = c0 + c1*x through 4 noisy points
x = np.array([0., 1., 2., 3.])
y = np.array([1.1, 2.9, 4.8, 7.2])
A = np.column_stack([np.ones_like(x), x]) # columns: [1, x]
b = y
x_hat = np.linalg.solve(A.T @ A, A.T @ b)
print("normal equations solution: c0={:.4f}, c1={:.4f}".format(*x_hat))
x_hat_direct = np.linalg.lstsq(A, b, rcond=None)[0]
print("np.linalg.lstsq: c0={:.4f}, c1={:.4f}".format(*x_hat_direct))
normal equations solution: c0=0.9700, c1=2.0200
np.linalg.lstsq: c0=0.9700, c1=2.0200
Confirming the residual is orthogonal to \operatorname{Col}(A):
import numpy as np
x = np.array([0., 1., 2., 3.])
y = np.array([1.1, 2.9, 4.8, 7.2])
A = np.column_stack([np.ones_like(x), x])
x_hat = np.linalg.solve(A.T @ A, A.T @ y)
residual = y - A @ x_hat
print("residual =", [round(r, 4) for r in residual.tolist()])
print("A^T @ residual =", [round(v, 10) or 0.0 for v in (A.T @ residual).tolist()], " (should be zero)")
print("residual . column0 =", round(float(residual @ A[:, 0]), 10))
print("residual . column1 =", round(float(residual @ A[:, 1]), 10))
residual = [0.13, -0.09, -0.21, 0.17]
A^T @ residual = [0.0, 0.0] (should be zero)
residual . column0 = 0.0
residual . column1 = 0.0
Worked example
Find the least-squares solution to A\vec x=\vec b for A=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}, \vec b=(1,2,2).
A^TA=\begin{pmatrix}1&1&1\\0&1&2\end{pmatrix}\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}=\begin{pmatrix}3&3\\3&5\end{pmatrix}
A^T\vec b=\begin{pmatrix}1&1&1\\0&1&2\end{pmatrix}\begin{pmatrix}1\\2\\2\end{pmatrix}=\begin{pmatrix}5\\6\end{pmatrix}
Solve \begin{pmatrix}3&3\\3&5\end{pmatrix}\hat{\vec x}=\begin{pmatrix}5\\6\end{pmatrix}: R_2\leftarrow R_2-R_1 gives \begin{pmatrix}3&3\\0&2\end{pmatrix}\hat{\vec x}=\begin{pmatrix}5\\1\end{pmatrix}\Rightarrow\hat x_2=\tfrac12; 3\hat x_1+3(\tfrac12)=5\Rightarrow\hat x_1=\tfrac76.
\boxed{\hat{\vec x}=\left(\tfrac76,\tfrac12\right)}
Sanity check. A\hat{\vec x}=\begin{pmatrix}7/6\\7/6+1/2\\7/6+1\end{pmatrix}=\begin{pmatrix}7/6\\5/3\\13/6\end{pmatrix}, residual \vec r=\vec b-A\hat{\vec x}=(1-\tfrac76,\,2-\tfrac53,\,2-\tfrac{13}6)=(-\tfrac16,\tfrac13,-\tfrac16). Check A^T\vec r=\vec0: column 1 of A is all 1s, so row 1 of A^T\vec r is just \vec r's entries summed: -\tfrac16+\tfrac13-\tfrac16=0 ✓. Column 2 of A is (0,1,2): 0(-\tfrac16)+1(\tfrac13)+2(-\tfrac16)=\tfrac13-\tfrac13=0 ✓ — the residual is genuinely orthogonal to both columns of A, exactly as the normal equations guarantee.
Your turn
1. Why does an inconsistent system have no exact solution but always have a least-squares solution (when A's columns are independent)?
2. What geometric object is A\hat{\vec x}, in relation to \vec b and \operatorname{Col}(A)?
3. True or false: if A\vec x=\vec b actually has an exact solution, the least-squares solution agrees with it.
Solutions
1. "No exact solution" means \vec b\notin\operatorname{Col}(A) — but \operatorname{proj}_{\operatorname{Col}(A)}\vec b (§20.4) is always defined, for any \vec b, since projection onto a subspace never requires the target to already be in it. The normal equations are just the algebraic route to computing that projection's coordinates, so a least-squares solution exists unconditionally (given independent columns), even though an exact one generally does not.
2. A\hat{\vec x}=\operatorname{proj}_{\operatorname{Col}(A)}\vec b — the closest point in \operatorname{Col}(A) to \vec b, i.e. exactly §20.4's orthogonal projection of \vec b onto the column space.
3. True. If \vec b\in\operatorname{Col}(A) exactly, then \vec b is already its own closest point in \operatorname{Col}(A) (distance 0 is unbeatable), so \operatorname{proj}_{\operatorname{Col}(A)}\vec b=\vec b and the residual is \vec0 — the normal equations reduce to solving A\vec x=\vec b exactly, and least squares is a strict generalization of ordinary solving, never a different answer when an exact one exists.
Check yourself in code
Find the least-squares solution to A\vec x=\vec b for A=\begin{pmatrix}1&1\\1&2\\1&3\end{pmatrix}, \vec b=(2,3,5), via the normal equations.
Print exactly this:
A^T A =
[3, 6]
[6, 14]
A^T b = [10, 23]
x_hat = [0.3333, 1.5]
import numpy as np
A = np.array([[1., 1.], [1., 2.], [1., 3.]])
b = np.array([2., 3., 5.])
AtA = A.T @ A
Atb = A.T @ b
print("A^T A =")
for row in AtA.tolist():
print([int(v) for v in row])
print("A^T b =", [int(v) for v in Atb.tolist()])
# solve AtA x_hat = Atb and print x_hat rounded to 4 decimal places
import numpy as np
A = np.array([[1., 1.], [1., 2.], [1., 3.]])
b = np.array([2., 3., 5.])
AtA = A.T @ A
Atb = A.T @ b
print("A^T A =")
for row in AtA.tolist():
print([int(v) for v in row])
print("A^T b =", [int(v) for v in Atb.tolist()])
x_hat = np.linalg.solve(AtA, Atb)
print("x_hat =", [round(v, 4) for v in x_hat.tolist()])
When A\vec x=\vec b has no solution, the least-squares solution \hat{\vec x} minimizes \|\vec b-A\vec x\| by making the residual orthogonal to \operatorname{Col}(A), giving the normal equations A^TA\hat{\vec x}=A^T\vec b — solvable uniquely whenever A's columns are independent, and exact whenever \vec b already was reachable.
Next: a numerically better way to compute this same projection — QR decomposition, which sidesteps forming A^TA at all.