{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "backpropagation.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# backpropagation \u2014 Python demo\n\nNumerical companion to the entry [backpropagation](https://dictionaryofml.org/terms/backpropagation.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 claims on a small artificial neural network (ANN) with one hidden layer, the network of the entry's figure: backpropagation delivers the partial derivative of the loss with respect to every weight in one forward and one backward pass, the result agrees with perturbing each weight separately, and the cost of the latter grows with the number of weights while the cost of backpropagation stays a small multiple of one forward pass. 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/backpropagation.py`](https://dictionaryofml.org/terms/backpropagation.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(), \"backpropagation.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nbackpropagation.py \u2014 numerical companion to the glossary entry\n'backpropagation'.\n\nPurpose\n-------\nBacks the entry's claims on a small artificial neural network (ANN) with\none hidden layer, the network of the entry's figure: backpropagation\ndelivers the partial derivative of the loss with respect to every weight\nin one forward and one backward pass, the result agrees with perturbing\neach weight separately, and the cost of the latter grows with the number\nof weights while the cost of backpropagation stays a small multiple of\none forward pass.  Self-contained (numpy/matplotlib only), fixed seed.\n\nSetup\n-----\nData points with d = 8 features and a numeric label, generated from a\nfixed network with added noise.  The ANN has one hidden layer of q\nneurons with the activation function tanh and a single output neuron;\nits weights are collected into one vector w.  The loss function is the\nsquared error, the objective function f(w) the average loss on a batch\nof 64 data points.  \"One pass\" counts the multiplications of a forward\npass through the network for the whole batch.\n\nBlocks\n------\n[B-net]      The forward pass: the batch travels through the hidden\n             layer and the output neuron, and the loss function compares\n             the predictions to the labels.\n[B-backward] The backward pass applies the chain rule from the output\n             back to the input and yields the whole gradient.  Perturbing\n             each weight by 1e-6 and evaluating f again gives the same\n             partial derivatives, to within 1e-6.\n[B-cost]     Multiplications counted: perturbing each weight costs one\n             forward pass per weight, so the cost grows in proportion to\n             the number of weights; backpropagation costs between two and\n             three forward passes whatever the number of weights.\n[B-gd]       Backpropagation is not itself an iteration: gradient\n             descent consumes the gradient in a gradient step and calls\n             backpropagation again.  Twenty steps lower f at every step.\n\nOutputs\n-------\nbackpropagation_cost.csv  : weights, forward passes needed by perturbation\n                            and by backpropagation, for growing hidden layers\nbackpropagation_grad.csv  : for every weight, the partial derivative from\n                            backpropagation and from perturbation\nbackpropagation_gd.csv    : iteration, objective value along gradient descent\nbackpropagation.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)\nD, M = 8, 64                                      # features, batch size\nX = rng.standard_normal((M, D))\nW_true = rng.standard_normal((D, 6)); v_true = rng.standard_normal(6)\ny = np.tanh(X @ W_true) @ v_true + 0.1 * rng.standard_normal(M)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-net]** The forward pass: the batch travels through the hidden layer and the output neuron, and the loss function compares the predictions to the labels."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def unpack(w, q):\n    \"\"\"The weight vector w of a network with q hidden neurons, as its layers.\"\"\"\n    W1 = w[:D * q].reshape(D, q); b1 = w[D * q:D * q + q]\n    W2 = w[D * q + q:D * q + 2 * q]; b2 = w[-1]\n    return W1, b1, W2, b2\n\n\ndef n_weights(q):\n    return D * q + q + q + 1\n\n\ndef forward(w, q):\n    \"\"\"Forward pass on the batch: hidden activations, predictions, average loss.\"\"\"\n    W1, b1, W2, b2 = unpack(w, q)\n    Z = X @ W1 + b1; H = np.tanh(Z)                  # hidden layer\n    pred = H @ W2 + b2                               # output neuron\n    f = float(np.mean((pred - y) ** 2))              # average loss: squared error of the predictions\n    return H, pred, f\n\n\ndef forward_mults(q):\n    \"\"\"Multiplications of one forward pass for the batch.\"\"\"\n    return M * (D * q + q)\n\n\nQ = 6\nw0 = 0.3 * rng.standard_normal(n_weights(Q))\nH0, pred0, f0 = forward(w0, Q)\ncheck(f\"[B-net]      forward pass of a batch of {M} data points through {Q} hidden \"\n      f\"neurons: {n_weights(Q)} weights, average loss {f0:.3f}\",\n      pred0.shape == (M,) and np.isfinite(f0))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-backward]** The backward pass applies the chain rule from the output back to the input and yields the whole gradient. Perturbing each weight by 1e-6 and evaluating f again gives the same partial derivatives, to within 1e-6."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def backward(w, q):\n    \"\"\"Backward pass: the gradient of f, layer by layer from the output back.\"\"\"\n    W1, b1, W2, b2 = unpack(w, q)\n    H, pred, f = forward(w, q)\n    g_pred = 2.0 * (pred - y) / M                    # d f / d prediction\n    gW2 = H.T @ g_pred; gb2 = float(g_pred.sum())    # output layer\n    gH = np.outer(g_pred, W2)                        # chain rule into the hidden layer\n    gZ = gH * (1.0 - H ** 2)                         # through tanh\n    gW1 = X.T @ gZ; gb1 = gZ.sum(axis=0)\n    return np.r_[gW1.ravel(), gb1, gW2, gb2], f\n\n\ndef backward_mults(q):\n    \"\"\"Multiplications of one backward pass (the products above).\"\"\"\n    return M * (q + q + D * q) + M * q\n\n\ndef perturbation_gradient(w, q, eps=1e-6):\n    \"\"\"One extra forward pass per weight.\"\"\"\n    base = forward(w, q)[2]; g = np.zeros_like(w)\n    for k in range(len(w)):\n        wk = w.copy(); wk[k] += eps\n        g[k] = (forward(wk, q)[2] - base) / eps\n    return g\n\n\ng_bp, _ = backward(w0, Q)\ng_pt = perturbation_gradient(w0, Q)\ngap = float(np.max(np.abs(g_bp - g_pt)))\ncheck(f\"[B-backward] the backward pass yields all {len(w0)} partial derivatives at once, \"\n      f\"and they agree with perturbing each weight separately to within {gap:.1e}\",\n      gap < 1e-5)\nwith open(OUT_DIR / \"backpropagation_grad.csv\", \"w\") as fh:\n    fh.write(\"weight,backprop,perturbation\\n\")\n    for k, (a, b) in enumerate(zip(g_bp, g_pt), 1):\n        fh.write(f\"{k},{a:.6f},{b:.6f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-cost]** Multiplications counted: perturbing each weight costs one forward pass per weight, so the cost grows in proportion to the number of weights; backpropagation costs between two and three forward passes whatever the number of weights."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rows = []\nfor q in (2, 4, 8, 16, 32, 64, 128):\n    n = n_weights(q); fwd = forward_mults(q)\n    passes_pt = n + 1                                # one forward pass per weight, plus the base\n    passes_bp = (fwd + backward_mults(q)) / fwd      # in units of one forward pass\n    rows.append((q, n, passes_pt, passes_bp))\ncheck(\"[B-cost]     perturbing each weight costs one forward pass per weight: \"\n      + \", \".join(f\"{n} weights -> {p} passes\" for _, n, p, _ in rows[::3]),\n      all(p == n + 1 for _, n, p, _ in rows))\ncheck(\"[B-cost]     backpropagation costs between two and three forward passes at every \"\n      f\"width ({min(r[3] for r in rows):.2f} to {max(r[3] for r in rows):.2f})\",\n      all(2.0 <= r[3] <= 3.0 for r in rows))\nwith open(OUT_DIR / \"backpropagation_cost.csv\", \"w\") as fh:\n    fh.write(\"hidden,weights,passes_perturbation,passes_backprop\\n\")\n    for q, n, p, b in rows:\n        fh.write(f\"{q},{n},{p},{b:.3f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-gd]** Backpropagation is not itself an iteration: gradient descent consumes the gradient in a gradient step and calls backpropagation again. Twenty steps lower f at every step."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "LRATE, STEPS = 0.05, 20\nw = w0.copy(); fs = [f0]\nfor t in range(STEPS):\n    g, _ = backward(w, Q)                            # backpropagation at the current weights\n    w = w - LRATE * g                                # the gradient step consumes it\n    fs.append(forward(w, Q)[2])\ncheck(f\"[B-gd]       {STEPS} gradient steps, each calling backpropagation once, lower \"\n      f\"the objective at every step: {fs[0]:.3f} -> {fs[-1]:.3f}\",\n      all(b < a for a, b in zip(fs, fs[1:])))\nwith open(OUT_DIR / \"backpropagation_gd.csv\", \"w\") as fh:\n    fh.write(\"iteration,objective\\n\")\n    for t, v in enumerate(fs):\n        fh.write(f\"{t},{v:.5f}\\n\")\n\n# -------------------------------------------------------------- preview\nfig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12.5, 3.7))\nax1.plot(g_pt, g_bp, \"k.\", ms=4)\nlim = [min(g_pt.min(), g_bp.min()), max(g_pt.max(), g_bp.max())]\nax1.plot(lim, lim, \"k--\", lw=0.8, label=\"equal\")\nax1.set_xlabel(\"partial derivative by perturbing one weight\")\nax1.set_ylabel(\"partial derivative by backpropagation\")\nax1.set_title(f\"the {len(w0)} partial derivatives agree\"); ax1.legend(frameon=False)\nax2.plot([r[1] for r in rows], [r[2] for r in rows], \"k-o\", ms=4, label=\"perturbing each weight\")\nax2.plot([r[1] for r in rows], [r[3] for r in rows], \"k--s\", ms=4, mfc=\"none\", label=\"backpropagation\")\nax2.set_xscale(\"log\"); ax2.set_yscale(\"log\")\nax2.set_xlabel(\"number of weights\"); ax2.set_ylabel(\"cost in forward passes\")\nax2.set_title(\"cost of all partial derivatives\"); ax2.legend(frameon=False)\nax3.plot(range(STEPS + 1), fs, \"k-o\", ms=3)\nax3.set_xlabel(\"gradient step\"); ax3.set_ylabel(\"objective $f(\\\\mathbf{w})$\")\nax3.set_title(\"gradient descent using backpropagation\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"backpropagation.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 / 'backpropagation_cost.csv'}, {OUT_DIR / 'backpropagation_grad.csv'}, \"\n      f\"{OUT_DIR / 'backpropagation_gd.csv'}, {OUT_DIR / 'backpropagation.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}