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