{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "training.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# training \u2014 Python demo\n\nNumerical companion to the entry [training](https://dictionaryofml.org/terms/training.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nBacks the entry's view of training as an iteration w^(t+1) = T(w^(t)) whose operator T is settled by the ML setting: offline ERM by gradient descent (GD), SGD when the training set is large, online GD when data points arrive one at a time, and alternating optimization in unsupervised settings (k-means, and the EM algorithm for a GMM). Every method stops at, or near, a fixed point of its operator. Self-contained (numpy/matplotlib only), fixed seeds.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/training.py`](https://dictionaryofml.org/terms/training.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(), \"training.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ntraining.py \u2014 numerical companion to the glossary entry 'training'.\n\nPurpose\n-------\nBacks the entry's view of training as an iteration w^(t+1) = T(w^(t))\nwhose operator T is settled by the ML setting: offline ERM by gradient\ndescent (GD), SGD when the training set is\nlarge, online GD when data points arrive one at a time, and alternating\noptimization in unsupervised settings (k-means, and the EM algorithm for\na GMM).  Every method stops at, or near, a fixed point of its operator.\nSelf-contained (numpy/matplotlib only), fixed seeds.\n\nSetup\n-----\nThe entry's shop: daily bread demand forecast from the weather with a\nlinear model, demand = w1 * temperature + w0.  A year of 365 past records\nis generated from a fixed line with noise; a second year with a slowly\nchanging clientele (the intercept drifts) supplies the online setting.\nFor the unsupervised settings, 300 customer profiles in the plane form\nthree groups.\n\nBlocks\n------\n[B-gd]      Offline ERM on the year of records: T(w) = w - eta grad f(w).\n            The change |w^(t+1) - w^(t)| shrinks at every iteration, and\n            the iterates approach the fixed point grad f(w) = 0, the\n            model parameters at which the average loss is smallest.\n[B-sgd]     SGD replaces grad f by the gradient of the average loss over\n            a few records; its iterates approach the same fixed point but\n            the change per iteration never settles below the noise of\n            the subset gradients.\n[B-online]  Online GD takes one gradient step per arriving day on that\n            day's loss alone; with a drifting clientele its model\n            parameters follow the drift, while the offline line, fixed\n            after its year, falls behind.\n[B-kmeans]  Alternating optimization: assign each profile to its nearest\n            cluster mean, recompute the means.  The loss never increases\n            and the assignments stop changing: a fixed point.\n[B-em]      The EM algorithm for a GMM alternates between the\n            responsibilities of the components for each profile and the\n            component parameters; the average loss (negative\n            log-likelihood) never increases and the parameters settle.\n\nOutputs\n-------\ntraining_steps.csv  : iteration, |w^(t+1) - w^(t)| for GD, SGD and online GD\ntraining_online.csv : day of the second year, true intercept, intercept of\n                      the offline line and of the online line\ntraining_kmeans.csv : iteration, k-means loss; and training_em.csv for EM\ntraining.png        : matplotlib preview (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\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\nrng = np.random.default_rng(0)\nN_DAYS = 365\ntemp = rng.uniform(-5.0, 30.0, N_DAYS)\nW_TRUE = np.array([120.0, -1.5])                 # intercept, slope: colder days, more bread\ndemand = W_TRUE[0] + W_TRUE[1] * temp + 6.0 * rng.standard_normal(N_DAYS)\nT_MEAN, T_SD = temp.mean(), temp.std()\nA = np.c_[np.ones(N_DAYS), (temp - T_MEAN) / T_SD]  # standardized temperature, for a well-conditioned f\n\n\ndef grad_f(w, idx=None):\n    \"\"\"Gradient of the average squared error loss over the records idx.\"\"\"\n    Ai = A if idx is None else A[idx]; yi = demand if idx is None else demand[idx]\n    return 2.0 * Ai.T @ (Ai @ w - yi) / len(yi)\n\n\ndef loss(w):\n    return float(np.mean((A @ w - demand) ** 2))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-gd]** Offline ERM on the year of records: T(w) = w - eta grad f(w). The change |w^(t+1) - w^(t)| shrinks at every iteration, and the iterates approach the fixed point grad f(w) = 0, the model parameters at which the average loss is smallest."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "L = 2.0 * np.linalg.eigvalsh(A.T @ A / N_DAYS).max()\nETA, T_GD = 0.5 / L, 40\nw = np.zeros(2); steps_gd = []; w_gd_path = [w.copy()]\nfor t in range(T_GD):\n    w_new = w - ETA * grad_f(w)                  # the operator T of offline ERM\n    steps_gd.append(float(np.linalg.norm(w_new - w))); w = w_new; w_gd_path.append(w.copy())\nw_gd = w\nw_star = np.linalg.solve(A.T @ A, A.T @ demand)  # the fixed point: grad f = 0\ncheck(f\"[B-gd]      the change per iteration of GD never grows and shrinks from {steps_gd[0]:.2f} \"\n      f\"to {steps_gd[-1]:.1e}\", all(b <= a for a, b in zip(steps_gd, steps_gd[1:])) and steps_gd[-1] < 1e-6)\ncheck(f\"[B-gd]      the iterates approach the fixed point of T, where the gradient vanishes: \"\n      f\"|grad f| = {np.linalg.norm(grad_f(w_gd)):.1e}, distance to it {np.linalg.norm(w_gd - w_star):.1e}\",\n      np.linalg.norm(w_gd - w_star) < 1e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-sgd]** SGD replaces grad f by the gradient of the average loss over a few records; its iterates approach the same fixed point but the change per iteration never settles below the noise of the subset gradients."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "BATCH = 8\nw = np.zeros(2); steps_sgd = []\nfor t in range(T_GD):\n    idx = rng.choice(N_DAYS, BATCH, replace=False)\n    w_new = w - ETA * grad_f(w, idx)\n    steps_sgd.append(float(np.linalg.norm(w_new - w))); w = w_new\nw_sgd = w\ncheck(f\"[B-sgd]     SGD with {BATCH} records per step reaches a distance {np.linalg.norm(w_sgd - w_star):.2f} \"\n      f\"from the fixed point (loss {loss(w_sgd):.1f} against {loss(w_star):.1f} at the fixed point)\",\n      loss(w_sgd) < 1.5 * loss(w_star))\ncheck(f\"[B-sgd]     its change per iteration does not settle: the last ten steps average \"\n      f\"{np.mean(steps_sgd[-10:]):.2f}, against {np.mean(steps_gd[-10:]):.1e} for GD\",\n      np.mean(steps_sgd[-10:]) > 100 * np.mean(steps_gd[-10:]))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-online]** Online GD takes one gradient step per arriving day on that day's loss alone; with a drifting clientele its model parameters follow the drift, while the offline line, fixed after its year, falls behind."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# a second year: the clientele changes, the true intercept drifts upward\ndrift = W_TRUE[0] + np.linspace(0.0, 40.0, N_DAYS)\ntemp2 = rng.uniform(-5.0, 30.0, N_DAYS)\ndemand2 = drift + W_TRUE[1] * temp2 + 6.0 * rng.standard_normal(N_DAYS)\nETA_ON = 0.05\nw = w_gd.copy(); on_path = []; steps_on = []; pred_on = []; pred_off = []\nfor d in range(N_DAYS):\n    a_d = np.array([1.0, (temp2[d] - T_MEAN) / T_SD])\n    pred_on.append(float(a_d @ w)); pred_off.append(float(a_d @ w_gd))   # forecasts before the day's sales arrive\n    g = 2.0 * (a_d @ w - demand2[d]) * a_d         # gradient of that day's loss alone\n    w_new = w - ETA_ON * g\n    steps_on.append(float(np.linalg.norm(w_new - w))); w = w_new; on_path.append(w[0])\non_path = np.array(on_path)\ntrue_w0 = drift + W_TRUE[1] * T_MEAN              # the intercept at the mean temperature, day by day\nerr_offline = float(np.mean((demand2 - np.array(pred_off)) ** 2))\nerr_online = float(np.mean((demand2 - np.array(pred_on)) ** 2))\ncheck(f\"[B-online]  with a drifting clientele, the online line follows: its intercept ends at \"\n      f\"{on_path[-1]:.1f} against a true {true_w0[-1]:.1f}, while the offline line stays at {w_gd[0]:.1f}\",\n      abs(on_path[-1] - true_w0[-1]) < abs(w_gd[0] - true_w0[-1]) / 3)\ncheck(f\"[B-online]  average loss over the second year: online {err_online:.0f}, offline {err_offline:.0f}\",\n      err_online < err_offline)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-kmeans]** Alternating optimization: assign each profile to its nearest cluster mean, recompute the means. The loss never increases and the assignments stop changing: a fixed point."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "K = 3\ncenters = np.array([[0.0, 0.0], [6.0, 1.0], [3.0, 5.0]])\nprofiles = np.vstack([c + rng.standard_normal((100, 2)) for c in centers])\nmeans = profiles[rng.choice(len(profiles), K, replace=False)]\nkm_loss = []; assign_prev = None\nfor t in range(50):\n    d2 = ((profiles[:, None, :] - means[None]) ** 2).sum(-1)\n    assign = d2.argmin(1)                         # assignment step\n    km_loss.append(float(d2[np.arange(len(profiles)), assign].mean()))\n    if assign_prev is not None and np.array_equal(assign, assign_prev):\n        break\n    means = np.array([profiles[assign == k].mean(0) for k in range(K)])  # update step\n    assign_prev = assign\ncheck(f\"[B-kmeans]  the k-means loss never increases over {len(km_loss)} iterations \"\n      f\"({km_loss[0]:.2f} -> {km_loss[-1]:.2f}) and the assignments stop changing: a fixed point\",\n      all(b <= a + 1e-12 for a, b in zip(km_loss, km_loss[1:])) and len(km_loss) < 50)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-em]** The EM algorithm for a GMM alternates between the responsibilities of the components for each profile and the component parameters; the average loss (negative log-likelihood) never increases and the parameters settle."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def gauss(x, mu, var):\n    d = x - mu\n    return np.exp(-0.5 * (d ** 2).sum(1) / var) / (2 * np.pi * var)\n\n\nmu = profiles[rng.choice(len(profiles), K, replace=False)]; var = np.full(K, 2.0); pi = np.full(K, 1 / K)\nem_loss = []\nfor t in range(100):\n    dens = np.stack([pi[k] * gauss(profiles, mu[k], var[k]) for k in range(K)], 1)\n    em_loss.append(float(-np.mean(np.log(dens.sum(1)))))      # average negative log-likelihood\n    resp = dens / dens.sum(1, keepdims=True)                   # responsibilities of the components\n    nk = resp.sum(0)\n    mu_new = (resp.T @ profiles) / nk[:, None]                 # component parameters\n    var_new = np.array([(resp[:, k] * ((profiles - mu_new[k]) ** 2).sum(1)).sum() / (2 * nk[k]) for k in range(K)])\n    pi_new = nk / len(profiles)\n    moved = float(np.abs(mu_new - mu).max()); mu, var, pi = mu_new, var_new, pi_new\n    if moved < 1e-6:\n        break\ncheck(f\"[B-em]      the EM algorithm's average loss never increases over {len(em_loss)} iterations \"\n      f\"({em_loss[0]:.3f} -> {em_loss[-1]:.3f}) and the component parameters settle\",\n      all(b <= a + 1e-9 for a, b in zip(em_loss, em_loss[1:])) and moved < 1e-6)\n\n# ----------------------------------------------------------------- CSV\nwith open(OUT_DIR / \"training_steps.csv\", \"w\") as fh:\n    fh.write(\"iteration,gd,sgd,online\\n\")\n    for t in range(T_GD):\n        fh.write(f\"{t + 1},{steps_gd[t]:.5f},{steps_sgd[t]:.5f},{steps_on[t]:.5f}\\n\")\nwith open(OUT_DIR / \"training_online.csv\", \"w\") as fh:\n    fh.write(\"day,true_intercept,offline,online\\n\")\n    for d in range(0, N_DAYS, 5):\n        fh.write(f\"{d + 1},{true_w0[d]:.3f},{w_gd[0]:.3f},{on_path[d]:.3f}\\n\")\nwith open(OUT_DIR / \"training_kmeans.csv\", \"w\") as fh:\n    fh.write(\"iteration,loss\\n\")\n    for t, v in enumerate(km_loss):\n        fh.write(f\"{t},{v:.5f}\\n\")\nwith open(OUT_DIR / \"training_em.csv\", \"w\") as fh:\n    fh.write(\"iteration,loss\\n\")\n    for t, v in enumerate(em_loss):\n        fh.write(f\"{t},{v:.5f}\\n\")\n\n# -------------------------------------------------------------- preview\nfig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12.5, 3.7))\nit = np.arange(1, T_GD + 1)\nax1.semilogy(it, steps_gd, \"k-\", label=\"GD\")\nax1.semilogy(it, steps_sgd, \"k--\", label=\"SGD, 8 records per step\")\nax1.semilogy(it, steps_on[:T_GD], \"k:\", label=\"online GD, one day per step\")\nax1.set_xlabel(\"iteration $t$\"); ax1.set_ylabel(\"$\\\\|\\\\mathbf{w}^{(t+1)} - \\\\mathbf{w}^{(t)}\\\\|$\")\nax1.set_title(\"change per iteration of three operators\"); ax1.legend(frameon=False, fontsize=8)\ndays = np.arange(1, N_DAYS + 1)\nax2.plot(days, true_w0, \"k-\", lw=2.5, alpha=0.35, label=\"true intercept (drifting clientele)\")\nax2.plot(days, np.full(N_DAYS, w_gd[0]), \"k--\", label=\"offline line, trained on year one\")\nax2.plot(days, on_path, \"k-\", label=\"online GD, one step per day\")\nax2.set_xlabel(\"day of the second year\"); ax2.set_ylabel(\"intercept $w_0$\")\nax2.set_title(\"online training follows a drift\"); ax2.legend(frameon=False, fontsize=8)\nax3.plot(range(len(km_loss)), km_loss, \"k-o\", ms=3, label=\"k-means loss\")\nax3b = ax3.twinx(); ax3b.plot(range(len(em_loss)), em_loss, \"k--s\", ms=3, mfc=\"none\", label=\"EM average loss\")\nax3.set_xlabel(\"iteration\"); ax3.set_ylabel(\"k-means loss\"); ax3b.set_ylabel(\"EM average loss\")\nax3b.set_xlabel(\"iteration\")\nax3.set_title(\"alternating optimization reaches a fixed point\"); ax3b.set_title(ax3.get_title())\nh1, l1 = ax3.get_legend_handles_labels(); h2, l2 = ax3b.get_legend_handles_labels()\nax3.legend(h1 + h2, l1 + l2, frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"training.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(f\"wrote {OUT_DIR / 'training_steps.csv'}, {OUT_DIR / 'training_online.csv'}, \"\n      f\"{OUT_DIR / 'training_kmeans.csv'}, {OUT_DIR / 'training_em.csv'}, {OUT_DIR / 'training.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}