Lecture 06: Linear Regression (1) – CS 189, Fall 2026
In this lecture we will explore the formulation of linear regression, basis functions, the design matrix, the error function and its minimization, the geometry of least squares, how to evaluate a fit, and regularization.
Reference: Bishop & Bishop, Deep Learning: Foundations and Concepts, §4.1.1–4.1.6, §4.2.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import pandas as pd
from plotly.subplots import make_subplots
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.datasets import load_diabetes, fetch_california_housing
import warnings
from sklearn.exceptions import ConvergenceWarning
warnings.filterwarnings("ignore", category=ConvergenceWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
import plotly.io as pio
pio.renderers.default = "notebook_connected"
np.random.seed(42)
# One colour vocabulary for the whole notebook, so the slides and the demo match.
C_DATA, C_FIT, C_ALT, C_RESID, C_SPAN = "#003262", "#FDB515", "#C4820E", "#D55E00", "#00553A"
plt.rcParams.update({"figure.figsize": (8, 5), "axes.grid": True, "grid.alpha": 0.3,
"font.size": 13, "axes.titlesize": 15, "axes.labelsize": 14})
1. The Model¶
A linear regression model predicts a scalar target $t$ as a linear combination of the input features. With $\mathbf{x} = (x_1, \dots, x_D)^\top \in \mathbb{R}^D$:
$$ y(\mathbf{x}, \mathbf{w}) = w_0 + w_1 x_1 + w_2 x_2 + \dots + w_D x_D $$- $w_0$ is the intercept (bias term),
- $w_1, \dots, w_D$ are the weights.
It is convenient to fold $w_0$ into the dot product by augmenting the input with a constant 1:
$$ \tilde{\mathbf{x}} = (1, x_1, \dots, x_D)^\top \qquad\Longrightarrow\qquad y(\mathbf{x}, \mathbf{w}) = \tilde{\mathbf{x}}^\top \mathbf{w} $$That augmentation is the thing to keep straight: $\mathbf{x}^\top\mathbf{w}$ has no intercept, $\tilde{\mathbf{x}}^\top\mathbf{w}$ does.
1.1 One input: slope and intercept¶
In [2]:
N = 50
x = np.random.rand(N) * 10
t = 2 * x + 1 + np.random.randn(N) * 2 # ground truth: w0 = 1, w1 = 2
w0_true, w1_true = 1.0, 2.0
x_grid = np.linspace(0, 10, 100)
y_grid = w0_true + w1_true * x_grid
fig, ax = plt.subplots()
ax.scatter(x, t, alpha=0.75, s=55, color=C_DATA, label=r"data $(x_n, t_n)$")
ax.plot(x_grid, y_grid, color=C_ALT, lw=3.5, ls="--",
label=rf"$y(x,\mathbf{{w}}) = {w0_true:.0f} + {w1_true:.0f}x$")
# Mark the intercept and draw the rise-over-run triangle.
ax.plot(0, w0_true, "o", ms=13, color=C_SPAN, zorder=5, label=rf"intercept $w_0 = {w0_true:.0f}$")
xa, xb = 2, 4
ya, yb = w0_true + w1_true * xa, w0_true + w1_true * xb
ax.plot([xa, xb], [ya, ya], "k--", lw=2.5)
ax.plot([xb, xb], [ya, yb], "k--", lw=2.5)
ax.text((xa + xb) / 2, ya - 1.6, rf"$\Delta x = {xb-xa}$", ha="center", va="top", fontsize=13)
ax.text(xb + 0.25, (ya + yb) / 2, rf"$\Delta y = {yb-ya:.0f}$", ha="left", va="center", fontsize=13)
ax.text(9.6, 2.0, rf"slope $w_1 = \Delta y / \Delta x = {w1_true:.0f}$", ha="right", fontsize=13)
ax.set(xlabel="$x$", ylabel="$t$ / $y(x,\\mathbf{w})$", title="Intercept and slope")
ax.legend(fontsize=11, loc="upper left")
plt.tight_layout(); plt.show()
1.2 More than one input: the hyperplane¶
With $D$ inputs the model traces out a hyperplane in $\mathbb{R}^{D+1}$. Rotate the figure below: every point on the red surface is the prediction for one $(x_1, x_2)$ pair.
In [3]:
N = 150
x1 = np.random.rand(N) * 10
x2 = np.random.rand(N) * 10
t2 = 5 + 2 * x1 + 3 * x2 + np.random.randn(N) * 1.5 # truth: [5, 2, 3]
X2 = np.column_stack([x1, x2])
model = LinearRegression().fit(X2, t2)
w1_hat, w2_hat = model.coef_
w0_hat = model.intercept_
print(f"true : t = 5.00 + 2.00*x1 + 3.00*x2")
print(f"fitted : t = {w0_hat:.2f} + {w1_hat:.2f}*x1 + {w2_hat:.2f}*x2")
In [4]:
g1, g2 = np.meshgrid(np.linspace(x1.min(), x1.max(), 12),
np.linspace(x2.min(), x2.max(), 12))
t_surf = w0_hat + w1_hat * g1 + w2_hat * g2
fig = go.Figure([
go.Scatter3d(x=x1, y=x2, z=t2, mode="markers", name="data",
marker=dict(size=5, color=C_DATA, opacity=0.85)),
go.Surface(x=g1, y=g2, z=t_surf, name="fitted hyperplane", showscale=False,
opacity=0.5, colorscale=[[0, C_RESID], [1, C_RESID]]),
])
fig.update_layout(
title="Least squares fit is a hyperplane in (x1, x2, t) space",
scene=dict(xaxis_title="x1", yaxis_title="x2", zaxis_title="t",
xaxis=dict(backgroundcolor="white", gridcolor="lightgray"),
yaxis=dict(backgroundcolor="white", gridcolor="lightgray"),
zaxis=dict(backgroundcolor="white", gridcolor="lightgray"),
bgcolor="white"),
margin=dict(l=0, r=0, b=0, t=40), font=dict(size=13), height=520)
fig.show()
2. Basis Functions¶
2.1 First, watch the straight line fail¶
Before introducing machinery, it is worth seeing the thing the machinery fixes. Below, the data comes from $t = \sin(5x) + \varepsilon$. A straight line cannot represent it — and no amount of fitting will help, because the problem is the model class, not the optimizer.
In [5]:
_rng_course = np.random.default_rng(189)
n = 200
x_s = np.sort(_rng_course.random(n) * 2 - 1)
t_s = np.sin(5 * x_s) + 0.1 * _rng_course.standard_normal(n)
x_dense = np.linspace(-1, 1, 400)
lin = LinearRegression().fit(x_s[:, None], t_s)
mse_lin = np.mean((t_s - lin.predict(x_s[:, None])) ** 2)
fig, ax = plt.subplots()
ax.scatter(x_s, t_s, s=45, alpha=0.7, color=C_DATA, label="data")
ax.plot(x_s, lin.predict(x_s[:, None]), color=C_ALT, lw=3,
label=f"best straight line (MSE = {mse_lin:.3f})")
ax.set(xlabel="$x$", ylabel="$t$", title="This is the best a straight line can do")
ax.legend(); plt.tight_layout(); plt.show()
print(f"MSE of the best straight line: {mse_lin:.4f}")
print(f"Variance of t : {t_s.var():.4f} <- the line explains almost nothing")
2.2 The fix: change the representation, not the model¶
We keep the model linear and replace the raw input by $M$ fixed non-linear functions of it:
$$ y(\mathbf{x}, \mathbf{w}) = w_0 + \sum_{j=1}^{M-1} w_j \phi_j(\mathbf{x}) = \mathbf{w}^\top \boldsymbol{\phi}(\mathbf{x}), \qquad \phi_0(\mathbf{x}) \equiv 1 $$Common families (Bishop §4.1.1):
| Family | $\phi_j(x)$ | Behaviour |
|---|---|---|
| Polynomial | $x^j$ | global — changing one $w_j$ moves the whole curve |
| Gaussian (RBF) | $\exp\!\big(-\tfrac{(x-\mu_j)^2}{2s^2}\big)$ | local — each $w_j$ affects a neighbourhood of $\mu_j$ |
| Sigmoidal | $\sigma\!\big(\tfrac{x-\mu_j}{s}\big)$, $\sigma(a)=\tfrac{1}{1+e^{-a}}$ | local step at $\mu_j$ |
| Fourier | $\sin(jx),\ \cos(jx)$ | periodic |
Note $s$ (the width) is a separate symbol from $\sigma$ (the sigmoid). Bishop uses $s$; mixing the two up is a common source of confusion.
In [6]:
xg = np.linspace(-1, 1, 400)
def phi_poly(x, j): return x ** j
def phi_rbf(x, mu, s): return np.exp(-((x - mu) ** 2) / (2 * s ** 2))
def phi_sigmoid(x, mu, s): return 1.0 / (1.0 + np.exp(-(x - mu) / s))
def phi_fourier(x, j): return np.sin(j * np.pi * x)
fig, axes = plt.subplots(1, 4, figsize=(15, 3.4), sharey=True)
for j in range(1, 5):
axes[0].plot(xg, phi_poly(xg, j), lw=2.5, label=f"$j={j}$")
for mu in np.linspace(-0.8, 0.8, 5):
axes[1].plot(xg, phi_rbf(xg, mu, 0.2), lw=2.5)
for mu in np.linspace(-0.8, 0.8, 5):
axes[2].plot(xg, phi_sigmoid(xg, mu, 0.1), lw=2.5)
for j in range(1, 5):
axes[3].plot(xg, phi_fourier(xg, j), lw=2.5, label=f"$j={j}$")
for ax, name in zip(axes, ["Polynomial $x^j$", "Gaussian (RBF)", "Sigmoidal", "Fourier $\\sin(j\\pi x)$"]):
ax.set(title=name, xlabel="$x$")
axes[0].set_ylabel(r"$\phi_j(x)$"); axes[0].legend(fontsize=9); axes[3].legend(fontsize=9)
plt.tight_layout(); plt.show()
2.3 Four models, one dataset¶
Each panel below fits the same linear least-squares machinery to the same data. The only thing that changes is $\Phi$.
Pay attention to what "linear" is doing here: LinearRegression is called four times, unchanged.
In [7]:
def design_matrix(x, kind, M=9):
"Build Phi for one basis family. Columns exclude the bias; sklearn fits the intercept."
x = np.asarray(x).ravel()
if kind == "identity":
return x[:, None]
if kind == "polynomial":
return np.column_stack([x ** j for j in range(1, M)])
if kind == "rbf":
mus = np.linspace(x.min(), x.max(), M - 1)
s = (mus[1] - mus[0])
return np.column_stack([np.exp(-((x - mu) ** 2) / (2 * s ** 2)) for mu in mus])
if kind == "fourier":
cols = []
for j in range(1, (M - 1) // 2 + 1):
cols += [np.sin(j * np.pi * x), np.cos(j * np.pi * x)]
return np.column_stack(cols)
raise ValueError(kind)
fig, axes = plt.subplots(2, 2, figsize=(12, 7.5))
for ax, kind in zip(axes.ravel(), ["identity", "polynomial", "rbf", "fourier"]):
Phi = design_matrix(x_s, kind)
Phi_dens = design_matrix(x_dense, kind)
m = LinearRegression().fit(Phi, t_s)
mse = np.mean((t_s - m.predict(Phi)) ** 2)
ax.scatter(x_s, t_s, s=25, alpha=0.5, color=C_DATA)
ax.plot(x_dense, m.predict(Phi_dens), color=C_ALT, lw=3)
ax.set(title=f"{kind} (M = {Phi.shape[1] + 1}, MSE = {mse:.4f})", xlabel="$x$", ylabel="$t$",
ylim=(-1.6, 1.6))
print(f"{kind:>11}: M = {Phi.shape[1] + 1:>2}, MSE = {mse:.4f}")
plt.tight_layout(); plt.show()
2.4 Linear in the parameters, not in the features¶
This is the point students most often get wrong, so let us make it a numerical claim rather than a slogan.
A model is linear in $\mathbf{w}$ if $y(\mathbf{x}, a\mathbf{w} + b\mathbf{v}) = a\,y(\mathbf{x}, \mathbf{w}) + b\,y(\mathbf{x}, \mathbf{v})$ for all $a, b$. Let's test it on the degree-8 polynomial model, whose features are wildly non-linear.
In [8]:
Phi = design_matrix(x_s, "polynomial") # highly non-linear in x
rng = np.random.default_rng(0)
w, v = rng.normal(size=Phi.shape[1]), rng.normal(size=Phi.shape[1])
a, b = 2.7, -0.4
lhs = Phi @ (a * w + b * v)
rhs = a * (Phi @ w) + b * (Phi @ v)
print("linear in the PARAMETERS w? max |lhs - rhs| =", np.abs(lhs - rhs).max())
# Now the same test in the INPUT x, at fixed w.
x_a, x_b = 0.3, -0.7
lhs_x = design_matrix([a * x_a + b * x_b], "polynomial") @ w
rhs_x = a * (design_matrix([x_a], "polynomial") @ w) + b * (design_matrix([x_b], "polynomial") @ w)
print("linear in the INPUT x? |lhs - rhs| =", float(np.abs(lhs_x - rhs_x)[0]))
3. The Design Matrix and the Error Function¶
Stack the $N$ training inputs row-wise into the $N \times (D+1)$ design matrix
$$ \Phi = \begin{pmatrix} \phi_0(\mathbf{x}_1) & \phi_1(\mathbf{x}_1) & \cdots & \phi_{D}(\mathbf{x}_1)\\ \phi_0(\mathbf{x}_2) & \phi_1(\mathbf{x}_2) & \cdots & \phi_{D}(\mathbf{x}_2)\\ \vdots & \vdots & \ddots & \vdots\\ \phi_0(\mathbf{x}_N) & \phi_1(\mathbf{x}_N) & \cdots & \phi_{D}(\mathbf{x}_N) \end{pmatrix} $$so that all $N$ predictions are the single matrix-vector product $\mathbf{y} = \Phi\mathbf{w}$, and the sum-of-squares error is
$$ E(\mathbf{w}) = \tfrac{1}{2}\sum_{n=1}^{N}\big(t_n - \mathbf{w}^\top\boldsymbol{\phi}(\mathbf{x}_n)\big)^2 = \tfrac{1}{2}\,\lVert \mathbf{t} - \Phi\mathbf{w} \rVert^2 $$$E(\mathbf{w})$ is non-negative, and zero only when every prediction hits its target exactly.
In [9]:
def build_Phi(x, D, bias=True):
"Polynomial design matrix WITH the bias column, written out explicitly."
x = np.asarray(x).ravel()
cols = [np.ones_like(x)] if bias else []
cols += [x ** j for j in range(1, D)]
return np.column_stack(cols)
Phi3 = build_Phi(x_s, D=4)
print("Phi shape (N x D+1):", Phi3.shape)
print(np.array2string(Phi3[:4], precision=3, suppress_small=True))
def E(w, Phi, t):
r = t - Phi @ w
return 0.5 * float(r @ r)
print("\nE(w) at w = 0 :", round(E(np.zeros(4), Phi3, t_s), 3))
print("E(w) at a random w:", round(E(np.array([0., 1., 0., 0.]), Phi3, t_s), 3))
3.1 What the error surface looks like¶
For the straight-line model there are only two parameters, so we can draw $E(w_0, w_1)$ directly. It is a quadratic bowl — one minimum, no local traps. That is exactly why a closed-form solution exists.
In [10]:
Phi_lin = build_Phi(x, D=2) # the section-1 data: N x 2, columns [1, x]
w0g = np.linspace(-6, 8, 200)
w1g = np.linspace(-1, 5, 200)
W0, W1 = np.meshgrid(w0g, w1g)
# E(w) = 0.5 ||t - Phi w||^2, evaluated on the whole grid at once.
R = t[None, None, :] - (W0[..., None] * Phi_lin[:, 0] + W1[..., None] * Phi_lin[:, 1])
Egrid = 0.5 * np.sum(R ** 2, axis=-1)
w_star = np.linalg.lstsq(Phi_lin, t, rcond=None)[0]
fig, ax = plt.subplots(figsize=(7.5, 5.5))
cs = ax.contourf(W0, W1, Egrid, levels=40, cmap="Blues_r")
ax.contour(W0, W1, Egrid, levels=20, colors="white", linewidths=0.6, alpha=0.6)
ax.plot(*w_star, marker="*", ms=22, color=C_ALT, mec="k", mew=1, zorder=5,
label=rf"$\hat{{\mathbf{{w}}}} = ({w_star[0]:.2f}, {w_star[1]:.2f})$")
ax.plot(w0_true, w1_true, marker="o", ms=10, color=C_RESID, mec="k", zorder=5,
label=rf"truth $= ({w0_true:.0f}, {w1_true:.0f})$")
ax.set(xlabel="$w_0$ (intercept)", ylabel="$w_1$ (slope)",
title="$E(\\mathbf{w})$ is a quadratic bowl")
ax.legend(); plt.colorbar(cs, ax=ax, label="$E(\\mathbf{w})$")
plt.tight_layout(); plt.show()
4. Minimizing the Error: the Normal Equations¶
Setting the gradient to zero,
$$ \nabla_{\mathbf{w}} E(\mathbf{w}) = -\Phi^\top(\mathbf{t} - \Phi\mathbf{w}) = \mathbf{0} \qquad\Longrightarrow\qquad \boxed{\;\Phi^\top\Phi\,\hat{\mathbf{w}} = \Phi^\top\mathbf{t}\;} $$these are the normal equations. When $\Phi^\top\Phi$ is invertible, $\hat{\mathbf{w}} = (\Phi^\top\Phi)^{-1}\Phi^\top\mathbf{t} = \Phi^{\dagger}\mathbf{t}$, where $\Phi^\dagger$ is the Moore–Penrose pseudo-inverse.
Three ways to compute it, in decreasing order of numerical virtue:
In [11]:
Phi5 = build_Phi(x_s, D=6)
w_solve = np.linalg.solve(Phi5.T @ Phi5, Phi5.T @ t_s) # form the normal equations
w_lstsq = np.linalg.lstsq(Phi5, t_s, rcond=None)[0] # QR / SVD on Phi directly
w_pinv = np.linalg.pinv(Phi5) @ t_s # explicit pseudo-inverse
print("solve(PhiT Phi, PhiT t):", np.array2string(w_solve, precision=4))
print("lstsq(Phi, t) :", np.array2string(w_lstsq, precision=4))
print("pinv(Phi) @ t :", np.array2string(w_pinv, precision=4))
print("\nE(w) at each:",
round(E(w_solve, Phi5, t_s), 6), round(E(w_lstsq, Phi5, t_s), 6), round(E(w_pinv, Phi5, t_s), 6))
4.1 The case where solve falls over¶
np.linalg.solve needs $\Phi^\top\Phi$ to be genuinely invertible. Duplicate a column — say a student accidentally includes the same feature twice — and it is not. Let's break it on purpose.
In [12]:
Phi_bad = np.column_stack([Phi5, Phi5[:, 1]]) # column 1 repeated: rank deficient
A, b = Phi_bad.T @ Phi_bad, Phi_bad.T @ t_s
print("Phi_bad shape :", Phi_bad.shape)
print("rank(Phi_bad) :", np.linalg.matrix_rank(Phi_bad), " <- should be 7, it is not")
print("cond(Phi^T Phi):", f"{np.linalg.cond(A):.3e}")
try:
w_bad = np.linalg.solve(A, b)
# It may not even raise. Check whether the answer actually solves the system.
print("solve -> did NOT raise. Residual of the normal equations:",
f"{np.linalg.norm(A @ w_bad - b):.3e}")
print(" returned w =", np.array2string(w_bad, precision=1))
except np.linalg.LinAlgError as err:
print("solve -> LinAlgError:", err)
w_ok = np.linalg.lstsq(Phi_bad, t_s, rcond=None)[0]
print("\nlstsq -> minimum-norm solution:", np.array2string(w_ok, precision=3))
print(" E(w) =", round(E(w_ok, Phi_bad, t_s), 6),
" (identical to the full-rank fit:", round(E(w_lstsq, Phi5, t_s), 6), ")")
Two things to notice. First, solve may not even complain — depending on rounding it either raises or quietly hands back a number, and the condition number is the only warning you get. Second, the fit is fine while the parameters are not unique: lstsq silently picks the shortest $\mathbf{w}$ among the infinitely many achieving that error. Hold on to both — they are the seed of the ill-conditioning story in Lecture 07.
5. The Geometry of Least Squares¶
Read $\Phi\mathbf{w}$ column-wise instead of row-wise:
$$ \Phi\mathbf{w} = w_0\,\Phi_{:,0} + w_1\,\Phi_{:,1} + \dots + w_{M-1}\,\Phi_{:,D} $$So the prediction vector $\mathbf{y} \in \mathbb{R}^N$ is a linear combination of the columns of $\Phi$. Whatever $\mathbf{w}$ we choose, $\mathbf{y}$ is trapped inside the $D+1$-dimensional subspace $\mathcal{S} = \text{span}(\Phi) \subset \mathbb{R}^N$ — but the target $\mathbf{t}$ generally is not.
Minimizing $\lVert \mathbf{t} - \Phi\mathbf{w}\rVert$ therefore means: find the point of $\mathcal{S}$ closest to $\mathbf{t}$. That point is the orthogonal projection of $\mathbf{t}$ onto $\mathcal{S}$, and the residual $\mathbf{e} = \mathbf{t} - \Phi\hat{\mathbf{w}}$ is perpendicular to every column of $\Phi$:
$$ \Phi^\top \mathbf{e} = \mathbf{0} $$which is the normal equations again, arrived at without any calculus.
5.1 Verify it on a tiny example you can hold in your head¶
In [13]:
# N = 3 observations, M = 2 basis functions -> a plane inside R^3.
Phi_t = np.array([[1.0, 0.0],
[1.0, 1.0],
[1.0, 2.0]])
t_t = np.array([1.0, 3.0, 2.0])
w_hat = np.linalg.lstsq(Phi_t, t_t, rcond=None)[0]
y_hat = Phi_t @ w_hat
e = t_t - y_hat
print("w_hat =", np.array2string(w_hat, precision=4))
print("y_hat =", np.array2string(y_hat, precision=4))
print("residual =", np.array2string(e, precision=4))
print("\nPhi^T e =", np.array2string(Phi_t.T @ e, precision=12), " <- zero: e is orthogonal to span(Phi)")
print("angle between e and column 0:", round(np.degrees(np.arccos(
Phi_t[:, 0] @ e / (np.linalg.norm(Phi_t[:, 0]) * np.linalg.norm(e)))), 4), "degrees")
# Pythagoras: ||t||^2 = ||y_hat||^2 + ||e||^2, because the two pieces are orthogonal.
print("\n||t||^2 =", round(t_t @ t_t, 6))
print("||y_hat||^2 + ||e||^2 =", round(y_hat @ y_hat + e @ e, 6))
5.2 The same picture, in 3D¶
In [14]:
SCALE, LINE_W, COEFF_RANGE, GRID_LINES, PAD_RATIO, BRACKET_FRAC = 2.0, 14, 3.0, 11, 0.20, 0.22
c1 = np.array([3.0, 0.4, 0.6]) * SCALE # first column of Phi
c2 = np.array([1.2, 2.4, -0.5]) * SCALE # second column of Phi
a1, a2, residual_height = 1.5, 0.9, 2.2 * SCALE
O = np.zeros(3)
Y = a1 * c1 + a2 * c2
def unit(v):
n = np.linalg.norm(v)
return v / n if n else v
u1 = unit(c1)
u2 = unit(c2 - (c2 @ u1) * u1)
n_hat = unit(np.cross(u1, u2))
t_vec = Y + residual_height * n_hat # t sits off the plane by construction
U, V = np.meshgrid(np.linspace(-COEFF_RANGE, COEFF_RANGE, 45),
np.linspace(-COEFF_RANGE, COEFF_RANGE, 45))
P = U[..., None] * c1 + V[..., None] * c2
plane = go.Surface(x=P[..., 0], y=P[..., 1], z=P[..., 2], opacity=0.22, showscale=False,
name="span(Phi)", surfacecolor=np.zeros_like(P[..., 0]),
colorscale=[[0, C_SPAN], [1, C_SPAN]], hoverinfo="skip")
grid = []
for s in np.linspace(-COEFF_RANGE, COEFF_RANGE, GRID_LINES):
for a, b in [(s * c1 - COEFF_RANGE * c2, s * c1 + COEFF_RANGE * c2),
(-COEFF_RANGE * c1 + s * c2, COEFF_RANGE * c1 + s * c2)]:
grid.append(go.Scatter3d(x=[a[0], b[0]], y=[a[1], b[1]], z=[a[2], b[2]], mode="lines",
line=dict(width=2, color="rgba(0,85,58,0.6)"),
showlegend=False, hoverinfo="skip"))
lines, cones = [], []
def arrow(start, end, name, color, width=LINE_W, head=1.6):
vec = end - start
L = np.linalg.norm(vec)
if L == 0:
return
tip_base = end - max(0.9 * head, 0.02 * L) * (vec / L)
lines.append(go.Scatter3d(x=[start[0], tip_base[0]], y=[start[1], tip_base[1]],
z=[start[2], tip_base[2]], mode="lines",
line=dict(width=width, color=color), name=name, hoverinfo="skip"))
cones.append(go.Cone(x=[end[0]], y=[end[1]], z=[end[2]], u=[vec[0]], v=[vec[1]], w=[vec[2]],
sizemode="absolute", sizeref=head, anchor="tip", showscale=False,
colorscale=[[0, color], [1, color]], showlegend=False))
arrow(O, c1, "Phi[:,0]", C_SPAN)
arrow(O, c2, "Phi[:,1]", C_SPAN)
arrow(O, Y, "y = Phi w", "black")
arrow(O, t_vec, "t", C_RESID, head=1.9)
arrow(Y, t_vec, "e = t - Phi w", C_ALT, width=LINE_W - 2, head=1.5)
# Right-angle bracket at the foot of the residual.
y_dir, e_dir = unit(Y - O), unit(t_vec - Y)
tick = BRACKET_FRAC * min(np.linalg.norm(c1), np.linalg.norm(c2)) * SCALE
base = Y - tick * y_dir
bracket = go.Scatter3d(x=[base[0], (base + tick * e_dir)[0], (Y + tick * e_dir)[0]],
y=[base[1], (base + tick * e_dir)[1], (Y + tick * e_dir)[1]],
z=[base[2], (base + tick * e_dir)[2], (Y + tick * e_dir)[2]],
mode="lines", line=dict(width=10, color="royalblue"),
showlegend=False, hoverinfo="skip")
pts = np.vstack([O, c1, c2, Y, t_vec])
mins, maxs = pts.min(axis=0), pts.max(axis=0)
pad = PAD_RATIO * float(np.max(maxs - mins))
# One shared cube so the right angle actually looks like a right angle.
lo, hi = float((mins - pad).min()), float((maxs + pad).max())
def label(pt, text, d=0.12 * SCALE):
return dict(x=pt[0] + d, y=pt[1] + d, z=pt[2] + d, text=text, showarrow=False,
bgcolor="rgba(255,255,255,0.85)", bordercolor="black")
fig = go.Figure(data=[plane, *grid, bracket, *lines, *cones])
fig.update_layout(
title="Least squares = orthogonal projection of t onto span(Phi)",
width=900, height=760,
scene=dict(xaxis=dict(title="", range=[lo, hi], showbackground=False, zeroline=False),
yaxis=dict(title="", range=[lo, hi], showbackground=False, zeroline=False),
zaxis=dict(title="", range=[lo, hi], showbackground=False, zeroline=False),
annotations=[label(c1, "Φ<sub>:,0</sub>"), label(c2, "Φ<sub>:,1</sub>"),
label(Y, "y = Φw"), label(t_vec, "t"), label((Y + t_vec) / 2, "e")],
aspectmode="cube", camera=dict(eye=dict(x=1, y=-1, z=1.4))),
legend=dict(x=0.02, y=0.98, bgcolor="rgba(255,255,255,0.75)"))
fig.show()
In [15]:
USE_CALIFORNIA = False
if USE_CALIFORNIA:
frame = fetch_california_housing(as_frame=True).frame
data = frame[["MedInc", "MedHouseVal"]].rename(columns={"MedInc": "x", "MedHouseVal": "y"})
data = data[(data["x"] < 10) & (data["y"] < 5)] # drop the censored top-coded block
rng_s = np.random.default_rng(7)
data = data.iloc[rng_s.choice(len(data), 450, replace=False)]
XLAB, YLAB, UNIT, SOURCE = ("Median income", "Median house value",
"$100k", "California housing")
else:
d = load_diabetes(as_frame=True, scaled=False)
data = pd.DataFrame({"x": d.data["bmi"].to_numpy(), "y": d.target.to_numpy()})
XLAB, YLAB, UNIT, SOURCE = ("Body mass index", "Disease progression after one year",
"progression units", "diabetes (bundled with scikit-learn)")
Xh = data[["x"]].to_numpy()
yh = data["y"].to_numpy()
Xh_tr, Xh_te, yh_tr, yh_te = train_test_split(Xh, yh, test_size=0.3, random_state=42)
house = LinearRegression().fit(Xh_tr, yh_tr)
pred_tr, pred_te = house.predict(Xh_tr), house.predict(Xh_te)
resid_tr = yh_tr - pred_tr
print(f"source : {SOURCE}")
print(f"n_train = {len(yh_tr)}, n_test = {len(yh_te)}")
print(f"fitted : y = {house.intercept_:.2f} + {house.coef_[0]:.3f} * x")
In [16]:
fig, axs = plt.subplots(1, 2, figsize=(13, 5))
axs[0].scatter(Xh_tr, yh_tr, s=40, alpha=0.6, color=C_DATA, edgecolor="k", linewidth=0.4)
xl = np.linspace(Xh_tr.min(), Xh_tr.max(), 200).reshape(-1, 1)
axs[0].plot(xl, house.predict(xl), lw=3, color=C_ALT, zorder=3, label="fitted line")
axs[0].set(xlabel=XLAB, ylabel=YLAB, title="Least-squares fit")
axs[0].legend(frameon=False)
axs[1].scatter(pred_tr, resid_tr, s=40, alpha=0.6, color=C_DATA, edgecolor="k", linewidth=0.4)
axs[1].axhline(0, color="k", lw=2, ls="--")
axs[1].set(xlabel="Fitted value", ylabel="Residual", title="Residuals vs fitted")
plt.tight_layout(); plt.show()
# Quantify the fan rather than eyeballing it: residual spread by fitted-value tercile.
edges = np.quantile(pred_tr, [0, 1/3, 2/3, 1.0])
for lo, hi, name in zip(edges[:-1], edges[1:], ["low", "mid", "high"]):
m = (pred_tr >= lo) & (pred_tr <= hi)
print(f"{name:>5} fitted values: residual std = {resid_tr[m].std():.1f}")
6.2 Metrics¶
$$ \text{MSE} = \frac{1}{N}\sum_n (t_n - y_n)^2 \qquad \text{RMSE} = \sqrt{\text{MSE}} \qquad R^2 = 1 - \frac{\sum_n (t_n - y_n)^2}{\sum_n (t_n - \bar{t})^2} $$- MSE is in squared units of the target.
- RMSE puts it back in the target's units.
- $R^2$ is unitless: it compares the model to the constant model $y = \bar{t}$. $R^2 = 0$ means no better than predicting the mean; $R^2 < 0$ means worse.
In [17]:
def metrics(t, y, baseline=None):
"Computed from the definitions, so the algebra on the slide is visible."
ss_res = np.sum((t - y) ** 2)
ss_tot = np.sum((t - (np.mean(t) if baseline is None else baseline)) ** 2)
mse = ss_res / len(t)
return mse, np.sqrt(mse), 1 - ss_res / ss_tot
mse, rmse, r2 = metrics(yh_te, pred_te)
print(f"Test MSE = {mse:.2f} (units: [{UNIT}]^2)")
print(f"Test RMSE = {rmse:.2f} (units: {UNIT} -> a typical miss)")
print(f"Test R^2 = {r2:.4f}")
print("\nsklearn agrees:", round(mean_squared_error(yh_te, pred_te), 2),
round(r2_score(yh_te, pred_te), 4))
bad = np.full_like(yh_te, yh_te.mean() + 1.5 * yh_te.std())
print(f"\nA constant-but-wrong predictor: R^2 = {metrics(yh_te, bad)[2]:.4f}")
6.3 $R^2$ as a ratio of areas¶
Square each residual and you get a literal square. $R^2$ is the fraction of the constant model's total square area that the regression model removes.
In [18]:
sub = np.random.default_rng(7).choice(len(Xh_tr), 12, replace=False)
xs, ts = Xh_tr[sub].ravel(), yh_tr[sub]
ys_model = house.predict(xs.reshape(-1, 1))
ys_const = np.full_like(ts, yh_tr.mean())
fig, axs = plt.subplots(1, 2, figsize=(13, 5.5), sharey=True)
for ax, ys, name in [(axs[0], ys_model, "regression model"),
(axs[1], ys_const, "constant model $\\bar{t}$")]:
area = 0.0
for xi, ti, yi in zip(xs, ts, ys):
r = ti - yi
ax.add_patch(plt.Rectangle((xi, min(ti, yi)), abs(r), abs(r),
facecolor=C_ALT, alpha=0.35, edgecolor=C_ALT))
ax.plot([xi, xi], [ti, yi], color=C_RESID, lw=1.5)
area += r ** 2
ax.scatter(xs, ts, s=70, color=C_DATA, zorder=5)
ax.plot(np.sort(xs), ys[np.argsort(xs)], color=C_SPAN, lw=3, zorder=4)
ax.set(xlabel=XLAB, title=f"{name}\ntotal area = {area:.0f}", aspect="equal")
axs[0].set_ylabel(YLAB)
plt.tight_layout(); plt.show()
ss_res = np.sum((ts - ys_model) ** 2); ss_tot = np.sum((ts - ys_const) ** 2)
print(f"R^2 = 1 - {ss_res:.0f}/{ss_tot:.0f} = {1 - ss_res/ss_tot:.3f}")
7. Model Complexity¶
Now turn the dial on model complexity, with a held-out test set this time.
In [19]:
i_tr, i_te = train_test_split(np.arange(n), test_size=0.55, random_state=0)
print(f"{'degree':>7} {'train MSE':>12} {'TEST MSE':>12} {'max |w_j|':>12}")
for d in [1, 3, 5, 9, 16]:
print(d)
P_tr, P_te = build_Phi(x_s[i_tr], d + 1), build_Phi(x_s[i_te], d + 1)
w = np.linalg.lstsq(P_tr, t_s[i_tr], rcond=None)[0]
print(f"{d:>7} {np.mean((t_s[i_tr] - P_tr @ w)**2):>12.5f} "
f"{np.mean((t_s[i_te] - P_te @ w)**2):>12.5f} {np.abs(w).max():>12.1f}")
In [20]:
ds = list(range(1, 17))
TR = []; TE = []
for d in ds:
P_tr, P_te = build_Phi(x_s[i_tr], d + 1), build_Phi(x_s[i_te], d + 1)
w = np.linalg.lstsq(P_tr, t_s[i_tr], rcond=None)[0]
TR.append(np.mean((t_s[i_tr] - P_tr @ w) ** 2))
TE.append(np.mean((t_s[i_te] - P_te @ w) ** 2))
fig, ax = plt.subplots(figsize=(8.5, 5))
ax.semilogy(ds, TR, "o-", lw=3, ms=8, color=C_DATA, label="training MSE")
ax.semilogy(ds, TE, "s-", lw=3, ms=8, color=C_RESID, label="test MSE")
ax.axvline(ds[int(np.argmin(TE))], color=C_ALT, ls="--", lw=2.5,
label=f"lowest test error: degree {ds[int(np.argmin(TE))]}")
ax.set(xlabel="polynomial degree", ylabel="MSE", xticks=ds[::2])
ax.legend(); plt.tight_layout(); plt.show()
print(f"lowest test MSE {min(TE):.5f} at degree {ds[int(np.argmin(TE))]}")
print(f"training MSE at degree 16: {TR[-1]:.5f} (test: {TE[-1]:.5f})")
8. Regularized Least Squares¶
Training error falls, test error turns around, and the coefficients explode. Regularization attacks the coefficients directly:
$$ E(\mathbf{w}) = \underbrace{\tfrac{1}{2}\lVert \mathbf{t} - \Phi\mathbf{w}\rVert^2}_{E_D(\mathbf{w})} \;+\;\lambda\, \underbrace{E_W(\mathbf{w})}_{\text{penalty}} $$| penalty | closed form? | effect on $\mathbf{w}$ | |
|---|---|---|---|
| Ridge (L2) | $\tfrac12\lVert\mathbf{w}\rVert_2^2$ | yes: $(\Phi^\top\Phi + \lambda I)^{-1}\Phi^\top\mathbf{t}$ | shrinks all coefficients smoothly |
| Lasso (L1) | $\lVert\mathbf{w}\rVert_1$ | no | drives some coefficients exactly to zero |
8.1 Coefficient paths¶
In [21]:
DEG = 10
i_tr2, i_te2 = train_test_split(np.arange(n), test_size=0.85, random_state=189)
# Standardize the polynomial columns: penalties are not scale-invariant, and x^10 is
# numerically tiny next to x^1. Skipping this is the most common bug in ridge/lasso demos.
Ptr_raw = build_Phi(x_s[i_tr2], DEG + 1, bias=False)
mu, sd = Ptr_raw.mean(0), Ptr_raw.std(0)
Ptr = (Ptr_raw - mu) / sd
Pte = (build_Phi(x_s[i_te2], DEG + 1, bias=False) - mu) / sd
ttr, tte = t_s[i_tr2], t_s[i_te2]
lambdas = np.logspace(-6, 2, 40)
def path(Model, penalty):
coefs, tr, te = [], [], []
for lam in lambdas:
m = Model(alpha=lam, fit_intercept=True, max_iter=500_000, tol=1e-8).fit(Ptr, ttr)
coefs.append(m.coef_.ravel())
tr.append(mean_squared_error(ttr, m.predict(Ptr)))
te.append(mean_squared_error(tte, m.predict(Pte)))
return np.array(coefs), np.array(tr), np.array(te)
ridge_c, ridge_tr, ridge_te = path(Ridge, None)
lasso_c, lasso_tr, lasso_te = path(Lasso, None)
print(f"n_train = {len(i_tr2)} points, degree {DEG} -> {DEG + 1} parameters.")
print(f"unregularized (lambda -> 0) test MSE: {ridge_te[0]:.4f}")
print(f"ridge: best test MSE {ridge_te.min():.4f} at lambda = {lambdas[ridge_te.argmin()]:.4g}"
f" ({ridge_te[0]/ridge_te.min():.0f}x better)")
print(f"lasso: best test MSE {lasso_te.min():.4f} at lambda = {lambdas[lasso_te.argmin()]:.4g}")
In [22]:
fig, axs = plt.subplots(2, 2, figsize=(13, 9))
for row, (coefs, tr, te, name) in enumerate(
[(ridge_c, ridge_tr, ridge_te, "Ridge (L2)"),
(lasso_c, lasso_tr, lasso_te, "Lasso (L1)")]):
ax = axs[row, 0]
for j in range(coefs.shape[1]):
ax.plot(lambdas, coefs[:, j], lw=2, label=f"degree {j+1}")
ax.axhline(0, color="k", lw=1)
ax.set(xscale="log", xlabel=r"$\lambda$", ylabel="coefficient",
title=f"{name}: coefficient paths")
ax.legend(fontsize=8, ncol=2)
ax = axs[row, 1]
ax.plot(lambdas, tr, lw=3, color=C_DATA, label="train MSE")
ax.plot(lambdas, te, lw=3, color=C_RESID, label="test MSE")
ax.axvline(lambdas[te.argmin()], color=C_ALT, ls="--", lw=2.5,
label=rf"best $\lambda$ = {lambdas[te.argmin()]:.2g}")
ax.set(xscale="log", yscale="log", xlabel=r"$\lambda$", ylabel="MSE",
title=f"{name}: train vs test")
ax.legend(fontsize=10)
plt.tight_layout(); plt.show()
8.2 Lasso sets coefficients to zero; ridge does not¶
Index convention: with the bias column excluded, column $j$ holds $x^{j+1}$, so the first row is degree 1, not degree 0.
In [23]:
show = [1e-6, 7e-4, 1e-2, 1e-1, 1e0]
for name, C in (("LASSO", lasso_c), ("RIDGE", ridge_c)):
rows = {}
for lam in show:
j = int(np.argmin(np.abs(lambdas - lam)))
rows[f"lam={lambdas[j]:.1e}"] = C[j]
tbl = pd.DataFrame(rows, index=[f"degree {d}" for d in range(1, DEG + 1)])
print(f"{name} coefficients\n"); print(tbl.round(3).to_string())
print("exact zeros per lambda:", {c: int((tbl[c] == 0).sum()) for c in tbl.columns}, "\n")
8.3 The constraint picture¶
The penalized problem is equivalent to minimizing $E_D(\mathbf{w})$ subject to $E_W(\mathbf{w}) \le c$ for some $c(\lambda)$. The contours of $E_D$ grow until they first touch the constraint region. The L2 region is a circle, smooth everywhere; the L1 region is a diamond whose corners lie on the axes, and a corner is a point where one coordinate is exactly zero.
In [24]:
W0 = np.array([3.2, 0.9])
th = np.deg2rad(25)
Rot = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
A_MAT = Rot @ np.diag([6.0, 1.0]) @ Rot.T * 2.5
def ridge_solution(lam, A=A_MAT, w0=W0):
w = np.linalg.solve(A + lam * np.eye(2), A @ w0)
return w, 0.5 * (w - w0) @ A @ (w - w0), 0.5 * (w @ w)
def lasso_solution(lam, A=A_MAT, w0=W0, iters=4000):
"Coordinate descent with soft-thresholding: reaches exact zeros, unlike Nelder-Mead."
w, b = np.zeros(2), A @ w0
for _ in range(iters):
for j in range(2):
rho = b[j] - A[j] @ w + A[j, j] * w[j]
w[j] = np.sign(rho) * max(abs(rho) - lam, 0.0) / A[j, j]
return w, 0.5 * (w - w0) @ A @ (w - w0), np.sum(np.abs(w))
for lam in [0.0, 2.0, 6.0, 12.0]:
wr, wl = ridge_solution(lam)[0], lasso_solution(lam)[0]
print(f"lambda={lam:5.1f} ridge w = [{wr[0]:6.3f} {wr[1]:6.3f}] "
f"lasso w = [{wl[0]:6.3f} {wl[1]:6.3f}]"
+ (" <- w2 is exactly 0" if wl[1] == 0.0 else ""))
In [25]:
def regularization_figure(solver, kind):
gx = np.linspace(W0[0] - 10, W0[0] + 10, 401)
gy = np.linspace(W0[1] - 10, W0[1] + 10, 401)
GX, GY = np.meshgrid(gx, gy)
U = np.stack([GX - W0[0], GY - W0[1]], axis=-1)
Z = 0.5 * np.einsum("...i,...i", U, U @ A_MAT)
zmax = float(np.percentile(Z, 95))
lams_curve = np.linspace(0.0, 15.0, 400)
ED, EW = np.array([[solver(l)[1], solver(l)[2]] for l in lams_curve]).T
Etot = ED + lams_curve * EW
ymax = float(Etot.max()) * 1.05
def shape(lam):
w = solver(lam)[0]
if kind == "ridge":
c = np.linalg.norm(w)
a = np.linspace(0, 2 * np.pi, 400)
return c * np.cos(a), c * np.sin(a)
c = np.linalg.norm(w, ord=1)
return c * np.array([1, 0, -1, 0, 1]), c * np.array([0, 1, 0, -1, 0])
norm_lbl = "||w||₂ = c(λ)" if kind == "ridge" else "||w||₁ = c(λ)"
fig = make_subplots(rows=1, cols=2, column_widths=[0.55, 0.45],
subplot_titles=(f"{kind.capitalize()}: data contours and constraint region",
"Error decomposition vs λ"))
fig.add_trace(go.Contour(x=gx, y=gy, z=np.clip(Z, 0, zmax), zmin=0, zmax=zmax,
colorscale="Blues", reversescale=True, showscale=False, opacity=0.96,
contours=dict(start=0.01 * zmax, end=0.99 * zmax,
size=0.98 * zmax / 20, showlines=False)), row=1, col=1)
fig.add_trace(go.Scatter(x=[W0[0]], y=[W0[1]], mode="markers",
marker=dict(symbol="star", size=16, color="crimson"),
name="unregularized optimum"), row=1, col=1)
for y, nm in [(ED, "E_D"), (lams_curve * EW, "λ·E_W"), (Etot, "E")]:
fig.add_trace(go.Scatter(x=lams_curve, y=y, mode="lines", line=dict(width=3), name=nm),
row=1, col=2)
cx, cy = shape(0.0)
dyn = [go.Scatter(x=cx, y=cy, mode="lines",
line=dict(width=5, color="darkmagenta"), name=norm_lbl),
go.Scatter(x=[solver(0.0)[0][0]], y=[solver(0.0)[0][1]], mode="markers",
marker=dict(size=14, symbol="x", color="teal"), name="ŵ(λ)"),
go.Scatter(x=[0, 0], y=[0, ymax], mode="lines",
line=dict(width=2, dash="dot", color="teal"), showlegend=False)]
for tr, col in zip(dyn, [1, 1, 2]):
fig.add_trace(tr, row=1, col=col)
dyn_ix = list(range(len(fig.data) - 3, len(fig.data)))
lams = np.linspace(0.0, 15.0, 16)
fig.frames = [go.Frame(name=f"{l:.2f}", traces=dyn_ix, data=[
go.Scatter(x=shape(l)[0], y=shape(l)[1]),
go.Scatter(x=[solver(l)[0][0]], y=[solver(l)[0][1]]),
go.Scatter(x=[l, l], y=[0, ymax])]) for l in lams]
fig.update_layout(
template="plotly_white", height=560,
sliders=[dict(active=0, pad=dict(l=100, t=55), steps=[
{"label": f"λ = {l:.1f}", "method": "animate",
"args": [[f"{l:.2f}"], {"mode": "immediate",
"frame": {"duration": 0, "redraw": True},
"transition": {"duration": 0}}]} for l in lams])],
updatemenus=[dict(type="buttons", x=0.02, y=0, xanchor="left", yanchor="bottom",
direction="left", buttons=[
dict(label="▶", method="animate", args=[None, {"fromcurrent": True,
"frame": {"duration": 600, "redraw": True}, "transition": {"duration": 50}}]),
dict(label="⏸", method="animate", args=[[None], {"mode": "immediate",
"frame": {"duration": 0, "redraw": False}}])])])
fig.update_xaxes(title_text="w₁", range=[W0[0] - 10, W0[0] + 10], row=1, col=1)
fig.update_yaxes(title_text="w₂", range=[W0[1] - 10, W0[1] + 10],
scaleanchor="x", scaleratio=1, row=1, col=1)
fig.update_xaxes(title_text="λ", row=1, col=2)
return fig
regularization_figure(ridge_solution, "ridge").show()
In [26]:
regularization_figure(lasso_solution, "lasso").show()