Phase A gave you the normal equations. This lecture turns them into a model you can actually use — basis functions, feature design, coefficient interpretation, residual diagnostics — and is honest about the four ways a linear fit lies to you.
Phase A Lecture 10 derived, twice, the solution to
\[\hat{\vw} = \argmin_{\vw}\ \|\mX\vw-\vy\|^2 \qquad\Longrightarrow\qquad \mX^\top\mX\,\hat{\vw} = \mX^\top\vy,\]
geometrically as the projection of \(\vy\) onto the column space of \(\mX\), and by calculus as the unique minimiser of a strictly convex quadratic. Nothing about that derivation is repeated here. What we do now is make it a model: decide what goes into \(\mX\), what the numbers coming out mean, and how to tell when the whole thing is inappropriate.
The single most common misconception about linear regression is that it can only fit straight lines. It cannot fit anything but a straight line in \(\vw\) — and that is a completely different restriction.
Choose functions \(\phi_1,\dots,\phi_m:\mathcal{X}\to\RR\) and model
\[ f(\vx) = \sum_{j=1}^{m} w_j\,\phi_j(\vx) = \vw^\top\boldsymbol\phi(\vx). \]The \(\phi_j\) are basis functions or features. They may be arbitrarily nonlinear in \(\vx\). Build the design matrix \(\mX_{ij}=\phi_j(\vx_i)\) and everything from Phase A Lecture 10 applies unchanged, because the objective is still quadratic in \(\vw\).
So a "linear" model can fit a sine wave, a step function, or a 12th-degree polynomial. What it cannot do is learn the basis — you choose \(\boldsymbol\phi\) up front, and that choice is the inductive bias. (Learning the basis is precisely what a neural network does, and it is the one structural difference between a linear model and a deep one.)
Four constructions cover most of what you will need.
inv fails or returns garbage; the coefficients are not
identifiable (you can add \(c\) to the intercept and subtract \(c\) from all three colour coefficients with
no change to any prediction); and interpreting any individual coefficient is meaningless.drop_first=True in pandas, drop='first' in scikit-learn's
OneHotEncoder). With all \(K\) and an intercept, the design is rank-deficient. Try it: build
the matrix and print np.linalg.matrix_rank.Both solve the same problem. The choice is about cost and conditioning.
\[\begin{array}{lll} \textbf{method} & \textbf{cost} & \textbf{use when} \\[3pt] \text{normal equations, } (\mX^\top\mX)^{-1}\mX^\top\vy & O(nd^2 + d^3) & \text{never — see below} \\[2pt] \texttt{np.linalg.solve} \text{ on } \mX^\top\mX & O(nd^2 + d^3) & d \text{ small, } \mX \text{ well conditioned} \\[2pt] \text{QR or SVD (} \texttt{lstsq} \text{)} & O(nd^2) & \textbf{the default} \\[2pt] \text{gradient descent} & O(nd) \text{ per step} & d \text{ large, or } n \text{ too large for memory} \\[2pt] \text{SGD} & O(Bd) \text{ per step} & n \text{ enormous, or streaming} \end{array}\]
The rule from Phase A Lecture 10 stands: never form
\((\mX^\top\mX)^{-1}\) explicitly. Forming \(\mX^\top\mX\) squares the condition number, so a design
matrix with \(\kappa=10^{8}\) — very ordinary for polynomial features — produces a Gram matrix with
\(\kappa=10^{16}\), which is all of float64's precision. lstsq works on \(\mX\) directly via the
SVD and never builds that product. The extra cost is a small constant factor; the accuracy difference can be
total.
The crossover to iterative methods is roughly at \(d\sim10^4\): above that, \(O(d^3)\) becomes unaffordable and \(O(nd)\) per gradient step wins. It also wins whenever \(\mX\) does not fit in memory, which is the usual situation at scale. Lecture 03 is about doing that well.
A fitted \(\hat w_j\) says: holding every other feature fixed, a one-unit increase in feature \(j\) is associated with a \(\hat w_j\) change in the prediction. Every clause of that sentence is load-bearing.
age and years_of_experience, no intervention increases one while
holding the other fixed, and the coefficient describes a comparison that never occurs.Least squares is the MLE under four assumptions (Phase A Lecture 14, and Lecture 04 here): the mean is linear in the features, the noise has constant variance, the noise is independent across observations, and — for inference, though not for the point estimate — the noise is Gaussian. Each assumption leaves a distinct fingerprint in the residuals \(\vr = \vy - \mX\hat{\vw}\) when it fails.
Note the last line of the "curvature" case. \(R^2\) stays respectable while the model is systematically wrong. Three things worth knowing about \(R^2 = 1 - \frac{\|\vr\|^2}{\|\vy-\bar y\|^2}\):
You fit \(\text{price} = w_0 + w_1\cdot\text{area}_{m^2} + w_2\cdot\text{bedrooms}\) and obtain \(\hat w_1 = 1800\), \(\hat w_2 = 12000\), with price in pounds.
1. \(\hat w_1 = 1800\) £/m²: holding bedroom count fixed, one additional square metre is associated with £1800 more. \(\hat w_2 = 12000\) £/bedroom: holding area fixed, one additional bedroom is associated with £12 000 more. Note that second one is already suspicious — adding a bedroom while holding total area constant means subdividing the space, which is a different intervention from what most people picture.
2. You cannot tell from these numbers. They are in different units, so their magnitudes are not comparable — 12000 > 1800 says nothing. The comparable quantity is the change per standard deviation: \(\hat w_j \times \text{sd}(x_j)\). If area has sd 50 m² and bedrooms has sd 1.1, then area contributes \(1800\times50 = £90{,}000\) per sd and bedrooms \(12000\times1.1 = £13{,}200\) per sd — area matters roughly seven times more. Standardising the features before fitting produces these numbers directly.
3. One square foot is \(0.0929\) m², so \(\hat w_1^{\text{new}} = 1800\times0.0929 = 167.2\) £/ft². \(R^2\), residuals, fitted values and predictions are all completely unchanged, because rescaling a column is an invertible linear change of basis that leaves the column space — and therefore the projection — identical (Phase A Lecture 9).
This is worth internalising: OLS is equivariant to rescaling features. Ridge and LASSO are not, because the penalty \(\|\vw\|^2\) depends on the units. That is precisely why standardisation is mandatory before regularising and merely convenient before plain OLS.
Show that if \(\mX\) contains a column of ones (an intercept), then the least-squares residuals satisfy \(\sum_i r_i = 0\), and hence the fitted line passes through \((\bar{\vx}, \bar y)\).
Then show \(R^2 = \operatorname{corr}(\vy,\hat{\vy})^2\) for a model with an intercept, and explain why that identity makes \(R^2\) blind to a systematic bias in the predictions.
Residuals sum to zero. The normal equations are \(\mX^\top\vr = \vzero\), i.e. the residual is orthogonal to every column. Applying this to the column of ones \(\mathbf{1}\):
\[\mathbf{1}^\top\vr = \sum_i r_i = 0. \qquad\blacksquare\]Dividing by \(n\): \(\bar y - \overline{\hat y} = 0\), so the mean fitted value equals the mean observed value, and by linearity \(\hat y(\bar{\vx}) = \bar{\vx}^\top\hat{\vw} = \bar y\). ∎
Note this holds only because of the intercept column. Drop it and residuals need not sum to zero — and \(R^2\) can then go negative.
\(R^2 = \operatorname{corr}(\vy,\hat{\vy})^2\). Write \(\tilde{\vy}=\vy-\bar y\mathbf{1}\) and \(\tilde{\hat{\vy}}=\hat{\vy}-\bar y\mathbf{1}\) (legitimate since the means agree). Since \(\vr\perp\mathcal{C}(\mX)\) and \(\hat{\vy}\in\mathcal{C}(\mX)\), Pythagoras gives \(\|\tilde{\vy}\|^2 = \|\tilde{\hat{\vy}}\|^2 + \|\vr\|^2\). Also \(\tilde{\vy}^\top\tilde{\hat{\vy}} = (\tilde{\hat{\vy}}+\vr)^\top\tilde{\hat{\vy}} = \|\tilde{\hat{\vy}}\|^2\), using \(\vr^\top\hat{\vy}=0\). Therefore
\[\operatorname{corr}(\vy,\hat{\vy})^2 = \frac{(\tilde{\vy}^\top\tilde{\hat{\vy}})^2} {\|\tilde{\vy}\|^2\|\tilde{\hat{\vy}}\|^2} = \frac{\|\tilde{\hat{\vy}}\|^4}{\|\tilde{\vy}\|^2\|\tilde{\hat{\vy}}\|^2} = \frac{\|\tilde{\hat{\vy}}\|^2}{\|\tilde{\vy}\|^2} = 1-\frac{\|\vr\|^2}{\|\tilde{\vy}\|^2} = R^2. \qquad\blacksquare\]Why this makes \(R^2\) blind to bias. Correlation is invariant to affine transformation of either argument: \(\operatorname{corr}(\vy, a\hat{\vy}+b) = \operatorname{corr}(\vy,\hat{\vy})\) for any \(a>0\). So a model that predicts half the true value plus a constant offset has the same correlation — and, on a fresh dataset where the intercept identity no longer holds, can report a high correlation-squared while every prediction is badly miscalibrated.
The practical warning: on training data with an intercept the two definitions of \(R^2\) coincide, so nothing goes wrong. On test data they diverge, and reporting \(\operatorname{corr}(y_{\text{test}},\hat y_{\text{test}})^2\) — which many people do, sometimes without realising — hides systematic over- or under-prediction. Report test MSE or the proper \(1-\text{SS}_{\text{res}}/\text{SS}_{\text{tot}}\) computed with the training mean, and plot predicted against actual with a \(y=x\) line.
Fit a Gaussian RBF model \(f(x)=\sum_{j=1}^{m}w_j\exp(-(x-c_j)^2/2s^2)\).
1. Grid over \(m\in\{5,10,20,40\}\) and \(s\in\{0.05,0.2,0.5,2\}\) with centres on an even grid, evaluating held-out MSE for each cell. You will find test error varies by orders of magnitude across \(s\) and comparatively little across \(m\) once \(m\) is large enough to tile the range. The width dominates.
The reason is that \(s\) controls the effective capacity — how smooth the achievable functions are — while \(m\) mostly controls resolution. Adding more narrow bumps than the data can support just produces near-duplicate columns, which regularisation handles. Adding a wrong width cannot be fixed by anything downstream.
2. \(s\) much smaller than the spacing: the bumps do not overlap, so between centres every basis function is essentially zero and the fit collapses to zero there — a row of spikes on a flat floor. The columns of \(\mX\) become nearly orthogonal and the model interpolates each centre independently, i.e. it memorises. \(s\) much larger than the range: every bump is nearly constant across the data, so all columns are nearly identical, \(\mX^\top\mX\) is horribly ill-conditioned, and the model can only express something close to a constant. Underfitting, plus numerical trouble.
A reasonable default is \(s\) equal to one to two times the centre spacing, then tuned by cross-validation.
3. Because the basis is local and bounded. Outside the range of the centres, every \(\phi_j(x)\to0\), so \(f(x)\to0\) — the model returns to its baseline and says, in effect, "no information here". A polynomial basis is global and unbounded: \(x^{12}\) grows without limit, so the coefficients fitted to interior data produce enormous extrapolated values, as you saw in Phase A Lecture 14's capstone (a fit reaching \(9\times10^4\) on data of magnitude 2).
The general principle: a model extrapolates according to the shape of its basis, not according to the data. Reverting to zero is at least a predictable failure. This is also why Gaussian-process regression — the Bayesian relative of RBF models — reverts to the prior mean with growing uncertainty away from the data, which is exactly the honest behaviour you want.
Coding. Build a complete regression study on a dataset with a nonlinear truth and a categorical feature.
lstsq,
solve, and the coefficients.import numpy as np
rng = np.random.default_rng(1)
n = 500
x1 = rng.uniform(-3, 3, n); x2 = rng.uniform(-2, 2, n)
g = rng.integers(0, 4, n)
geff = np.array([0., 2.5, -1.5, 4.0])[g]
y = 3*np.sin(x1) + 0.5*x2**2 + geff + 0.6*rng.normal(size=n)
def dummies(g, K=4): return np.column_stack([(g==k).astype(float) for k in range(1,K)])
def rbf(x, cs, s): return np.column_stack([np.exp(-(x-c)**2/(2*s*s)) for c in cs])
cs = np.linspace(-3, 3, 12); s = 0.55
models = {
"linear": lambda: np.column_stack([np.ones(n), x1, x2, dummies(g)]),
"+ quadratic": lambda: np.column_stack([np.ones(n), x1, x2, x1**2, x2**2, dummies(g)]),
"+ RBF on x1": lambda: np.column_stack([np.ones(n), rbf(x1,cs,s), x2, x2**2, dummies(g)]),
"+ interaction": lambda: np.column_stack([np.ones(n), rbf(x1,cs,s), x2, x2**2,
dummies(g), (x1*x2)[:,None]]),
}
perm = rng.permutation(n); tr, te = perm[:375], perm[375:]
best = None
for name, build in models.items():
X = build()
w, *_ = np.linalg.lstsq(X[tr], y[tr], rcond=None)
a = np.mean((X[tr]@w - y[tr])**2); b = np.mean((X[te]@w - y[te])**2)
print(f"{name:15s} cols {X.shape[1]:3d} train {a:.4f} test {b:.4f}")
if best is None or b < best[0]: best = (b, name, X, w)
_, name, X, w = best
r = y[tr] - X[tr]@w; f = X[tr]@w
print(f"\\nbest: {name}")
print(f" residual mean {r.mean():.2e} (should be ~0 with an intercept)")
print(f" corr(|r|, fitted) {np.corrcoef(np.abs(r), f)[0,1]:+.4f} (funnel test)")
lo, hi = f < np.median(f), f >= np.median(f)
print(f" mean r low/high {r[lo].mean():+.4f} / {r[hi].mean():+.4f} (curvature test)")
print(f" residuals > 3 sd {(np.abs(r) > 3*r.std()).sum()} of {len(r)} (outlier test)")
Xg = np.column_stack([np.ones(n), x1, x2] + [(g==k).astype(float) for k in range(4)])
print(f"\\nall 4 dummies + intercept: shape {Xg.shape}, rank {np.linalg.matrix_rank(Xg)}")
wl, *_ = np.linalg.lstsq(Xg[tr], y[tr], rcond=None)
print("lstsq ->", np.round(wl, 3), f" ||w|| {np.linalg.norm(wl):.3f}")
try:
ws = np.linalg.solve(Xg[tr].T@Xg[tr], Xg[tr].T@y[tr])
print("solve ->", np.round(ws, 3), f" ||w|| {np.linalg.norm(ws):.3e}")
except np.linalg.LinAlgError as e:
print("solve -> LinAlgError:", e)2. Typical test MSE: linear ≈ 5.0, + quadratic ≈ 4.6, + RBF ≈ 0.40, + interaction ≈ 0.40. The jump happens at the RBF step, because \(\sin(x_1)\) is the dominant nonlinearity and no polynomial of degree 2 can capture it. The interaction term adds nothing — correctly, since the truth is additive. Note that its training MSE is slightly lower (an extra column always helps training error) while test MSE is flat: that is the smallest possible instance of overfitting, and exactly why you compare on test.
3. The best model should show residual mean at machine zero, \(\operatorname{corr}(|r|, \hat y)\) near zero (no funnel), similar residual means in the low and high halves of the fitted range (no curvature), and roughly \(0.3\%\) of points beyond 3 sd, as a Gaussian predicts. Any of those failing points at the corresponding pathology in the figure above.
4. With all four dummies plus an intercept, \(\rank = 5\) but there are 6 columns.
lstsq succeeds — it uses the SVD and silently returns the minimum-norm solution among
the infinitely many, which is a perfectly reasonable but arbitrary choice. solve either raises
LinAlgError or returns coefficients of magnitude \(10^{13}\) or more, built entirely from
round-off.
The important observation: the predictions are identical in both cases. Only the
coefficients are broken. If you are predicting, you may never notice; if you are interpreting, you will report
nonsense. That is the difference between a model that works and a model you understand, and it is worth
checking matrix_rank whenever you intend to read coefficients.
Next: the closed form runs out at around \(d\sim10^4\). Lecture 03 replaces it with gradient descent done properly — convergence rates from the condition number, stochastic gradients, momentum, and Adam — which is the optimiser every model from here on will use.