{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "execute:\n",
    "  enabled: false\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-1 fw=all -->\n",
    "\n",
    "# Divergences and Distances Between Distributions\n",
    "\n",
    "Strip away the architecture diagrams and every modern generative model is the\n",
    "same optimization problem: *make my model's distribution close to the data's*.\n",
    "A language model pushes its next-token distribution toward the corpus; a VAE\n",
    "pushes an approximate posterior toward the true one; a GAN pushes generated\n",
    "images toward real ones; a diffusion model pushes a learned vector field toward\n",
    "the data's score. What distinguishes these model families is the **notion of\n",
    "\"close\"**: that choice silently determines the training objective, the\n",
    "gradients it produces, and the characteristic ways it fails.\n",
    "\n",
    "In that section we built one such notion, the\n",
    "Kullback--Leibler divergence, and saw that minimizing it is maximum likelihood.\n",
    "This section builds the rest of the map. There are three great families:\n",
    "*f-divergences*, which average a convex function of the density ratio $p/q$\n",
    "(KL, reverse KL, $\\chi^2$, Hellinger, Jensen--Shannon); *integral probability\n",
    "metrics*, which measure the largest gap in expectation any test function from a\n",
    "chosen class can detect (total variation, MMD, Wasserstein-1); and *optimal\n",
    "transport* distances, which measure how far probability mass must move. We\n",
    "prove the short, structural results that organize the families: Jensen for\n",
    "non-negativity, Fenchel duality for the adversarial (f-GAN) view, Pinsker's\n",
    "inequality tying total variation to KL, and Kantorovich--Rubinstein duality\n",
    "with the one-dimensional closed form for Wasserstein. We close with the *score*:\n",
    "the gradient $\\nabla_{\\mathbf{x}} \\log p(\\mathbf{x})$ of the log-density with\n",
    "respect to the data, the object on which all of\n",
    "that section runs. The capstone is a table\n",
    "that reads off, for each modern generative objective, which divergence it\n",
    "secretly minimizes.\n",
    "\n",
    "As in that section, **every logarithmic quantity is in\n",
    "nats** (natural logarithms); bits are a fixed $\\ln 2$ rescaling. Not every\n",
    "divergence carries units, though: KL, reverse KL, and Jensen--Shannon are in\n",
    "nats, but total variation, $\\chi^2$, and squared Hellinger contain no logarithm\n",
    "and are dimensionless, and transport distances carry the units of the sample\n",
    "space itself. After the imports, everything below is plain NumPy and SciPy."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-imports"
   },
   "source": [
    "%matplotlib inline\n",
    "from d2l import torch as d2l\n",
    "import numpy as np\n",
    "from scipy import integrate, optimize"
   ],
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-2 fw=all -->\n",
    "\n",
    "## What Is a Divergence? The f-Divergence Family\n",
    "\n",
    "### Axioms, Metrics, and Three Families\n",
    "\n",
    "What is the weakest useful notion of \"how far apart\" two distributions are? A\n",
    "**divergence** on a space of distributions is a function $D(P, Q) \\geq 0$ with\n",
    "\n",
    "$$\n",
    "D(P, Q) = 0 \\quad \\textrm{if and only if} \\quad P = Q,\n",
    "$$\n",
    "\n",
    "and nothing more. Compare this with the norms of\n",
    "that section: a *metric* additionally demands symmetry,\n",
    "$D(P, Q) = D(Q, P)$, and the triangle inequality,\n",
    "$D(P, R) \\leq D(P, Q) + D(Q, R)$. We drop both on purpose, because the most\n",
    "useful divergence in machine learning fails both. The KL divergence of\n",
    "that section is asymmetric: for\n",
    "$P = \\mathcal{N}(0, 1)$ and $Q = \\mathcal{N}(0, 4)$ the closed form\n",
    "the equation gives $D_{\\textrm{KL}}(P\\|Q) \\approx 0.318$ nats\n",
    "but $D_{\\textrm{KL}}(Q\\|P) \\approx 0.807$ nats. And KL fails the triangle\n",
    "inequality too, not *because* it is asymmetric (asymmetric functions can\n",
    "perfectly well satisfy a directed triangle inequality; such *quasimetrics*\n",
    "include one-way travel times), but because it scales like a *squared*\n",
    "distance: for unit Gaussians centered at $0$, $1$, and $2$, the same closed\n",
    "form gives $D_{\\textrm{KL}} = \\tfrac{1}{2}$ for each adjacent pair but $2$\n",
    "for the outer pair, and $2 > \\tfrac{1}{2} + \\tfrac{1}{2}$ in either argument\n",
    "order. The asymmetry is informative: that section shows\n",
    "that the *direction* of KL is a modeling decision that changes what the\n",
    "fitted model does. Still, some divergences *are* genuine metrics\n",
    "(total variation, the Hellinger distance, and the Wasserstein\n",
    "distances among them), and when a problem needs the triangle inequality\n",
    "(coupling arguments, convergence proofs), those are the ones to reach for.\n",
    "\n",
    "The divergences that matter in deep learning organize into three families,\n",
    "sketched in the figure:\n",
    "\n",
    "1. **f-divergences** average a convex function of the *density ratio*\n",
    "   $p(x)/q(x)$. They are the information-theoretic family: KL, reverse KL,\n",
    "   $\\chi^2$, squared Hellinger, total variation, Jensen--Shannon.\n",
    "2. **Integral probability metrics (IPMs)** report the largest gap in\n",
    "   expectation $E_P[f] - E_Q[f]$ over a class of test functions $f$. They are\n",
    "   the statistician's family, computable from samples alone, and include\n",
    "   total variation, maximum mean discrepancy (MMD), and Wasserstein-1.\n",
    "3. **Optimal transport distances** measure the minimum cost of physically\n",
    "   moving the mass of $P$ onto $Q$. They metrize the *geometry* of the sample\n",
    "   space, and they alone stay informative when $P$ and $Q$ do not overlap.\n",
    "\n",
    "The families intersect: total variation is both an f-divergence and an IPM,\n",
    "and Wasserstein-1 is both an IPM and a transport distance. The\n",
    "intersections are where the most useful theorems live.\n",
    "\n",
    "![The three families of divergences and where the common ones live. f-divergences average a convex function of the density ratio; integral probability metrics take a supremum of mean differences over a function class; optimal-transport distances move mass. Total variation sits in the first two families; Wasserstein-1 sits in the last two.](../img/mdl-it-divergence-taxonomy.svg)\n",
    "\n",
    "### The f-Divergence Template\n",
    "\n",
    "Much of this zoo is one definition. Let\n",
    "$f : (0, \\infty) \\to \\mathbb{R}$ be convex with $f(1) = 0$, and let $P$ and\n",
    "$Q$ have densities (or p.m.f.s) $p$ and $q$ with $p(x) = 0$ wherever\n",
    "$q(x) = 0$. The **f-divergence** with *generator* $f$ is\n",
    "[@Csiszar.1967]\n",
    "\n",
    "$$\n",
    "D_f(P\\|Q) = E_{x \\sim Q}\\!\\left[ f\\!\\left( \\frac{p(x)}{q(x)} \\right) \\right].\n",
    "$$\n",
    "\n",
    "(When the support condition fails, the definition extends by the standard\n",
    "conventions $0 \\cdot f(0/0) = 0$ and, for outcomes with $q(x) = 0 < p(x)$, a\n",
    "contribution $p(x)\\, f'(\\infty)$ with\n",
    "$f'(\\infty) = \\lim_{t \\to \\infty} f(t)/t$, possibly infinite; the\n",
    "disjoint-support evaluations of that section use this\n",
    "extension.)\n",
    "\n",
    "The density ratio $u = p/q$ says, point by point, how $P$ over- or\n",
    "under-represents $x$ relative to $Q$; the generator $f$ decides how to *score*\n",
    "that discrepancy ($f(1) = 0$: no discrepancy, no score); and the expectation\n",
    "under $Q$ aggregates. Convexity of $f$ is exactly what makes the score sound:\n",
    "\n",
    "**Proposition (non-negativity).** *For any convex $f$ with $f(1) = 0$,*\n",
    "\n",
    "$$\n",
    "D_f(P\\|Q) \\geq 0,\n",
    "$$\n",
    "\n",
    "*and if $f$ is strictly convex at $1$, then $D_f(P\\|Q) = 0$ if and only if\n",
    "$P = Q$.*\n",
    "\n",
    "**Proof.** Jensen's inequality (that section) applied to the\n",
    "convex $f$ and the random variable $u(x) = p(x)/q(x)$ under $x \\sim Q$:\n",
    "\n",
    "$$\n",
    "D_f(P\\|Q) = E_Q[f(u)] \\geq f(E_Q[u])\n",
    "= f\\!\\left( \\sum_x q(x) \\frac{p(x)}{q(x)} \\right) = f(1) = 0.\n",
    "$$\n",
    "\n",
    "If $f$ is strictly convex at $1$, equality forces $u$ to be constant\n",
    "$Q$-almost surely; a constant density ratio between two normalized\n",
    "distributions must equal $1$, i.e., $P = Q$. $\\blacksquare$\n",
    "\n",
    "The proof is Gibbs' inequality from that section, run\n",
    "for every generator at once: taking $f(u) = u \\log u$ recovers\n",
    "$E_Q[(p/q)\\log(p/q)] = E_P[\\log(p/q)] = D_{\\textrm{KL}}(P\\|Q)$ and the proof\n",
    "specializes to the one we gave there.\n",
    "\n",
    "### A Gallery of Generators\n",
    "\n",
    "Each row of the following table is one choice of $f$; the curves are plotted\n",
    "in the figure.\n",
    "\n",
    "| Divergence | Generator $f(u)$ | Symmetric? | Metric? |\n",
    "|:---|:---|:---|:---|\n",
    "| Kullback--Leibler (forward) | $u \\log u$ | no | no |\n",
    "| reverse KL | $-\\log u$ | no | no |\n",
    "| Pearson $\\chi^2$ | $(u - 1)^2$ | no | no |\n",
    "| squared Hellinger $H^2$ | $(\\sqrt{u} - 1)^2$ | yes | $H$ is a metric |\n",
    "| total variation | $\\tfrac{1}{2}\\lvert u - 1 \\rvert$ | yes | yes |\n",
    "| Jensen--Shannon | $\\tfrac{u}{2} \\log u - \\tfrac{u+1}{2} \\log \\tfrac{u+1}{2}$ | yes | $\\sqrt{\\textrm{JS}}$ is a metric |\n",
    "| $\\alpha$-divergence ($\\alpha \\neq 0, 1$) | $\\dfrac{u^\\alpha-\\alpha u+\\alpha-1}{\\alpha(\\alpha-1)}$ | no | no |\n",
    "\n",
    "![Six generators on the same axes. Every curve is convex and touches zero at the no-discrepancy point of the density ratio, but they penalize over- and under-representation very differently: the KL generator grows slowly for small ratios and superlinearly for large ones, the reverse-KL generator diverges as the ratio approaches zero, the chi-squared generator is the symmetric parabola, and the total-variation generator grows only linearly.](../img/mdl-it-f-div-generators.svg)\n",
    "\n",
    "*Reverse KL* is just KL with its arguments swapped: $D_f(P\\|Q)$ with\n",
    "$f(u) = -\\log u$ equals $E_Q[\\log(q/p)] = D_{\\textrm{KL}}(Q\\|P)$, so the\n",
    "asymmetry of KL becomes a *choice of generator* rather than a quirk\n",
    "(Exercise 1). *Total variation* with\n",
    "$f(u) = \\tfrac{1}{2}|u - 1|$ unwinds to\n",
    "$\\tfrac{1}{2}\\sum_x |p(x) - q(x)|$, half the $\\ell_1$ distance between the\n",
    "probability vectors; we study it in that section. And the\n",
    "**Jensen--Shannon divergence** [@Lin.1991], whose generator looks\n",
    "unguessable, is the natural symmetrization of KL: writing\n",
    "$M = \\tfrac{1}{2}(P + Q)$ for the even mixture,\n",
    "\n",
    "$$\n",
    "\\textrm{JS}(P, Q)\n",
    "= \\frac{1}{2} D_{\\textrm{KL}}(P \\,\\|\\, M)\n",
    "+ \\frac{1}{2} D_{\\textrm{KL}}(Q \\,\\|\\, M).\n",
    "$$\n",
    "\n",
    "Each distribution is compared with the mixture $M$, which dominates both;\n",
    "comparing them with each other directly would risk the infinities of KL when\n",
    "supports differ. As a result JS is symmetric, always finite, and bounded:\n",
    "$0 \\leq \\textrm{JS}(P, Q) \\leq \\log 2$, with the upper bound attained exactly\n",
    "when $P$ and $Q$ have disjoint supports (each $\\log(p/m)$ becomes $\\log 2$ on\n",
    "its own support); its square root is moreover a genuine metric\n",
    "[@Endres.Schindelin.2003], as the gallery table records. That\n",
    "boundedness will return as a *bug* in\n",
    "that section: on disjoint supports JS is constant at\n",
    "$\\log 2$, so it provides no gradient. Expanding the equation in\n",
    "terms of the ratio $u = p/q$ produces the generator in the table (Exercise 1\n",
    "checks a similar unwinding).\n",
    "\n",
    "The table's last row is a *family within the family*\n",
    "[@Amari.2016]. For each $\\alpha \\neq 0, 1$ use the normalized generator\n",
    "\n",
    "$$\n",
    "f_\\alpha(u)=\\frac{u^\\alpha-\\alpha u+\\alpha-1}{\\alpha(\\alpha-1)}.\n",
    "$$\n",
    "\n",
    "It is convex with $f_\\alpha(1)=0$, and its pointwise limits are\n",
    "$u\\log u-u+1$ as $\\alpha\\to1$ and $-\\log u+u-1$ as\n",
    "$\\alpha\\to0$. The added linear terms integrate to zero because\n",
    "$E_Q[p/q-1]=0$, so the resulting divergences are forward and reverse KL.\n",
    "At $\\alpha=\\tfrac12$ it gives twice the squared Hellinger divergence under the\n",
    "table's convention. The shorter generator\n",
    "$(u^\\alpha-1)/(\\alpha(\\alpha-1))$ defines the same divergence for fixed\n",
    "$\\alpha$—the two differ only by a multiple of $u-1$—but it does **not** have\n",
    "the claimed pointwise limits. Closely related, and the form that appears in\n",
    "applications, is the **Rényi divergence** [@Renyi.1961],\n",
    "\n",
    "$$\n",
    "D_\\alpha(P\\|Q) = \\frac{1}{\\alpha - 1}\n",
    "\\log \\sum_x p(x)^\\alpha\\, q(x)^{1-\\alpha}.\n",
    "$$\n",
    "\n",
    "The outer logarithm breaks the template, so Rényi is *not* an f-divergence.\n",
    "But it is a monotone increasing function of one: the sum inside is\n",
    "$E_Q[(p/q)^\\alpha] = 1 + \\alpha(\\alpha - 1)\\, D_{f_\\alpha}(P\\|Q)$, and\n",
    "$t \\mapsto \\frac{1}{\\alpha-1}\\log\\big(1 + \\alpha(\\alpha-1)\\,t\\big)$ is\n",
    "increasing for $\\alpha > 0$, so non-negativity and the data-processing\n",
    "inequality (that section) transfer verbatim. One application:\n",
    "*Rényi differential privacy* [@Mironov.2017] measures a randomized\n",
    "algorithm's privacy loss in $D_\\alpha$, which adds exactly across independent\n",
    "compositions. For the full atlas of the family (limits, orderings, and the\n",
    "fact that $D_\\alpha$ is nondecreasing in $\\alpha$) see\n",
    "@Van-Erven.Harremoes.2014.\n",
    "\n",
    "The whole gallery is six lines of code. We evaluate every divergence on the\n",
    "same pair of categorical distributions used in\n",
    "that section, $P = (0.6, 0.3, 0.1)$ and\n",
    "$Q = (0.2, 0.5, 0.3)$, in both argument orders."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-f-table"
   },
   "source": [
    "def f_divergence(f, p, q):\n",
    "    \"\"\"D_f(P||Q) = sum_x q(x) f(p(x)/q(x)); nats for KL/JS, else unitless.\"\"\"\n",
    "    return np.sum(q * f(p / q))\n",
    "\n",
    "generators = {\n",
    "    'KL':          lambda u: u * np.log(u),\n",
    "    'reverse KL':  lambda u: -np.log(u),\n",
    "    'chi^2':       lambda u: (u - 1) ** 2,\n",
    "    'Hellinger^2': lambda u: (np.sqrt(u) - 1) ** 2,\n",
    "    'TV':          lambda u: 0.5 * np.abs(u - 1),\n",
    "    'JS':          lambda u: u / 2 * np.log(u) - (u + 1) / 2 * np.log((u + 1) / 2),\n",
    "}\n",
    "\n",
    "p = np.array([0.6, 0.3, 0.1])\n",
    "q = np.array([0.2, 0.5, 0.3])\n",
    "for name, f in generators.items():\n",
    "    print(f'{name:12s} D_f(P||Q) = {f_divergence(f, p, q):.4f}   '\n",
    "          f'D_f(Q||P) = {f_divergence(f, q, p):.4f}')\n",
    "\n",
    "m = (p + q) / 2\n",
    "js_mix = 0.5 * f_divergence(generators['KL'], p, m) \\\n",
    "    + 0.5 * f_divergence(generators['KL'], q, m)\n",
    "print(f'JS via the mixture formula: {js_mix:.4f}')"
   ],
   "execution_count": 2,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "KL           D_f(P||Q) = 0.3961   D_f(Q||P) = 0.3653\nreverse KL   D_f(P||Q) = 0.3653   D_f(Q||P) = 0.3961\nchi^2        D_f(P||Q) = 1.0133   D_f(Q||P) = 0.8000\nHellinger^2  D_f(P||Q) = 0.1862   D_f(Q||P) = 0.1862\nTV           D_f(P||Q) = 0.4000   D_f(Q||P) = 0.4000\nJS           D_f(P||Q) = 0.0911   D_f(Q||P) = 0.0911\nJS via the mixture formula: 0.0911\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-3 fw=all -->\n",
    "\n",
    "The KL row is asymmetric ($0.3961$ vs.\n",
    "$0.3653$ nats, the same two numbers as in\n",
    "that section), and the reverse-KL row is the KL row\n",
    "with its columns swapped, exactly as the generator algebra predicts. The\n",
    "$\\chi^2$ row is also asymmetric ($1.0133$ vs. $0.8000$). Hellinger, TV, and JS\n",
    "are symmetric, JS comfortably under its $\\log 2 \\approx 0.693$ ceiling at\n",
    "$0.0911$ nats, and the mixture formula the equation reproduces the\n",
    "generator's value exactly.\n",
    "\n",
    "One more structural fact ties the family together. Adding $c\\,(u - 1)$ to any\n",
    "generator changes nothing, since $E_Q[c\\,(p/q - 1)] = c\\,(1 - 1) = 0$ (Exercise\n",
    "1), so each $f$ is really an equivalence class, and what distinguishes\n",
    "divergences *locally* is the curvature $f''(1)$. For $f$ twice continuously\n",
    "differentiable near $1$ and density ratios uniformly near $1$, a second-order\n",
    "expansion of the equation around $u = 1$ makes every such\n",
    "f-divergence the same quadratic to leading order,\n",
    "\n",
    "$$\n",
    "D_f(P\\|Q) \\approx \\frac{f''(1)}{2}\\, \\chi^2(P\\|Q),\n",
    "$$\n",
    "\n",
    "with error third order in the deviation of $p/q$ from $1$\n",
    "[@Van-Erven.Harremoes.2014]: $f''(1) = 1$ for KL and reverse KL,\n",
    "$2$ for $\\chi^2$, $\\tfrac{1}{2}$ for $H^2$, $\\tfrac{1}{4}$ for JS. Near\n",
    "equality all smooth f-divergences agree. This shared local quadratic is the\n",
    "Fisher-information geometry behind natural-gradient methods\n",
    "[@Amari.1998]: for a parametric family $q_{\\boldsymbol{\\theta}}$, the\n",
    "quadratic's Hessian in $\\boldsymbol{\\theta}$ is $f''(1)$ times the Fisher\n",
    "information matrix that that section develops as the\n",
    "curvature of the log-likelihood. The families disagree only about *distant*\n",
    "distributions, which is precisely the regime early in training.\n",
    "\n",
    "## Duality: The Variational View\n",
    "\n",
    "The definition the equation has a practical flaw: it needs the\n",
    "densities. A generative model can *sample*, and the data are samples, but\n",
    "neither side hands you $p(x)/q(x)$. The fix is convex duality, which rewrites\n",
    "every f-divergence as the value of a *game* that only ever evaluates\n",
    "expectations, and expectations can be estimated from samples.\n",
    "\n",
    "### The Fenchel Conjugate and the f-GAN Bound\n",
    "\n",
    "Recall the convex conjugate from that section: for a\n",
    "convex $f$,\n",
    "\n",
    "$$\n",
    "f^*(t) = \\sup_{u > 0}\\, \\big( ut - f(u) \\big),\n",
    "$$\n",
    "\n",
    "and for closed convex $f$ the biconjugation theorem gives back\n",
    "$f(u) = \\sup_t \\big( ut - f^*(t) \\big)$. Geometrically, $-f^*(t)$ is the\n",
    "intercept of the tangent line to $f$ with slope $t$, and biconjugation says a\n",
    "convex function is the upper envelope of its tangent lines\n",
    "(the figure). Keeping *one* tangent instead of the\n",
    "envelope gives the Fenchel--Young inequality\n",
    "\n",
    "$$\n",
    "f(u) \\geq ut - f^*(t) \\quad \\textrm{for all } t,\n",
    "$$\n",
    "\n",
    "with equality when $t$ is the slope of $f$ at $u$, i.e., $t = f'(u)$.\n",
    "\n",
    "![A convex generator is the upper envelope of its tangent lines. Each tangent of slope t has intercept given by the negative conjugate, so keeping one tangent yields the Fenchel--Young lower bound and taking the supremum over slopes recovers the function exactly. Applied pointwise to the density ratio, this picture is the f-GAN variational bound.](../img/mdl-it-f-gan-tangent-bound.svg)\n",
    "\n",
    "Now do this *pointwise inside the divergence*, letting the slope vary with\n",
    "$x$: choose any function $T(x)$ (the **critic**) and apply\n",
    "the equation at $u = p(x)/q(x)$, $t = T(x)$.\n",
    "\n",
    "**Proposition (f-GAN variational bound).** *For every function $T$,*\n",
    "\n",
    "$$\n",
    "D_f(P\\|Q) \\;\\geq\\; E_{x \\sim P}[T(x)] - E_{x \\sim Q}[f^*(T(x))],\n",
    "$$\n",
    "\n",
    "*and the bound is attained at the optimal critic $T^\\star(x) = f'(p(x)/q(x))$,\n",
    "so the supremum over $T$ equals $D_f(P\\|Q)$* [@Nowozin.Cseke.Tomioka.2016].\n",
    "\n",
    "**Proof.** Multiply the equation at $u = p(x)/q(x)$,\n",
    "$t = T(x)$ by $q(x) \\geq 0$:\n",
    "\n",
    "$$\n",
    "q(x)\\, f\\!\\left(\\frac{p(x)}{q(x)}\\right)\n",
    "\\;\\geq\\; p(x)\\, T(x) - q(x)\\, f^*(T(x)).\n",
    "$$\n",
    "\n",
    "Summing (or integrating) over $x$ gives the equation. Since\n",
    "Fenchel--Young holds with equality at $t = f'(u)$, the choice\n",
    "$T^\\star(x) = f'(p(x)/q(x))$ makes the inequality an equality pointwise, hence\n",
    "in expectation. $\\blacksquare$\n",
    "\n",
    "Look at what the right-hand side of the equation asks for:\n",
    "an average of $T$ over samples from $P$ and an average of $f^*(T)$ over\n",
    "samples from $Q$. *No densities anywhere.* Parameterize $T$ by a neural\n",
    "network and you can estimate, and ascend, this bound from minibatches. That\n",
    "is the **f-GAN** recipe: an inner player (the critic) pushes the bound up\n",
    "toward the true divergence, an outer player (the generator, which controls\n",
    "$Q$) pushes it down. The original GAN of\n",
    "@Goodfellow.Pouget-Abadie.Mirza.ea.2014 is the special case\n",
    "corresponding to the Jensen--Shannon generator: at the optimal discriminator,\n",
    "the classic GAN value function equals $2\\,\\textrm{JS}(P, Q) - \\log 4$. For the\n",
    "KL generator $f(u) = u \\log u$ the conjugate is $f^*(t) = e^{t-1}$, giving\n",
    "$D_{\\textrm{KL}}(P\\|Q) \\geq E_P[T] - E_Q[e^{T-1}]$, a bound we will meet\n",
    "again, tightened into Donsker--Varadhan form, when\n",
    "that section estimates mutual information\n",
    "variationally.\n",
    "\n",
    "Two practical caveats. The bound is tight only at the *optimal* critic, so\n",
    "with an undertrained critic the game systematically *underestimates* the\n",
    "divergence: adversarial losses are biased low. And the critic that attains\n",
    "the bound depends on the density ratio, so on disjoint supports (where the\n",
    "ratio is $0$ or $\\infty$) optimal critics saturate, previewing the gradient\n",
    "problems of that section.\n",
    "\n",
    "Let's verify the proposition where we can see everything: on the categorical\n",
    "pair from before, with the $\\chi^2$ generator $f(u) = (u-1)^2$. Its conjugate\n",
    "is $f^*(t) = t + t^2/4$ and the optimal critic is $T^\\star = 2(p/q - 1)$\n",
    "(Exercise 2 derives both). One fine point: $t + t^2/4$ is the conjugate taken\n",
    "over all $u \\in \\mathbb{R}$; over the generator's true domain $u \\in (0, \\infty)$\n",
    "the supremum flattens to $f^*(t) = -1$ for $t \\leq -2$. Using the larger\n",
    "$\\mathbb{R}$-conjugate is safe (a bigger $f^*$ only *lowers* the bound\n",
    "the equation) and costs nothing at the optimum, where\n",
    "$T^\\star = 2(p/q - 1) > -2$ automatically. On three outcomes a critic is just\n",
    "a vector of three numbers."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-f-gan-bound"
   },
   "source": [
    "chi_sq = generators['chi^2']            # f(u) = (u - 1)^2\n",
    "f_star = lambda t: t + t ** 2 / 4       # its convex conjugate\n",
    "\n",
    "def fgan_bound(T, p, q):\n",
    "    \"\"\"E_P[T] - E_Q[f*(T)]: a lower bound on D_f for any critic T.\"\"\"\n",
    "    return np.sum(p * T) - np.sum(q * f_star(T))\n",
    "\n",
    "T_star = 2 * (p / q - 1)                # the optimal critic f'(p/q)\n",
    "print(f'exact chi^2(P||Q)           = {f_divergence(chi_sq, p, q):.4f}')\n",
    "print(f'bound at the optimal critic = {fgan_bound(T_star, p, q):.4f}')\n",
    "rng = np.random.default_rng(42)\n",
    "for scale in (0.5, 1.0, 2.0):\n",
    "    T = T_star + scale * rng.standard_normal(3)\n",
    "    print(f'bound at a perturbed critic (scale {scale}): '\n",
    "          f'{fgan_bound(T, p, q):.4f}')"
   ],
   "execution_count": 3,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "exact chi^2(P||Q)           = 1.0133\nbound at the optimal critic = 1.0133\nbound at a perturbed critic (scale 0.5): 0.9678\nbound at a perturbed critic (scale 1.0): 0.3661\nbound at a perturbed critic (scale 2.0): 0.9600\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-4 fw=all -->\n",
    "\n",
    "The optimal critic reproduces the exact divergence, $1.0133$, to every\n",
    "printed digit, and every random perturbation of it ($0.9678$, $0.3661$,\n",
    "$0.9600$) lands strictly below: the variational objective is a floor that\n",
    "only the right critic touches.\n",
    "\n",
    "### Forward vs. Reverse KL: Mode-Covering vs. Mode-Seeking\n",
    "\n",
    "When fitting a model $Q_\\theta$ to a target $P$, the two directions of KL ask\n",
    "for different things:\n",
    "\n",
    "* **Forward KL**, $D_{\\textrm{KL}}(P \\,\\|\\, Q_\\theta) = E_P[\\log(p/q_\\theta)]$,\n",
    "  samples from the *truth*. Wherever $P$ puts mass, $q_\\theta$ appears in a\n",
    "  denominator: if $q_\\theta(x) \\to 0$ while $p(x) > 0$, the divergence blows\n",
    "  up. Forward KL is **zero-avoiding**: the model must cover every mode of the\n",
    "  data, even at the price of smearing mass over regions $P$ never visits.\n",
    "* **Reverse KL**, $D_{\\textrm{KL}}(Q_\\theta \\,\\|\\, P) = E_{Q_\\theta}[\\log(q_\\theta/p)]$,\n",
    "  samples from the *model*. Now $p$ is the denominator: the model is punished\n",
    "  for putting mass where the truth has none, while modes it never visits\n",
    "  leave the objective untouched. Reverse KL is **zero-forcing**: the model hugs the\n",
    "  mass it can explain and confidently drops the rest.\n",
    "\n",
    "The two directions correspond to the two great fitting paradigms. Maximum\n",
    "likelihood *is* forward KL minimization: we proved in\n",
    "that section that the average negative log-likelihood\n",
    "is the cross-entropy from the empirical distribution to the model. Variational\n",
    "inference and the ELBO (that section) minimize the\n",
    "*reverse* KL from the approximate posterior to the true one, which is why\n",
    "variational posteriors are characteristically too narrow and why VAEs can drop\n",
    "modes.\n",
    "\n",
    "When the model family contains the target, both directions agree on the\n",
    "answer. The interesting case is a *misspecified* family, and the simplest\n",
    "instance is fitting a single Gaussian to a bimodal target. For the forward\n",
    "direction the optimum is fully characterized:\n",
    "\n",
    "**Proposition (forward KL fits moments).** *Over the Gaussian family\n",
    "$Q = \\mathcal{N}(\\mu, \\sigma^2)$, the forward divergence\n",
    "$D_{\\textrm{KL}}(P \\,\\|\\, Q)$ is minimized at*\n",
    "\n",
    "$$\n",
    "\\mu^* = E_P[X], \\qquad (\\sigma^*)^2 = \\mathrm{Var}_P(X).\n",
    "$$\n",
    "\n",
    "**Proof.** For $P$ with a density and finite differential entropy $h(P)$,\n",
    "$D_{\\textrm{KL}}(P\\|Q) = -h(P) - E_P[\\log q(X)]$ with $h(P)$ fixed, so we\n",
    "maximize\n",
    "$E_P[\\log q(X)] = -\\tfrac{1}{2}\\log(2\\pi\\sigma^2) - E_P[(X-\\mu)^2]/(2\\sigma^2)$.\n",
    "For any $\\sigma$, $E_P[(X-\\mu)^2] = \\mathrm{Var}_P(X) + (E_P[X] - \\mu)^2$ is\n",
    "minimized at $\\mu = E_P[X]$; substituting $v = \\mathrm{Var}_P(X)$ and setting\n",
    "the $\\sigma^2$-derivative of $-\\tfrac{1}{2}\\log\\sigma^2 - v/(2\\sigma^2)$ to\n",
    "zero gives $\\sigma^2 = v$. $\\blacksquare$\n",
    "\n",
    "(The same argument runs for any exponential family: the forward-KL projection\n",
    "matches expected sufficient statistics, the *M-projection* of information\n",
    "geometry [@Amari.2016].) The reverse direction has no such closed form, and that is the\n",
    "point: its optima are local, one per mode. Let's compute both fits for the\n",
    "mixture $P = 0.7\\,\\mathcal{N}(-2, 0.6^2) + 0.3\\,\\mathcal{N}(2, 0.6^2)$,\n",
    "evaluating each KL by quadrature on a grid and minimizing over\n",
    "$(\\mu, \\log\\sigma)$ with a derivative-free optimizer."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-fwd-rev-kl"
   },
   "source": [
    "w, mus, sigmas = np.array([0.7, 0.3]), np.array([-2.0, 2.0]), np.array([0.6, 0.6])\n",
    "x = np.linspace(-8.0, 8.0, 4001)\n",
    "\n",
    "def log_p(x):\n",
    "    comps = np.stack([np.log(wi) - 0.5 * np.log(2 * np.pi * si ** 2)\n",
    "                      - (x - mi) ** 2 / (2 * si ** 2)\n",
    "                      for wi, mi, si in zip(w, mus, sigmas)])\n",
    "    cmax = comps.max(axis=0)\n",
    "    return cmax + np.log(np.exp(comps - cmax).sum(axis=0))\n",
    "\n",
    "def log_q(x, theta):\n",
    "    mu, log_sigma = theta\n",
    "    return (-0.5 * np.log(2 * np.pi) - log_sigma\n",
    "            - (x - mu) ** 2 / (2 * np.exp(2 * log_sigma)))\n",
    "\n",
    "def kl(log_a, log_b):\n",
    "    \"\"\"D_KL(A||B) on the grid, by quadrature; in nats.\"\"\"\n",
    "    return integrate.trapezoid(np.exp(log_a) * (log_a - log_b), x)\n",
    "\n",
    "fwd = lambda th: kl(log_p(x), log_q(x, th))     # D_KL(P || Q_theta)\n",
    "rev = lambda th: kl(log_q(x, th), log_p(x))     # D_KL(Q_theta || P)\n",
    "\n",
    "th_fwd = optimize.minimize(fwd, x0=[0.0, 0.0], method='Nelder-Mead').x\n",
    "th_rev = optimize.minimize(rev, x0=[0.0, 0.0], method='Nelder-Mead').x\n",
    "th_rev2 = optimize.minimize(rev, x0=[2.5, -0.5], method='Nelder-Mead').x\n",
    "\n",
    "for name, th, val in [('forward KL', th_fwd, fwd(th_fwd)),\n",
    "                      ('reverse KL', th_rev, rev(th_rev)),\n",
    "                      ('reverse KL, 2nd start', th_rev2, rev(th_rev2))]:\n",
    "    print(f'{name:22s} mu = {th[0]:+.3f}, sigma = {np.exp(th[1]):.3f}, '\n",
    "          f'KL = {val:.3f} nats')\n",
    "\n",
    "d2l.plot(x, [np.exp(log_p(x)), np.exp(log_q(x, th_fwd)),\n",
    "             np.exp(log_q(x, th_rev))], 'x', 'density',\n",
    "         legend=['target P (mixture)', 'forward-KL fit', 'reverse-KL fit'],\n",
    "         figsize=(6, 3))"
   ],
   "execution_count": 4,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "forward KL             mu = -0.800, sigma = 1.929, KL = 0.558 nats\nreverse KL             mu = -1.998, sigma = 0.603, KL = 0.356 nats\nreverse KL, 2nd start  mu = +1.995, sigma = 0.607, KL = 1.202 nats\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"385.78125pt\" height=\"211.07625pt\" viewBox=\"0 0 385.78125 211.07625\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n <metadata>\n  <rdf:RDF xmlns:dc=\"http://purl.org/dc/elements/1.1/\" xmlns:cc=\"http://creativecommons.org/ns#\" xmlns:rdf=\"http://www.w3.org/1999/02/22-rdf-syntax-ns#\">\n   <cc:Work>\n    <dc:type rdf:resource=\"http://purl.org/dc/dcmitype/StillImage\"/>\n    <dc:date>2026-07-19T13:41:28.488412</dc:date>\n    <dc:format>image/svg+xml</dc:format>\n    <dc:creator>\n     <cc:Agent>\n      <dc:title>Matplotlib v3.10.8, https://matplotlib.org/</dc:title>\n     </cc:Agent>\n    </dc:creator>\n   </cc:Work>\n  </rdf:RDF>\n </metadata>\n <defs>\n  <style type=\"text/css\">*{stroke-linejoin: round; stroke-linecap: butt}</style>\n </defs>\n <g id=\"figure_1\">\n  <g id=\"patch_1\">\n   <path d=\"M 0 211.07625 \nL 385.78125 211.07625 \nL 385.78125 0 \nL 0 0 \nz\n\" style=\"fill: #ffffff\"/>\n  </g>\n  <g id=\"axes_1\">\n   <g id=\"patch_2\">\n    <path d=\"M 43.78125 173.52 \nL 378.58125 173.52 \nL 378.58125 7.2 \nL 43.78125 7.2 \nz\n\" style=\"fill: #ffffff\"/>\n   </g>\n   <g id=\"matplotlib.axis_1\">\n    <g id=\"xtick_1\">\n     <g id=\"line2d_1\">\n      <path d=\"M 58.999432 173.52 \nL 58.999432 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_2\">\n      <defs>\n       <path id=\"m5c7fefe946\" d=\"M 0 0 \nL 0 3.5 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"58.999432\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_1\">\n      <!-- −8 -->\n      <g transform=\"translate(51.628338 188.118438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-2212\" d=\"M 678 2272 \nL 4684 2272 \nL 4684 1741 \nL 678 1741 \nL 678 2272 \nz\n\" transform=\"scale(0.015625)\"/>\n        <path id=\"DejaVuSans-38\" d=\"M 2034 2216 \nQ 1584 2216 1326 1975 \nQ 1069 1734 1069 1313 \nQ 1069 891 1326 650 \nQ 1584 409 2034 409 \nQ 2484 409 2743 651 \nQ 3003 894 3003 1313 \nQ 3003 1734 2745 1975 \nQ 2488 2216 2034 2216 \nz\nM 1403 2484 \nQ 997 2584 770 2862 \nQ 544 3141 544 3541 \nQ 544 4100 942 4425 \nQ 1341 4750 2034 4750 \nQ 2731 4750 3128 4425 \nQ 3525 4100 3525 3541 \nQ 3525 3141 3298 2862 \nQ 3072 2584 2669 2484 \nQ 3125 2378 3379 2068 \nQ 3634 1759 3634 1313 \nQ 3634 634 3220 271 \nQ 2806 -91 2034 -91 \nQ 1263 -91 848 271 \nQ 434 634 434 1313 \nQ 434 1759 690 2068 \nQ 947 2378 1403 2484 \nz\nM 1172 3481 \nQ 1172 3119 1398 2916 \nQ 1625 2713 2034 2713 \nQ 2441 2713 2670 2916 \nQ 2900 3119 2900 3481 \nQ 2900 3844 2670 4047 \nQ 2441 4250 2034 4250 \nQ 1625 4250 1398 4047 \nQ 1172 3844 1172 3481 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-2212\"/>\n       <use xlink:href=\"#DejaVuSans-38\" transform=\"translate(83.789062 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_2\">\n     <g id=\"line2d_3\">\n      <path d=\"M 97.044886 173.52 \nL 97.044886 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_4\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"97.044886\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_2\">\n      <!-- −6 -->\n      <g transform=\"translate(89.673793 188.118438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-36\" d=\"M 2113 2584 \nQ 1688 2584 1439 2293 \nQ 1191 2003 1191 1497 \nQ 1191 994 1439 701 \nQ 1688 409 2113 409 \nQ 2538 409 2786 701 \nQ 3034 994 3034 1497 \nQ 3034 2003 2786 2293 \nQ 2538 2584 2113 2584 \nz\nM 3366 4563 \nL 3366 3988 \nQ 3128 4100 2886 4159 \nQ 2644 4219 2406 4219 \nQ 1781 4219 1451 3797 \nQ 1122 3375 1075 2522 \nQ 1259 2794 1537 2939 \nQ 1816 3084 2150 3084 \nQ 2853 3084 3261 2657 \nQ 3669 2231 3669 1497 \nQ 3669 778 3244 343 \nQ 2819 -91 2113 -91 \nQ 1303 -91 875 529 \nQ 447 1150 447 2328 \nQ 447 3434 972 4092 \nQ 1497 4750 2381 4750 \nQ 2619 4750 2861 4703 \nQ 3103 4656 3366 4563 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-2212\"/>\n       <use xlink:href=\"#DejaVuSans-36\" transform=\"translate(83.789062 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_3\">\n     <g id=\"line2d_5\">\n      <path d=\"M 135.090341 173.52 \nL 135.090341 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_6\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"135.090341\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_3\">\n      <!-- −4 -->\n      <g transform=\"translate(127.719247 188.118438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-34\" d=\"M 2419 4116 \nL 825 1625 \nL 2419 1625 \nL 2419 4116 \nz\nM 2253 4666 \nL 3047 4666 \nL 3047 1625 \nL 3713 1625 \nL 3713 1100 \nL 3047 1100 \nL 3047 0 \nL 2419 0 \nL 2419 1100 \nL 313 1100 \nL 313 1709 \nL 2253 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-2212\"/>\n       <use xlink:href=\"#DejaVuSans-34\" transform=\"translate(83.789062 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_4\">\n     <g id=\"line2d_7\">\n      <path d=\"M 173.135795 173.52 \nL 173.135795 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_8\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"173.135795\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_4\">\n      <!-- −2 -->\n      <g transform=\"translate(165.764702 188.118438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-32\" d=\"M 1228 531 \nL 3431 531 \nL 3431 0 \nL 469 0 \nL 469 531 \nQ 828 903 1448 1529 \nQ 2069 2156 2228 2338 \nQ 2531 2678 2651 2914 \nQ 2772 3150 2772 3378 \nQ 2772 3750 2511 3984 \nQ 2250 4219 1831 4219 \nQ 1534 4219 1204 4116 \nQ 875 4013 500 3803 \nL 500 4441 \nQ 881 4594 1212 4672 \nQ 1544 4750 1819 4750 \nQ 2544 4750 2975 4387 \nQ 3406 4025 3406 3419 \nQ 3406 3131 3298 2873 \nQ 3191 2616 2906 2266 \nQ 2828 2175 2409 1742 \nQ 1991 1309 1228 531 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-2212\"/>\n       <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(83.789062 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_5\">\n     <g id=\"line2d_9\">\n      <path d=\"M 211.18125 173.52 \nL 211.18125 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_10\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"211.18125\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_5\">\n      <!-- 0 -->\n      <g transform=\"translate(208 188.118438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-30\" d=\"M 2034 4250 \nQ 1547 4250 1301 3770 \nQ 1056 3291 1056 2328 \nQ 1056 1369 1301 889 \nQ 1547 409 2034 409 \nQ 2525 409 2770 889 \nQ 3016 1369 3016 2328 \nQ 3016 3291 2770 3770 \nQ 2525 4250 2034 4250 \nz\nM 2034 4750 \nQ 2819 4750 3233 4129 \nQ 3647 3509 3647 2328 \nQ 3647 1150 3233 529 \nQ 2819 -91 2034 -91 \nQ 1250 -91 836 529 \nQ 422 1150 422 2328 \nQ 422 3509 836 4129 \nQ 1250 4750 2034 4750 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_6\">\n     <g id=\"line2d_11\">\n      <path d=\"M 249.226705 173.52 \nL 249.226705 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_12\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"249.226705\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_6\">\n      <!-- 2 -->\n      <g transform=\"translate(246.045455 188.118438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-32\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_7\">\n     <g id=\"line2d_13\">\n      <path d=\"M 287.272159 173.52 \nL 287.272159 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_14\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"287.272159\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_7\">\n      <!-- 4 -->\n      <g transform=\"translate(284.090909 188.118438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-34\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_8\">\n     <g id=\"line2d_15\">\n      <path d=\"M 325.317614 173.52 \nL 325.317614 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_16\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"325.317614\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_8\">\n      <!-- 6 -->\n      <g transform=\"translate(322.136364 188.118438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-36\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_9\">\n     <g id=\"line2d_17\">\n      <path d=\"M 363.363068 173.52 \nL 363.363068 7.2 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_18\">\n      <g>\n       <use xlink:href=\"#m5c7fefe946\" x=\"363.363068\" y=\"173.52\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_9\">\n      <!-- 8 -->\n      <g transform=\"translate(360.181818 188.118438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-38\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_10\">\n     <!-- x -->\n     <g transform=\"translate(208.221875 201.796563) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-78\" d=\"M 3513 3500 \nL 2247 1797 \nL 3578 0 \nL 2900 0 \nL 1881 1375 \nL 863 0 \nL 184 0 \nL 1544 1831 \nL 300 3500 \nL 978 3500 \nL 1906 2253 \nL 2834 3500 \nL 3513 3500 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-78\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"matplotlib.axis_2\">\n    <g id=\"ytick_1\">\n     <g id=\"line2d_19\">\n      <path d=\"M 43.78125 165.96 \nL 378.58125 165.96 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_20\">\n      <defs>\n       <path id=\"mfc859cc296\" d=\"M 0 0 \nL -3.5 0 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#mfc859cc296\" x=\"43.78125\" y=\"165.96\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_11\">\n      <!-- 0.0 -->\n      <g transform=\"translate(20.878125 169.759219) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-2e\" d=\"M 684 794 \nL 1344 794 \nL 1344 0 \nL 684 0 \nL 684 794 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(95.410156 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_2\">\n     <g id=\"line2d_21\">\n      <path d=\"M 43.78125 143.099036 \nL 378.58125 143.099036 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_22\">\n      <g>\n       <use xlink:href=\"#mfc859cc296\" x=\"43.78125\" y=\"143.099036\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_12\">\n      <!-- 0.1 -->\n      <g transform=\"translate(20.878125 146.898254) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-31\" d=\"M 794 531 \nL 1825 531 \nL 1825 4091 \nL 703 3866 \nL 703 4441 \nL 1819 4666 \nL 2450 4666 \nL 2450 531 \nL 3481 531 \nL 3481 0 \nL 794 0 \nL 794 531 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(95.410156 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_3\">\n     <g id=\"line2d_23\">\n      <path d=\"M 43.78125 120.238071 \nL 378.58125 120.238071 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_24\">\n      <g>\n       <use xlink:href=\"#mfc859cc296\" x=\"43.78125\" y=\"120.238071\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_13\">\n      <!-- 0.2 -->\n      <g transform=\"translate(20.878125 124.03729) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(95.410156 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_4\">\n     <g id=\"line2d_25\">\n      <path d=\"M 43.78125 97.377107 \nL 378.58125 97.377107 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_26\">\n      <g>\n       <use xlink:href=\"#mfc859cc296\" x=\"43.78125\" y=\"97.377107\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_14\">\n      <!-- 0.3 -->\n      <g transform=\"translate(20.878125 101.176326) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-33\" d=\"M 2597 2516 \nQ 3050 2419 3304 2112 \nQ 3559 1806 3559 1356 \nQ 3559 666 3084 287 \nQ 2609 -91 1734 -91 \nQ 1441 -91 1130 -33 \nQ 819 25 488 141 \nL 488 750 \nQ 750 597 1062 519 \nQ 1375 441 1716 441 \nQ 2309 441 2620 675 \nQ 2931 909 2931 1356 \nQ 2931 1769 2642 2001 \nQ 2353 2234 1838 2234 \nL 1294 2234 \nL 1294 2753 \nL 1863 2753 \nQ 2328 2753 2575 2939 \nQ 2822 3125 2822 3475 \nQ 2822 3834 2567 4026 \nQ 2313 4219 1838 4219 \nQ 1578 4219 1281 4162 \nQ 984 4106 628 3988 \nL 628 4550 \nQ 988 4650 1302 4700 \nQ 1616 4750 1894 4750 \nQ 2613 4750 3031 4423 \nQ 3450 4097 3450 3541 \nQ 3450 3153 3228 2886 \nQ 3006 2619 2597 2516 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-33\" transform=\"translate(95.410156 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_5\">\n     <g id=\"line2d_27\">\n      <path d=\"M 43.78125 74.516143 \nL 378.58125 74.516143 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_28\">\n      <g>\n       <use xlink:href=\"#mfc859cc296\" x=\"43.78125\" y=\"74.516143\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_15\">\n      <!-- 0.4 -->\n      <g transform=\"translate(20.878125 78.315361) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-34\" transform=\"translate(95.410156 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_6\">\n     <g id=\"line2d_29\">\n      <path d=\"M 43.78125 51.655178 \nL 378.58125 51.655178 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_30\">\n      <g>\n       <use xlink:href=\"#mfc859cc296\" x=\"43.78125\" y=\"51.655178\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_16\">\n      <!-- 0.5 -->\n      <g transform=\"translate(20.878125 55.454397) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-35\" d=\"M 691 4666 \nL 3169 4666 \nL 3169 4134 \nL 1269 4134 \nL 1269 2991 \nQ 1406 3038 1543 3061 \nQ 1681 3084 1819 3084 \nQ 2600 3084 3056 2656 \nQ 3513 2228 3513 1497 \nQ 3513 744 3044 326 \nQ 2575 -91 1722 -91 \nQ 1428 -91 1123 -41 \nQ 819 9 494 109 \nL 494 744 \nQ 775 591 1075 516 \nQ 1375 441 1709 441 \nQ 2250 441 2565 725 \nQ 2881 1009 2881 1497 \nQ 2881 1984 2565 2268 \nQ 2250 2553 1709 2553 \nQ 1456 2553 1204 2497 \nQ 953 2441 691 2322 \nL 691 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-35\" transform=\"translate(95.410156 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_7\">\n     <g id=\"line2d_31\">\n      <path d=\"M 43.78125 28.794214 \nL 378.58125 28.794214 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_32\">\n      <g>\n       <use xlink:href=\"#mfc859cc296\" x=\"43.78125\" y=\"28.794214\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_17\">\n      <!-- 0.6 -->\n      <g transform=\"translate(20.878125 32.593433) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-36\" transform=\"translate(95.410156 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_18\">\n     <!-- density -->\n     <g transform=\"translate(14.798438 108.692813) rotate(-90) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-64\" d=\"M 2906 2969 \nL 2906 4863 \nL 3481 4863 \nL 3481 0 \nL 2906 0 \nL 2906 525 \nQ 2725 213 2448 61 \nQ 2172 -91 1784 -91 \nQ 1150 -91 751 415 \nQ 353 922 353 1747 \nQ 353 2572 751 3078 \nQ 1150 3584 1784 3584 \nQ 2172 3584 2448 3432 \nQ 2725 3281 2906 2969 \nz\nM 947 1747 \nQ 947 1113 1208 752 \nQ 1469 391 1925 391 \nQ 2381 391 2643 752 \nQ 2906 1113 2906 1747 \nQ 2906 2381 2643 2742 \nQ 2381 3103 1925 3103 \nQ 1469 3103 1208 2742 \nQ 947 2381 947 1747 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-65\" d=\"M 3597 1894 \nL 3597 1613 \nL 953 1613 \nQ 991 1019 1311 708 \nQ 1631 397 2203 397 \nQ 2534 397 2845 478 \nQ 3156 559 3463 722 \nL 3463 178 \nQ 3153 47 2828 -22 \nQ 2503 -91 2169 -91 \nQ 1331 -91 842 396 \nQ 353 884 353 1716 \nQ 353 2575 817 3079 \nQ 1281 3584 2069 3584 \nQ 2775 3584 3186 3129 \nQ 3597 2675 3597 1894 \nz\nM 3022 2063 \nQ 3016 2534 2758 2815 \nQ 2500 3097 2075 3097 \nQ 1594 3097 1305 2825 \nQ 1016 2553 972 2059 \nL 3022 2063 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6e\" d=\"M 3513 2113 \nL 3513 0 \nL 2938 0 \nL 2938 2094 \nQ 2938 2591 2744 2837 \nQ 2550 3084 2163 3084 \nQ 1697 3084 1428 2787 \nQ 1159 2491 1159 1978 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1366 3272 1645 3428 \nQ 1925 3584 2291 3584 \nQ 2894 3584 3203 3211 \nQ 3513 2838 3513 2113 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-73\" d=\"M 2834 3397 \nL 2834 2853 \nQ 2591 2978 2328 3040 \nQ 2066 3103 1784 3103 \nQ 1356 3103 1142 2972 \nQ 928 2841 928 2578 \nQ 928 2378 1081 2264 \nQ 1234 2150 1697 2047 \nL 1894 2003 \nQ 2506 1872 2764 1633 \nQ 3022 1394 3022 966 \nQ 3022 478 2636 193 \nQ 2250 -91 1575 -91 \nQ 1294 -91 989 -36 \nQ 684 19 347 128 \nL 347 722 \nQ 666 556 975 473 \nQ 1284 391 1588 391 \nQ 1994 391 2212 530 \nQ 2431 669 2431 922 \nQ 2431 1156 2273 1281 \nQ 2116 1406 1581 1522 \nL 1381 1569 \nQ 847 1681 609 1914 \nQ 372 2147 372 2553 \nQ 372 3047 722 3315 \nQ 1072 3584 1716 3584 \nQ 2034 3584 2315 3537 \nQ 2597 3491 2834 3397 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-69\" d=\"M 603 3500 \nL 1178 3500 \nL 1178 0 \nL 603 0 \nL 603 3500 \nz\nM 603 4863 \nL 1178 4863 \nL 1178 4134 \nL 603 4134 \nL 603 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-74\" d=\"M 1172 4494 \nL 1172 3500 \nL 2356 3500 \nL 2356 3053 \nL 1172 3053 \nL 1172 1153 \nQ 1172 725 1289 603 \nQ 1406 481 1766 481 \nL 2356 481 \nL 2356 0 \nL 1766 0 \nQ 1100 0 847 248 \nQ 594 497 594 1153 \nL 594 3053 \nL 172 3053 \nL 172 3500 \nL 594 3500 \nL 594 4494 \nL 1172 4494 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-79\" d=\"M 2059 -325 \nQ 1816 -950 1584 -1140 \nQ 1353 -1331 966 -1331 \nL 506 -1331 \nL 506 -850 \nL 844 -850 \nQ 1081 -850 1212 -737 \nQ 1344 -625 1503 -206 \nL 1606 56 \nL 191 3500 \nL 800 3500 \nL 1894 763 \nL 2988 3500 \nL 3597 3500 \nL 2059 -325 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-64\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(63.476562 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(125 0)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(188.378906 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(240.478516 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(268.261719 0)\"/>\n      <use xlink:href=\"#DejaVuSans-79\" transform=\"translate(307.470703 0)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"line2d_33\">\n    <path d=\"M 58.999432 165.96 \nL 130.82925 165.849469 \nL 134.633795 165.600291 \nL 137.144795 165.222528 \nL 139.123159 164.705135 \nL 140.797159 164.038052 \nL 142.242886 163.230261 \nL 143.612523 162.209812 \nL 144.906068 160.964492 \nL 146.199614 159.390563 \nL 147.417068 157.557977 \nL 148.634523 155.335749 \nL 149.928068 152.496183 \nL 151.221614 149.115435 \nL 152.59125 144.903105 \nL 154.036977 139.722076 \nL 155.634886 133.119171 \nL 157.384977 124.90026 \nL 159.515523 113.75356 \nL 166.591977 75.68392 \nL 168.037705 69.659595 \nL 169.179068 65.762905 \nL 170.16825 63.093906 \nL 171.00525 61.395317 \nL 171.766159 60.320951 \nL 172.374886 59.793793 \nL 172.907523 59.578884 \nL 173.440159 59.595431 \nL 173.972795 59.843325 \nL 174.581523 60.407774 \nL 175.266341 61.395317 \nL 176.02725 62.91774 \nL 176.940341 65.307647 \nL 177.929523 68.540336 \nL 179.146977 73.336842 \nL 180.592705 80.00655 \nL 182.494977 89.937595 \nL 189.799705 129.309404 \nL 191.549795 137.003047 \nL 193.147705 143.082575 \nL 194.669523 148.011612 \nL 196.11525 151.939268 \nL 197.484886 155.026666 \nL 198.778432 157.428386 \nL 200.071977 159.386283 \nL 201.365523 160.957563 \nL 202.659068 162.19874 \nL 204.028705 163.212324 \nL 205.474432 164.008668 \nL 207.072341 164.630218 \nL 208.822432 165.071468 \nL 210.800795 165.342908 \nL 213.007432 165.425073 \nL 215.290159 165.295301 \nL 217.496795 164.947497 \nL 219.399068 164.43164 \nL 221.149159 163.732198 \nL 222.747068 162.859872 \nL 224.268886 161.780085 \nL 225.790705 160.418305 \nL 227.312523 158.739538 \nL 228.834341 156.716515 \nL 230.43225 154.205864 \nL 232.10625 151.155344 \nL 234.008523 147.212787 \nL 236.443432 141.605147 \nL 241.845886 128.962832 \nL 243.519886 125.717235 \nL 244.889523 123.535322 \nL 246.030886 122.111951 \nL 247.020068 121.203292 \nL 247.933159 120.650895 \nL 248.770159 120.39544 \nL 249.607159 120.384301 \nL 250.444159 120.617656 \nL 251.281159 121.091759 \nL 252.19425 121.874531 \nL 253.183432 123.018388 \nL 254.324795 124.688398 \nL 255.694432 127.122951 \nL 257.368432 130.60259 \nL 259.72725 136.093716 \nL 264.597068 147.544832 \nL 266.651523 151.741063 \nL 268.401614 154.840168 \nL 270.075614 157.361613 \nL 271.673523 159.366783 \nL 273.195341 160.932463 \nL 274.717159 162.193895 \nL 276.315068 163.231984 \nL 277.989068 164.054527 \nL 279.81525 164.702997 \nL 281.869705 165.196521 \nL 284.380705 165.562782 \nL 287.652614 165.802336 \nL 292.750705 165.928281 \nL 306.370977 165.959836 \nL 363.363068 165.96 \nL 363.363068 165.96 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square\"/>\n   </g>\n   <g id=\"line2d_34\">\n    <path d=\"M 58.999432 165.915466 \nL 76.42425 165.725751 \nL 86.620432 165.402643 \nL 94.229523 164.947918 \nL 100.545068 164.352806 \nL 106.023614 163.616427 \nL 110.969523 162.728191 \nL 115.534977 161.681691 \nL 119.796068 160.47837 \nL 123.904977 159.086857 \nL 127.861705 157.514718 \nL 131.818432 155.70296 \nL 135.775159 153.646571 \nL 139.884068 151.255777 \nL 144.22125 148.46648 \nL 148.938886 145.161185 \nL 154.645705 140.87795 \nL 170.700886 128.652666 \nL 174.733705 125.962011 \nL 178.233886 123.884014 \nL 181.429705 122.241447 \nL 184.39725 120.965559 \nL 187.212614 119.999479 \nL 189.875795 119.319983 \nL 192.462886 118.888368 \nL 194.973886 118.690922 \nL 197.484886 118.714461 \nL 199.995886 118.958655 \nL 202.582977 119.437399 \nL 205.246159 120.163598 \nL 207.985432 121.145735 \nL 210.876886 122.423671 \nL 213.996614 124.054574 \nL 217.420705 126.107412 \nL 221.301341 128.706973 \nL 226.018977 132.153154 \nL 233.475886 137.923432 \nL 241.921977 144.382766 \nL 247.096159 148.055773 \nL 251.585523 150.975474 \nL 255.770523 153.436781 \nL 259.803341 155.554173 \nL 263.760068 157.384579 \nL 267.792886 159.003017 \nL 271.901795 160.40766 \nL 276.086795 161.603696 \nL 280.500068 162.634229 \nL 285.217705 163.507364 \nL 290.391886 164.236869 \nL 296.250886 164.832065 \nL 303.022977 165.290469 \nL 311.392977 165.624753 \nL 322.730523 165.839105 \nL 341.524977 165.941938 \nL 363.363068 165.958574 \nL 363.363068 165.958574 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke-dasharray: 5.55,2.4; stroke-dashoffset: 0; stroke: #bf00bf; stroke-width: 1.5\"/>\n   </g>\n   <g id=\"line2d_35\">\n    <path d=\"M 58.999432 165.96 \nL 129.535705 165.850536 \nL 133.264159 165.602726 \nL 135.699068 165.229288 \nL 137.601341 164.721467 \nL 139.19925 164.071232 \nL 140.644977 163.238943 \nL 141.938523 162.238123 \nL 143.155977 161.019817 \nL 144.373432 159.476123 \nL 145.590886 157.545328 \nL 146.808341 155.161823 \nL 148.025795 152.258302 \nL 149.24325 148.768686 \nL 150.536795 144.350484 \nL 151.830341 139.139833 \nL 153.199977 132.710856 \nL 154.721795 124.444456 \nL 156.395795 114.024961 \nL 158.298068 100.684371 \nL 160.809068 81.327484 \nL 165.374523 45.922255 \nL 167.048523 34.817419 \nL 168.418159 27.186167 \nL 169.559523 22.063393 \nL 170.472614 18.880906 \nL 171.309614 16.734287 \nL 171.994432 15.551158 \nL 172.527068 14.996695 \nL 172.983614 14.779324 \nL 173.440159 14.801135 \nL 173.896705 15.062026 \nL 174.429341 15.666855 \nL 175.038068 16.749344 \nL 175.722886 18.454898 \nL 176.559886 21.214895 \nL 177.472977 25.024833 \nL 178.53825 30.432272 \nL 179.831795 38.209246 \nL 181.429705 49.268377 \nL 183.712432 66.819212 \nL 188.58225 104.613371 \nL 190.560614 118.013408 \nL 192.234614 127.931666 \nL 193.756432 135.719744 \nL 195.202159 142.028844 \nL 196.571795 147.065711 \nL 197.865341 151.041765 \nL 199.158886 154.329864 \nL 200.452432 157.007741 \nL 201.669886 159.042903 \nL 202.887341 160.675233 \nL 204.180886 162.038389 \nL 205.474432 163.086682 \nL 206.844068 163.921358 \nL 208.365886 164.590758 \nL 210.115977 165.112292 \nL 212.246523 165.501745 \nL 214.985795 165.762466 \nL 219.170795 165.911079 \nL 228.529977 165.958667 \nL 363.363068 165.96 \nL 363.363068 165.96 \n\" clip-path=\"url(#p37a0adeb25)\" style=\"fill: none; stroke-dasharray: 9.6,2.4,1.5,2.4; stroke-dashoffset: 0; stroke: #008000; stroke-width: 1.5\"/>\n   </g>\n   <g id=\"patch_3\">\n    <path d=\"M 43.78125 173.52 \nL 43.78125 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_4\">\n    <path d=\"M 378.58125 173.52 \nL 378.58125 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_5\">\n    <path d=\"M 43.78125 173.52 \nL 378.58125 173.52 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_6\">\n    <path d=\"M 43.78125 7.2 \nL 378.58125 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"legend_1\">\n    <g id=\"patch_7\">\n     <path d=\"M 250.251563 59.234375 \nL 371.58125 59.234375 \nQ 373.58125 59.234375 373.58125 57.234375 \nL 373.58125 14.2 \nQ 373.58125 12.2 371.58125 12.2 \nL 250.251563 12.2 \nQ 248.251563 12.2 248.251563 14.2 \nL 248.251563 57.234375 \nQ 248.251563 59.234375 250.251563 59.234375 \nz\n\" style=\"fill: #ffffff; opacity: 0.8; stroke: #cccccc; stroke-linejoin: miter\"/>\n    </g>\n    <g id=\"line2d_36\">\n     <path d=\"M 252.251563 20.298438 \nL 262.251562 20.298438 \nL 272.251562 20.298438 \n\" style=\"fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square\"/>\n    </g>\n    <g id=\"text_19\">\n     <!-- target P (mixture) -->\n     <g transform=\"translate(280.251562 23.798438) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-61\" d=\"M 2194 1759 \nQ 1497 1759 1228 1600 \nQ 959 1441 959 1056 \nQ 959 750 1161 570 \nQ 1363 391 1709 391 \nQ 2188 391 2477 730 \nQ 2766 1069 2766 1631 \nL 2766 1759 \nL 2194 1759 \nz\nM 3341 1997 \nL 3341 0 \nL 2766 0 \nL 2766 531 \nQ 2569 213 2275 61 \nQ 1981 -91 1556 -91 \nQ 1019 -91 701 211 \nQ 384 513 384 1019 \nQ 384 1609 779 1909 \nQ 1175 2209 1959 2209 \nL 2766 2209 \nL 2766 2266 \nQ 2766 2663 2505 2880 \nQ 2244 3097 1772 3097 \nQ 1472 3097 1187 3025 \nQ 903 2953 641 2809 \nL 641 3341 \nQ 956 3463 1253 3523 \nQ 1550 3584 1831 3584 \nQ 2591 3584 2966 3190 \nQ 3341 2797 3341 1997 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-72\" d=\"M 2631 2963 \nQ 2534 3019 2420 3045 \nQ 2306 3072 2169 3072 \nQ 1681 3072 1420 2755 \nQ 1159 2438 1159 1844 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1341 3275 1631 3429 \nQ 1922 3584 2338 3584 \nQ 2397 3584 2469 3576 \nQ 2541 3569 2628 3553 \nL 2631 2963 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-67\" d=\"M 2906 1791 \nQ 2906 2416 2648 2759 \nQ 2391 3103 1925 3103 \nQ 1463 3103 1205 2759 \nQ 947 2416 947 1791 \nQ 947 1169 1205 825 \nQ 1463 481 1925 481 \nQ 2391 481 2648 825 \nQ 2906 1169 2906 1791 \nz\nM 3481 434 \nQ 3481 -459 3084 -895 \nQ 2688 -1331 1869 -1331 \nQ 1566 -1331 1297 -1286 \nQ 1028 -1241 775 -1147 \nL 775 -588 \nQ 1028 -725 1275 -790 \nQ 1522 -856 1778 -856 \nQ 2344 -856 2625 -561 \nQ 2906 -266 2906 331 \nL 2906 616 \nQ 2728 306 2450 153 \nQ 2172 0 1784 0 \nQ 1141 0 747 490 \nQ 353 981 353 1791 \nQ 353 2603 747 3093 \nQ 1141 3584 1784 3584 \nQ 2172 3584 2450 3431 \nQ 2728 3278 2906 2969 \nL 2906 3500 \nL 3481 3500 \nL 3481 434 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-20\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-50\" d=\"M 1259 4147 \nL 1259 2394 \nL 2053 2394 \nQ 2494 2394 2734 2622 \nQ 2975 2850 2975 3272 \nQ 2975 3691 2734 3919 \nQ 2494 4147 2053 4147 \nL 1259 4147 \nz\nM 628 4666 \nL 2053 4666 \nQ 2838 4666 3239 4311 \nQ 3641 3956 3641 3272 \nQ 3641 2581 3239 2228 \nQ 2838 1875 2053 1875 \nL 1259 1875 \nL 1259 0 \nL 628 0 \nL 628 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-28\" d=\"M 1984 4856 \nQ 1566 4138 1362 3434 \nQ 1159 2731 1159 2009 \nQ 1159 1288 1364 580 \nQ 1569 -128 1984 -844 \nL 1484 -844 \nQ 1016 -109 783 600 \nQ 550 1309 550 2009 \nQ 550 2706 781 3412 \nQ 1013 4119 1484 4856 \nL 1984 4856 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6d\" d=\"M 3328 2828 \nQ 3544 3216 3844 3400 \nQ 4144 3584 4550 3584 \nQ 5097 3584 5394 3201 \nQ 5691 2819 5691 2113 \nL 5691 0 \nL 5113 0 \nL 5113 2094 \nQ 5113 2597 4934 2840 \nQ 4756 3084 4391 3084 \nQ 3944 3084 3684 2787 \nQ 3425 2491 3425 1978 \nL 3425 0 \nL 2847 0 \nL 2847 2094 \nQ 2847 2600 2669 2842 \nQ 2491 3084 2119 3084 \nQ 1678 3084 1418 2786 \nQ 1159 2488 1159 1978 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1356 3278 1631 3431 \nQ 1906 3584 2284 3584 \nQ 2666 3584 2933 3390 \nQ 3200 3197 3328 2828 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-75\" d=\"M 544 1381 \nL 544 3500 \nL 1119 3500 \nL 1119 1403 \nQ 1119 906 1312 657 \nQ 1506 409 1894 409 \nQ 2359 409 2629 706 \nQ 2900 1003 2900 1516 \nL 2900 3500 \nL 3475 3500 \nL 3475 0 \nL 2900 0 \nL 2900 538 \nQ 2691 219 2414 64 \nQ 2138 -91 1772 -91 \nQ 1169 -91 856 284 \nQ 544 659 544 1381 \nz\nM 1991 3584 \nL 1991 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-29\" d=\"M 513 4856 \nL 1013 4856 \nQ 1481 4119 1714 3412 \nQ 1947 2706 1947 2009 \nQ 1947 1309 1714 600 \nQ 1481 -109 1013 -844 \nL 513 -844 \nQ 928 -128 1133 580 \nQ 1338 1288 1338 2009 \nQ 1338 2731 1133 3434 \nQ 928 4138 513 4856 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-74\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(39.208984 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(100.488281 0)\"/>\n      <use xlink:href=\"#DejaVuSans-67\" transform=\"translate(139.851562 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(203.328125 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(264.851562 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(304.060547 0)\"/>\n      <use xlink:href=\"#DejaVuSans-50\" transform=\"translate(335.847656 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(396.150391 0)\"/>\n      <use xlink:href=\"#DejaVuSans-28\" transform=\"translate(427.9375 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6d\" transform=\"translate(466.951172 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(564.363281 0)\"/>\n      <use xlink:href=\"#DejaVuSans-78\" transform=\"translate(592.146484 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(651.326172 0)\"/>\n      <use xlink:href=\"#DejaVuSans-75\" transform=\"translate(690.535156 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(753.914062 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(792.777344 0)\"/>\n      <use xlink:href=\"#DejaVuSans-29\" transform=\"translate(854.300781 0)\"/>\n     </g>\n    </g>\n    <g id=\"line2d_37\">\n     <path d=\"M 252.251563 34.976563 \nL 262.251562 34.976563 \nL 272.251562 34.976563 \n\" style=\"fill: none; stroke-dasharray: 5.55,2.4; stroke-dashoffset: 0; stroke: #bf00bf; stroke-width: 1.5\"/>\n    </g>\n    <g id=\"text_20\">\n     <!-- forward-KL fit -->\n     <g transform=\"translate(280.251562 38.476563) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-66\" d=\"M 2375 4863 \nL 2375 4384 \nL 1825 4384 \nQ 1516 4384 1395 4259 \nQ 1275 4134 1275 3809 \nL 1275 3500 \nL 2222 3500 \nL 2222 3053 \nL 1275 3053 \nL 1275 0 \nL 697 0 \nL 697 3053 \nL 147 3053 \nL 147 3500 \nL 697 3500 \nL 697 3744 \nQ 697 4328 969 4595 \nQ 1241 4863 1831 4863 \nL 2375 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6f\" d=\"M 1959 3097 \nQ 1497 3097 1228 2736 \nQ 959 2375 959 1747 \nQ 959 1119 1226 758 \nQ 1494 397 1959 397 \nQ 2419 397 2687 759 \nQ 2956 1122 2956 1747 \nQ 2956 2369 2687 2733 \nQ 2419 3097 1959 3097 \nz\nM 1959 3584 \nQ 2709 3584 3137 3096 \nQ 3566 2609 3566 1747 \nQ 3566 888 3137 398 \nQ 2709 -91 1959 -91 \nQ 1206 -91 779 398 \nQ 353 888 353 1747 \nQ 353 2609 779 3096 \nQ 1206 3584 1959 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-77\" d=\"M 269 3500 \nL 844 3500 \nL 1563 769 \nL 2278 3500 \nL 2956 3500 \nL 3675 769 \nL 4391 3500 \nL 4966 3500 \nL 4050 0 \nL 3372 0 \nL 2619 2869 \nL 1863 0 \nL 1184 0 \nL 269 3500 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-2d\" d=\"M 313 2009 \nL 1997 2009 \nL 1997 1497 \nL 313 1497 \nL 313 2009 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-4b\" d=\"M 628 4666 \nL 1259 4666 \nL 1259 2694 \nL 3353 4666 \nL 4166 4666 \nL 1850 2491 \nL 4331 0 \nL 3500 0 \nL 1259 2247 \nL 1259 0 \nL 628 0 \nL 628 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-4c\" d=\"M 628 4666 \nL 1259 4666 \nL 1259 531 \nL 3531 531 \nL 3531 0 \nL 628 0 \nL 628 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-66\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(35.205078 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(96.386719 0)\"/>\n      <use xlink:href=\"#DejaVuSans-77\" transform=\"translate(137.5 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(219.287109 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(280.566406 0)\"/>\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(319.929688 0)\"/>\n      <use xlink:href=\"#DejaVuSans-2d\" transform=\"translate(383.40625 0)\"/>\n      <use xlink:href=\"#DejaVuSans-4b\" transform=\"translate(419.490234 0)\"/>\n      <use xlink:href=\"#DejaVuSans-4c\" transform=\"translate(485.066406 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(540.779297 0)\"/>\n      <use xlink:href=\"#DejaVuSans-66\" transform=\"translate(572.566406 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(607.771484 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(635.554688 0)\"/>\n     </g>\n    </g>\n    <g id=\"line2d_38\">\n     <path d=\"M 252.251563 49.654688 \nL 262.251562 49.654688 \nL 272.251562 49.654688 \n\" style=\"fill: none; stroke-dasharray: 9.6,2.4,1.5,2.4; stroke-dashoffset: 0; stroke: #008000; stroke-width: 1.5\"/>\n    </g>\n    <g id=\"text_21\">\n     <!-- reverse-KL fit -->\n     <g transform=\"translate(280.251562 53.154688) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-76\" d=\"M 191 3500 \nL 800 3500 \nL 1894 563 \nL 2988 3500 \nL 3597 3500 \nL 2284 0 \nL 1503 0 \nL 191 3500 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-72\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(38.863281 0)\"/>\n      <use xlink:href=\"#DejaVuSans-76\" transform=\"translate(100.386719 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(159.566406 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(221.089844 0)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(262.203125 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(314.302734 0)\"/>\n      <use xlink:href=\"#DejaVuSans-2d\" transform=\"translate(375.826172 0)\"/>\n      <use xlink:href=\"#DejaVuSans-4b\" transform=\"translate(411.910156 0)\"/>\n      <use xlink:href=\"#DejaVuSans-4c\" transform=\"translate(477.486328 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(533.199219 0)\"/>\n      <use xlink:href=\"#DejaVuSans-66\" transform=\"translate(564.986328 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(600.191406 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(627.974609 0)\"/>\n     </g>\n    </g>\n   </g>\n  </g>\n </g>\n <defs>\n  <clipPath id=\"p37a0adeb25\">\n   <rect x=\"43.78125\" y=\"7.2\" width=\"334.8\" height=\"166.32\"/>\n  </clipPath>\n </defs>\n</svg>\n"
     },
     "metadata": {}
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-5 fw=all -->\n",
    "\n",
    "The two objectives choose *different Gaussians for the same target*. The\n",
    "forward fit lands at $\\mu = -0.800$, $\\sigma = 1.929$, precisely the\n",
    "mixture's mean $0.7(-2) + 0.3(2) = -0.8$ and standard deviation\n",
    "$\\sqrt{3.72} \\approx 1.929$, as the proposition demands: a broad Gaussian\n",
    "draped across both modes, with substantial mass in the valley between them\n",
    "where $P$ has almost none. The reverse fit lands at $\\mu = -1.998$,\n",
    "$\\sigma = 0.603$: it *is* the dominant component, and the minor mode is\n",
    "simply dropped. Its divergence, $0.356$ nats, is essentially\n",
    "$\\log(1/0.7) \\approx 0.357$, the divergence incurred by treating the $70\\%$\n",
    "component as the whole distribution. And reverse KL is genuinely multimodal\n",
    "as an objective: restarting the optimizer near the minor mode converges to a\n",
    "*second* local optimum at $\\mu = +1.995$ with KL\n",
    "$\\approx 1.202 \\approx \\log(1/0.3)$ nats. Which local optimum a variational\n",
    "method finds depends on initialization, a failure mode every practitioner of\n",
    "variational inference eventually meets.\n",
    "\n",
    "The generative-modeling translation: maximum-likelihood families\n",
    "(autoregressive models, normalizing flows) inherit forward KL's\n",
    "mass-covering: they rarely miss a mode but can produce blurry,\n",
    "over-dispersed samples; adversarial and variational objectives built on\n",
    "reverse-type divergences produce sharp samples but can silently drop modes.\n",
    "The capstone table in that section files every\n",
    "major objective under exactly this dichotomy.\n",
    "\n",
    "## Metrics: Total Variation, MMD, and Optimal Transport\n",
    "\n",
    "The f-divergence family compares densities *pointwise* through the ratio\n",
    "$p/q$. This section develops the complementary view: divergences defined\n",
    "through *test functions* and *transport*, which see the geometry of the\n",
    "sample space and remain estimable and informative when density ratios are\n",
    "degenerate or unavailable.\n",
    "\n",
    "### Total Variation and Pinsker's Inequality\n",
    "\n",
    "The most interpretable distance between distributions answers the question:\n",
    "*what is the largest disagreement in probability that $P$ and $Q$ assign to\n",
    "any event?*\n",
    "\n",
    "$$\n",
    "\\textrm{TV}(P, Q) = \\sup_{A} \\,\\lvert P(A) - Q(A) \\rvert.\n",
    "$$\n",
    "\n",
    "If $\\textrm{TV}(P, Q) = 0.03$, then *no event you can test* will see the two\n",
    "distributions differ by more than three percentage points. TV is a genuine metric (symmetry is\n",
    "visible in the equation; the triangle inequality is Exercise 4),\n",
    "and the supremum has a closed form.\n",
    "\n",
    "**Proposition (TV is half the $\\ell_1$ distance).** *For discrete $P$, $Q$,*\n",
    "\n",
    "$$\n",
    "\\textrm{TV}(P, Q) = \\frac{1}{2} \\sum_x \\lvert p(x) - q(x) \\rvert,\n",
    "$$\n",
    "\n",
    "*and the supremum in the equation is attained at\n",
    "$A^\\star = \\{x : p(x) > q(x)\\}$.*\n",
    "\n",
    "**Proof.** For any event $A$,\n",
    "$P(A) - Q(A) = \\sum_{x \\in A} (p(x) - q(x)) \\leq \\sum_{x \\in A^\\star} (p(x) - q(x))$:\n",
    "enlarging $A$ to include every $x$ with $p(x) > q(x)$ and discarding the rest\n",
    "only adds non-negative terms and removes non-positive ones. Since\n",
    "$\\sum_x (p(x) - q(x)) = 0$, the positive part equals the negative part in\n",
    "magnitude, so\n",
    "\n",
    "$$\n",
    "\\sum_{x \\in A^\\star} (p(x) - q(x))\n",
    "= \\frac{1}{2} \\sum_x \\lvert p(x) - q(x) \\rvert.\n",
    "$$\n",
    "\n",
    "The same argument bounds $Q(A) - P(A)$ by the same quantity. $\\blacksquare$\n",
    "\n",
    "the figure shows the picture: TV is half the total area where\n",
    "the two densities disagree, and the optimal distinguishing event is simply\n",
    "\"the region where $P$ is the better explanation\". This gives TV its\n",
    "operational meaning. Hand a tester one sample, drawn from $P$ or $Q$ with\n",
    "equal probability, and ask which distribution produced it: the best possible\n",
    "test (guess $P$ exactly on $A^\\star$) succeeds with probability\n",
    "$\\tfrac{1}{2}\\big(1 + \\textrm{TV}(P, Q)\\big)$, an excess of\n",
    "$\\textrm{TV}/2$ over coin-flipping, which cryptographers double and call the\n",
    "*advantage*, so that under their convention the best achievable advantage is\n",
    "exactly $\\textrm{TV}(P, Q)$.\n",
    "\n",
    "![Two densities and the region where they disagree. The total variation distance is half the shaded area, and the event that attains the supremum is the set where the first density exceeds the second: betting on P exactly there is the best possible single-sample test.](../img/mdl-it-tv-area.svg)\n",
    "\n",
    "So TV bounds what any test can detect, which makes the next inequality\n",
    "consequential, because it says *KL bounds TV*. Train a model to small KL\n",
    "divergence and you have certified indistinguishability against every test.\n",
    "\n",
    "**Proposition (Pinsker's inequality).** [@Pinsker.1964]\n",
    "\n",
    "$$\n",
    "\\textrm{TV}(P, Q) \\;\\leq\\; \\sqrt{ \\tfrac{1}{2}\\, D_{\\textrm{KL}}(P\\|Q) }.\n",
    "$$\n",
    "\n",
    "**Proof.** *Step 1: two outcomes.* For $a, b \\in (0, 1)$ let\n",
    "$d(a\\|b) = a \\log\\frac{a}{b} + (1-a) \\log\\frac{1-a}{1-b}$ be the KL divergence\n",
    "between coins with heads-probabilities $a$ and $b$; their TV distance is\n",
    "$|a - b|$, so we must show $d(a\\|b) \\geq 2(a-b)^2$. Fix $b$ and let\n",
    "$h(a) = d(a\\|b) - 2(a - b)^2$. Then $h(b) = 0$, $h'(b) = 0$, and\n",
    "\n",
    "$$\n",
    "h''(a) = \\frac{1}{a(1-a)} - 4 \\geq 0,\n",
    "$$\n",
    "\n",
    "since $a(1-a) \\leq \\tfrac{1}{4}$. A convex function with value and slope zero\n",
    "at $a = b$ is non-negative everywhere.\n",
    "\n",
    "*Step 2: reduction to two outcomes.* Let $A^\\star = \\{p > q\\}$, and set\n",
    "$a = P(A^\\star)$, $b = Q(A^\\star)$, so that\n",
    "$\\textrm{TV}(P, Q) = a - b$ by the previous proposition. Merging the outcomes\n",
    "inside $A^\\star$ and inside its complement can only *decrease* KL: by the\n",
    "log-sum inequality (for non-negative numbers,\n",
    "$\\sum_i a_i \\log\\frac{a_i}{b_i} \\geq \\big(\\sum_i a_i\\big) \\log \\frac{\\sum_i a_i}{\\sum_i b_i}$,\n",
    "itself one application of Jensen to $t \\mapsto t\\log t$), applied separately\n",
    "to the terms in $A^\\star$ and in its complement,\n",
    "\n",
    "$$\n",
    "D_{\\textrm{KL}}(P\\|Q) \\;\\geq\\; d(a\\|b) \\;\\geq\\; 2(a - b)^2\n",
    "= 2\\,\\textrm{TV}(P, Q)^2,\n",
    "$$\n",
    "\n",
    "which rearranges to the equation. $\\blacksquare$\n",
    "\n",
    "The merging step, that coarsening the outcome space cannot increase an\n",
    "f-divergence, holds in full generality, and the proof above already contains\n",
    "the general argument.\n",
    "\n",
    "**Remark (data-processing for f-divergences).** *Passing $P$ and $Q$ through\n",
    "any channel $K$ (any deterministic or random map from $x$ to $y$, with\n",
    "output distributions $(pK)(y) = \\sum_x K(y \\mid x)\\, p(x)$ and likewise\n",
    "$qK$) can only lose distinguishability:*\n",
    "\n",
    "$$\n",
    "D_f(PK \\,\\|\\, QK) \\;\\leq\\; D_f(P\\|Q).\n",
    "$$\n",
    "\n",
    "**Proof.** Since $\\sum_y K(y \\mid x) = 1$,\n",
    "$D_f(P\\|Q) = \\sum_y \\sum_x K(y \\mid x)\\, q(x)\\, f\\big(p(x)/q(x)\\big)$. For\n",
    "each fixed $y$, Jensen's inequality with weights proportional to\n",
    "$K(y \\mid x)\\, q(x)$, applied at the points $p(x)/q(x)$, moves $f$ outside\n",
    "the inner sum and leaves exactly the $y$-th term of $D_f(PK\\|QK)$; summing\n",
    "over $y$ with $(qK)(y) > 0$ finishes, since terms with $(qK)(y) = 0$\n",
    "contribute nothing. $\\blacksquare$\n",
    "\n",
    "Merging outcomes is the deterministic special case used above, and the\n",
    "data-processing inequality for *mutual information*, which\n",
    "that section states and proves, is the same\n",
    "principle in its best-known form. Note also what Pinsker does *not* say: it has no useful converse.\n",
    "TV is bounded by $1$ while KL is unbounded, so the bound goes slack for\n",
    "distant pairs (two unit-variance Gaussians $50$ apart have\n",
    "$\\textrm{TV} \\approx 1$ but KL $= 1250$ nats), and small TV does *not* imply\n",
    "small KL (a model can assign $q = 0$ to a rare event and have infinite KL at\n",
    "tiny TV). Let's check both the bound and its tightness numerically: first\n",
    "over $10{,}000$ random pairs of distributions on five outcomes, then on pairs\n",
    "of coins approaching each other, where the binary-case analysis says the\n",
    "ratio should approach $1$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-pinsker"
   },
   "source": [
    "rng = np.random.default_rng(0)\n",
    "worst = 0.0\n",
    "for _ in range(10000):\n",
    "    pr, qr = rng.dirichlet(np.ones(5)), rng.dirichlet(np.ones(5))\n",
    "    tv = 0.5 * np.abs(pr - qr).sum()\n",
    "    kl_pq = np.sum(pr * np.log(pr / qr))\n",
    "    worst = max(worst, tv / np.sqrt(0.5 * kl_pq))\n",
    "print(f'max TV / sqrt(KL/2) over 10,000 random pairs: {worst:.4f}')\n",
    "\n",
    "for eps in (0.1, 0.01, 0.001):\n",
    "    pr, qr = np.array([0.5, 0.5]), np.array([0.5 + eps, 0.5 - eps])\n",
    "    tv = 0.5 * np.abs(pr - qr).sum()\n",
    "    kl_pq = np.sum(pr * np.log(pr / qr))\n",
    "    print(f'coins 1/2 vs 1/2+{eps}: TV / sqrt(KL/2) = '\n",
    "          f'{tv / np.sqrt(0.5 * kl_pq):.6f}')"
   ],
   "execution_count": 5,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "max TV / sqrt(KL/2) over 10,000 random pairs: 0.9926\ncoins 1/2 vs 1/2+0.1: TV / sqrt(KL/2) = 0.989881\ncoins 1/2 vs 1/2+0.01: TV / sqrt(KL/2) = 0.999900\ncoins 1/2 vs 1/2+0.001: TV / sqrt(KL/2) = 0.999999\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-6 fw=all -->\n",
    "\n",
    "The ratio $\\textrm{TV}/\\sqrt{\\textrm{KL}/2}$ never exceeds $1$ (the worst of\n",
    "$10{,}000$ random pairs reaches $0.9926$), and on nearly-fair coins it climbs\n",
    "to $0.999999$: the constant $\\tfrac{1}{2}$ is the best possible (Exercise 5),\n",
    "a sharpening due to Csiszár [@Csiszar.1967] rather than to Pinsker's\n",
    "original argument.\n",
    "\n",
    "### Integral Probability Metrics and MMD\n",
    "\n",
    "Total variation's definition, a supremum of differences over events,\n",
    "generalizes wholesale. Replace indicator functions of events by an arbitrary\n",
    "class $\\mathcal{F}$ of test functions and you get the **integral probability\n",
    "metrics** [@Muller.1997]:\n",
    "\n",
    "$$\n",
    "\\textrm{IPM}_{\\mathcal{F}}(P, Q)\n",
    "= \\sup_{f \\in \\mathcal{F}} \\,\\big( E_{x \\sim P}[f(x)] - E_{x \\sim Q}[f(x)] \\big).\n",
    "$$\n",
    "\n",
    "The class $\\mathcal{F}$ is a panel of auditors; the IPM reports the largest\n",
    "discrepancy any auditor in the panel can certify. Three choices of panel give\n",
    "three famous distances: bounded functions\n",
    "$\\mathcal{F} = \\{f : \\|f\\|_\\infty \\leq \\tfrac{1}{2}\\}$ recover total\n",
    "variation; 1-Lipschitz functions give the Wasserstein-1 distance (next\n",
    "subsection); and the unit ball of a *reproducing kernel Hilbert space* gives\n",
    "the **maximum mean discrepancy** [@Gretton.Borgwardt.Rasch.ea.2012].\n",
    "Contrast the structure with the equation:\n",
    "f-divergences integrate a function of the *density ratio* and need densities;\n",
    "IPMs compare *expectations* and need only samples.\n",
    "\n",
    "MMD is the member built for computation. A **kernel** $k(x, y)$ is a\n",
    "symmetric, positive-definite similarity function; the standard example is\n",
    "the RBF kernel $k(x,y) = \\exp(-\\|x - y\\|^2 / (2\\ell^2))$, which scores two\n",
    "points by how close they are on the length scale $\\ell$. Every such kernel\n",
    "generates a *reproducing kernel Hilbert space* (RKHS) $\\mathcal{H}$, a space\n",
    "of functions in which $k(x, \\cdot)$ evaluates: $f(x) = \\langle f,\n",
    "k(x, \\cdot) \\rangle_{\\mathcal{H}}$, the *reproducing property*. (Kernels and\n",
    "their Hilbert spaces are developed at length in [@Scholkopf.Smola.2002];\n",
    "these two facts are all we need.) Every distribution gets a *mean embedding*\n",
    "$\\mu_P = E_{x \\sim P}[k(x, \\cdot)] \\in \\mathcal{H}$, its average feature\n",
    "function; for bounded kernels the embedding exists and expectations commute\n",
    "with inner products, so for $f$ in the unit\n",
    "ball the reproducing property turns expectations into inner products,\n",
    "$E_P[f] - E_Q[f] = \\langle f, \\mu_P - \\mu_Q \\rangle_{\\mathcal{H}}$. The\n",
    "supremum over the unit ball of an inner product against a fixed vector is\n",
    "that vector's norm, so the IPM collapses to\n",
    "$\\textrm{MMD}(P, Q) = \\|\\mu_P - \\mu_Q\\|_{\\mathcal{H}}$, and squaring expands\n",
    "the norm into three kernel expectations:\n",
    "\n",
    "$$\n",
    "\\textrm{MMD}^2(P, Q)\n",
    "= E_{x, x' \\sim P}[k(x, x')] + E_{y, y' \\sim Q}[k(y, y')]\n",
    "- 2\\, E_{x \\sim P, y \\sim Q}[k(x, y)].\n",
    "$$\n",
    "\n",
    "Within-sample similarity under $P$, plus within-sample similarity under $Q$,\n",
    "minus twice the across-sample similarity: if the two samples interleave, the\n",
    "three terms cancel; if they form separate clumps, the within terms beat the\n",
    "across term. Replacing expectations by sample averages, *excluding* the\n",
    "diagonal terms $k(x_i, x_i)$ whose inclusion would bias the within-sample\n",
    "terms upward, gives the standard unbiased estimator\n",
    "[@Gretton.Borgwardt.Rasch.ea.2012], computable in a few lines with no\n",
    "optimization and no densities. For *characteristic* kernels (the\n",
    "RBF kernel is one) the embedding $P \\mapsto \\mu_P$ is injective, so MMD is a\n",
    "genuine metric: zero only at equality."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-mmd"
   },
   "source": [
    "def mmd2_unbiased(x, y, ell=1.0):\n",
    "    \"\"\"Unbiased MMD^2 with the RBF kernel exp(-(a-b)^2 / (2 ell^2)).\"\"\"\n",
    "    k = lambda a, b: np.exp(-(a[:, None] - b[None, :]) ** 2 / (2 * ell ** 2))\n",
    "    kxx, kyy, kxy = k(x, x), k(y, y), k(x, y)\n",
    "    n, m = len(x), len(y)\n",
    "    return ((kxx.sum() - np.trace(kxx)) / (n * (n - 1))\n",
    "            + (kyy.sum() - np.trace(kyy)) / (m * (m - 1))\n",
    "            - 2 * kxy.mean())\n",
    "\n",
    "rng = np.random.default_rng(3)\n",
    "n = 250\n",
    "x1, y_same = rng.standard_normal(n), rng.standard_normal(n)\n",
    "y_shift = rng.standard_normal(n) + 0.5\n",
    "print(f'MMD^2, same distribution : {mmd2_unbiased(x1, y_same):+.5f}')\n",
    "print(f'MMD^2, mean shifted 0.5  : {mmd2_unbiased(x1, y_shift):+.5f}')"
   ],
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "MMD^2, same distribution : +0.00054\nMMD^2, mean shifted 0.5  : +0.05890\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-7 fw=all -->\n",
    "\n",
    "Two samples from the same standard Gaussian give\n",
    "$\\textrm{MMD}^2 \\approx 0.0005$ (noise around zero; the unbiased\n",
    "estimator is even allowed to go slightly negative), while shifting one\n",
    "sample's mean by half a standard deviation produces $\\approx 0.059$, two\n",
    "orders of magnitude larger and consistent, up to sampling noise at\n",
    "$n = 250$, with the population value $\\approx 0.047$ that Exercise 8\n",
    "derives in closed form. This\n",
    "sample-only, optimization-free property is why MMD powers kernel two-sample\n",
    "tests and adversary-free generative training (MMD-GANs and generative moment\n",
    "matching [@Li.Chang.Cheng.ea.2017; @Li.Swersky.Zemel.2015]): the \"critic\" is the whole RKHS ball\n",
    "at once, and the equation evaluates its supremum in closed form.\n",
    "\n",
    "### Optimal Transport and the Wasserstein Distance\n",
    "\n",
    "Every divergence so far shares a blind spot. Let $P = \\delta_0$ be a point\n",
    "mass at the origin and $Q_d = \\delta_d$ a point mass at distance $d$. For any\n",
    "$d \\neq 0$ the supports are disjoint, so the density ratio is degenerate\n",
    "everywhere and *every f-divergence is a constant*:\n",
    "$D_{\\textrm{KL}} = \\infty$, $\\textrm{TV} = 1$, $\\textrm{JS} = \\log 2$,\n",
    "whether $d = 0.01$ or $d = 100$. A generator early in training is in exactly\n",
    "this situation (its samples and the data occupy disjoint slivers of image\n",
    "space), and an objective that is *constant* in the generator's parameters\n",
    "supplies zero gradient. This is the vanishing-gradient pathology of GAN\n",
    "training, and it is a property of the divergence, not of the optimizer\n",
    "[@Arjovsky.Chintala.Bottou.2017].\n",
    "\n",
    "What survives disjoint supports is geometry: $\\delta_{0.01}$ is *near*\n",
    "$\\delta_0$ because mass need only move $0.01$. The **Wasserstein-1\n",
    "distance** (earth-mover's distance) makes this precise. A **coupling**\n",
    "$\\gamma$ of $P$ and $Q$ is a joint distribution over pairs $(x, y)$ with\n",
    "marginals $P$ and $Q$ (a transport plan specifying how much mass travels\n",
    "from each source to each destination, as in\n",
    "the figure), and\n",
    "\n",
    "$$\n",
    "W_1(P, Q) = \\inf_{\\gamma \\in \\Pi(P, Q)} E_{(x, y) \\sim \\gamma}\\big[ \\|x - y\\| \\big],\n",
    "$$\n",
    "\n",
    "the cheapest total mass-times-distance over all plans, for $P$ and $Q$ with\n",
    "finite first moments (without that hypothesis the infimum may be infinite).\n",
    "On the point masses,\n",
    "$W_1(\\delta_0, \\delta_d) = |d|$: smooth in $d$, gradient $\\pm 1$, exactly\n",
    "the training signal the f-divergences withheld.\n",
    "\n",
    "![Optimal transport between two distributions on the line. Left: the transport plan moves each parcel of mass from source to destination, and the Wasserstein-1 distance is the minimum total mass-times-distance; the shaded area between the two CDFs computes it exactly. Right: for two point masses at separation d, the Jensen--Shannon divergence is frozen at log 2 for every nonzero d while the Wasserstein distance grows linearly: only transport gives a gradient.](../img/mdl-it-ot-transport-plan.svg)\n",
    "\n",
    "The primal the equation is an optimization over joint\n",
    "distributions, hard to even parameterize from samples. Its dual is the\n",
    "reason $W_1$ matters to deep learning.\n",
    "\n",
    "**Proposition (Kantorovich--Rubinstein duality).**\n",
    "\n",
    "$$\n",
    "W_1(P, Q) = \\sup_{\\|f\\|_{\\textrm{Lip}} \\leq 1}\n",
    "\\big( E_{x \\sim P}[f(x)] - E_{y \\sim Q}[f(y)] \\big),\n",
    "$$\n",
    "\n",
    "*where the supremum runs over all 1-Lipschitz functions\n",
    "($|f(x) - f(y)| \\leq \\|x - y\\|$ for all $x, y$).*\n",
    "\n",
    "We will not prove the general statement (see\n",
    "@Peyre.Cuturi.2019 for a proof and for optimal transport at large),\n",
    "but here is why it is plausible, in one economic metaphor that goes back to\n",
    "Kantorovich. Interpret $f(x)$ as a *price* paid per unit of mass at location\n",
    "$x$. A logistics contractor offers: \"rather than moving the mass yourself,\n",
    "sell it to me at your sources and buy it back at your destinations; I pocket\n",
    "$E_P[f] - E_Q[f]$.\" You would accept any offer where the price gradient is\n",
    "nowhere steeper than the transport cost (the 1-Lipschitz\n",
    "condition), because then no parcel's price difference exceeds its shipping\n",
    "cost, term by term, under any plan you might have used; this gives weak\n",
    "duality ($\\textrm{sup} \\leq \\textrm{inf}$, compare\n",
    "that section). The content of the theorem is that\n",
    "the *best* contractor extracts the *entire* transport cost. So $W_1$ is the\n",
    "1-Lipschitz IPM: the test class in the equation becomes the\n",
    "functions with slope at most $1$. The **WGAN** [@Arjovsky.Chintala.Bottou.2017]\n",
    "trains exactly this dual: a neural critic plays the role of $f$, constrained\n",
    "to be (approximately) 1-Lipschitz, by weight clipping originally, then by\n",
    "gradient penalties [@Gulrajani.Ahmed.Arjovsky.ea.2017] (penalizing\n",
    "$\\|\\nabla f\\| \\neq 1$, since an optimal Kantorovich potential has unit slope\n",
    "along transport rays) and by spectral normalization\n",
    "[@Miyato.Kataoka.Koyama.ea.2018], which caps each layer's Lipschitz\n",
    "constant through its largest singular value\n",
    "(that section).\n",
    "\n",
    "In one dimension, optimal transport collapses to a formula with no\n",
    "optimization left in it.\n",
    "\n",
    "**Proposition (1-D Wasserstein in closed form).** *For distributions on\n",
    "$\\mathbb{R}$ with CDFs $F_P$ and $F_Q$ and finite first moments, granting\n",
    "the duality the equation,*\n",
    "\n",
    "$$\n",
    "W_1(P, Q) = \\int_{-\\infty}^{\\infty} \\lvert F_P(t) - F_Q(t) \\rvert \\, dt.\n",
    "$$\n",
    "\n",
    "**Proof.** Let $f$ be 1-Lipschitz; then $f$ is absolutely continuous with\n",
    "$|f'(t)| \\leq 1$ wherever the derivative exists, and\n",
    "$f(x) = f(0) + \\int_0^x f'(t)\\, dt$, a standard fact of real analysis\n",
    "[@Folland.1999] that we grant. The constant $f(0)$ cancels in\n",
    "$E_P[f] - E_Q[f]$. For $X \\sim P$, writing the inner integral with\n",
    "indicators, $\\int_0^X f'(t)\\,dt = \\int_0^{\\infty} f'(t)\\, \\mathbf{1}[X > t]\\, dt\n",
    "- \\int_{-\\infty}^0 f'(t)\\, \\mathbf{1}[X \\leq t]\\, dt$, and swapping\n",
    "expectation and integral (Fubini's theorem from\n",
    "that section, justified by the moment assumption),\n",
    "\n",
    "$$\n",
    "E_P\\!\\left[ \\int_0^X f'(t)\\, dt \\right]\n",
    "= \\int_0^{\\infty} f'(t)\\,(1 - F_P(t))\\, dt - \\int_{-\\infty}^0 f'(t)\\, F_P(t)\\, dt.\n",
    "$$\n",
    "\n",
    "Subtracting the same expression for $Q$, the constants cancel and both pieces\n",
    "merge into one integral over all of $\\mathbb{R}$:\n",
    "\n",
    "$$\n",
    "E_P[f] - E_Q[f] = \\int_{-\\infty}^{\\infty} f'(t)\\, \\big( F_Q(t) - F_P(t) \\big)\\, dt\n",
    "\\;\\leq\\; \\int_{-\\infty}^{\\infty} \\lvert F_P(t) - F_Q(t) \\rvert\\, dt,\n",
    "$$\n",
    "\n",
    "using $|f'| \\leq 1$. Equality is attained by integrating the choice\n",
    "$f'(t) = \\mathrm{sign}\\big(F_Q(t) - F_P(t)\\big)$, which defines a 1-Lipschitz\n",
    "function. By Kantorovich--Rubinstein the equation the supremum of\n",
    "the left side over 1-Lipschitz $f$ is $W_1$, proving\n",
    "the equation. $\\blacksquare$\n",
    "\n",
    "Geometrically, $W_1$ is the *area between the two CDFs* (the shaded region of\n",
    "the figure). Slicing that area horizontally instead\n",
    "of vertically gives the equivalent quantile form\n",
    "$W_1 = \\int_0^1 |F_P^{-1}(u) - F_Q^{-1}(u)|\\, du$, which for two equal-size\n",
    "empirical samples is just the mean absolute difference of their *sorted*\n",
    "values (Exercise 7). Let's verify the equation against the primal\n",
    "the equation, solved exactly as a linear program over transport\n",
    "plans (the plan $\\gamma$ is a matrix with row sums $p$ and column sums $q$,\n",
    "compare that section)."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-w1"
   },
   "source": [
    "atoms = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])\n",
    "p_w = np.array([0.30, 0.20, 0.25, 0.10, 0.10, 0.05])\n",
    "q_w = np.array([0.05, 0.10, 0.10, 0.25, 0.20, 0.30])\n",
    "\n",
    "F_p, F_q = np.cumsum(p_w), np.cumsum(q_w)\n",
    "w1_cdf = np.sum(np.abs(F_p - F_q)[:-1] * np.diff(atoms))\n",
    "\n",
    "C = np.abs(atoms[:, None] - atoms[None, :])      # cost matrix |x_i - y_j|\n",
    "k = len(atoms)\n",
    "A_eq = np.zeros((2 * k, k * k))\n",
    "for i in range(k):\n",
    "    A_eq[i, i * k:(i + 1) * k] = 1               # row sums    = p\n",
    "    A_eq[k + i, i::k] = 1                        # column sums = q\n",
    "res = optimize.linprog(C.ravel(), A_eq=A_eq,\n",
    "                       b_eq=np.concatenate([p_w, q_w]), method='highs')\n",
    "print(f'W1 via the CDF formula : {w1_cdf:.10f}')\n",
    "print(f'W1 via the primal LP   : {res.fun:.10f}')"
   ],
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "W1 via the CDF formula : 1.7000000000\nW1 via the primal LP   : 1.7000000000\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-8 fw=all -->\n",
    "\n",
    "Both routes give $W_1 = 1.7$ to ten decimal places: a one-line integral of\n",
    "CDFs replaces a $36$-variable linear program.\n",
    "\n",
    "Beyond one dimension no such formula exists, and the LP scales poorly: with\n",
    "$n$ atoms a side it has $n^2$ variables. The standard fix at scale is **entropic\n",
    "regularization** [@Cuturi.2013]: add $-\\varepsilon H(\\gamma)$ to the\n",
    "primal objective, penalizing low-entropy (overly deterministic) plans. The\n",
    "regularized problem is strictly convex and its solution has the form\n",
    "$\\gamma = \\mathrm{diag}(u)\\, K\\, \\mathrm{diag}(v)$ with\n",
    "$K = e^{-C/\\varepsilon}$, where the **Sinkhorn iterations** (alternately\n",
    "rescaling rows to match $p$ and columns to match $q$) converge fast, run on\n",
    "matrix-vector products (GPU-friendly), and are differentiable end to end. As\n",
    "$\\varepsilon \\to 0$ the entropic cost recovers the unregularized optimum\n",
    "[@Peyre.Cuturi.2019]."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-sinkhorn"
   },
   "source": [
    "def sinkhorn(p, q, C, eps, iters=5000):\n",
    "    \"\"\"Entropic OT by Sinkhorn iterations; returns plan and cost.\"\"\"\n",
    "    K = np.exp(-C / eps)\n",
    "    u = np.ones_like(p)\n",
    "    for _ in range(iters):\n",
    "        v = q / (K.T @ u)\n",
    "        u = p / (K @ v)\n",
    "    plan = u[:, None] * K * v[None, :]\n",
    "    return plan, np.sum(plan * C)\n",
    "\n",
    "fig, axes = d2l.plt.subplots(1, 3, figsize=(9, 2.8))\n",
    "axes[0].set_ylabel('source $i$')\n",
    "for ax, eps in zip(axes, (1.0, 0.1, 0.02)):\n",
    "    plan, cost = sinkhorn(p_w, q_w, C, eps)\n",
    "    print(f'epsilon = {eps:4.2f}: entropic cost = {cost:.4f}'\n",
    "          f'   (unregularized LP: {res.fun:.4f})')\n",
    "    ax.imshow(plan, cmap='Blues', origin='lower')\n",
    "    ax.set_title(f'$\\\\gamma$ at $\\\\varepsilon = {eps}$')\n",
    "    ax.set_xlabel('destination $j$')"
   ],
   "execution_count": 8,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "epsilon = 1.00: entropic cost = 1.7700   (unregularized LP: 1.7000)\nepsilon = 0.10: entropic cost = 1.7000   (unregularized LP: 1.7000)\nepsilon = 0.02: entropic cost = 1.7000   (unregularized LP: 1.7000)\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"543.640625pt\" height=\"207.577757pt\" viewBox=\"0 0 543.640625 207.577757\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n <metadata>\n  <rdf:RDF xmlns:dc=\"http://purl.org/dc/elements/1.1/\" xmlns:cc=\"http://creativecommons.org/ns#\" xmlns:rdf=\"http://www.w3.org/1999/02/22-rdf-syntax-ns#\">\n   <cc:Work>\n    <dc:type rdf:resource=\"http://purl.org/dc/dcmitype/StillImage\"/>\n    <dc:date>2026-07-19T13:41:29.261096</dc:date>\n    <dc:format>image/svg+xml</dc:format>\n    <dc:creator>\n     <cc:Agent>\n      <dc:title>Matplotlib v3.10.8, https://matplotlib.org/</dc:title>\n     </cc:Agent>\n    </dc:creator>\n   </cc:Work>\n  </rdf:RDF>\n </metadata>\n <defs>\n  <style type=\"text/css\">*{stroke-linejoin: round; stroke-linecap: butt}</style>\n </defs>\n <g id=\"figure_1\">\n  <g id=\"patch_1\">\n   <path d=\"M 0 207.577757 \nL 543.640625 207.577757 \nL 543.640625 0 \nL 0 0 \nz\n\" style=\"fill: #ffffff\"/>\n  </g>\n  <g id=\"axes_1\">\n   <g id=\"patch_2\">\n    <path d=\"M 34.240625 169.999632 \nL 181.946507 169.999632 \nL 181.946507 22.29375 \nL 34.240625 22.29375 \nz\n\" style=\"fill: #ffffff\"/>\n   </g>\n   <g clip-path=\"url(#p6122c0cef4)\">\n    <image xlink:href=\"data:image/png;base64,\niVBORw0KGgoAAAANSUhEUgAAAM4AAADOCAYAAAB2Hz3EAAADRUlEQVR4nO3dsarXZRzH8a95kGgIWwPBzhyF4Bh4ASGEk3gRDe4t0RzU0A3k4OguuFTbQQnEIfScxUkThDDDw/G4yP8C3svzC16vK/jAw5vv+Jy58euD09mAp89frZ6wc+3yp6snzMzMd7/8tnrCzMxc/fqL1RN2Hjx6tnrCzMx8sHoA/B8JBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBDs/fzN56s3zMzMD/cer56ws3/+o9UTZmbmy8v7qyfMzMylCx+vnrA5Lg4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFg78NzZ1dvmJmZm199tnrCzk9/HK2eMDMzB7dur57w3vXVA3b+PDhcPWFmXBxIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIzrw+Pj1dPWJm5nQbM2Zm5vhkG1u+v/vX6gkzM3Pl4ierJ+wcvvx39YSZcXEgEQ4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAgEA4EwoFAOBAIBwLhQCAcCIQDgXAg2MyPbFuyld/h/v7nzeoJMzPz4+9Hqyfs3H/yYvWEmXFxIBEOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBL4y3LCTt9t4mv/enKyesPPtnYerJ8yMiwOJcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAPBO/KoQsmLJhhDAAAAAElFTkSuQmCC\" id=\"image1d287c4f57\" transform=\"scale(1 -1) translate(0 -148.32)\" x=\"34.240625\" y=\"-21.679632\" width=\"148.32\" height=\"148.32\"/>\n   </g>\n   <g id=\"matplotlib.axis_1\">\n    <g id=\"xtick_1\">\n     <g id=\"line2d_1\">\n      <defs>\n       <path id=\"m419e685706\" d=\"M 0 0 \nL 0 3.5 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"46.549449\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_1\">\n      <!-- 0 -->\n      <g transform=\"translate(43.368199 184.59807) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-30\" d=\"M 2034 4250 \nQ 1547 4250 1301 3770 \nQ 1056 3291 1056 2328 \nQ 1056 1369 1301 889 \nQ 1547 409 2034 409 \nQ 2525 409 2770 889 \nQ 3016 1369 3016 2328 \nQ 3016 3291 2770 3770 \nQ 2525 4250 2034 4250 \nz\nM 2034 4750 \nQ 2819 4750 3233 4129 \nQ 3647 3509 3647 2328 \nQ 3647 1150 3233 529 \nQ 2819 -91 2034 -91 \nQ 1250 -91 836 529 \nQ 422 1150 422 2328 \nQ 422 3509 836 4129 \nQ 1250 4750 2034 4750 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_2\">\n     <g id=\"line2d_2\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"95.784743\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_2\">\n      <!-- 2 -->\n      <g transform=\"translate(92.603493 184.59807) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-32\" d=\"M 1228 531 \nL 3431 531 \nL 3431 0 \nL 469 0 \nL 469 531 \nQ 828 903 1448 1529 \nQ 2069 2156 2228 2338 \nQ 2531 2678 2651 2914 \nQ 2772 3150 2772 3378 \nQ 2772 3750 2511 3984 \nQ 2250 4219 1831 4219 \nQ 1534 4219 1204 4116 \nQ 875 4013 500 3803 \nL 500 4441 \nQ 881 4594 1212 4672 \nQ 1544 4750 1819 4750 \nQ 2544 4750 2975 4387 \nQ 3406 4025 3406 3419 \nQ 3406 3131 3298 2873 \nQ 3191 2616 2906 2266 \nQ 2828 2175 2409 1742 \nQ 1991 1309 1228 531 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-32\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_3\">\n     <g id=\"line2d_3\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"145.020037\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_3\">\n      <!-- 4 -->\n      <g transform=\"translate(141.838787 184.59807) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-34\" d=\"M 2419 4116 \nL 825 1625 \nL 2419 1625 \nL 2419 4116 \nz\nM 2253 4666 \nL 3047 4666 \nL 3047 1625 \nL 3713 1625 \nL 3713 1100 \nL 3047 1100 \nL 3047 0 \nL 2419 0 \nL 2419 1100 \nL 313 1100 \nL 313 1709 \nL 2253 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-34\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_4\">\n     <!-- destination $j$ -->\n     <g transform=\"translate(77.093566 198.277757) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-64\" d=\"M 2906 2969 \nL 2906 4863 \nL 3481 4863 \nL 3481 0 \nL 2906 0 \nL 2906 525 \nQ 2725 213 2448 61 \nQ 2172 -91 1784 -91 \nQ 1150 -91 751 415 \nQ 353 922 353 1747 \nQ 353 2572 751 3078 \nQ 1150 3584 1784 3584 \nQ 2172 3584 2448 3432 \nQ 2725 3281 2906 2969 \nz\nM 947 1747 \nQ 947 1113 1208 752 \nQ 1469 391 1925 391 \nQ 2381 391 2643 752 \nQ 2906 1113 2906 1747 \nQ 2906 2381 2643 2742 \nQ 2381 3103 1925 3103 \nQ 1469 3103 1208 2742 \nQ 947 2381 947 1747 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-65\" d=\"M 3597 1894 \nL 3597 1613 \nL 953 1613 \nQ 991 1019 1311 708 \nQ 1631 397 2203 397 \nQ 2534 397 2845 478 \nQ 3156 559 3463 722 \nL 3463 178 \nQ 3153 47 2828 -22 \nQ 2503 -91 2169 -91 \nQ 1331 -91 842 396 \nQ 353 884 353 1716 \nQ 353 2575 817 3079 \nQ 1281 3584 2069 3584 \nQ 2775 3584 3186 3129 \nQ 3597 2675 3597 1894 \nz\nM 3022 2063 \nQ 3016 2534 2758 2815 \nQ 2500 3097 2075 3097 \nQ 1594 3097 1305 2825 \nQ 1016 2553 972 2059 \nL 3022 2063 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-73\" d=\"M 2834 3397 \nL 2834 2853 \nQ 2591 2978 2328 3040 \nQ 2066 3103 1784 3103 \nQ 1356 3103 1142 2972 \nQ 928 2841 928 2578 \nQ 928 2378 1081 2264 \nQ 1234 2150 1697 2047 \nL 1894 2003 \nQ 2506 1872 2764 1633 \nQ 3022 1394 3022 966 \nQ 3022 478 2636 193 \nQ 2250 -91 1575 -91 \nQ 1294 -91 989 -36 \nQ 684 19 347 128 \nL 347 722 \nQ 666 556 975 473 \nQ 1284 391 1588 391 \nQ 1994 391 2212 530 \nQ 2431 669 2431 922 \nQ 2431 1156 2273 1281 \nQ 2116 1406 1581 1522 \nL 1381 1569 \nQ 847 1681 609 1914 \nQ 372 2147 372 2553 \nQ 372 3047 722 3315 \nQ 1072 3584 1716 3584 \nQ 2034 3584 2315 3537 \nQ 2597 3491 2834 3397 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-74\" d=\"M 1172 4494 \nL 1172 3500 \nL 2356 3500 \nL 2356 3053 \nL 1172 3053 \nL 1172 1153 \nQ 1172 725 1289 603 \nQ 1406 481 1766 481 \nL 2356 481 \nL 2356 0 \nL 1766 0 \nQ 1100 0 847 248 \nQ 594 497 594 1153 \nL 594 3053 \nL 172 3053 \nL 172 3500 \nL 594 3500 \nL 594 4494 \nL 1172 4494 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-69\" d=\"M 603 3500 \nL 1178 3500 \nL 1178 0 \nL 603 0 \nL 603 3500 \nz\nM 603 4863 \nL 1178 4863 \nL 1178 4134 \nL 603 4134 \nL 603 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6e\" d=\"M 3513 2113 \nL 3513 0 \nL 2938 0 \nL 2938 2094 \nQ 2938 2591 2744 2837 \nQ 2550 3084 2163 3084 \nQ 1697 3084 1428 2787 \nQ 1159 2491 1159 1978 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1366 3272 1645 3428 \nQ 1925 3584 2291 3584 \nQ 2894 3584 3203 3211 \nQ 3513 2838 3513 2113 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-61\" d=\"M 2194 1759 \nQ 1497 1759 1228 1600 \nQ 959 1441 959 1056 \nQ 959 750 1161 570 \nQ 1363 391 1709 391 \nQ 2188 391 2477 730 \nQ 2766 1069 2766 1631 \nL 2766 1759 \nL 2194 1759 \nz\nM 3341 1997 \nL 3341 0 \nL 2766 0 \nL 2766 531 \nQ 2569 213 2275 61 \nQ 1981 -91 1556 -91 \nQ 1019 -91 701 211 \nQ 384 513 384 1019 \nQ 384 1609 779 1909 \nQ 1175 2209 1959 2209 \nL 2766 2209 \nL 2766 2266 \nQ 2766 2663 2505 2880 \nQ 2244 3097 1772 3097 \nQ 1472 3097 1187 3025 \nQ 903 2953 641 2809 \nL 641 3341 \nQ 956 3463 1253 3523 \nQ 1550 3584 1831 3584 \nQ 2591 3584 2966 3190 \nQ 3341 2797 3341 1997 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6f\" d=\"M 1959 3097 \nQ 1497 3097 1228 2736 \nQ 959 2375 959 1747 \nQ 959 1119 1226 758 \nQ 1494 397 1959 397 \nQ 2419 397 2687 759 \nQ 2956 1122 2956 1747 \nQ 2956 2369 2687 2733 \nQ 2419 3097 1959 3097 \nz\nM 1959 3584 \nQ 2709 3584 3137 3096 \nQ 3566 2609 3566 1747 \nQ 3566 888 3137 398 \nQ 2709 -91 1959 -91 \nQ 1206 -91 779 398 \nQ 353 888 353 1747 \nQ 353 2609 779 3096 \nQ 1206 3584 1959 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-20\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-Oblique-6a\" d=\"M 928 3500 \nL 1503 3500 \nL 813 -63 \nL 809 -78 \nQ 694 -675 544 -897 \nQ 403 -1106 132 -1218 \nQ -138 -1331 -506 -1331 \nL -722 -1331 \nL -628 -844 \nL -481 -844 \nQ -144 -844 -1 -703 \nQ 141 -563 238 -63 \nL 928 3500 \nz\nM 1197 4863 \nL 1772 4863 \nL 1631 4134 \nL 1056 4134 \nL 1197 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(0 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(63.476562 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(125 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(177.099609 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(216.308594 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(244.091797 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(307.470703 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(368.75 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(407.958984 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(435.742188 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(496.923828 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(560.302734 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-Oblique-6a\" transform=\"translate(592.089844 0.015625)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"matplotlib.axis_2\">\n    <g id=\"ytick_1\">\n     <g id=\"line2d_4\">\n      <defs>\n       <path id=\"m0be083e38b\" d=\"M 0 0 \nL -3.5 0 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"34.240625\" y=\"157.690809\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_5\">\n      <!-- 0 -->\n      <g transform=\"translate(20.878125 161.490028) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_2\">\n     <g id=\"line2d_5\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"34.240625\" y=\"133.073162\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_6\">\n      <!-- 1 -->\n      <g transform=\"translate(20.878125 136.872381) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-31\" d=\"M 794 531 \nL 1825 531 \nL 1825 4091 \nL 703 3866 \nL 703 4441 \nL 1819 4666 \nL 2450 4666 \nL 2450 531 \nL 3481 531 \nL 3481 0 \nL 794 0 \nL 794 531 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-31\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_3\">\n     <g id=\"line2d_6\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"34.240625\" y=\"108.455515\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_7\">\n      <!-- 2 -->\n      <g transform=\"translate(20.878125 112.254733) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-32\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_4\">\n     <g id=\"line2d_7\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"34.240625\" y=\"83.837868\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_8\">\n      <!-- 3 -->\n      <g transform=\"translate(20.878125 87.637086) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-33\" d=\"M 2597 2516 \nQ 3050 2419 3304 2112 \nQ 3559 1806 3559 1356 \nQ 3559 666 3084 287 \nQ 2609 -91 1734 -91 \nQ 1441 -91 1130 -33 \nQ 819 25 488 141 \nL 488 750 \nQ 750 597 1062 519 \nQ 1375 441 1716 441 \nQ 2309 441 2620 675 \nQ 2931 909 2931 1356 \nQ 2931 1769 2642 2001 \nQ 2353 2234 1838 2234 \nL 1294 2234 \nL 1294 2753 \nL 1863 2753 \nQ 2328 2753 2575 2939 \nQ 2822 3125 2822 3475 \nQ 2822 3834 2567 4026 \nQ 2313 4219 1838 4219 \nQ 1578 4219 1281 4162 \nQ 984 4106 628 3988 \nL 628 4550 \nQ 988 4650 1302 4700 \nQ 1616 4750 1894 4750 \nQ 2613 4750 3031 4423 \nQ 3450 4097 3450 3541 \nQ 3450 3153 3228 2886 \nQ 3006 2619 2597 2516 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-33\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_5\">\n     <g id=\"line2d_8\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"34.240625\" y=\"59.220221\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_9\">\n      <!-- 4 -->\n      <g transform=\"translate(20.878125 63.019439) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-34\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_6\">\n     <g id=\"line2d_9\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"34.240625\" y=\"34.602574\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_10\">\n      <!-- 5 -->\n      <g transform=\"translate(20.878125 38.401792) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-35\" d=\"M 691 4666 \nL 3169 4666 \nL 3169 4134 \nL 1269 4134 \nL 1269 2991 \nQ 1406 3038 1543 3061 \nQ 1681 3084 1819 3084 \nQ 2600 3084 3056 2656 \nQ 3513 2228 3513 1497 \nQ 3513 744 3044 326 \nQ 2575 -91 1722 -91 \nQ 1428 -91 1123 -41 \nQ 819 9 494 109 \nL 494 744 \nQ 775 591 1075 516 \nQ 1375 441 1709 441 \nQ 2250 441 2565 725 \nQ 2881 1009 2881 1497 \nQ 2881 1984 2565 2268 \nQ 2250 2553 1709 2553 \nQ 1456 2553 1204 2497 \nQ 953 2441 691 2322 \nL 691 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-35\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_11\">\n     <!-- source $i$ -->\n     <g transform=\"translate(14.798437 115.846691) rotate(-90) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-75\" d=\"M 544 1381 \nL 544 3500 \nL 1119 3500 \nL 1119 1403 \nQ 1119 906 1312 657 \nQ 1506 409 1894 409 \nQ 2359 409 2629 706 \nQ 2900 1003 2900 1516 \nL 2900 3500 \nL 3475 3500 \nL 3475 0 \nL 2900 0 \nL 2900 538 \nQ 2691 219 2414 64 \nQ 2138 -91 1772 -91 \nQ 1169 -91 856 284 \nQ 544 659 544 1381 \nz\nM 1991 3584 \nL 1991 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-72\" d=\"M 2631 2963 \nQ 2534 3019 2420 3045 \nQ 2306 3072 2169 3072 \nQ 1681 3072 1420 2755 \nQ 1159 2438 1159 1844 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1341 3275 1631 3429 \nQ 1922 3584 2338 3584 \nQ 2397 3584 2469 3576 \nQ 2541 3569 2628 3553 \nL 2631 2963 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-63\" d=\"M 3122 3366 \nL 3122 2828 \nQ 2878 2963 2633 3030 \nQ 2388 3097 2138 3097 \nQ 1578 3097 1268 2742 \nQ 959 2388 959 1747 \nQ 959 1106 1268 751 \nQ 1578 397 2138 397 \nQ 2388 397 2633 464 \nQ 2878 531 3122 666 \nL 3122 134 \nQ 2881 22 2623 -34 \nQ 2366 -91 2075 -91 \nQ 1284 -91 818 406 \nQ 353 903 353 1747 \nQ 353 2603 823 3093 \nQ 1294 3584 2113 3584 \nQ 2378 3584 2631 3529 \nQ 2884 3475 3122 3366 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-Oblique-69\" d=\"M 1172 4863 \nL 1747 4863 \nL 1606 4134 \nL 1031 4134 \nL 1172 4863 \nz\nM 909 3500 \nL 1484 3500 \nL 800 0 \nL 225 0 \nL 909 3500 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(0 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(52.099609 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-75\" transform=\"translate(113.28125 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(176.660156 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-63\" transform=\"translate(217.773438 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(272.753906 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(334.277344 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-Oblique-69\" transform=\"translate(366.064453 0.015625)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"patch_3\">\n    <path d=\"M 34.240625 169.999632 \nL 34.240625 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_4\">\n    <path d=\"M 181.946507 169.999632 \nL 181.946507 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_5\">\n    <path d=\"M 34.240625 169.999632 \nL 181.946507 169.999632 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_6\">\n    <path d=\"M 34.240625 22.29375 \nL 181.946507 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"text_12\">\n    <!-- $\\gamma$ at $\\varepsilon = 1.0$ -->\n    <g transform=\"translate(74.493566 16.29375) scale(0.12 -0.12)\">\n     <defs>\n      <path id=\"DejaVuSans-Oblique-3b3\" d=\"M 1491 2950 \nL 1838 788 \nL 3456 3500 \nL 4066 3500 \nL 1972 0 \nL 1713 -1331 \nL 1138 -1331 \nL 1397 0 \nL 988 2613 \nQ 925 3006 628 3006 \nL 475 3006 \nL 569 3500 \nL 788 3500 \nQ 1403 3500 1491 2950 \nz\n\" transform=\"scale(0.015625)\"/>\n      <path id=\"DejaVuSans-Oblique-3b5\" d=\"M 1263 1888 \nQ 925 1963 778 2166 \nQ 669 2309 669 2503 \nQ 669 3034 1175 3344 \nQ 1575 3588 2188 3588 \nQ 2425 3588 2684 3550 \nQ 2944 3513 3234 3438 \nL 3134 2916 \nQ 2847 3006 2606 3047 \nQ 2359 3088 2138 3088 \nQ 1766 3088 1522 2944 \nQ 1228 2772 1228 2522 \nQ 1228 2356 1381 2241 \nQ 1563 2103 1925 2103 \nL 2409 2103 \nL 2319 1628 \nL 1856 1628 \nQ 1425 1628 1172 1469 \nQ 828 1253 828 916 \nQ 828 703 1013 563 \nQ 1244 388 1716 388 \nQ 2006 388 2284 444 \nQ 2563 503 2806 619 \nL 2700 84 \nQ 2403 -3 2131 -47 \nQ 1859 -91 1609 -91 \nQ 866 -91 516 194 \nQ 250 413 250 781 \nQ 250 1278 600 1584 \nQ 859 1813 1263 1888 \nz\n\" transform=\"scale(0.015625)\"/>\n      <path id=\"DejaVuSans-3d\" d=\"M 678 2906 \nL 4684 2906 \nL 4684 2381 \nL 678 2381 \nL 678 2906 \nz\nM 678 1631 \nL 4684 1631 \nL 4684 1100 \nL 678 1100 \nL 678 1631 \nz\n\" transform=\"scale(0.015625)\"/>\n      <path id=\"DejaVuSans-2e\" d=\"M 684 794 \nL 1344 794 \nL 1344 0 \nL 684 0 \nL 684 794 \nz\n\" transform=\"scale(0.015625)\"/>\n     </defs>\n     <use xlink:href=\"#DejaVuSans-Oblique-3b3\" transform=\"translate(0 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(59.179688 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(90.966797 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(152.246094 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(191.455078 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-Oblique-3b5\" transform=\"translate(223.242188 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-3d\" transform=\"translate(296.777344 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(400.048828 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(463.671875 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(495.458984 0.78125)\"/>\n    </g>\n   </g>\n  </g>\n  <g id=\"axes_2\">\n   <g id=\"patch_7\">\n    <path d=\"M 211.487684 169.999632 \nL 359.193566 169.999632 \nL 359.193566 22.29375 \nL 211.487684 22.29375 \nz\n\" style=\"fill: #ffffff\"/>\n   </g>\n   <g clip-path=\"url(#pb5a7a00e21)\">\n    <image xlink:href=\"data:image/png;base64,\niVBORw0KGgoAAAANSUhEUgAAAM4AAADOCAYAAAB2Hz3EAAADCklEQVR4nO3doWoXUBiH4W+yIEaDrFgM4jBabRbrmmjQO9gN7Cbsohi0LnoHJhenuwCFFUENOjDMov8kCG85E57nCn4cePni2dp79u58LoDP385WT9h4cvf66gkzM3Pw4mj1hJmZeXB/d/WEjfefvq6eMDMzl1YPgP+RcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAPB9uvHd1ZvmJmZ/cPj1RM2dq5cXj1hZmZu7e6snjAzMzevXYz3uEhcHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA8H26gF/PN27vXrCxsGbk9UTZmbm7fNXqyf89mj1gI2TD6erJ8yMiwOJcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAPB1o+f5+erR/B3F+VnuHs3rq6esHH6/Wz1hJlxcSARDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCDwIxv/tH94vHrCxsnHL6snzIyLA4lwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIPCVIf+Vhy+PVk+YGRcHEuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHgl/TxTDmYLeE5QAAAABJRU5ErkJggg==\" id=\"image3deb632da8\" transform=\"scale(1 -1) translate(0 -148.32)\" x=\"211.487684\" y=\"-21.679632\" width=\"148.32\" height=\"148.32\"/>\n   </g>\n   <g id=\"matplotlib.axis_3\">\n    <g id=\"xtick_4\">\n     <g id=\"line2d_10\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"223.796507\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_13\">\n      <!-- 0 -->\n      <g transform=\"translate(220.615257 184.59807) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_5\">\n     <g id=\"line2d_11\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"273.031801\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_14\">\n      <!-- 2 -->\n      <g transform=\"translate(269.850551 184.59807) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-32\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_6\">\n     <g id=\"line2d_12\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"322.267096\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_15\">\n      <!-- 4 -->\n      <g transform=\"translate(319.085846 184.59807) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-34\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_16\">\n     <!-- destination $j$ -->\n     <g transform=\"translate(254.340625 198.277757) scale(0.1 -0.1)\">\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(0 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(63.476562 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(125 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(177.099609 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(216.308594 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(244.091797 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(307.470703 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(368.75 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(407.958984 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(435.742188 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(496.923828 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(560.302734 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-Oblique-6a\" transform=\"translate(592.089844 0.015625)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"matplotlib.axis_4\">\n    <g id=\"ytick_7\">\n     <g id=\"line2d_13\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"211.487684\" y=\"157.690809\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_17\">\n      <!-- 0 -->\n      <g transform=\"translate(198.125184 161.490028) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_8\">\n     <g id=\"line2d_14\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"211.487684\" y=\"133.073162\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_18\">\n      <!-- 1 -->\n      <g transform=\"translate(198.125184 136.872381) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_9\">\n     <g id=\"line2d_15\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"211.487684\" y=\"108.455515\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_19\">\n      <!-- 2 -->\n      <g transform=\"translate(198.125184 112.254733) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-32\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_10\">\n     <g id=\"line2d_16\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"211.487684\" y=\"83.837868\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_20\">\n      <!-- 3 -->\n      <g transform=\"translate(198.125184 87.637086) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-33\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_11\">\n     <g id=\"line2d_17\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"211.487684\" y=\"59.220221\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_21\">\n      <!-- 4 -->\n      <g transform=\"translate(198.125184 63.019439) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-34\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_12\">\n     <g id=\"line2d_18\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"211.487684\" y=\"34.602574\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_22\">\n      <!-- 5 -->\n      <g transform=\"translate(198.125184 38.401792) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-35\"/>\n      </g>\n     </g>\n    </g>\n   </g>\n   <g id=\"patch_8\">\n    <path d=\"M 211.487684 169.999632 \nL 211.487684 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_9\">\n    <path d=\"M 359.193566 169.999632 \nL 359.193566 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_10\">\n    <path d=\"M 211.487684 169.999632 \nL 359.193566 169.999632 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_11\">\n    <path d=\"M 211.487684 22.29375 \nL 359.193566 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"text_23\">\n    <!-- $\\gamma$ at $\\varepsilon = 0.1$ -->\n    <g transform=\"translate(251.740625 16.29375) scale(0.12 -0.12)\">\n     <use xlink:href=\"#DejaVuSans-Oblique-3b3\" transform=\"translate(0 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(59.179688 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(90.966797 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(152.246094 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(191.455078 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-Oblique-3b5\" transform=\"translate(223.242188 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-3d\" transform=\"translate(296.777344 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(400.048828 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(463.671875 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(495.458984 0.78125)\"/>\n    </g>\n   </g>\n  </g>\n  <g id=\"axes_3\">\n   <g id=\"patch_12\">\n    <path d=\"M 388.734743 169.999632 \nL 536.440625 169.999632 \nL 536.440625 22.29375 \nL 388.734743 22.29375 \nz\n\" style=\"fill: #ffffff\"/>\n   </g>\n   <g clip-path=\"url(#p9272133272)\">\n    <image xlink:href=\"data:image/png;base64,\niVBORw0KGgoAAAANSUhEUgAAAM4AAADOCAYAAAB2Hz3EAAADCklEQVR4nO3doWoXUBiH4W+yIEaDrFgM4jBabRbrmmjQO9gN7Cbsohi0LnoHJhenuwCFFUENOjDMov8kCG85E57nCn4cePni2dp79u58LoDP385WT9h4cvf66gkzM3Pw4mj1hJmZeXB/d/WEjfefvq6eMDMzl1YPgP+RcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAPB9uvHd1ZvmJmZ/cPj1RM2dq5cXj1hZmZu7e6snjAzMzevXYz3uEhcHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA8H26gF/PN27vXrCxsGbk9UTZmbm7fNXqyf89mj1gI2TD6erJ8yMiwOJcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAPB1o+f5+erR/B3F+VnuHs3rq6esHH6/Wz1hJlxcSARDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCDwIxv/tH94vHrCxsnHL6snzIyLA4lwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIBAOBMKBQDgQCAcC4UAgHAiEA4FwIPCVIf+Vhy+PVk+YGRcHEuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHAuFAIBwIhAOBcCAQDgTCgUA4EAgHgl/TxTDmYLeE5QAAAABJRU5ErkJggg==\" id=\"image573d0c1a48\" transform=\"scale(1 -1) translate(0 -148.32)\" x=\"388.734743\" y=\"-21.679632\" width=\"148.32\" height=\"148.32\"/>\n   </g>\n   <g id=\"matplotlib.axis_5\">\n    <g id=\"xtick_7\">\n     <g id=\"line2d_19\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"401.043566\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_24\">\n      <!-- 0 -->\n      <g transform=\"translate(397.862316 184.59807) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_8\">\n     <g id=\"line2d_20\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"450.27886\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_25\">\n      <!-- 2 -->\n      <g transform=\"translate(447.09761 184.59807) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-32\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_9\">\n     <g id=\"line2d_21\">\n      <g>\n       <use xlink:href=\"#m419e685706\" x=\"499.514154\" y=\"169.999632\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_26\">\n      <!-- 4 -->\n      <g transform=\"translate(496.332904 184.59807) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-34\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_27\">\n     <!-- destination $j$ -->\n     <g transform=\"translate(431.587684 198.277757) scale(0.1 -0.1)\">\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(0 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(63.476562 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(125 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(177.099609 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(216.308594 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(244.091797 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(307.470703 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(368.75 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(407.958984 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(435.742188 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(496.923828 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(560.302734 0.015625)\"/>\n      <use xlink:href=\"#DejaVuSans-Oblique-6a\" transform=\"translate(592.089844 0.015625)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"matplotlib.axis_6\">\n    <g id=\"ytick_13\">\n     <g id=\"line2d_22\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"388.734743\" y=\"157.690809\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_28\">\n      <!-- 0 -->\n      <g transform=\"translate(375.372243 161.490028) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_14\">\n     <g id=\"line2d_23\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"388.734743\" y=\"133.073162\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_29\">\n      <!-- 1 -->\n      <g transform=\"translate(375.372243 136.872381) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_15\">\n     <g id=\"line2d_24\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"388.734743\" y=\"108.455515\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_30\">\n      <!-- 2 -->\n      <g transform=\"translate(375.372243 112.254733) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-32\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_16\">\n     <g id=\"line2d_25\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"388.734743\" y=\"83.837868\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_31\">\n      <!-- 3 -->\n      <g transform=\"translate(375.372243 87.637086) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-33\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_17\">\n     <g id=\"line2d_26\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"388.734743\" y=\"59.220221\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_32\">\n      <!-- 4 -->\n      <g transform=\"translate(375.372243 63.019439) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-34\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_18\">\n     <g id=\"line2d_27\">\n      <g>\n       <use xlink:href=\"#m0be083e38b\" x=\"388.734743\" y=\"34.602574\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_33\">\n      <!-- 5 -->\n      <g transform=\"translate(375.372243 38.401792) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-35\"/>\n      </g>\n     </g>\n    </g>\n   </g>\n   <g id=\"patch_13\">\n    <path d=\"M 388.734743 169.999632 \nL 388.734743 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_14\">\n    <path d=\"M 536.440625 169.999632 \nL 536.440625 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_15\">\n    <path d=\"M 388.734743 169.999632 \nL 536.440625 169.999632 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_16\">\n    <path d=\"M 388.734743 22.29375 \nL 536.440625 22.29375 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"text_34\">\n    <!-- $\\gamma$ at $\\varepsilon = 0.02$ -->\n    <g transform=\"translate(425.207684 16.29375) scale(0.12 -0.12)\">\n     <use xlink:href=\"#DejaVuSans-Oblique-3b3\" transform=\"translate(0 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(59.179688 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(90.966797 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(152.246094 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(191.455078 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-Oblique-3b5\" transform=\"translate(223.242188 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-3d\" transform=\"translate(296.777344 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(400.048828 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(463.671875 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(495.458984 0.78125)\"/>\n     <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(559.082031 0.78125)\"/>\n    </g>\n   </g>\n  </g>\n </g>\n <defs>\n  <clipPath id=\"p6122c0cef4\">\n   <rect x=\"34.240625\" y=\"22.29375\" width=\"147.705882\" height=\"147.705882\"/>\n  </clipPath>\n  <clipPath id=\"pb5a7a00e21\">\n   <rect x=\"211.487684\" y=\"22.29375\" width=\"147.705882\" height=\"147.705882\"/>\n  </clipPath>\n  <clipPath id=\"p9272133272\">\n   <rect x=\"388.734743\" y=\"22.29375\" width=\"147.705882\" height=\"147.705882\"/>\n  </clipPath>\n </defs>\n</svg>\n"
     },
     "metadata": {}
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-9 fw=all -->\n",
    "\n",
    "At $\\varepsilon = 1$ the blurred plan costs more ($1.77$ vs. $1.70$); by\n",
    "$\\varepsilon = 0.1$ the entropic cost matches the LP to four decimals. The\n",
    "heatmaps show why. At $\\varepsilon = 1$ entropy dominates and the plan is a\n",
    "haze: every source hedges its mass across many destinations, close to the\n",
    "independent coupling $p\\, q^\\top$, so parcels take detours and the cost runs\n",
    "high. Shrinking $\\varepsilon$ anneals the haze away: the plan sharpens\n",
    "toward a vertex of the transport polytope (an extreme point of the set of\n",
    "feasible plans), and at $\\varepsilon = 0.02$ only\n",
    "a thin monotone staircase of routes survives, the never-crossing assignment\n",
    "that one-dimensional optimal transport always produces and the LP finds\n",
    "exactly. With squared cost $\\|x - y\\|^2$ in\n",
    "the equation one obtains the Wasserstein-2 distance, whose\n",
    "dynamical (Benamou--Brenier) formulation as a minimum-kinetic-energy flow is\n",
    "the natural language for diffusion models and flow matching; we develop it\n",
    "where it is needed, in that section.\n",
    "\n",
    "## Scores: Fisher Divergence, Stein Discrepancy, and the Objective Map\n",
    "\n",
    "One last family compares distributions through *derivatives* of their\n",
    "log-densities. It looks exotic until you notice that it is the only family\n",
    "that never asks for a normalizing constant, which is exactly the term modern\n",
    "energy-based and diffusion models cannot compute.\n",
    "\n",
    "### The Score and the Fisher Divergence\n",
    "\n",
    "The **score** of a distribution $P$ with differentiable density $p$ is the\n",
    "gradient of its log-density *with respect to the data point*:\n",
    "\n",
    "$$\n",
    "s_P(\\mathbf{x}) = \\nabla_{\\mathbf{x}} \\log p(\\mathbf{x}).\n",
    "$$\n",
    "\n",
    "(A terminology hazard: in classical statistics, \"score\" means the gradient\n",
    "with respect to *parameters*, $\\nabla_{\\boldsymbol{\\theta}} \\log\n",
    "p_{\\boldsymbol{\\theta}}(\\mathbf{x})$, the object behind maximum likelihood\n",
    "and Fisher information in that section. Here the\n",
    "gradient is in $\\mathbf{x}$; this *data score* is the one that\n",
    "that section builds diffusion models\n",
    "from.) The score is a vector field on the sample space, pointing in the\n",
    "direction in which the density increases fastest: \"uphill, toward where the\n",
    "mass is\", as the figure shows. For a Gaussian\n",
    "$\\mathcal{N}(\\mu, \\sigma^2)$ it is\n",
    "\n",
    "$$\n",
    "s(x) = -\\frac{x - \\mu}{\\sigma^2},\n",
    "$$\n",
    "\n",
    "a spring pulling toward the mean with stiffness $1/\\sigma^2$; for a mixture\n",
    "$p = \\sum_k w_k\\, p_k$ the chain rule gives\n",
    "$s(x) = \\sum_k r_k(x)\\, s_k(x)$ with\n",
    "$r_k(x) = w_k p_k(x) / p(x)$: each component's spring, weighted by the\n",
    "posterior probability (\"responsibility\") that $x$ came from it.\n",
    "\n",
    "![The score field of a two-component Gaussian mixture. Arrows show the gradient of the log-density: near each mode the field is a spring pulling toward that mode's mean, and along the valley between modes the responsibilities hand the point off from one attractor to the other. Density contours are underlaid for reference.](../img/mdl-it-score-field.svg)\n",
    "\n",
    "The property that earns the score its central role is one line long. Suppose\n",
    "we can only write the density up to a constant,\n",
    "$p(\\mathbf{x}) = \\tilde{p}(\\mathbf{x}) / Z$ with\n",
    "$Z = \\int \\tilde{p}$ intractable, the universal situation for energy-based\n",
    "models. Then\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf{x}} \\log p(\\mathbf{x})\n",
    "= \\nabla_{\\mathbf{x}} \\log \\tilde{p}(\\mathbf{x})\n",
    "- \\nabla_{\\mathbf{x}} \\log Z\n",
    "= \\nabla_{\\mathbf{x}} \\log \\tilde{p}(\\mathbf{x}),\n",
    "$$\n",
    "\n",
    "since $Z$ does not depend on $\\mathbf{x}$. *The score never sees the\n",
    "normalizer.* A divergence built from scores therefore lets us fit\n",
    "unnormalized models, and the natural choice is the mean squared mismatch of\n",
    "the two vector fields under the data distribution, the **Fisher\n",
    "divergence**:\n",
    "\n",
    "$$\n",
    "D_{\\textrm{F}}(P\\|Q) = \\frac{1}{2}\\, E_{\\mathbf{x} \\sim P}\\!\\left[\n",
    "\\big\\| \\nabla_{\\mathbf{x}} \\log p(\\mathbf{x})\n",
    "- \\nabla_{\\mathbf{x}} \\log q(\\mathbf{x}) \\big\\|^2 \\right].\n",
    "$$\n",
    "\n",
    "It is non-negative. If $p$ and $q$ are positive on the same connected support,\n",
    "score equality $P$-almost everywhere forces $\\log p-\\log q$ to be constant\n",
    "there, and normalization then gives $P=Q$. Connected support of $P$ alone is\n",
    "not enough: $Q$ could agree with $P$ on that support while assigning additional\n",
    "mass elsewhere. For two\n",
    "equal-variance Gaussians the scores differ by the constant\n",
    "$(\\mu_2 - \\mu_1)/\\sigma^2$, so\n",
    "$D_{\\textrm{F}} = (\\mu_1 - \\mu_2)^2 / (2\\sigma^4)$; compare KL's\n",
    "$(\\mu_1 - \\mu_2)^2/(2\\sigma^2)$ from the equation; the general\n",
    "unequal-variance form is Exercise 9. Let's verify the mixture score formula,\n",
    "the normalizer-blindness, and the Gaussian closed form numerically."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-score"
   },
   "source": [
    "def score_p(x):\n",
    "    \"\"\"Score of the mixture P: responsibility-weighted component scores.\"\"\"\n",
    "    log_comps = np.stack([np.log(wi) - 0.5 * np.log(2 * np.pi * si ** 2)\n",
    "                          - (x - mi) ** 2 / (2 * si ** 2)\n",
    "                          for wi, mi, si in zip(w, mus, sigmas)])\n",
    "    r = np.exp(log_comps - log_p(x))             # responsibilities r_k(x)\n",
    "    comp_scores = np.stack([-(x - mi) / si ** 2 for mi, si in zip(mus, sigmas)])\n",
    "    return np.sum(r * comp_scores, axis=0)\n",
    "\n",
    "numeric = np.gradient(log_p(x), x)               # numerical d/dx log p\n",
    "gap = np.abs(score_p(x) - numeric)[10:-10].max()\n",
    "print(f'mixture score, max |analytic - numerical|: {gap:.2e}')\n",
    "shifted = np.gradient(np.log(2.7) + log_p(x), x)  # unnormalized: 2.7 * p\n",
    "print(f'score change from rescaling p by 2.7    : '\n",
    "      f'{np.abs(shifted - numeric).max():.2e}')\n",
    "\n",
    "mu1, s1, mu2, s2 = 0.0, 1.0, 1.0, 1.0\n",
    "xs = np.linspace(-10.0, 12.0, 4001)\n",
    "pdf1 = np.exp(-(xs - mu1) ** 2 / (2 * s1 ** 2)) / np.sqrt(2 * np.pi * s1 ** 2)\n",
    "sc1, sc2 = -(xs - mu1) / s1 ** 2, -(xs - mu2) / s2 ** 2\n",
    "fisher = 0.5 * integrate.trapezoid(pdf1 * (sc1 - sc2) ** 2, xs)\n",
    "print(f'Fisher divergence N(0,1)||N(1,1): quadrature {fisher:.6f}, '\n",
    "      f'closed form {(mu1 - mu2) ** 2 / (2 * s1 ** 4):.6f}')"
   ],
   "execution_count": 9,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "mixture score, max |analytic - numerical|: 3.52e-04\nscore change from rescaling p by 2.7    : 1.82e-12\nFisher divergence N(0,1)||N(1,1): quadrature 0.500000, closed form 0.500000\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-10 fw=all -->\n",
    "\n",
    "The analytic mixture score matches the numerical gradient of $\\log p$ to the\n",
    "grid's finite-difference accuracy ($\\approx 3.5 \\times 10^{-4}$), rescaling\n",
    "the density by $2.7$ changes the score by nothing but floating-point dust\n",
    "($\\approx 2 \\times 10^{-12}$), and the quadrature reproduces the Gaussian\n",
    "closed form $0.5$ to six decimals.\n",
    "\n",
    "One apparent obstacle remains: the equation is an\n",
    "expectation involving $\\nabla \\log p$ of the *data* distribution, which we do\n",
    "not know either. The resolution, Hyvärinen's **score matching** identity\n",
    "[@Hyvarinen.2005], is an integration by parts showing that, up to a\n",
    "constant independent of the model,\n",
    "\n",
    "$$\n",
    "D_{\\textrm{F}}(P\\|Q_{\\boldsymbol{\\theta}}) = E_{\\mathbf{x} \\sim P}\\!\\left[\n",
    "\\tfrac{1}{2} \\|\\mathbf{s}_{\\boldsymbol{\\theta}}(\\mathbf{x})\\|^2\n",
    "+ \\nabla_{\\mathbf{x}} \\cdot \\mathbf{s}_{\\boldsymbol{\\theta}}(\\mathbf{x})\n",
    "\\right] + \\textrm{const},\n",
    "$$\n",
    "\n",
    "an objective containing *only the model's score*\n",
    "$\\mathbf{s}_{\\boldsymbol{\\theta}}$, estimable from data samples alone. We\n",
    "state it here and prove it, together with its denoising variant (the actual\n",
    "training loss of diffusion models), in\n",
    "that section.\n",
    "\n",
    "### Stein's Identity\n",
    "\n",
    "The score can also *test* a sample. The starting point is a fact\n",
    "about expectations under a distribution whose score we know.\n",
    "\n",
    "**Proposition (Stein's identity).** [@Stein.1981] *Let $P$ have a differentiable density\n",
    "$p > 0$ on $\\mathbb{R}$ with score $s_P = (\\log p)'$, and let $f$ be\n",
    "differentiable with $f(x)\\, p(x) \\to 0$ as $x \\to \\pm\\infty$. Then*\n",
    "\n",
    "$$\n",
    "E_{x \\sim P}\\big[ f'(x) + f(x)\\, s_P(x) \\big] = 0.\n",
    "$$\n",
    "\n",
    "**Proof.** Since $p' = s_P\\, p$, integration by parts\n",
    "(that section) gives\n",
    "\n",
    "$$\n",
    "\\int f'(x)\\, p(x)\\, dx\n",
    "= \\big[ f(x)\\, p(x) \\big]_{-\\infty}^{\\infty} - \\int f(x)\\, p'(x)\\, dx\n",
    "= 0 - \\int f(x)\\, s_P(x)\\, p(x)\\, dx,\n",
    "$$\n",
    "\n",
    "and moving the right-hand side over is the equation.\n",
    "$\\blacksquare$\n",
    "\n",
    "For $P = \\mathcal{N}(0, 1)$ the score is $s_P(x) = -x$ and the identity reads\n",
    "$E[f'(X)] = E[X f(X)]$; with $f(x) = x$ it recovers $E[X^2] = 1$. The point\n",
    "is that the identity holds for *huge classes of $f$ simultaneously*, and only\n",
    "for $P$ itself: if a sample's averages violate\n",
    "the equation for some $f$, the sample is not from $P$.\n",
    "A Monte-Carlo check of the equation:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-stein"
   },
   "source": [
    "rng = np.random.default_rng(7)\n",
    "z = rng.standard_normal(1_000_000)\n",
    "# Stein operator for P = N(0,1):  (A_P f)(x) = f'(x) - x f(x)\n",
    "for name, f, fprime in [('x^3', lambda t: t ** 3, lambda t: 3 * t ** 2),\n",
    "                        ('sin x', np.sin, np.cos)]:\n",
    "    val = (fprime(z) - z * f(z)).mean()\n",
    "    print(f\"E[ f'(Z) - Z f(Z) ] for f(x) = {name}: {val:+.4f}\")"
   ],
   "execution_count": 10,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "E[ f'(Z) - Z f(Z) ] for f(x) = x^3: +0.0020\nE[ f'(Z) - Z f(Z) ] for f(x) = sin x: +0.0005\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-11 fw=all -->\n",
    "\n",
    "Both averages ($+0.0020$ and $+0.0005$ on a million samples) are zero to\n",
    "within Monte-Carlo error, for two quite different test functions, a glimpse\n",
    "of the infinite family of constraints the identity imposes.\n",
    "\n",
    "### The Kernel Stein Discrepancy\n",
    "\n",
    "To turn the identity into a divergence, run the IPM construction of\n",
    "that section on it: apply the **Stein operator**\n",
    "$(\\mathcal{A}_P f)(x) = f'(x) + f(x)\\, s_P(x)$ to every $f$ in the unit ball\n",
    "of an RKHS and take the largest violation,\n",
    "$\\sup_{\\|f\\|_{\\mathcal{H}} \\leq 1} E_{x \\sim Q}[(\\mathcal{A}_P f)(x)]$. By\n",
    "Stein's identity the supremum is zero when $Q = P$; when $Q \\neq P$, some\n",
    "test function witnesses the mismatch. Exactly as with MMD, the supremum over\n",
    "an RKHS ball has a closed form: it is the square root of an expected kernel.\n",
    "For a positive-definite kernel $k$ and score\n",
    "$s_P = \\nabla \\log p$, define the **Stein kernel** of $P$,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "u_P(\\mathbf{x}, \\mathbf{x}') = {}& s_P(\\mathbf{x})^\\top k(\\mathbf{x}, \\mathbf{x}')\\, s_P(\\mathbf{x}')\n",
    "+ s_P(\\mathbf{x})^\\top \\nabla_{\\mathbf{x}'} k(\\mathbf{x}, \\mathbf{x}') \\\\\n",
    "&+ \\nabla_{\\mathbf{x}} k(\\mathbf{x}, \\mathbf{x}')^\\top s_P(\\mathbf{x}')\n",
    "+ \\nabla_{\\mathbf{x}} \\cdot \\nabla_{\\mathbf{x}'} k(\\mathbf{x}, \\mathbf{x}'),\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where the last term is the sum of mixed partials\n",
    "$\\sum_i \\partial_{x_i} \\partial_{x'_i} k$; in one dimension it is just\n",
    "$\\partial_x \\partial_{x'} k$. The squared **kernel Stein discrepancy (KSD)**\n",
    "is the expected Stein kernel under two independent draws from $Q$:\n",
    "\n",
    "$$\n",
    "\\mathrm{KSD}^2(Q, P) = E_{\\mathbf{x}, \\mathbf{x}' \\sim Q}\\big[ u_P(\\mathbf{x}, \\mathbf{x}') \\big].\n",
    "$$\n",
    "\n",
    "Three properties make this the right tool for modern models. First,\n",
    "the equation involves $P$ only through its score, so by the\n",
    "normalizer-blindness of that section an unnormalized\n",
    "model works exactly as well as a normalized one, and no samples from $P$ are\n",
    "ever drawn. Second, for suitable kernels (the RBF kernel among them, under\n",
    "mild conditions on the score) $\\mathrm{KSD}^2(Q, P) = 0$ if and only if\n",
    "$Q = P$ [@Liu.Lee.Jordan.2016; @Chwialkowski.Strathmann.Gretton.2016].\n",
    "Third, the equation is a double expectation under $Q$ alone, so a\n",
    "sample $x_1, \\ldots, x_n \\sim Q$ gives the unbiased U-statistic estimator\n",
    "\n",
    "$$\n",
    "\\widehat{\\mathrm{KSD}}^2 = \\frac{1}{n(n-1)} \\sum_{i \\neq j} u_P(x_i, x_j),\n",
    "$$\n",
    "\n",
    "the same diagonal-excluding average as the MMD estimator. Let's compute it.\n",
    "For the 1-D RBF kernel $k(x, y) = e^{-(x-y)^2/(2\\ell^2)}$ the derivatives in\n",
    "the equation are closed-form:\n",
    "$\\partial_y k = \\frac{x-y}{\\ell^2}\\, k$,\n",
    "$\\partial_x k = -\\frac{x-y}{\\ell^2}\\, k$, and\n",
    "$\\partial_x \\partial_y k = \\big( \\frac{1}{\\ell^2} - \\frac{(x-y)^2}{\\ell^4} \\big)\\, k$.\n",
    "We test one sample from $Q = \\mathcal{N}(0, 1)$ against two models: the true\n",
    "$P = \\mathcal{N}(0, 1)$, whose score is $s(x) = -x$, and the wrong\n",
    "$P' = \\mathcal{N}(1, 1)$, whose score is $s(x) = -(x - 1)$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "divergences-distances-ksd"
   },
   "source": [
    "def ksd2_ustat(x, score, ell=1.0):\n",
    "    \"\"\"U-statistic KSD^2 with the RBF kernel exp(-(a-b)^2 / (2 ell^2)).\"\"\"\n",
    "    d = x[:, None] - x[None, :]\n",
    "    k = np.exp(-d ** 2 / (2 * ell ** 2))\n",
    "    dk_dx = -d / ell ** 2 * k\n",
    "    dk_dy = d / ell ** 2 * k\n",
    "    d2k = (1 / ell ** 2 - d ** 2 / ell ** 4) * k\n",
    "    s = score(x)\n",
    "    u = (s[:, None] * s[None, :] * k + s[:, None] * dk_dy\n",
    "         + dk_dx * s[None, :] + d2k)\n",
    "    n = len(x)\n",
    "    return (u.sum() - np.trace(u)) / (n * (n - 1))\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "x_q = rng.standard_normal(1000)          # the sample: Q = N(0, 1)\n",
    "for name, s in [('true model  N(0,1)', lambda t: -t),\n",
    "                ('wrong model N(1,1)', lambda t: -(t - 1))]:\n",
    "    print(f'KSD^2 vs the {name}: {ksd2_ustat(x_q, s):+.5f}')"
   ],
   "execution_count": 11,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "KSD^2 vs the true model  N(0,1): +0.00019\n"
     ]
    },
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "KSD^2 vs the wrong model N(1,1): +0.64425\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-12 fw=all -->\n",
    "\n",
    "Against the true model the U-statistic is $\\approx 0.0002$, consistent with\n",
    "zero (like the unbiased MMD estimator, it may even dip slightly negative);\n",
    "against the model whose mean is off by one it is $\\approx 0.64$, clearly\n",
    "positive and three orders of magnitude larger, although the *sample* never\n",
    "changed and $P'$ was never sampled at all: the discrepancy reads the mismatch\n",
    "straight off the score. This sample-versus-model comparison makes KSD the\n",
    "natural goodness-of-fit test for unnormalized models\n",
    "[@Liu.Lee.Jordan.2016], and the descent direction it induces on a\n",
    "particle set is Stein variational gradient descent (SVGD)\n",
    "[@Liu.Wang.2016], which transports particles toward $P$ using only its\n",
    "score.\n",
    "\n",
    "### The Divergence-to-Objective Map\n",
    "\n",
    "We can now deliver the payoff promised at the start: a single table that\n",
    "reads the modern generative objectives as choices of divergence. Each\n",
    "row names a training objective, the divergence its loss minimizes (often\n",
    "implicitly), where this section treats it, and the behavior the divergence\n",
    "makes inevitable.\n",
    "\n",
    "| Training objective | Divergence minimized | Treated in | Characteristic behavior |\n",
    "|:---|:---|:---|:---|\n",
    "| maximum likelihood: autoregressive models, normalizing flows | forward KL | that section | zero-avoiding and often mass-covering; may over-spread |\n",
    "| variational inference, VAE posterior (ELBO) | reverse KL | that section | mode-seeking; sharp but can drop modes; local optima |\n",
    "| original GAN | Jensen--Shannon | that section | sharp samples; no gradient on disjoint supports |\n",
    "| f-GAN | any chosen f-divergence | that section | inherits the chosen generator's behavior; critic-limited |\n",
    "| WGAN | Wasserstein-1 | that section | gradients survive disjoint supports; needs Lipschitz critic |\n",
    "| MMD-GAN, two-sample tests | MMD | that section | closed-form, adversary-free; kernel choice sets sensitivity |\n",
    "| score matching, diffusion models | Fisher divergence | that section | normalizer-free; trains on the score field |\n",
    "| SVGD, model criticism | kernel Stein discrepancy | that section | needs only the model's score; no model samples |\n",
    "\n",
    "*Maximum likelihood*\n",
    "(row 1) is forward KL by the NLL--cross-entropy equivalence of\n",
    "that section, so a language model trained on\n",
    "next-token prediction is mass-covering by construction: it would rather\n",
    "assign some probability to text it will never produce than risk assigning\n",
    "none to text that occurs. *The GAN family* (rows 3--5) is one design axis:\n",
    "the original GAN minimizes a bounded f-divergence and inherits its\n",
    "flat-on-disjoint-supports gradient pathology; WGAN swaps families to optimal\n",
    "transport, which supplies gradients everywhere but requires a constrained critic.\n",
    "*Diffusion models* (row 7) sidestep densities entirely: Fisher divergence\n",
    "compares score fields, the normalizer cancels, and the training loss becomes\n",
    "a regression onto $\\nabla_{\\mathbf{x}} \\log p$, which is where\n",
    "that section picks up the story.\n",
    "\n",
    "The thesis of the table, and of this section: **choose your divergence,\n",
    "inherit its failure modes.** Mode-dropping, blurriness, vanishing gradients,\n",
    "and training instability are theorems about the objective, visible in the\n",
    "geometry of $f$, the direction of the KL, the boundedness of JS, or the\n",
    "function class of the IPM, before a single parameter is trained.\n",
    "\n",
    "## Summary\n",
    "\n",
    "* A divergence demands only $D(P, Q) \\geq 0$ with equality iff $P = Q$;\n",
    "  metrics add symmetry and the triangle inequality. KL is a divergence but\n",
    "  not a metric; TV, Hellinger, and Wasserstein are metrics.\n",
    "* f-divergences $D_f(P\\|Q) = E_Q[f(p/q)]$ unify KL, reverse KL, $\\chi^2$,\n",
    "  Hellinger, TV, and Jensen--Shannon; non-negativity is Jensen's inequality,\n",
    "  and near $P = Q$ all smooth members agree up to the factor $f''(1)$. The\n",
    "  Rényi/$\\alpha$ family sweeps between the two KL directions with a single\n",
    "  knob and appears in differential-privacy accounting.\n",
    "* Fenchel duality turns any f-divergence into an adversarial game,\n",
    "  $D_f = \\sup_T \\{ E_P[T] - E_Q[f^*(T)] \\}$, estimable from samples; the\n",
    "  original GAN is the Jensen--Shannon case, and an undertrained critic\n",
    "  biases the estimate low.\n",
    "* The direction of KL is a modeling decision: forward KL (maximum\n",
    "  likelihood) is zero-avoiding and mass-covering; reverse KL (variational\n",
    "  inference) is zero-forcing and often mode-seeking on multimodal targets;\n",
    "  the number and location of local optima depend on the model family.\n",
    "* Total variation is the largest probability any event can disagree by, and\n",
    "  Pinsker's inequality $\\textrm{TV} \\leq \\sqrt{D_{\\textrm{KL}}/2}$ means\n",
    "  small KL certifies indistinguishability under every test.\n",
    "* IPMs replace events by a test-function class; the RKHS ball gives MMD,\n",
    "  with a closed-form unbiased estimator from samples alone.\n",
    "* Wasserstein distances measure mass transport and remain informative and\n",
    "  continuous when supports are disjoint (they need not be differentiable or\n",
    "  smooth there), equal an integral of CDF differences in 1-D, and\n",
    "  are computed at scale by entropic regularization and Sinkhorn iterations.\n",
    "* The score $\\nabla_{\\mathbf{x}} \\log p$ is blind to the normalizing\n",
    "  constant; the Fisher divergence compares score fields and underlies score\n",
    "  matching and diffusion; Stein's identity and the kernel Stein discrepancy\n",
    "  turn the score into goodness-of-fit tests.\n",
    "\n",
    "## Exercises\n",
    "\n",
    "1. Recover KL and reverse KL from the template the equation\n",
    "   with the generators $f(u) = u \\log u$ and $f(u) = -\\log u$. Then show\n",
    "   that for any constant $c$, the generators $f(u)$ and $f(u) + c\\,(u - 1)$\n",
    "   define the same divergence, and use this freedom to find a generator for\n",
    "   reverse KL that is non-negative everywhere.\n",
    "2. Derive the convex conjugate $f^*(t) = t + t^2/4$ of the $\\chi^2$\n",
    "   generator $f(u) = (u - 1)^2$ over $u \\in \\mathbb{R}$, then redo the\n",
    "   computation over the generator's true domain $u \\in (0, \\infty)$ and show\n",
    "   that the supremum is $t + t^2/4$ for $t \\geq -2$ but $-1$ for $t \\leq -2$.\n",
    "   Explain why plugging the (larger) $\\mathbb{R}$-conjugate into\n",
    "   the equation still yields a valid lower bound, merely a\n",
    "   weaker one for critics that dip below $-2$. Finally, write the explicit\n",
    "   f-GAN objective for this generator and verify that the critic\n",
    "   $T^\\star = f'(p/q) = 2(p/q - 1)$ attains the bound with equality.\n",
    "3. Two point masses at distance $d$: show that $D_{\\textrm{KL}}$, TV, and JS\n",
    "   are constant in $d$ (for $d \\neq 0$) while $W_1 = |d|$. What does this\n",
    "   imply about the gradient each objective supplies to a generator whose\n",
    "   samples are far from the data?\n",
    "4. Prove that total variation satisfies the triangle inequality, and decide\n",
    "   which of KL, reverse KL, squared Hellinger, and $W_1$ are metrics. (For\n",
    "   the Hellinger *distance* $H = \\sqrt{H^2}$, relate it to an $\\ell_2$ norm\n",
    "   of $\\sqrt{p} - \\sqrt{q}$.)\n",
    "5. Sharpness of Pinsker: for coins with biases $\\tfrac{1}{2}$ and\n",
    "   $\\tfrac{1}{2} + \\epsilon$, expand the KL divergence to second order in\n",
    "   $\\epsilon$ to show $D_{\\textrm{KL}} = 2\\epsilon^2 + O(\\epsilon^4)$ and\n",
    "   conclude that the ratio $\\textrm{TV}/\\sqrt{D_{\\textrm{KL}}/2}$ tends to\n",
    "   $1$, matching the experiment. Why is the bound loose for two\n",
    "   unit-variance Gaussians with distant means?\n",
    "6. Hellinger and total variation control each other. With the squared\n",
    "   Hellinger divergence of the gallery,\n",
    "   $H^2(P, Q) = \\sum_x \\big(\\sqrt{p(x)} - \\sqrt{q(x)}\\big)^2$, and\n",
    "   $H = \\sqrt{H^2}$ the Hellinger distance, prove the two-sided bound\n",
    "\n",
    "   $$\n",
    "   \\tfrac{1}{2}\\, H^2(P, Q) \\;\\leq\\; \\textrm{TV}(P, Q)\n",
    "   \\;\\leq\\; H(P, Q)\\, \\sqrt{ 1 - \\tfrac{1}{4} H^2(P, Q) }.\n",
    "   $$\n",
    "\n",
    "   (Hint: write $|p - q| = |\\sqrt{p} - \\sqrt{q}|\\,(\\sqrt{p} + \\sqrt{q})$.\n",
    "   The lower bound is $\\sqrt{p} + \\sqrt{q} \\geq |\\sqrt{p} - \\sqrt{q}|$; the\n",
    "   upper bound is Cauchy--Schwarz\n",
    "   (that section) together with\n",
    "   $\\sum_x (\\sqrt{p} + \\sqrt{q})^2 = 4 - H^2$.) Then verify both\n",
    "   inequalities numerically over $10{,}000$ random Dirichlet pairs on five\n",
    "   outcomes, in the style of the Pinsker experiment, and record how close\n",
    "   each side comes to equality. Unlike Pinsker, this sandwich is two-sided:\n",
    "   Hellinger and TV agree about which sequences of distributions converge,\n",
    "   whereas KL can be infinite at arbitrarily small TV.\n",
    "7. From the equation, derive the quantile form\n",
    "   $W_1 = \\int_0^1 |F_P^{-1}(u) - F_Q^{-1}(u)|\\, du$, and show that for two\n",
    "   empirical distributions on samples $x_1, \\ldots, x_n$ and\n",
    "   $y_1, \\ldots, y_n$ it equals $\\tfrac{1}{n}\\sum_i |x_{(i)} - y_{(i)}|$,\n",
    "   the mean absolute difference of sorted values. Verify numerically against\n",
    "   the linear program on a small example.\n",
    "8. For $P = \\mathcal{N}(0, 1)$ and $Q = \\mathcal{N}(\\delta, 1)$ with the RBF\n",
    "   kernel of bandwidth $\\ell = 1$, use the Gaussian integral\n",
    "   $E[e^{-Z^2/2}] = (1 + s^2)^{-1/2} e^{-m^2/(2(1+s^2))}$ for\n",
    "   $Z \\sim \\mathcal{N}(m, s^2)$ to derive\n",
    "   $\\textrm{MMD}^2 = \\tfrac{2}{\\sqrt{3}}\\big(1 - e^{-\\delta^2/6}\\big)$.\n",
    "   Evaluate it at $\\delta = 0.5$ and compare with the estimate from the code\n",
    "   cell.\n",
    "9. Compute the Fisher divergence between $\\mathcal{N}(\\mu_1, \\sigma_1^2)$\n",
    "   and $\\mathcal{N}(\\mu_2, \\sigma_2^2)$ in closed form (the score difference\n",
    "   is affine in $x$, so only Gaussian first and second moments are needed),\n",
    "   and check it reduces to $(\\mu_1 - \\mu_2)^2/(2\\sigma^4)$ for equal\n",
    "   variances. Then verify Stein's identity the equation\n",
    "   for $\\mathcal{N}(0, 1)$ with $f(x) = x$ by hand. What classical fact\n",
    "   about the standard Gaussian do you recover?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-divergences-distances-md-13 fw=pytorch -->\n",
    "\n",
    "[Discussions](https://d2l.discourse.group/t/divergences)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "d2l-pytorch",
   "display_name": "d2l (PyTorch)",
   "language": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
