{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "gd.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# gradient descent (GD) \u2014 Python demo\n\nNumerical companion to the entry [gradient descent (GD)](https://dictionaryofml.org/terms/gd.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 claim of the entry (marked [P...]): each block verifies numerically what the corresponding statement asserts, so the entry's claims 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/gd.py`](https://dictionaryofml.org/terms/gd.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(), \"gd.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ngd.py \u2014 numerical companion to the glossary entry 'gradient descent (GD)'.\n\nPurpose\n-------\nOne block per claim of the entry (marked [P...]): each block verifies\nnumerically what the corresponding statement asserts, so the entry's\nclaims are backed by a small reproducible experiment. Self-contained\n(numpy/matplotlib only), fixed seed.\n\nBlocks\n------\n[P-descent]  On a strongly convex quadratic f(w) = (1/m)||Xw - y||^2\n             (the ERM objective of least-squares linear regression), the\n             GD step w <- w - eta grad f(w) with eta = 1/L does not\n             increase f at any iteration and converges to the unique\n             minimizer w_hat = argmin f.\n[P-rate]     For the convex L-smooth objective, the constant step eta=1/L\n             obeys the suboptimality bound\n                 f(w^(T)) - f* <= L ||w^(0) - w_hat||^2 / (2 T),\n             i.e. the suboptimality falls in proportion to 1/T.\n[P-momentum] Polyak's heavy-ball momentum\n                 w^(t+1) = w^(t) - eta grad f(w^(t)) + beta (w^(t) - w^(t-1)),\n             which combines gradients over several steps, reaches a given\n             suboptimality in far fewer iterations than plain GD on the\n             same strongly convex quadratic.\n[P-euler]    GD is exactly the explicit Euler discretization of the\n             gradient flow dw/dt = -grad f(w): one GD step with step size\n             eta equals one explicit-Euler step of the ODE, and the\n             deviation of the GD iterate from the exact flow trajectory\n             at a fixed time horizon shrinks in proportion to eta\n             (properties of the flow transfer to GD for small eta).\n\nOutputs\n-------\ngd_convergence.csv : iter, gd, momentum  (suboptimality f(w^(t)) - f* per\n                     iteration) \u2014 read by the entry's pgfplots figure.\ngd_flow.csv        : t, w1, w2 \u2014 exact gradient-flow trajectory on a 2-d\n                     quadratic (read by the entry's flow figure).\ngd_flowsmall.csv   : k, w1, w2 \u2014 GD iterates, small step (eta = 0.02).\ngd_flowlarge.csv   : k, w1, w2 \u2014 GD iterates, large step (eta = 0.3).\ngd.png             : preview figure (checking only).\n\nData generated by pythondemos/gd.py.\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# ----------------------------------------------------------------------\n# Strongly convex quadratic: the ERM objective of least-squares\n# linear regression, f(w) = (1/m) || X w - y ||^2.\n# ----------------------------------------------------------------------\nm, d = 200, 20\nX = rng.standard_normal((m, d))\n# mild anisotropy so the ratio mu/L (and hence the gap between GD\n# and momentum) is visible but not pathological.\nX *= np.linspace(1.0, 3.0, d)\nw_true = rng.standard_normal(d)\ny = X @ w_true + 0.1 * rng.standard_normal(m)\n\nA = (2.0 / m) * (X.T @ X)            # curvature matrix of f; its eigenvalues give mu and L\neigs = np.linalg.eigvalsh(A)\nL, mu = float(eigs[-1]), float(eigs[0])   # smoothness / strong-convexity\nw_hat = np.linalg.solve(X.T @ X, X.T @ y)  # unique minimizer\n\n\ndef f(w):\n    return float(np.mean((X @ w - y) ** 2))\n\n\ndef grad(w):\n    return (2.0 / m) * (X.T @ (X @ w - y))\n\n\nf_star = f(w_hat)\nw0 = np.zeros(d)\nT = 60\n\n# ----------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-descent]** On a strongly convex quadratic f(w) = (1/m)||Xw - y||^2 (the ERM objective of least-squares linear regression), the GD step w <- w - eta grad f(w) with eta = 1/L does not increase f at any iteration and converges to the unique minimizer w_hat = argmin f."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# ----------------------------------------------------------------------\neta = 1.0 / L\nw = w0.copy()\ngd_gap = [f(w) - f_star]\ngd_vals = [f(w)]\nfor _ in range(T):\n    w = w - eta * grad(w)\n    gd_vals.append(f(w))\n    gd_gap.append(f(w) - f_star)\n\nmonotone = all(gd_vals[t + 1] <= gd_vals[t] + 1e-12 for t in range(T))\ncheck(\"[P-descent] GD is monotone non-increasing (eta = 1/L)\", monotone)\n# \"ideally converge to a minimum\": the iterate moves toward w_hat and the\n# suboptimality shrinks by orders of magnitude over the run.\ncloser = np.linalg.norm(w - w_hat) < np.linalg.norm(w0 - w_hat)\ncheck(\"[P-descent] GD iterate moves toward w_hat\", closer)\ncheck(\"[P-descent] GD reduces suboptimality >1e5x\", gd_gap[-1] < 1e-5 * gd_gap[0])\n\n# ----------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-rate]** For the convex L-smooth objective, the constant step eta=1/L obeys the suboptimality bound f(w^(T)) - f* <= L ||w^(0) - w_hat||^2 / (2 T), i.e. the suboptimality falls in proportion to 1/T."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# ----------------------------------------------------------------------\ndist0_sq = float(np.linalg.norm(w0 - w_hat) ** 2)\nbound_ok = all(\n    gd_gap[t] <= L * dist0_sq / (2.0 * t) + 1e-9 for t in range(1, T + 1)\n)\ncheck(\"[P-rate] f(w^T) - f* <= L||w0-w_hat||^2/(2T) for all T\", bound_ok)\n\n# ----------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-momentum]** Polyak's heavy-ball momentum w^(t+1) = w^(t) - eta grad f(w^(t)) + beta (w^(t) - w^(t-1)), which combines gradients over several steps, reaches a given suboptimality in far fewer iterations than plain GD on the same strongly convex quadratic."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# ----------------------------------------------------------------------\nsqL, sqmu = np.sqrt(L), np.sqrt(mu)\neta_hb = 4.0 / (sqL + sqmu) ** 2\nbeta_hb = ((sqL - sqmu) / (sqL + sqmu)) ** 2\nw_prev = w0.copy()\nw = w0.copy()\nhb_gap = [f(w) - f_star]\nfor _ in range(T):\n    w_next = w - eta_hb * grad(w) + beta_hb * (w - w_prev)\n    w_prev, w = w, w_next\n    hb_gap.append(f(w) - f_star)\n\ntol = 1e-6\ngd_hit = next((t for t, g in enumerate(gd_gap) if g <= tol), None)\nhb_hit = next((t for t, g in enumerate(hb_gap) if g <= tol), None)\nfaster = hb_hit is not None and (gd_hit is None or hb_hit < gd_hit)\ncheck(f\"[P-momentum] momentum hits {tol:g} first (hb={hb_hit}, gd={gd_hit})\", faster)\n\n# ----------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-fixedpoint]** the GD operator F^(eta): solutions are fixed points;"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# non-expansive for eta = 1/L; a contraction with factor 1 - mu*eta under\n# mu-strong convexity, so the iterates converge geometrically to the\n# unique fixed point (Banach fixed-point theorem).\n# ----------------------------------------------------------------------\nF = lambda wv: wv - eta * grad(wv)\ncheck(\"[P-fixedpoint] F(w_hat) = w_hat (solutions are fixed points)\",\n      np.linalg.norm(F(w_hat) - w_hat) < 1e-10)\npairs = rng.standard_normal((300, 2, d))\nne = all(np.linalg.norm(F(a) - F(b)) <= np.linalg.norm(a - b) + 1e-12\n         for a, b in pairs)\ncheck(\"[P-fixedpoint] F non-expansive for eta = 1/L (300 random pairs)\", ne)\nkappa = 1.0 - mu * eta\nctr = all(np.linalg.norm(F(a) - F(b)) <= kappa * np.linalg.norm(a - b) + 1e-12\n          for a, b in pairs)\ncheck(f\"[P-fixedpoint] contraction with factor 1 - mu*eta = {kappa:.4f}\", ctr)\nwfp = w0.copy()\nit_err = [np.linalg.norm(wfp - w_hat)]\nfor _ in range(T):\n    wfp = F(wfp)\n    it_err.append(np.linalg.norm(wfp - w_hat))\ngeo_it = all(it_err[t + 1] <= kappa * it_err[t] + 1e-12 for t in range(T))\ncheck(\"[P-fixedpoint] ||w^(t+1) - w_hat|| <= (1 - mu*eta) ||w^(t) - w_hat||\",\n      geo_it)\n\n# ----------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-rate-geo]** mu-strongly convex geometric bound"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# f(w^(t)) - f* <= (1 - mu/L)^t (f(w^(0)) - f*)\n# ----------------------------------------------------------------------\ngeo_bound = all(gd_gap[t] <= (1.0 - mu / L) ** t * gd_gap[0] + 1e-12\n                for t in range(T + 1))\ncheck(\"[P-rate-geo] f(w^t) - f* <= (1 - mu/L)^t (f(w^0) - f*)\", geo_bound)\n\n# ----------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-nesterov]** Nesterov's accelerated variant of GD obeys the improved"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# bound f(w^(T)) - f* <= 2 L ||w^(0) - w_hat||^2 / (T+1)^2.\n# ----------------------------------------------------------------------\ny_seq = w0.copy()\nw_cur = w0.copy()\nnes_gap = [f(w_cur) - f_star]\nfor t in range(1, T + 1):\n    w_next = y_seq - (1.0 / L) * grad(y_seq)\n    y_seq = w_next + (t - 1.0) / (t + 2.0) * (w_next - w_cur)\n    w_cur = w_next\n    nes_gap.append(f(w_cur) - f_star)\nnes_bound = all(nes_gap[t] <= 2.0 * L * dist0_sq / (t + 1) ** 2 + 1e-9\n                for t in range(1, T + 1))\ncheck(\"[P-nesterov] Nesterov obeys 2L||w0 - w_hat||^2/(T+1)^2\", nes_bound)\ncheck(\"[P-nesterov] Nesterov reaches a smaller gap than plain GD at T\",\n      nes_gap[-1] < gd_gap[-1])\n\n# ----------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-euler]** GD is exactly the explicit Euler discretization of the gradient flow dw/dt = -grad f(w): one GD step with step size eta equals one explicit-Euler step of the ODE, and the deviation of the GD iterate from the exact flow trajectory at a fixed time horizon shrinks in proportion to eta (properties of the flow transfer to GD for small eta)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# dw/dt = -grad f(w), on the 2-d quadratic f(w) = (w1^2 + 5 w2^2)/2\n# (grad f(w) = A2 w with diagonal A2, so the flow w(t) = exp(-A2 t) w(0) is exact).\n# ----------------------------------------------------------------------\nA2 = np.diag([1.0, 5.0])\nw0_2d = np.array([2.0, 1.8])\n\n\ndef grad2(w):\n    return A2 @ w\n\n\ndef flow(t):\n    return np.exp(-np.diag(A2) * t) * w0_2d\n\n\ndef euler_step(w, eta):\n    return w + eta * (-grad2(w))         # explicit Euler for dw/dt = -grad f\n\n\ndef gd_step(w, eta):\n    return w - eta * grad2(w)            # the GD update\n\n\n# one GD step IS one explicit-Euler step (the identification is exact)\nw_probe = rng.standard_normal(2)\ncheck(\n    \"[P-euler] GD step equals explicit-Euler step of the flow\",\n    np.allclose(gd_step(w_probe, 0.1), euler_step(w_probe, 0.1)),\n)\n\n# global discretization error at fixed horizon tau shrinks like eta\ntau = 1.0\nerrs = {}\nfor eta2 in (0.1, 0.01):\n    k = int(round(tau / eta2))\n    w = w0_2d.copy()\n    for _ in range(k):\n        w = gd_step(w, eta2)\n    errs[eta2] = float(np.linalg.norm(w - flow(tau)))\nratio = errs[0.01] / errs[0.1]\ncheck(\n    f\"[P-euler] error at tau=1 falls ~ eta (ratio {ratio:.3f}, expect ~0.1)\",\n    ratio < 0.2,\n)\n\n# trajectories for the entry's flow figure\nwith open(OUT_DIR / \"gd_flow.csv\", \"w\") as fh:\n    fh.write(\"t,w1,w2\\n\")\n    for t in np.linspace(0.0, 6.0, 300):\n        w1, w2 = flow(t)\n        fh.write(f\"{t:.4f},{w1:.6e},{w2:.6e}\\n\")\n\nfor eta2, steps, name in ((0.02, 250, \"small\"), (0.3, 14, \"large\")):\n    w = w0_2d.copy()\n    with open(OUT_DIR / f\"gd_flow{name}.csv\", \"w\") as fh:\n        fh.write(\"k,w1,w2\\n\")\n        fh.write(f\"0,{w[0]:.6e},{w[1]:.6e}\\n\")\n        for k in range(1, steps + 1):\n            w = gd_step(w, eta2)\n            fh.write(f\"{k},{w[0]:.6e},{w[1]:.6e}\\n\")\n\n# ----------------------------------------------------------------------\n# CSV for the entry's pgfplots figure (clip tiny values so log axis is\n# well-defined).\n# ----------------------------------------------------------------------\nfloor = 1e-12\nwith open(OUT_DIR / \"gd_convergence.csv\", \"w\") as fh:\n    fh.write(\"iter,gd,momentum\\n\")\n    for t in range(T + 1):\n        fh.write(f\"{t},{max(gd_gap[t], floor):.8e},{max(hb_gap[t], floor):.8e}\\n\")\n\n# ----------------------------------------------------------------------\n# preview figure (checking only)\n# ----------------------------------------------------------------------\nfig, (ax, ax2) = plt.subplots(1, 2, figsize=(11, 4))\nit = np.arange(T + 1)\nax.semilogy(it, np.maximum(gd_gap, floor), \"b-o\", ms=3, label=\"GD ($\\\\eta=1/L$)\")\nax.semilogy(it, np.maximum(hb_gap, floor), \"r-s\", ms=3, label=\"heavy-ball momentum\")\nax.semilogy(\n    it[1:], L * dist0_sq / (2.0 * it[1:]), \"k--\", label=\"$L\\\\|w^{(0)}-\\\\hat w\\\\|^2/(2T)$\"\n)\nax.set_xlabel(\"iteration $t$\")\nax.set_ylabel(\"suboptimality $f(w^{(t)})-f^\\\\star$\")\nax.legend()\nax.set_title(\"GD vs momentum on a strongly convex quadratic\")\n\nts = np.linspace(0.0, 6.0, 300)\nfl = np.array([flow(t) for t in ts])\nax2.plot(fl[:, 0], fl[:, 1], \"k-\", label=\"gradient flow\")\nfor eta2, steps, style, lab in (\n    (0.02, 250, \"b--o\", \"GD, $\\\\eta=0.02$\"),\n    (0.3, 14, \"r:s\", \"GD, $\\\\eta=0.3$\"),\n):\n    w = w0_2d.copy()\n    traj = [w.copy()]\n    for _ in range(steps):\n        w = gd_step(w, eta2)\n        traj.append(w.copy())\n    traj = np.array(traj)\n    ax2.plot(traj[:, 0], traj[:, 1], style, ms=3, markevery=5, label=lab)\nax2.set_xlabel(\"$w_1$\")\nax2.set_ylabel(\"$w_2$\")\nax2.legend(frameon=False)\nax2.set_title(\"GD as explicit Euler of the gradient flow\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"gd.png\", dpi=110)\n\nprint(f\"\\nL={L:.3f}, mu={mu:.3f}, cond={L/mu:.1f}\")\nok = all(v for _, v in report)\nprint(\"ALL OK\" if ok else \"SOME CHECKS FAILED\")"
  }
 ]
}