Real data gives you more equations than unknowns, and no exact solution. Least squares is the principled answer: project the target onto what your model can actually reach. This is linear regression, derived twice — geometrically and by calculus — and shown to be the same theorem.
We have \(n\) data points and \(d\) features, and we want weights \(\vw\) such that predictions match targets exactly:
\[\mX\vw = \vy, \qquad \mX\in\RR^{n\times d},\ \vw\in\RR^{d},\ \vy\in\RR^{n}.\]
In any realistic setting \(n \gg d\): thousands of houses, a handful of features. That is thousands of equations in a handful of unknowns — an overdetermined system. By Lecture 9's column picture, \(\mX\vw\) ranges over the column space \(\mathcal{C}(\mX)\), a subspace of \(\RR^n\) of dimension at most \(d\). Meanwhile \(\vy\) is an arbitrary point of \(\RR^n\). A \(d\)-dimensional subspace inside an \(n\)-dimensional space is a vanishingly thin sliver, so \(\vy\) is essentially never in it.
Picture \(n=3\), \(d=2\): the reachable predictions form a plane through the origin in three-dimensional space, and \(\vy\) is a point floating somewhere off that plane. No \(\vw\) can hit it. So we stop asking for an exact hit and ask a better question: which point of the plane is closest to \(\vy\)? That is the projection question from Lecture 8, now with a subspace instead of a line.
\[ \hat{\vw} \;=\; \argmin_{\vw\in\RR^d}\ \|\mX\vw-\vy\|^2 \;=\; \argmin_{\vw}\ \sum_{i=1}^{n}\big(\vx_i^\top\vw - y_i\big)^2. \]
The vector \(\vr = \vy-\mX\hat{\vw}\) is the residual, and \(\hat{\vy}=\mX\hat{\vw}\) is the vector of fitted values.
This is Lecture 8's projection argument with a subspace in place of a line. The logic is identical; only the bookkeeping grows.
Now the same result by minimising, which is how you would get there without the geometry — and which reassures us that "closest point" and "minimum squared error" really are the same problem.
np.linalg.inv(X.T @ X) will either raise
LinAlgError or, worse, return an enormous garbage matrix built from round-off. Only routines
built on the SVD — np.linalg.lstsq, pinv — degrade gracefully, and they do so by
silently choosing the minimum-norm solution among the infinitely many. That is a reasonable choice, but it
is a choice, and you should know it is being made for you.Add a penalty on the size of the weights:
\[J_\lambda(\vw) = \|\mX\vw-\vy\|^2 + \lambda\|\vw\|^2, \qquad \lambda>0.\]
Differentiating gives \(2\mX^\top\mX\vw - 2\mX^\top\vy + 2\lambda\vw = \vzero\), so
\[\boxed{\;\hat{\vw}_{\text{ridge}} = \big(\mX^\top\mX + \lambda\mI\big)^{-1}\mX^\top\vy\;}\]
The matrix \(\mX^\top\mX+\lambda\mI\) is invertible for every \(\lambda>0\), even when \(\mX^\top\mX\) is singular. Here is the one-line reason, which we will make fully rigorous in Lecture 11: for any \(\vv\neq\vzero\),
\[\vv^\top\big(\mX^\top\mX+\lambda\mI\big)\vv = \|\mX\vv\|^2 + \lambda\|\vv\|^2 \ \ge\ \lambda\|\vv\|^2 \ > \ 0,\]
so no nonzero vector is annihilated, so the null space is trivial, so the matrix is invertible. The \(\lambda\mI\) term puts a floor under every eigenvalue.
Three things this figure is showing you at once, all of which recur throughout Phase B:
weight_decay parameter in every deep learning
optimiser is \(\lambda\) here, and adding \(\frac{\lambda}{2}\|\vw\|^2\) to the loss adds \(-\eta\lambda\vw\)
to each gradient step — literally decaying the weights.Why is inv worse? Two reasons. First, forming \(\mX^\top\mX\) squares the
condition number: if \(\mX\) has condition number \(10^{8}\) — quite ordinary — then \(\mX^\top\mX\) has
\(10^{16}\), which exhausts float64 entirely. lstsq works on \(\mX\) directly via the SVD and never
forms that product. Second, \(\mathbf{A}^{-1}\vb\) computed as "invert, then multiply" does more arithmetic and
accumulates more error than solving the system directly. The rule never invert, always solve is close
to universal in numerical linear algebra.
Fit \(y = w_1 x + w_0\) by least squares to the three points \((1,2),(2,2),(3,4)\), by hand.
1. With the bias column first, \(\mX=\begin{bmatrix}1&1\\1&2\\1&3\end{bmatrix}\), \(\vy=\begin{bmatrix}2\\2\\4\end{bmatrix}\), \(\vw=\begin{bmatrix}w_0\\w_1\end{bmatrix}\).
2. \(\mX^\top\mX=\begin{bmatrix}3&6\\6&14\end{bmatrix}\) (entries: \(n=3\), \(\sum x=6\), \(\sum x^2=1+4+9=14\)). \(\mX^\top\vy=\begin{bmatrix}8\\ 2+4+12\end{bmatrix}=\begin{bmatrix}8\\18\end{bmatrix}\).
3. \(3w_0+6w_1=8\) and \(6w_0+14w_1=18\). Multiply the first by 2 and subtract: \(2w_1=2\), so \(w_1=1\). Then \(3w_0=8-6=2\), so \(w_0=\tfrac23\).
Cross-check with Lecture 7's formula: \(\bar x=2\), \(\bar y=\tfrac83\); \(\sum(x-\bar x)(y-\bar y) = (-1)(-\tfrac23)+(0)(-\tfrac23)+(1)(\tfrac43)=2\); \(\sum(x-\bar x)^2=2\); so \(\hat w_1=1\) ✓ and \(\hat w_0=\tfrac83-2=\tfrac23\) ✓.
4. Fitted values: \(\tfrac23+1,\ \tfrac23+2,\ \tfrac23+3 = \tfrac53,\tfrac83,\tfrac{11}3\).
Residuals: \(2-\tfrac53=\tfrac13\), \(2-\tfrac83=-\tfrac23\), \(4-\tfrac{11}3=\tfrac13\).
\(\sum r_i = \tfrac13-\tfrac23+\tfrac13=0\) ✓ (orthogonality to the bias column).
\(\sum x_ir_i = \tfrac13-\tfrac43+1=0\) ✓ (orthogonality to the \(x\) column).
Note the first check is where "residuals sum to zero" comes from — it holds precisely because there is a column of ones. Drop the intercept and it stops being true.
Define \(\mP=\mX(\mX^\top\mX)^{-1}\mX^\top\), the projection matrix (or "hat matrix", since \(\hat{\vy}=\mP\vy\)). Prove:
1. Write \(\mS=(\mX^\top\mX)^{-1}\).
\[\mP^2 = \mX\mS\mX^\top\,\mX\mS\mX^\top = \mX\mS\,\underbrace{(\mX^\top\mX)\mS}_{=\,\mI}\,\mX^\top = \mX\mS\mX^\top = \mP.\]Geometrically: projecting a vector onto a subspace lands it in the subspace; projecting again moves nothing, because it is already there. Idempotence is "projection" in algebraic form.
2. \(\mP^\top = (\mX\mS\mX^\top)^\top = \mX\mS^\top\mX^\top\). And \(\mS^\top = ((\mX^\top\mX)^{-1})^\top = ((\mX^\top\mX)^\top)^{-1} = (\mX^\top\mX)^{-1}=\mS\), using that transpose and inverse commute and that \(\mX^\top\mX\) is symmetric. So \(\mP^\top=\mP\). ∎
(Symmetry is what makes it an orthogonal projection rather than an oblique one.)
3. \((\mI-\mP)\mX = \mX - \mX\mS\mX^\top\mX = \mX - \mX\mS(\mX^\top\mX) = \mX-\mX\mI
= \mathbf{0}\). ∎
Meaning: \(\mI-\mP\) is the projection onto the orthogonal complement, and it annihilates everything
in the column space. The residual is \(\vr=(\mI-\mP)\vy\), which is why \(\mX^\top\vr=\vzero\).
4. Use the cyclic property of the trace, \(\tr(\mA\mB)=\tr(\mB\mA)\):
\[\tr(\mP) = \tr\big(\mX\mS\mX^\top\big) = \tr\big(\mS\mX^\top\mX\big) = \tr\big((\mX^\top\mX)^{-1}(\mX^\top\mX)\big) = \tr(\mI_d) = d. \qquad\blacksquare\]Why anyone cares: \(\tr(\mP)\) counts the effective degrees of freedom used by the fit. For ordinary least squares it is exactly \(d\), one per parameter. For ridge it is \(\tr(\mX(\mX^\top\mX+\lambda\mI)^{-1}\mX^\top) = \sum_j \frac{\sigma_j^2}{\sigma_j^2+\lambda} < d\) — the penalty buys you a model that behaves as if it had fewer parameters, continuously tunable by \(\lambda\). This is the precise sense in which regularisation reduces model complexity, and it is the basis of AIC, generalised cross-validation, and the whole effective-capacity literature.
Show that least squares with squared error is exactly maximum likelihood under the model
\[ y_i = \vx_i^\top\vw + \epsilon_i, \qquad \epsilon_i \sim \mathcal{N}(0,\sigma^2) \text{ independent}. \]You may use that the Gaussian density is \(p(\epsilon)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\big(-\frac{\epsilon^2}{2\sigma^2}\big)\).
What does this tell you about when squared error is the wrong loss?
The likelihood of the whole dataset, using independence:
\[ L(\vw) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\left(-\frac{(y_i-\vx_i^\top\vw)^2}{2\sigma^2}\right).\]Take logs — legal and argmax-preserving, by Lecture 2's theorem:
\[ \log L(\vw) = -\frac{n}{2}\log(2\pi\sigma^2) \;-\; \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-\vx_i^\top\vw)^2.\]The first term does not involve \(\vw\). Maximising over \(\vw\) therefore means maximising \(-\frac{1}{2\sigma^2}\|\vy-\mX\vw\|^2\), i.e. minimising \(\|\vy-\mX\vw\|^2\) — the least-squares objective, with \(\frac{1}{2\sigma^2}\) an irrelevant positive constant. ∎
When squared error is wrong. The derivation used three assumptions, and squared error is the wrong loss exactly when one of them fails:
And for classification, \(y\in\{0,1\}\) is not Gaussian at all. The corresponding derivation with a Bernoulli likelihood produces cross-entropy, exactly as in Lecture 2. Every standard loss function is a negative log-likelihood for some noise model; choosing a loss is choosing a noise model, whether or not you realise it.
Coding. Investigate the bias–variance effect of ridge empirically.
import numpy as np
rng = np.random.default_rng(0)
n, d, sigma, reps = 30, 20, 1.0, 200
w_true = rng.normal(size=d)
X = rng.normal(size=(n, d)) # fixed design
Xt = rng.normal(size=(2000, d)) # test design
yt = Xt @ w_true # noiseless test target
lams = np.logspace(-3, 3, 40)
rows = []
for lam in lams:
W = np.empty((reps, d))
for r in range(reps):
y = X @ w_true + sigma*rng.normal(size=n)
W[r] = np.linalg.solve(X.T@X + lam*np.eye(d), X.T@y)
bias2 = np.sum((W.mean(0) - w_true)**2)
var = np.sum(W.var(0))
test = np.mean([(np.mean((Xt@w - yt)**2)) for w in W])
rows.append((lam, bias2, var, test))
print(f"lam={lam:9.4f} bias^2={bias2:9.4f} var={var:9.4f} test MSE={test:9.4f}")
best = min(rows, key=lambda r: r[3])
print("\\nbest lambda:", best[0], " test MSE:", best[3])What you will see. With \(n=30\) and \(d=20\), the design is badly under-determined and the unregularised estimator is terrible.
The answer to part 4: the optimum coincides with neither. Minimum bias is at \(\lambda=0\); minimum variance is at \(\lambda=\infty\) (where the estimator is the constant \(\vzero\), with zero variance and maximal bias). The best model deliberately accepts some bias to buy a larger reduction in variance. That trade is the single most important idea in statistical learning, and it is why an unbiased estimator is not automatically a good one.
Try re-running with \(n=500\): the variance term shrinks by roughly \(n\), the U-shape flattens, and the optimal \(\lambda\) moves towards zero. More data reduces the need for regularisation — which is the same statement as "regularisation substitutes for data you do not have".
Next: we have been saying "\(\mX^\top\mX\) is positive semi-definite" and "\(\lambda\mI\) puts a floor under the eigenvalues" on credit. Lecture 11 introduces eigenvalues properly, and with them the spectral theorem, quadratic forms, condition numbers, and the multi-dimensional learning-rate bound.