{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "gradient.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# gradient \u2014 Python demo\n\nNumerical companion to the entry [gradient](https://dictionaryofml.org/terms/gradient.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nVerifies, for the quadratic f(w) = (1/2) w^T Q w with Q = [[2, 0.6], [0.6, 1]], the defining and geometric properties of the gradient stated in the entry, and generates the level-set data for the entry's figure. Self-contained (numpy/matplotlib only), deterministic.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/gradient.py`](https://dictionaryofml.org/terms/gradient.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(), \"gradient.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ngradient.py \u2014 numerical companion to the glossary entry 'gradient'.\n\nPurpose\n-------\nVerifies, for the quadratic f(w) = (1/2) w^T Q w with\nQ = [[2, 0.6], [0.6, 1]], the defining and geometric properties of the\ngradient stated in the entry, and generates the level-set data for the\nentry's figure.  Self-contained (numpy/matplotlib only), deterministic.\n\nBlocks\n------\n[B-fd]       The analytic gradient grad f(w') = Q w' matches central finite\n             differences at w' = (1.4, 1.1) to 1e-8.\n[B-taylor]   Local linear approximation: the error\n             |f(w) - f(w') - grad^T (w - w')| / ||w - w'|| vanishes as\n             w -> w' (ratio decreases by ~10x per 10x step shrink).\n[B-steepest] Among 720 unit directions d, the directional derivative\n             grad^T d is maximized by d = grad/||grad|| (within 0.5 deg).\n[B-orth]     The directional derivative along the level-set tangent at w'\n             is zero (up to 1e-12): the gradient is orthogonal to the\n             level set.\n[B-partials] Partial derivatives alone do not guarantee a gradient for a\n             non-convex function: g(w) = w1 w2^2 / (w1^2 + w2^4) has both\n             partials equal to 0 at the origin, yet the linear-\n             approximation error ratio DIVERGES along w = (t^2, t) \u2014 no\n             gradient exists there (the entry's convexity assumption is\n             not dispensable).\n[B-hilbert]  Inner-product dependence of the gradient: w.r.t. the inner\n             product <u,v>_M = u^T M v (an SPD M), the gradient of f at\n             w' is M^{-1} Q w'. It represents the directional derivative\n             (<M^{-1} Q w', d>_M = (Q w')^T d for all d) and is the\n             steepest-ascent direction among 720 directions of unit\n             M-norm.\n[B-min]      At the minimizer w-hat = 0 of f, the gradient vanishes.\n[B-erm]      ML relevance: for a synthetic training set (m = 30) and a\n             linear regression, the analytic ERM\n             gradient -(2/m) sum (y - w^T x) x matches finite\n             differences, vanishes at the least-squares solution, and\n             one GD step decreases the objective.\n[B-backprop] For a one-hidden-layer network with tanh activation, the\n             gradient of the ERM objective computed by backpropagation\n             (chain rule, layer by layer) matches finite differences on\n             all weights.\n\nOutputs\n-------\ngradient_levelsets.csv : polylines of three level sets f(w) = const.\n                         through and inside w' (columns x,y; nan rows\n                         separate the level sets) for the entry's pgfplots\n                         figure.\ngradient.png           : matplotlib preview of the figure (checking only).\n\nThe figure's gradient arrow at w' is the unit vector of\ngrad f(w') = Q w' = (3.46, 1.94); its coordinates are printed below and\nused verbatim in the entry's TikZ code.\n\"\"\"\n\nimport numpy as np                  # the only numerical dependency\nimport matplotlib                   # imported before pyplot to set the backend\n\nmatplotlib.use(\"Agg\")               # non-interactive backend: no display needed\nimport matplotlib.pyplot as plt    # plotting API for the preview figure\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nreport = []                         # collects (check name, pass/fail) pairs\n\n\ndef check(name, ok):                # records and prints one verification\n    report.append((name, bool(ok)))               # store the verdict\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")  # echo it immediately\n\n\nQ = np.array([[2.0, 0.6], [0.6, 1.0]])  # SPD matrix defining f(w) = w^T Q w / 2\n\n\ndef f(w):                           # the quadratic objective f(w) = (1/2) w^T Q w\n    # einsum contracts w_i Q_ij w_j over the LAST axis, so w may be a\n    # single point (shape (2,)) or a whole grid (shape (..., 2))\n    return 0.5 * np.einsum(\"...i,ij,...j->...\", w, Q, w)\n\n\ndef grad(w):                        # analytic gradient of f: grad f(w) = Q w\n    return Q @ w                    # matrix-vector product\n\n\nwp = np.array([1.4, 1.1])           # the point w' marked in the entry's figure\ng = grad(wp)                        # gradient at w' (drawn as the solid arrow)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fd]** The analytic gradient grad f(w') = Q w' matches central finite differences at w' = (1.4, 1.1) to 1e-8."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "eps = 1e-6                          # half-width of the finite-difference stencil\n# central difference (f(w'+eps e_j) - f(w'-eps e_j)) / (2 eps) per coordinate\ng_fd = np.array([(f(wp + eps * e) - f(wp - eps * e)) / (2 * eps)\n                 for e in np.eye(2)])\ncheck(\"[B-fd]       analytic gradient matches finite differences\",\n      np.linalg.norm(g - g_fd) < 1e-8)           # entries agree to 1e-8"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-taylor]** Local linear approximation: the error |f(w) - f(w') - grad^T (w - w')| / ||w - w'|| vanishes as w -> w' (ratio decreases by ~10x per 10x step shrink)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rng = np.random.default_rng(0)      # fixed seed: deterministic direction\nd = rng.standard_normal(2)          # a random approach direction\nd /= np.linalg.norm(d)              # normalized so ||w - w'|| = h below\nratios = []                         # error ratios for shrinking displacements h\nfor h in (1e-1, 1e-2, 1e-3):        # three decades of displacement h\n    w = wp + h * d                  # approach point w = w' + h d\n    # the entry's defining ratio |f(w) - f(w') - g^T (w - w')| / ||w - w'||\n    ratios.append(abs(f(w) - f(wp) - g @ (w - wp)) / h)\ncheck(\"[B-taylor]   linear-approximation error ratio vanishes\",\n      # for a quadratic the ratio is O(h): each 10x shrink of h must shrink\n      # the ratio ~10x (the 1.5 slack absorbs rounding)\n      ratios[0] > 9 * ratios[1] > 81 * ratios[2] / 1.5)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-steepest]** Among 720 unit directions d, the directional derivative grad^T d is maximized by d = grad/||grad|| (within 0.5 deg)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "angles = np.linspace(0.0, 2.0 * np.pi, 720, endpoint=False)  # 0.5-degree grid\ndirs = np.c_[np.cos(angles), np.sin(angles)]  # 720 unit direction vectors\nbest = dirs[np.argmax(dirs @ g)]    # direction maximizing the directional derivative g^T d\n# angle (in degrees) between that maximizer and the normalized gradient\nang_err = np.degrees(np.arccos(np.clip(best @ (g / np.linalg.norm(g)),\n                                       -1.0, 1.0)))\ncheck(\"[B-steepest] gradient direction maximizes directional derivative\",\n      ang_err < 0.5)                # agreement within the grid resolution"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-orth]** The directional derivative along the level-set tangent at w' is zero (up to 1e-12): the gradient is orthogonal to the level set."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "t = np.array([-g[1], g[0]]) / np.linalg.norm(g)   # level-set tangent: rotate g by 90 deg\ncheck(\"[B-orth]     zero directional derivative along the level set\",\n      abs(g @ t) < 1e-12)           # g^T t = 0: gradient orthogonal to the level set"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-partials]** Partial derivatives alone do not guarantee a gradient for a non-convex function: g(w) = w1 w2^2 / (w1^2 + w2^4) has both partials equal to 0 at the origin, yet the linear- approximation error ratio DIVERGES along w = (t^2, t) \u2014 no gradient exists there (the entry's convexity assumption is not dispensable)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# g(w) = w1 w2^2 / (w1^2 + w2^4): both partials exist (= 0) at the origin,\n# but g is not differentiable there (g is not convex).\ndef g_cex(w1, w2):                  # the counterexample function\n    den = w1 ** 2 + w2 ** 4         # denominator, zero only at the origin\n    # define g(0,0) = 0 (its limit along both axes); elsewhere the formula\n    return np.where(den == 0.0, 0.0, w1 * w2 ** 2 / den)\n\n\neps = 1e-6                          # finite-difference half-width (as in [B-fd])\npx = (g_cex(eps, 0.0) - g_cex(-eps, 0.0)) / (2 * eps)  # partial w.r.t. w1 at 0\npy = (g_cex(0.0, eps) - g_cex(0.0, -eps)) / (2 * eps)  # partial w.r.t. w2 at 0\n# along w = (t^2, t), g = 1/2 while ||w|| -> 0: the error ratio\n# |g(w) - 0 - 0| / ||w|| blows up, so no gradient exists at 0.\nratios_cex = [abs(g_cex(t ** 2, t)) / np.hypot(t ** 2, t)\n              for t in (1e-1, 1e-2, 1e-3)]\ncheck(\"[B-partials] partials exist at 0 but the error ratio diverges\",\n      # partials vanish exactly, and the ratio GROWS ~10x per decade\n      abs(px) < 1e-12 and abs(py) < 1e-12\n      and ratios_cex[2] > 9 * ratios_cex[1] > 81 * ratios_cex[0] / 1.5)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-hilbert]** Inner-product dependence of the gradient: w.r.t. the inner product <u,v>_M = u^T M v (an SPD M), the gradient of f at w' is M^{-1} Q w'. It represents the directional derivative (<M^{-1} Q w', d>_M = (Q w')^T d for all d) and is the steepest-ascent direction among 720 directions of unit M-norm."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "M = np.array([[2.0, 0.5], [0.5, 1.0]])            # SPD inner-product matrix\ngM = np.linalg.solve(M, g)                        # gradient w.r.t. <.,.>_M: M^{-1} Q w'\nrng_h = np.random.default_rng(1)    # fixed seed: deterministic test directions\nD = rng_h.standard_normal((100, 2))  # 100 random directions d\ncheck(\"[B-hilbert]  <g_M, d>_M equals the directional derivative g^T d\",\n      # representation property: d^T M g_M = d^T g for every d\n      np.allclose(D @ M @ gM, D @ g, atol=1e-12))\n# rescale the 720 directions from [B-steepest] to unit M-norm ||d||_M = 1\ndirs_M = dirs / np.sqrt(np.einsum(\"ij,jk,ik->i\", dirs, M, dirs))[:, None]\nbest_M = dirs_M[np.argmax(dirs_M @ g)]  # maximizer of the directional derivative g^T d\ncheck(\"[B-hilbert]  steepest ascent w.r.t. the M-norm is along g_M\",\n      # angle between that maximizer and g_M is below the grid resolution\n      np.degrees(np.arccos(np.clip(\n          best_M @ gM / (np.linalg.norm(best_M) * np.linalg.norm(gM)),\n          -1.0, 1.0))) < 0.5)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-min]** At the minimizer w-hat = 0 of f, the gradient vanishes."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "check(\"[B-min]      gradient vanishes at the minimizer w-hat = 0\",\n      np.linalg.norm(grad(np.zeros(2))) == 0.0)   # grad f(0) = Q 0 = 0 exactly"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-erm]** ML relevance: for a synthetic training set (m = 30) and a linear regression, the analytic ERM gradient -(2/m) sum (y - w^T x) x matches finite differences, vanishes at the least-squares solution, and one GD step decreases the objective."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# linear regression: f(w) = (1/m) sum (y - w^T x)^2 with\n# gradient -(2/m) sum (y - w^T x) x, as displayed in the entry.\nrng_e = np.random.default_rng(2)    # fixed seed: reproducible training set\nm = 30                              # size of the training set\nX = rng_e.standard_normal((m, 2))   # rows are the vectors x^(r)\ny = X @ np.array([1.0, -0.5]) + 0.1 * rng_e.standard_normal(m)  # noisy labels\n\n\ndef f_erm(w):                       # the ERM objective: average squared error\n    return np.mean((y - X @ w) ** 2)\n\n\ndef grad_erm(w):                    # its analytic gradient: -(2/m) X^T (y - X w)\n    return -(2.0 / m) * X.T @ (y - X @ w)\n\n\nw0 = np.array([0.4, 0.8])           # an arbitrary (non-optimal) parameter vector\n# central finite differences of f_erm at w0, coordinate by coordinate\ng_erm_fd = np.array([(f_erm(w0 + eps * e) - f_erm(w0 - eps * e)) / (2 * eps)\n                     for e in np.eye(2)])\ncheck(\"[B-erm]      analytic ERM gradient matches finite differences\",\n      np.linalg.norm(grad_erm(w0) - g_erm_fd) < 1e-8)  # agree to 1e-8\nw_hat = np.linalg.lstsq(X, y, rcond=None)[0]  # least-squares solution (ERM minimizer)\ncheck(\"[B-erm]      ERM gradient vanishes at the least-squares solution\",\n      np.linalg.norm(grad_erm(w_hat)) < 1e-10)  # zero-gradient condition at w-hat\ncheck(\"[B-erm]      one GD step decreases the ERM objective\",\n      f_erm(w0 - 0.1 * grad_erm(w0)) < f_erm(w0))  # step along -grad with lrate 0.1"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-backprop]** For a one-hidden-layer network with tanh activation, the gradient of the ERM objective computed by backpropagation (chain rule, layer by layer) matches finite differences on all weights."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# one-hidden-layer network h(x) = w2^T tanh(W1 x); the gradient of the ERM\n# objective via the chain rule (backpropagation) vs finite differences.\nrng_b = np.random.default_rng(3)    # fixed seed: reproducible weights\nW1 = rng_b.standard_normal((3, 2)) * 0.5  # hidden-layer weights (3 units, 2 inputs)\nw2 = rng_b.standard_normal(3) * 0.5       # output-layer weights\n\n\ndef f_net(W1_, w2_):                # ERM objective of the network on (X, y)\n    return np.mean((y - np.tanh(X @ W1_.T) @ w2_) ** 2)\n\n\nA = np.tanh(X @ W1.T)                              # m x 3 activations (forward pass)\nres = y - A @ w2                                   # m residuals y - h(x)\ndpred = -(2.0 / m) * res                           # dL/dprediction, per sample\ngb_w2 = A.T @ dpred                                # backprop: output layer (dL/dw2)\n# hidden layer: chain rule through tanh' = 1 - A^2, then through W1 x\ngb_W1 = ((dpred[:, None] * (1.0 - A ** 2) * w2).T @ X)   # hidden layer\ntheta_fd = []                       # finite-difference gradient, same parameter order\nfor i in range(3):                  # loop over hidden units ...\n    for j in range(2):              # ... and input coordinates\n        P = np.zeros_like(W1); P[i, j] = eps  # perturbation of one W1 entry\n        # central difference w.r.t. that single entry\n        theta_fd.append((f_net(W1 + P, w2) - f_net(W1 - P, w2)) / (2 * eps))\nfor i in range(3):                  # loop over the output-layer weights\n    p = np.zeros(3); p[i] = eps     # perturbation of one w2 entry\n    theta_fd.append((f_net(W1, w2 + p) - f_net(W1, w2 - p)) / (2 * eps))\ncheck(\"[B-backprop] backprop gradient matches finite differences\",\n      # flatten (W1, w2) backprop gradients and compare to the stencil\n      np.linalg.norm(np.r_[gb_W1.ravel(), gb_w2] - theta_fd) < 1e-8)\n\n# --------------------------------------------------------------- level sets\nG = 500                             # grid resolution per axis\n# the grid must contain the outermost level set f(w) = f(w') entirely\n# (max |w_2| on it is sqrt(2 f(w') (Q^{-1})_{22}) ~ 2.92): a too-small grid\n# clips the contour and the closing segment below draws a straight chord.\ngx = np.linspace(-3.2, 3.6, G)      # horizontal grid coordinates\ngy = np.linspace(-3.2, 3.2, G)      # vertical grid coordinates\nGX, GY = np.meshgrid(gx, gy)        # full 2-D evaluation grid\nF = f(np.stack([GX, GY], axis=-1))  # f evaluated on the whole grid at once\nlevels = [0.8, 2.0, float(f(wp))]   # two inner levels plus the level through w'\n\nfig, ax = plt.subplots(figsize=(4.8, 3.8))  # preview canvas (checking only)\ncs = ax.contour(GX, GY, F, levels=levels, colors=\"k\",\n                linewidths=0.9)     # level sets; polylines reused for the CSV\n\nwith open(OUT_DIR / \"gradient_levelsets.csv\", \"w\") as fh:  # pgfplots data file\n    fh.write(\"x,y\\n\")               # column header expected by \\addplot table\n    first = True                    # tracks whether a nan separator is needed\n    for segs in cs.allsegs:         # one entry per contour level ...\n        for s in segs:              # ... each holding its polyline segments\n            if not first:           # between polylines ...\n                fh.write(\"nan,nan\\n\")  # ... a nan row breaks the pgfplots path\n            first = False           # subsequent polylines need the separator\n            for p in s[::5]:        # every 5th vertex suffices at print size\n                fh.write(f\"{p[0]:.4f},{p[1]:.4f}\\n\")  # one vertex per row\n            if np.allclose(s[0], s[-1]):          # close closed loops only\n                fh.write(f\"{s[0][0]:.4f},{s[0][1]:.4f}\\n\")  # repeat first vertex\n\nu = g / np.linalg.norm(g) * 1.2     # unit gradient arrow, scaled to length 1.2\n# print the arrow coordinates that the entry's TikZ code uses verbatim\nprint(f\"\\nw' = {wp} | grad f(w') = {g} | unit arrow (len 1.2) = \"\n      f\"({u[0]:.2f}, {u[1]:.2f})\")\n\n# -------------------------------------------------------------- preview\nax.plot(*wp, \"ko\", ms=4)            # mark the point w'\nax.annotate(\"\", xy=wp + u, xytext=wp,   # solid arrow: the gradient at w'\n            arrowprops=dict(arrowstyle=\"-|>\", lw=2))\nax.annotate(\"\", xy=wp - u, xytext=wp,   # dashed arrow: the negative gradient\n            arrowprops=dict(arrowstyle=\"-|>\", lw=1.2, linestyle=\"--\"))\nax.plot(0, 0, \"k.\", ms=6)           # mark the minimizer w-hat = 0\nax.set_aspect(\"equal\")              # equal axis scaling (as in the TikZ figure)\nax.set_xlabel(\"$w_1$\"); ax.set_ylabel(\"$w_2$\")\nax.set_title(\"level sets of $f$ and the gradient at $w'$\")\nfig.tight_layout()                  # trim whitespace\nfig.savefig(OUT_DIR / \"gradient.png\", dpi=110)  # write the preview PDF\n\nn_ok = sum(ok for _, ok in report)  # count the passed checks\nprint(f\"{n_ok}/{len(report)} checks pass\")  # summary line\nprint(f\"wrote {OUT_DIR / 'gradient_levelsets.csv'}, {OUT_DIR / 'gradient.png'}\")\nif n_ok != len(report):             # any failed check ...\n    raise SystemExit(1)             # ... makes the script exit non-zero"
  }
 ]
}