{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "linreg.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# linear regression \u2014 Python demo\n\nNumerical companion to the entry [linear regression](https://dictionaryofml.org/terms/linreg.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nOne block per paragraph of the entry (marked [P1], [P2], ...): each block verifies numerically what the corresponding paragraph claims, so the entry's statements are backed by a small reproducible experiment. Self-contained (numpy/matplotlib only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/linreg.py`](https://dictionaryofml.org/terms/linreg.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"linreg.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nlinreg.py \u2014 numerical companion to the glossary entry 'linear regression'.\n\nPurpose\n-------\nOne block per paragraph of the entry (marked [P1], [P2], ...): each block\nverifies numerically what the corresponding paragraph claims, so the entry's\nstatements are backed by a small reproducible experiment. Self-contained\n(numpy/matplotlib only), fixed seed.\n\nBlocks\n------\n[P1/P2]  Prediction of a numeric label via h(x) = w^T x; constant feature\n         as intercept.\n[P3]     Least-squares ERM; the entry's Fig. panel (a): with the constant\n         feature x = 1 and labels (2, 3, 4), the solution is the average 3.\n[P4-P6]  Matrix form f(w) = (1/m)||y - Xw||^2; normal equations\n         X^T X w = X^T y; unique closed-form solution under full column rank.\n[P7]     Underdetermined case m < d: two solutions with identical training\n         loss but different predictions (the generalization issue); ridge\n         (l2 penalty, closed form) and Lasso (l1 penalty, proximal GD /\n         ISTA) as regularized variants.\n[P8]     Statistical interpretation: for jointly Gaussian (x, y) the Bayes\n         estimator has w_star = C_x^{-1} c_xy; least squares recovers it\n         from a large sample.\n[P9-P12] GD step operator F^(eta)(w) = w - eta grad f(w): fixed points solve\n         the normal equations; contraction w.r.t. the Euclidean norm for\n         0 < eta < m / lambda_max; convergence speed governed by the\n         condition number lambda_max / lambda_min  ->  linreg_gdconv.csv.\n[P13]    Online GD / LMS: one GD step per arriving data point.\n[P14]    Stability: label-only perturbation, Delta w = X^+ Delta y, the\n         spectral-norm bound, and the entry's Fig. panel (b) numbers\n         (Delta y^(3) = 6 shifts the average from 3 to 5).\n[P15]    Perturbed online GD update = clean update + perturbation term.\n\nOutputs\n-------\nlinreg_gdconv.csv : GD suboptimality per iteration for a well-conditioned\n                    and an ill-conditioned feature matrix (columns:\n                    iter, wellcond, illcond).\nlinreg.png        : 2x2 preview figure (checking only).\n\"\"\"\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nrng = np.random.default_rng(42)\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P1/P2]** Linear regression predicts a numeric label from features via the"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# linear hypothesis map h(x) = w^T x. A constant feature makes it affine.\n# ===========================================================================\nprint(\"[P1/P2] linear hypothesis map with intercept via constant feature\")\nm, d = 50, 3\nw_true = np.array([0.8, -0.5, 2.0, 10.0])          # last entry: intercept\nX_raw = rng.normal(size=(m, d))                     # \"weather measurements\"\nX = np.hstack([X_raw, np.ones((m, 1))])             # append constant feature\ny = X @ w_true + 0.1 * rng.normal(size=m)           # tomorrow's temperature\nx_new = np.append(rng.normal(size=d), 1.0)\ncheck(\"h(x) = w^T x yields a numeric label\", np.isscalar((w_true @ x_new).item()))\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P3]** Least squares as ERM; entry Fig. panel (a): with the constant scalar"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# feature x = 1 and labels (2, 3, 4), the minimizer is the average 3.\n# ===========================================================================\nprint(\"[P3] constant-feature special case reduces to the average\")\ny_fig = np.array([2.0, 3.0, 4.0])\nX_fig = np.ones((3, 1))\nw_avg = np.linalg.lstsq(X_fig, y_fig, rcond=None)[0].item()\ncheck(\"average = 3 as in Fig. panel (a)\", np.isclose(w_avg, 3.0))\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P4-P6]** Matrix form, normal equations, unique closed-form solution when"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# X has full column rank (m >= d).\n# ===========================================================================\nprint(\"[P4-P6] normal equations and closed-form solution\")\nG = X.T @ X                                          # X^T X\nw_hat = np.linalg.solve(G, X.T @ y)                  # closed form\ncheck(\"normal equations X^T X w = X^T y hold\", np.allclose(G @ w_hat, X.T @ y))\ncheck(\"closed form matches lstsq\", np.allclose(w_hat, np.linalg.lstsq(X, y, rcond=None)[0]))\ncheck(\"full column rank (m >= d)\", np.linalg.matrix_rank(X) == X.shape[1])\n\n\ndef f_avg(Xm, ym, w):\n    return np.mean((ym - Xm @ w) ** 2)\n\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P7]** Underdetermined m < d: many solutions tie on the training set but"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# predict differently outside it; ridge and Lasso as regularized variants.\n# ===========================================================================\nprint(\"[P7] underdetermined case, ridge, Lasso\")\nm_u, d_u = 3, 6\nX_u = rng.normal(size=(m_u, d_u))\ny_u = rng.normal(size=m_u)\nw_min = np.linalg.pinv(X_u) @ y_u                    # minimum-norm solution\nnull_dir = np.linalg.svd(X_u)[2][-1]                 # a null-space direction\nw_alt = w_min + 5.0 * null_dir                       # second solution\ncheck(\"both solutions fit the trainset exactly\",\n      np.allclose(X_u @ w_min, y_u) and np.allclose(X_u @ w_alt, y_u))\nx_out = rng.normal(size=d_u)\ncheck(\"but they predict differently outside it\",\n      abs(x_out @ w_min - x_out @ w_alt) > 1e-3)\n\nalpha = 0.1\nw_ridge = np.linalg.solve(X_u.T @ X_u + alpha * m_u * np.eye(d_u), X_u.T @ y_u)\ncheck(\"ridge (l2 penalty) solution is unique/well-defined\",\n      np.all(np.isfinite(w_ridge)))\n\n\ndef ista(Xm, ym, alpha, iters=2000):\n    \"\"\"Proximal GD (ISTA) for the Lasso: average sqerrloss + alpha*||w||_1.\"\"\"\n    mm = Xm.shape[0]\n    L = 2.0 * np.linalg.eigvalsh(Xm.T @ Xm).max() / mm\n    w = np.zeros(Xm.shape[1])\n    for _ in range(iters):\n        g = w - (2.0 / (mm * L)) * Xm.T @ (Xm @ w - ym)\n        w = np.sign(g) * np.maximum(np.abs(g) - alpha / L, 0.0)\n    return w\n\n\nw_lasso = ista(X_u, y_u, alpha)\ncheck(\"Lasso (l1 penalty) drives entries to exactly zero\",\n      np.sum(np.isclose(w_lasso, 0.0)) > 0)\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P8]** Statistical interpretation: jointly Gaussian (x, y) with zero mean;"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Bayes estimator w_star = C_x^{-1} c_xy; least squares recovers it from a\n# large sample.\n# ===========================================================================\nprint(\"[P8] Bayes estimator from covariances vs sample-based least squares\")\nd_g = 3\nA = rng.normal(size=(d_g, d_g))\nC_x = A @ A.T + d_g * np.eye(d_g)                    # invertible covariance\nw_pop = np.array([1.0, -2.0, 0.5])\nm_big = 200_000\nX_g = rng.multivariate_normal(np.zeros(d_g), C_x, size=m_big)\ny_g = X_g @ w_pop + rng.normal(size=m_big)           # zero-mean jointly Gaussian\nc_xy = C_x @ w_pop                                   # E[x y]\nw_star = np.linalg.solve(C_x, c_xy)\nw_ls = np.linalg.solve(X_g.T @ X_g, X_g.T @ y_g)\ncheck(\"w_star = C_x^{-1} c_xy equals the population weights\",\n      np.allclose(w_star, w_pop))\ncheck(\"sample least squares approximates w_star (m = 2e5)\",\n      np.allclose(w_ls, w_star, atol=2e-2))\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P9-P12]** GD step operator: fixed points = normal-equation solutions;"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# contraction w.r.t. the Euclidean norm for 0 < eta < m/lambda_max; speed\n# governed by the condition number  ->  linreg_gdconv.csv.\n# ===========================================================================\nprint(\"[P9-P12] GD step operator, contraction, condition number\")\nlam = np.linalg.eigvalsh(G)\nlam_max, lam_min = lam.max(), lam.min()\neta = 0.9 * m / lam_max                              # 0 < eta < m/lambda_max\n\n\ndef gdstep(Xm, ym, eta_, w):\n    mm = Xm.shape[0]\n    return w + (2.0 * eta_ / mm) * Xm.T @ (ym - Xm @ w)\n\n\ncheck(\"w_hat is a fixed point of the GD step operator\",\n      np.allclose(gdstep(X, y, eta, w_hat), w_hat))\nwa, wb = rng.normal(size=d + 1), rng.normal(size=d + 1)\nq = np.max(np.abs(1.0 - 2.0 * eta * lam / m))        # contraction factor\ncheck(\"contraction w.r.t. the Euclidean norm with factor < 1\",\n      (np.linalg.norm(gdstep(X, y, eta, wa) - gdstep(X, y, eta, wb))\n       <= q * np.linalg.norm(wa - wb) + 1e-12) and q < 1)\n\n\ndef gd_curve(Xm, ym, iters=60):\n    Gm = Xm.T @ Xm\n    lmax = np.linalg.eigvalsh(Gm).max()\n    eta_ = 0.9 * Xm.shape[0] / lmax\n    wh = np.linalg.solve(Gm, Xm.T @ ym)\n    w = np.zeros(Xm.shape[1])\n    errs = []\n    for _ in range(iters + 1):\n        errs.append(np.linalg.norm(w - wh))\n        w = gdstep(Xm, ym, eta_, w)\n    return np.array(errs)\n\n\nX_well = rng.normal(size=(200, 2))                   # cond(X^T X) close to 1\nscales = np.array([1.0, 12.0])\nX_ill = X_well * scales                              # cond larger by ~144\nerr_well = gd_curve(X_well, X_well @ np.ones(2))\nerr_ill = gd_curve(X_ill, X_ill @ np.ones(2))\ncheck(\"larger condition number slows GD convergence\",\n      err_ill[-1] > err_well[-1])\niters = np.arange(len(err_well))\nnp.savetxt(OUT_DIR / \"linreg_gdconv.csv\",\n           np.column_stack([iters, err_well, err_ill]),\n           delimiter=\",\", header=\"iter,wellcond,illcond\", comments=\"\")\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P13]** Online GD / LMS: one GD step per arriving data point."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# ===========================================================================\nprint(\"[P13] online GD (LMS) on streaming data points\")\nw_on = np.zeros(d + 1)\neta_on = 0.02\nfor t in range(m):\n    x_t, y_t = X[t], y[t]\n    w_on = w_on + 2.0 * eta_on * (y_t - w_on @ x_t) * x_t\nerr_online_final = np.linalg.norm(w_on - w_hat)\ncheck(\"one pass of LMS approaches the least-squares solution\",\n      err_online_final < 0.5 * np.linalg.norm(w_hat))\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P14]** Stability under a label-only perturbation: Delta w = X^+ Delta y,"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# the spectral-norm bound, and the entry's Fig. panel (b) numbers.\n# ===========================================================================\nprint(\"[P14] label perturbation: pseudoinverse formula and bound\")\ndy = rng.normal(size=m)\nw_pert = np.linalg.solve(G, X.T @ (y + dy))\nX_pinv = np.linalg.pinv(X)\ncheck(\"Delta w = X^+ Delta y\", np.allclose(w_pert - w_hat, X_pinv @ dy))\ncheck(\"||Delta w|| <= ||X^+||_2 ||Delta y||\",\n      np.linalg.norm(w_pert - w_hat)\n      <= np.linalg.norm(X_pinv, 2) * np.linalg.norm(dy) + 1e-12)\ndy_fig = np.array([0.0, 0.0, 6.0])                   # Fig. panel (b)\nshift = (np.linalg.pinv(X_fig) @ dy_fig).item()\ncheck(\"Fig. panel (b): Delta y^(3) = 6 shifts the average by 2 (3 -> 5)\",\n      np.isclose(shift, 2.0) and np.isclose(w_avg + shift, 5.0))\n\n# ==========================================================================="
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P15]** Perturbed online GD update = clean update + perturbation term."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# ===========================================================================\nprint(\"[P15] perturbed online GD update decomposition\")\nt = 7\ndx_t, dy_t = 0.05 * rng.normal(size=d + 1), 0.3\nx_t, y_t = X[t], y[t]\nw_cur = rng.normal(size=d + 1)\nupd_pert = w_cur + 2 * eta_on * ((y_t + dy_t) - w_cur @ (x_t + dx_t)) * (x_t + dx_t)\nupd_clean = w_cur + 2 * eta_on * (y_t - w_cur @ x_t) * x_t\neps = upd_pert - upd_clean                           # perturbation term\ncheck(\"perturbed update = clean update + perturbation term (depends on \"\n      \"the data-point perturbation and the current w)\",\n      np.allclose(upd_pert, upd_clean + eps))\n\n# ===========================================================================\n# Preview figure (checking only)\n# ===========================================================================\nfig, ax = plt.subplots(2, 2, figsize=(9, 7))\nax[0, 0].scatter([1, 2, 3], y_fig, label=\"labels\")\nax[0, 0].axhline(w_avg, ls=\"--\", label=r\"$\\hat w = 3$\")\nax[0, 0].axhline(w_avg + shift, ls=\":\", label=r\"$\\tilde w = 5$\")\nax[0, 0].set_title(\"[P3/P14] constant feature: average and outlier shift\")\nax[0, 0].set_xlabel(\"data point index $r$\")\nax[0, 0].set_ylabel(\"label $y$\")\nax[0, 0].legend(frameon=False)\nax[0, 1].semilogy(iters, err_well, label=\"well-conditioned\")\nax[0, 1].semilogy(iters, err_ill, \"--\", label=\"ill-conditioned\")\nax[0, 1].set_title(r\"[P9-P12] GD error vs iteration $t$\")\nax[0, 1].set_xlabel(\"iteration\")\nax[0, 1].set_ylabel(\"error\")\nax[0, 1].legend(frameon=False)\nax[1, 0].stem(w_lasso)\nax[1, 0].set_title(\"[P7] Lasso coefficients (sparse)\")\nax[1, 0].set_xlabel(\"coefficient index $j$\")\nax[1, 0].set_ylabel(\"coefficient value\")\nax[1, 1].plot(np.abs(X @ w_on - y), \".\", ms=3)\nax[1, 1].set_title(\"[P13] residuals after one LMS pass\")\nax[1, 1].set_xlabel(\"data point index $r$\")\nax[1, 1].set_ylabel(\"absolute residual\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"linreg.png\", dpi=110)\n\nn_fail = sum(1 for _, ok in report if not ok)\nprint(f\"\\n{len(report)} checks, {n_fail} failed.\")\nraise SystemExit(1 if n_fail else 0)"
  }
 ]
}