{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "execute:\n",
    "  enabled: false\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-1 fw=all -->\n",
    "\n",
    "# Numerical Stability and Conditioning\n",
    "\n",
    "The math can be right and your loss can still go to `NaN`. Every result in this\n",
    "chapter so far was proved over the real numbers; your GPU computes over a\n",
    "finite, gappy imitation of them. This section explains the floating-point\n",
    "failure modes that bite real training runs (overflow, underflow, and\n",
    "catastrophic cancellation) and the handful of two-line fixes that keep\n",
    "softmax, cross-entropy, and ill-conditioned least squares alive:\n",
    "max-subtraction, log-space arithmetic, and ridge regularization. Two ideas\n",
    "organize everything. The first, due to numerical analysis and crystallized in\n",
    "@Higham.2002, is to split the blame for a wrong answer between the\n",
    "*algorithm* (did it solve a nearby problem?) and the *problem* (do nearby\n",
    "problems have wildly different answers?). The second is that the dividing line\n",
    "is a single number we have already met: the **condition number**\n",
    "$\\kappa = \\sigma_{\\max}/\\sigma_{\\min}$ of that section,\n",
    "the same $\\kappa$ that sets gradient descent's convergence rate in\n",
    "that section: one number, two consequences.\n",
    "The payoffs land downstream: the log-space trick rescues naive Bayes in\n",
    "that section from underflow, and the stable cross-entropy here\n",
    "is the same computation analyzed in that section.\n",
    "\n",
    "We proceed in four steps: what floating-point numbers are and where their\n",
    "cliffs lie; how to compute softmax, log-sum-exp, and cross-entropy safely;\n",
    "why subtracting nearly equal numbers destroys digits\n",
    "and how reformulation (not higher precision) repairs it; and finally\n",
    "conditioning: backward versus forward error, the Hilbert matrix,\n",
    "why normal equations square the pain, and how ridge regularization\n",
    "conditions the problem the way a preconditioner does. The standard references\n",
    "are @Goldberg.1991 for floating point and @Higham.2002 for\n",
    "everything else; @Goodfellow.Bengio.Courville.2016 (chapter 4) gives\n",
    "the deep-learning framing. Most code in this section is plain NumPy, since\n",
    "these phenomena belong to the arithmetic rather than to any library; the\n",
    "exception is the cross-entropy experiment, where library behavior genuinely\n",
    "differs."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-imports"
   },
   "source": [
    "%matplotlib inline\n",
    "from d2l import jax as d2l\n",
    "import jax\n",
    "from jax import numpy as jnp\n",
    "import optax\n",
    "import numpy as np"
   ],
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-2 fw=all -->\n",
    "\n",
    "## Floating-Point Arithmetic\n",
    "\n",
    "### A Number System with Gaps\n",
    "\n",
    "A floating-point number is scientific notation in base $2$ with a fixed budget\n",
    "of digits:\n",
    "\n",
    "$$\n",
    "x = (-1)^s \\cdot (1.m_1 m_2 \\ldots m_p)_2 \\cdot 2^{e},\n",
    "$$\n",
    "\n",
    "a sign bit $s$, a *mantissa* (significand) with $p$ stored bits, and an\n",
    "integer exponent $e$ from a fixed range. The exponent gives enormous *range*,\n",
    "the mantissa gives fixed *relative*\n",
    "precision, and between consecutive powers of two the representable values are\n",
    "evenly spaced, so the spacing *doubles* every time the magnitude does.\n",
    "the figure shows the resulting number line:\n",
    "representable values crowd near zero and thin out toward the overflow\n",
    "threshold,\n",
    "while the *relative* gap between neighbors stays essentially constant.\n",
    "\n",
    "![Floating-point numbers on the real line. Representable values are dense near zero and sparse far out: the absolute gap between neighbors doubles at every power of two while the relative gap stays near $\\varepsilon_{\\text{mach}}$. Each format ends at an overflow threshold (fp16 at $65504$, long before fp32 at about $3.4 \\times 10^{38}$) and in an underflow region below its smallest normal number.](../img/mdl-opt-fp-number-line.svg)\n",
    "\n",
    "That constant relative gap has a name. **Machine epsilon**\n",
    "$\\varepsilon_{\\text{mach}}$ is the distance from $1$ to the next representable\n",
    "number, $\\varepsilon_{\\text{mach}} = 2^{-p}$ for a $p$-bit mantissa, and it is\n",
    "the relative-error floor of the entire arithmetic: rounding any real $x$ to\n",
    "the nearest float $\\mathrm{fl}(x)$ obeys\n",
    "\n",
    "$$\n",
    "\\mathrm{fl}(x) = x\\,(1 + \\delta), \\qquad |\\delta| \\le u = \\tfrac12\\,\\varepsilon_{\\text{mach}},\n",
    "$$\n",
    "\n",
    "For a correctly rounded basic operation whose exact finite result lies in the\n",
    "normal range, IEEE arithmetic gives the analogous model: the computed\n",
    "$x \\oplus y$ equals $(x + y)(1 + \\delta)$ with $|\\delta| \\le u$\n",
    "([@IEEE.754.2019; @Goldberg.1991]). Overflow, division by zero, NaNs, and\n",
    "results in the subnormal/underflow region require separate absolute-error\n",
    "reasoning; those exceptions are central rather than incidental below. The quantity $u$ is the *unit\n",
    "roundoff*. Everything\n",
    "in this section is bookkeeping on these $(1+\\delta)$ factors: one of them is\n",
    "harmless; trouble starts when millions of them interact, or when a\n",
    "subtraction promotes one from relative to *catastrophic* (we get there in\n",
    "that section).\n",
    "\n",
    "Deep learning juggles three formats; the table below prints their parameters\n",
    "directly from the library:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-finfo"
   },
   "source": [
    "import numpy as onp\n",
    "header = f'{\"dtype\":>10} {\"eps\":>12} {\"smallest normal\":>17} {\"max\":>12}'\n",
    "print(header)\n",
    "for dt in [onp.float16, onp.float32]:\n",
    "    fi = onp.finfo(dt)\n",
    "    print(f'{onp.dtype(dt).name:>10} {fi.eps:12.3e} '\n",
    "          f'{fi.smallest_normal:17.3e} {fi.max:12.3e}')\n",
    "\n",
    "def to_bf16(x):\n",
    "    \"\"\"Round float32 to the nearest bfloat16 (round half to even).\"\"\"\n",
    "    bits = onp.atleast_1d(onp.asarray(x, onp.float32)).view(onp.uint32)\n",
    "    bits = (bits + 0x7FFF + ((bits >> 16) & 1)) & 0xFFFF0000\n",
    "    return bits.astype(onp.uint32).view(onp.float32)\n",
    "\n",
    "eps_bf16 = (to_bf16(1.0 + 2.0**-7) - 1.0).item()  # emulated: mxnet has no bf16\n",
    "print(f'{\"bfloat16\":>10} {eps_bf16:12.3e}   (exponent range = float32)')\n",
    "print('bfloat16 eps equals 2^-7:', eps_bf16 == 2.0**-7,\n",
    "      ' and 1 + 2^-8 rounds back to 1:', to_bf16(1.0 + 2.0**-8).item() == 1.0)"
   ],
   "execution_count": 2,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "     dtype          eps   smallest normal          max\n   float16    9.766e-04         6.104e-05    6.550e+04\n   float32    1.192e-07         1.175e-38    3.403e+38\n  bfloat16    7.812e-03   (exponent range = float32)\nbfloat16 eps equals 2^-7: True  and 1 + 2^-8 rounds back to 1: True\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-3 fw=all -->\n",
    "\n",
    "Read the three rows as three different bargains. **fp32** ($p = 23$ mantissa\n",
    "bits) has $\\varepsilon_{\\text{mach}} = 2^{-23} \\approx 1.19 \\times 10^{-7}$,\n",
    "about seven decimal digits, with range up to\n",
    "$3.4 \\times 10^{38}$. **fp16** ($p = 10$) keeps a respectable\n",
    "$\\varepsilon_{\\text{mach}} = 2^{-10} \\approx 9.8 \\times 10^{-4}$ but pays for\n",
    "it with a *tiny* exponent range: it overflows at $65504$ and its smallest\n",
    "normal number is about $6.1 \\times 10^{-5}$, so big activations overflow and\n",
    "small gradients underflow. **bfloat16** ($p = 7$) makes the opposite\n",
    "trade: it keeps fp32's full exponent range and sacrifices the mantissa,\n",
    "leaving $\\varepsilon_{\\text{mach}} = 2^{-7} = 0.0078$, between two and\n",
    "three decimal digits. The printout confirms the value that is easy to misquote:\n",
    "bfloat16's epsilon is $2^{-7}$, not $2^{-8}$; the eighth mantissa bit people\n",
    "sometimes count is the *implicit* leading $1$ in\n",
    "the equation, which contributes precision to no gap.\n",
    "\n",
    "Since 2022, hardware has pushed the same ladder one rung lower, to **fp8**,\n",
    "standardized in two flavors [@Micikevicius.Stosic.Burgess.ea.2022].\n",
    "**E4M3** ($p = 3$) keeps what resolution it can,\n",
    "$\\varepsilon_{\\text{mach}} = 2^{-3} = 0.125$ or roughly *one* decimal digit,\n",
    "inside a range that tops out at $448$; **E5M2** ($p = 2$) sacrifices another\n",
    "mantissa bit to buy fp16's full exponent range (maximum $57344$, smallest\n",
    "normal $6.1 \\times 10^{-5}$). The division of labor is visible in the numbers:\n",
    "E4M3 holds weights and activations, which need digits; E5M2 holds gradients,\n",
    "which need range. And at $\\varepsilon_{\\text{mach}} = 0.125$ there is no slack\n",
    "left in the format, so fp8 training must pair every tensor with a *scale\n",
    "factor* that recenters its values on the representable window. Every halving\n",
    "of the bit budget is paid for in digits, in range, or in bookkeeping."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-4 fw=all -->\n",
    "\n",
    "Two quick experiments make $\\varepsilon_{\\text{mach}}$ tangible: adding half\n",
    "an epsilon to $1$ vanishes without a trace, and the absolute gap between\n",
    "neighbors is a million times larger at $2^{20}$ than at $1$:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-spacing"
   },
   "source": [
    "eps = np.finfo(np.float32).eps\n",
    "one = np.float32(1.0)\n",
    "print('1 + eps   != 1 :', one + eps != one)\n",
    "print('1 + eps/2 == 1 :', one + np.float32(eps / 2) == one)\n",
    "print('gap between adjacent float32 values near 1    :',\n",
    "      np.spacing(np.float32(1.0)))\n",
    "print('gap between adjacent float32 values near 2^20 :',\n",
    "      np.spacing(np.float32(2.0**20)))\n",
    "for dt in [np.float16, np.float32]:\n",
    "    print(f'{np.dtype(dt).name}: exp(x) overflows for x >',\n",
    "          f'{np.log(np.finfo(dt).max):.2f}')"
   ],
   "execution_count": 3,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "1 + eps   != 1 : True\n1 + eps/2 == 1 : True\ngap between adjacent float32 values near 1    : 1.1920929e-07\ngap between adjacent float32 values near 2^20 : 0.125\nfloat16: exp(x) overflows for x > 11.09\nfloat32: exp(x) overflows for x > 88.72\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-5 fw=all -->\n",
    "\n",
    "### Overflow, Underflow, and Mixed Precision\n",
    "\n",
    "The last two printed lines locate the thresholds that matter most in practice.\n",
    "Because $e^x$ turns additive scale into multiplicative scale, the overflow\n",
    "threshold of each format translates into a modest *logit*:\n",
    "$e^x = \\infty$ in fp32 once $x > \\ln(3.4 \\times 10^{38}) \\approx 88.72$, and\n",
    "in fp16 once $x > \\ln(65504) \\approx 11.09$. Logits near $88.7$ are rare in\n",
    "healthy training, so fp32 softmax overflow is uncommon in practice; fp16's\n",
    "threshold of $11.09$ sits well inside the range of ordinary unnormalized\n",
    "scores, and it is the one that bites, which is why the mixed-precision recipe\n",
    "below exists. At the\n",
    "other end, $e^{-x}$ *underflows*: below the smallest normal number the format\n",
    "degrades gracefully through *subnormal* numbers with fewer and fewer\n",
    "significant bits, and then hits exactly $0$, at which point a subsequent\n",
    "$\\log$ returns $-\\infty$ and the backward pass turns to `NaN`.\n",
    "\n",
    "This is the arithmetic behind **mixed-precision training**\n",
    "[@Micikevicius.Narang.Alben.ea.2018]. In fp16, gradients routinely fall\n",
    "below $6 \\times 10^{-5}$ and vanish, so the loss is multiplied by a scale\n",
    "factor before the backward pass (and the gradients divided after) purely to\n",
    "shift them into representable territory: *loss scaling* is purely underflow\n",
    "management. bfloat16 makes that bookkeeping\n",
    "unnecessary: with fp32's exponent range nothing reasonable overflows or\n",
    "underflows, and the price, relative precision of only $2^{-7}$, is paid\n",
    "where deep learning is most tolerant, in the noise-dominated mantissa of each\n",
    "weight update. The reason a master copy of the weights is kept in fp32 is\n",
    "the equation again: a weight update of relative size\n",
    "below $\\varepsilon_{\\text{mach}}/2$ rounds to *no update at all*, and at\n",
    "$\\varepsilon_{\\text{mach}} = 2^{-7}$ that threshold is hit by perfectly\n",
    "healthy learning rates. A third escape, increasingly common in fp8 and\n",
    "integer training, is **stochastic rounding**\n",
    "[@Gupta.Agrawal.Gopalakrishnan.ea.2015]: round up or down at random with\n",
    "probabilities proportional to proximity, so the *expected* stored value is\n",
    "exact and updates too small to survive round-to-nearest still make progress\n",
    "on average.\n",
    "\n",
    "Both fp16 failure modes, and both rescues, fit in one cell. A gradient\n",
    "whose true value is $10^{-8}$ underflows to zero during an fp16 backward\n",
    "pass, and multiplying the *loss* by $2^{14}$ before differentiating shifts\n",
    "the whole gradient chain into representable territory; a perfectly healthy\n",
    "update of relative size $10^{-4}$ is lost to round-to-nearest in\n",
    "fp16, and an fp32 master copy of the weights accepts it without complaint:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-6 fw=all -->\n",
    "\n",
    "## Making Softmax and Cross-Entropy Safe\n",
    "\n",
    "### Softmax Overflows and the Shift That Fixes It\n",
    "\n",
    "The most common stability bug in machine learning is one line long. The\n",
    "softmax\n",
    "\n",
    "$$\n",
    "\\mathrm{softmax}(\\mathbf{z})_i = \\frac{e^{z_i}}{\\sum_{j=1}^n e^{z_j}}\n",
    "$$\n",
    "\n",
    "exponentiates its logits, and the table above says exactly where that goes\n",
    "wrong: fp32 overflows the moment any logit exceeds $88.72$, fp16 already at\n",
    "$11.09$: the numerator becomes `inf`, the ratio becomes `inf/inf = NaN`,\n",
    "and the model dies even though the *probabilities* it was computing are\n",
    "perfectly tame numbers in $[0, 1]$. The failure is entirely an artifact of\n",
    "the route. You have met the repair before:\n",
    "that section derived it when the fused\n",
    "cross-entropy loss was introduced, and we will not re-derive it here.\n",
    "Factoring the positive constant $e^{-c}$ out of numerator and denominator\n",
    "shows that softmax is *shift-invariant*,\n",
    "\n",
    "$$\n",
    "\\mathrm{softmax}(\\mathbf{z} - c\\mathbf{1}) = \\mathrm{softmax}(\\mathbf{z})\n",
    "\\qquad \\textrm{for every } c \\in \\mathbb{R},\n",
    "$$\n",
    "\n",
    "and the same factoring, applied under the logarithm of the denominator,\n",
    "rewrites the **log-sum-exp** $\\mathrm{lse}(\\mathbf{z}) = \\log \\sum_j e^{z_j}$\n",
    "exactly, for any shift $c$:\n",
    "\n",
    "$$\n",
    "\\mathrm{lse}(\\mathbf{z}) = c + \\log \\sum_{j=1}^n e^{z_j - c} .\n",
    "$$\n",
    "\n",
    "that section also gave the reason the\n",
    "shifted route is safe for finite logits: with $c = \\max_i z_i$ every\n",
    "exponent is at most $0$, so each $e^{z_i-c}$ lies in $(0,1]$ and exponential\n",
    "overflow is avoided; the largest term equals $1$, so the denominator, the sum\n",
    "in the equation, lies in $[1,n]$ and cannot underflow to\n",
    "zero. For any practically storable vector this sum is also representable.\n",
    "Non-finite inputs need separate handling (for example, subtracting an\n",
    "all-$-\\infty$ maximum is undefined). What the floating-point model of this\n",
    "section adds is the quantification: the overflow thresholds ($88.72$ in fp32,\n",
    "$11.09$ in fp16) say exactly when the naive route dies, while the shifted\n",
    "route sends only the nonpositive differences $z_i-c$ to the exponential. The cell below watches the naive route produce `NaN`\n",
    "on logits that the shifted route handles, and checks whether, where both\n",
    "routes work, the two agree:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-stable-softmax"
   },
   "source": [
    "def softmax_naive(z):\n",
    "    e = np.exp(z)\n",
    "    return e / e.sum()\n",
    "\n",
    "def softmax_stable(z):\n",
    "    e = np.exp(z - z.max())        # shift by the max: largest exponent is 0\n",
    "    return e / e.sum()\n",
    "\n",
    "z = np.array([1.0, 2.0, 3.0], dtype=np.float32)\n",
    "with np.errstate(over='ignore', invalid='ignore'):\n",
    "    print('naive,  logits z      :', softmax_naive(z))\n",
    "    print('naive,  logits z + 100:', softmax_naive(z + 100.0))\n",
    "print('stable, logits z + 100:', softmax_stable(z + 100.0))\n",
    "print('naive and stable agree where both work:',\n",
    "      bool((softmax_naive(z) == softmax_stable(z + 100.0)).all()))"
   ],
   "execution_count": 4,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "naive,  logits z      : [0.09003058 0.24472846 0.66524094]\nnaive,  logits z + 100: [nan nan nan]\nstable, logits z + 100: [0.09003057 0.24472846 0.66524094]\nnaive and stable agree where both work: False\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-7 fw=all -->\n",
    "\n",
    "The shifted logits give back $(0.090, 0.245, 0.665)$, the same probabilities\n",
    "as the small logits, while the naive route returns three `NaN`s. Whether the\n",
    "two routes agree *bit for bit* is a subtler question:\n",
    "the equation promises equality of real numbers, and two\n",
    "different exp/sum/divide routes may round differently along the way. Here\n",
    "they agree to the last printed digit; the final check prints `True` on most\n",
    "builds, but under one NumPy build in our environments it prints `False`, the\n",
    "naive route's first entry differing by a single ulp ($0.09003058$ versus\n",
    "$0.09003057$). That one-bit discrepancy is the section's lesson in\n",
    "miniature: the identity constrains the real values, not the rounded ones.\n",
    "The library's `softmax` does this max-subtraction internally; the trap is\n",
    "re-implementing it yourself, which is why the practical rule is to use a\n",
    "library's stable softmax or log-sum-exp rather than exponentiating raw logits.\n",
    "\n",
    "### The Log-Sum-Exp Sandwich\n",
    "\n",
    "The log-sum-exp recurs far beyond softmax: it is the normalizer of every\n",
    "exponential-family model (that section) and, as we prove\n",
    "in a moment, the backbone of cross-entropy. Beyond making it safe, the shift\n",
    "in the equation makes it *legible*, pinning lse between two\n",
    "bounds tight enough to reason with:\n",
    "\n",
    "**Proposition (log-sum-exp sandwich).** *For every\n",
    "$\\mathbf{z} \\in \\mathbb{R}^n$,*\n",
    "\n",
    "$$\n",
    "\\max_j z_j \\;\\le\\; \\mathrm{lse}(\\mathbf{z}) \\;\\le\\; \\max_j z_j + \\log n .\n",
    "$$\n",
    "\n",
    "**Proof.** Put $c = \\max_j z_j$ in the equation. Every\n",
    "term $e^{z_j - c} \\le 1$ and the maximizing term equals $1$, so the sum lies\n",
    "in $[1, n]$ and its logarithm in $[0, \\log n]$; adding $c$ gives both\n",
    "inequalities. $\\blacksquare$\n",
    "\n",
    "The sandwich says lse is a *soft maximum*, within $\\log n$ of the true max,\n",
    "which is also the intuition for why it is convex\n",
    "(that section; its gradient is exactly the softmax, a fact you\n",
    "will prove in the exercises). Numerically, the identity gives us log-space\n",
    "arithmetic for free: the **log-softmax** is\n",
    "\n",
    "$$\n",
    "\\log \\mathrm{softmax}(\\mathbf{z})_i = z_i - \\mathrm{lse}(\\mathbf{z}),\n",
    "$$\n",
    "\n",
    "a subtraction that does not materialize the probability itself. It therefore\n",
    "avoids the route in which a tiny probability first underflows to $0$ and its\n",
    "log then becomes $-\\infty$. Logits around $1000$ would overflow even float64;\n",
    "in log space they are effortless:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-logsumexp"
   },
   "source": [
    "def log_sum_exp(z):\n",
    "    c = z.max()\n",
    "    return c + np.log(np.exp(z - c).sum())\n",
    "\n",
    "z = np.array([1000.0, 1001.0, 1002.0], dtype=np.float32)\n",
    "with np.errstate(over='ignore'):\n",
    "    print('naive  log(sum(exp(z))) :', np.log(np.exp(z).sum()))\n",
    "print('stable log_sum_exp(z)   :', log_sum_exp(z))\n",
    "log_p = z - log_sum_exp(z)             # log softmax, eq. above\n",
    "print('log softmax             :', log_p)\n",
    "print('probabilities sum to 1  :', f'{np.exp(log_p).sum():.6f}')"
   ],
   "execution_count": 5,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "naive  log(sum(exp(z))) : inf\nstable log_sum_exp(z)   : 1002.4076\nlog softmax             : [-2.4075928  -1.4075928  -0.40759277]\nprobabilities sum to 1  : 1.000013\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-8 fw=all -->\n",
    "\n",
    "The naive route says `inf`; the stable route reports\n",
    "$\\mathrm{lse} = 1002.4076$, snugly inside the sandwich\n",
    "$[1002, 1002 + \\log 3]$, and the log-probabilities\n",
    "$(-2.408, -1.408, -0.408)$ exponentiate back to a distribution summing to\n",
    "$1.000013$: equal to $1$ up to the float32 spacing at magnitude $1000$,\n",
    "which is all one can ask of subtractions performed there. This identity is precisely how that section\n",
    "multiplies hundreds of per-pixel probabilities without underflowing to zero:\n",
    "sums of logs, never products of probabilities.\n",
    "\n",
    "### Pass Logits, Not Probabilities\n",
    "\n",
    "Cross-entropy is where the stakes are highest, because it is the loss\n",
    "*gradient* that dies. For a label $y$,\n",
    "the equation gives\n",
    "\n",
    "$$\n",
    "-\\log \\mathrm{softmax}(\\mathbf{z})_y = \\mathrm{lse}(\\mathbf{z}) - z_y,\n",
    "$$\n",
    "\n",
    "computable *directly from the logits* with one stable lse and one\n",
    "subtraction. This is what the library's \"from logits\" loss does (the\n",
    "fused implementation of that section),\n",
    "and it is why those APIs exist. The alternative, computing probabilities\n",
    "first and then taking the log, forces the loss through the representable\n",
    "range of probabilities: a true-class probability below about $10^{-45}$ (the\n",
    "smallest fp32 subnormal is $\\approx 1.4 \\times 10^{-45}$) underflows fp32 to\n",
    "exactly $0$, and the loss becomes $\\infty$. The cell below pits the\n",
    "two routes against each other on a two-class problem where the label is the\n",
    "*unlikely* class, with logit gap $t$, so the true loss is\n",
    "$\\log(1 + e^{t}) \\approx t$. This is the one computation in this section\n",
    "where library behavior genuinely differs; the from-probabilities route fails\n",
    "in one of three ways, depending on the library: subnormal noise followed by\n",
    "`inf`, `inf` outright, or a silent clip. Run the cell and read on:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-cross-entropy"
   },
   "source": [
    "print('gap    CE from logits    CE via probabilities')\n",
    "for t in [20.0, 60.0, 103.0, 104.0]:\n",
    "    logits = jnp.array([[0.0, t]])\n",
    "    y = jnp.array([0])                           # label = the unlikely class\n",
    "    from_logits = optax.softmax_cross_entropy_with_integer_labels(logits, y)\n",
    "    from_probs = -jnp.log(jax.nn.softmax(logits)[0, 0])\n",
    "    print(f'{t:5.0f}  {float(from_logits[0]):15.4f}  {float(from_probs):15.4f}')"
   ],
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "gap    CE from logits    CE via probabilities\n"
     ]
    },
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "   20          20.0000          20.0000\n   60          60.0000          60.0000\n  103         103.0000              inf\n  104         104.0000              inf\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-9 fw=jax -->\n",
    "\n",
    "The from-logits column reads $20$, $60$, $103$, $104$: exact at every gap.\n",
    "The from-probabilities column matches at gaps $20$ and $60$; at gap $103$,\n",
    "where $e^{-t}$ would survive only as a *subnormal* number, XLA does not\n",
    "linger in the subnormal range and the probability underflows to exactly $0$,\n",
    "and at gap $104$ the underflow is unconditional, so the loss reads `inf` at\n",
    "both gaps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-10 fw=all -->\n",
    "\n",
    "The lesson generalizes far beyond this toy: losses, likelihoods, and\n",
    "posteriors should live in log space from birth, and the conversion to\n",
    "probabilities, if it ever happens, should be the *last* step, for human\n",
    "eyes only. that section analyzes what cross-entropy\n",
    "*means*; the equation is how it is *computed*.\n",
    "\n",
    "## Catastrophic Cancellation\n",
    "\n",
    "### Subtraction Annihilates Digits\n",
    "\n",
    "Overflow announces itself with `inf`; cancellation gives no signal at all.\n",
    "Subtracting two nearly equal numbers is *exact* (no new rounding error is\n",
    "committed) but it strips away the leading digits on which both numbers\n",
    "agreed, leaving only their trailing digits, which is exactly where each\n",
    "number's *previous* rounding errors live. If $a$ and $b$ are correct to\n",
    "relative error $u$, their difference can carry relative error as large as\n",
    "\n",
    "$$\n",
    "\\frac{|a| + |b|}{|a - b|}\\; u,\n",
    "$$\n",
    "\n",
    "an amplification factor that blows up precisely when $a \\approx b$. The\n",
    "phenomenon is called **catastrophic cancellation**, and a two-line experiment\n",
    "shows both the disease and a cure. In float32, $1 + 10^{-8}$ rounds to\n",
    "exactly $1$ (the increment is below $\\varepsilon_{\\text{mach}}/2$), so the\n",
    "textbook expression $\\log(1 + x)$ returns $0$, a $100\\%$ relative error,\n",
    "while the library function `log1p`, which evaluates the *reformulated*\n",
    "series around $0$, is exact. Likewise two floats agreeing to seven digits\n",
    "leave a difference made of pure noise:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-log1p"
   },
   "source": [
    "x = np.float32(1e-8)\n",
    "print('float32 rounds 1 + x to     :', np.float32(1.0) + x)\n",
    "print('log(1 + x) =', np.log(np.float32(1.0) + x),\n",
    "      '   log1p(x) =', np.log1p(x))\n",
    "a, b = np.float32(1.0002344), np.float32(1.0002341)\n",
    "print('a - b in float32            :', a - b, '  (true value 3.0e-07)')\n",
    "print('amplification (|a|+|b|)/|a-b| ~', f'{(a + b) / abs(a - b):.1e}')"
   ],
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "float32 rounds 1 + x to     : 1.0\nlog(1 + x) = 0.0    log1p(x) = 1e-08\na - b in float32            : 2.3841858e-07   (true value 3.0e-07)\namplification (|a|+|b|)/|a-b| ~ 8.4e+06\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-11 fw=all -->\n",
    "\n",
    "The computed difference $2.384 \\times 10^{-7}$ misses the true\n",
    "$3.0 \\times 10^{-7}$ by twenty percent: with an amplification factor near\n",
    "$10^{7}$, float32's seven digits are gone in one subtraction. The catalogue\n",
    "of standard victims is short: $\\log(1+x)$ and $e^x - 1$\n",
    "near $0$ (use `log1p` and `expm1`), $1 - \\cos x$ near $0$ (use\n",
    "$2\\sin^2(x/2)$), the quadratic formula near a double root, and finite\n",
    "differences with too small a step, which is exactly the trade-off we\n",
    "quantified in that section. In every case the\n",
    "cure is the same: *reformulate so that the subtraction happens analytically*,\n",
    "where it costs nothing, rather than numerically, where it costs everything.\n",
    "Higher precision merely postpones the failure; reformulation removes it\n",
    "[@Higham.2002].\n",
    "\n",
    "### Case Study: Variance in One Pass\n",
    "\n",
    "The classic cancellation bug in data science is the \"computational formula\"\n",
    "for variance,\n",
    "\n",
    "$$\n",
    "\\mathrm{Var}(x) = \\mathbb{E}[x^2] - \\mathbb{E}[x]^2,\n",
    "$$\n",
    "\n",
    "beloved because it needs one pass over the data. As algebra it is flawless; as\n",
    "arithmetic it is a trap. For data with mean $\\mu$ and standard deviation\n",
    "$\\sigma \\ll |\\mu|$, both terms are about $\\mu^2$ while their difference is\n",
    "$\\sigma^2$, so the equation predicts an error\n",
    "amplification of about $\\mu^2/\\sigma^2$, and with $\\mu = 10^9$ and\n",
    "$\\sigma = 1$ that is $10^{18}$: more than every digit float64 has. The naive\n",
    "formula can even return a *negative* variance.\n",
    "\n",
    "The repair is a reformulation due to\n",
    "@Welford.1962 that keeps a running mean $m_k$ and a running sum of\n",
    "*centered* squares $M_k = \\sum_{i \\le k} (x_i - m_k)^2$, so no large numbers\n",
    "are ever subtracted:\n",
    "\n",
    "$$\n",
    "m_k = m_{k-1} + \\frac{x_k - m_{k-1}}{k},\n",
    "\\qquad\n",
    "M_k = M_{k-1} + (x_k - m_{k-1})(x_k - m_k).\n",
    "$$\n",
    "\n",
    "Note the two different factors in the $M_k$ update: the deviation from the\n",
    "*old* mean times the deviation from the *new* mean. That asymmetry is\n",
    "exactly what makes the recursion exact:\n",
    "\n",
    "**Proposition (Welford's recursion is exact).** *With $m_0 = M_0 = 0$, the\n",
    "recursions the equation satisfy, for every $k \\ge 1$ and in\n",
    "exact arithmetic,*\n",
    "\n",
    "$$\n",
    "m_k = \\frac{1}{k} \\sum_{i=1}^k x_i,\n",
    "\\qquad\n",
    "M_k = \\sum_{i=1}^k (x_i - m_k)^2 .\n",
    "$$\n",
    "\n",
    "**Proof.** The mean claim is the identity $k\\, m_k = (k-1)\\, m_{k-1} + x_k$,\n",
    "immediate from the first recursion. For the second claim, induct on $k$ and\n",
    "write $\\delta = x_k - m_{k-1}$, so that $m_k - m_{k-1} = \\delta/k$ and\n",
    "$x_k - m_k = \\delta\\,(k-1)/k$. Splitting the new sum of squares at its last\n",
    "term and re-centering the first $k-1$ terms around $m_{k-1}$,\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^{k} (x_i - m_k)^2\n",
    "= \\sum_{i=1}^{k-1} \\left( (x_i - m_{k-1}) - \\tfrac{\\delta}{k} \\right)^2 + (x_k - m_k)^2\n",
    "= M_{k-1} + (k-1)\\tfrac{\\delta^2}{k^2} + \\tfrac{(k-1)^2}{k^2}\\delta^2,\n",
    "$$\n",
    "\n",
    "where the cross term vanished because $\\sum_{i \\le k-1} (x_i - m_{k-1}) = 0$\n",
    "and the inductive hypothesis named the first sum $M_{k-1}$. The two correction\n",
    "terms combine to $\\frac{k-1}{k}\\,\\delta^2 = \\delta \\cdot \\delta \\frac{k-1}{k}\n",
    "= (x_k - m_{k-1})(x_k - m_k)$, which is precisely what the recursion adds.\n",
    "$\\blacksquare$\n",
    "\n",
    "Numerically, every quantity Welford touches is of size $\\sigma$, not $\\mu$,\n",
    "so the $\\mu^2/\\sigma^2$ amplification never materializes. The test:\n",
    "$10^5$ samples with mean $10^9$ and true variance $1$, all in float64:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-welford"
   },
   "source": [
    "rng = np.random.default_rng(0)\n",
    "x = 1e9 + rng.normal(0.0, 1.0, size=100_000)    # huge mean, unit variance\n",
    "\n",
    "naive = (x**2).mean() - x.mean()**2             # one pass, cancels\n",
    "two_pass = ((x - x.mean())**2).mean()           # subtract the mean first\n",
    "\n",
    "mean, m2 = 0.0, 0.0                             # Welford: one pass, stable\n",
    "for k, xk in enumerate(x, start=1):\n",
    "    delta = xk - mean\n",
    "    mean += delta / k                           # m_k\n",
    "    m2 += delta * (xk - mean)                   # M_k\n",
    "welford = m2 / len(x)\n",
    "\n",
    "print(f'naive E[x^2] - E[x]^2 : {naive:12.6f}')\n",
    "print(f'Welford, one pass     : {welford:12.6f}')\n",
    "print(f'two-pass reference    : {two_pass:12.6f}')"
   ],
   "execution_count": 8,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "naive E[x^2] - E[x]^2 :   384.000000\nWelford, one pass     :     1.000257\ntwo-pass reference    :     1.000257\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-12 fw=all -->\n",
    "\n",
    "The naive formula reports a variance in the *hundreds*: off by a factor of\n",
    "several hundred, in *double* precision, on a statistic every analyst computes\n",
    "daily. The answer is pure amplified rounding noise, so even its sign depends\n",
    "on the summation order of your NumPy build: the same cell printed $384$ under\n",
    "one build and $-256$, a negative variance, under another. Welford's\n",
    "one-pass answer $1.000257$ agrees with the two-pass reference to eight\n",
    "significant digits on every build. This recursion (and its batch-merging\n",
    "generalization, which you will derive in the exercises) is how `BatchNorm`\n",
    "layers (that section) and streaming-statistics utilities track\n",
    "running moments: one pass, bounded memory, no cancellation.\n",
    "\n",
    "The build-dependence is general, because it belongs to summation itself.\n",
    "Summing $n$ floats one after another commits one\n",
    "$(1 + \\delta)$ factor per addition, and the worst case compounds to a\n",
    "relative error of about $n u$ (at $n = 10^{5}$ in the cell above, some\n",
    "$10^{5}$ units of roundoff feeding the cancellation). **Pairwise summation**\n",
    "recursively sums halves, so each term passes through only $\\log_2 n$\n",
    "additions and the error growth drops to $O(u \\log n)$; this is what NumPy\n",
    "does inside `sum`, and its build-dependent blocking is why the naive\n",
    "formula's noise changed sign between builds. **Kahan (compensated)\n",
    "summation** carries each addition's rounding error explicitly in a second\n",
    "accumulator and drives the growth to $O(u)$, independent of $n$\n",
    "[@Kahan.1965; @Higham.2002]. Welford composes with either: the pairwise\n",
    "merge rule you will derive in Exercise 4 is precisely Welford in pairwise\n",
    "form, and it is how running moments are combined across devices.\n",
    "\n",
    "## Conditioning\n",
    "\n",
    "### Backward and Forward Error\n",
    "\n",
    "So far the *algorithms* were at fault, and rewriting them fixed everything.\n",
    "The deeper half of numerical analysis begins with a different question: when a\n",
    "computed answer is wrong, is the algorithm to blame, or the problem?\n",
    "@Higham.2002 makes the split precise with two definitions. The\n",
    "**forward error** is what you care about: the distance between the computed\n",
    "answer $\\hat{\\mathbf{x}}$ and the true answer $\\mathbf{x}$. The **backward\n",
    "error** is what the algorithm should be judged by: the size of the smallest\n",
    "perturbation of the *inputs* for which $\\hat{\\mathbf{x}}$ would be exactly\n",
    "correct: \"you got the right answer to a nearby question; how nearby?\" An\n",
    "algorithm with backward error at the rounding floor $u$ is called\n",
    "**backward stable**: it did everything that can be asked of finite-precision\n",
    "arithmetic, since merely *storing* the inputs already perturbs them by $u$.\n",
    "Gaussian elimination with pivoting (elimination with row swaps that keep the\n",
    "multipliers small; the algorithm inside `np.linalg.solve`) is backward\n",
    "stable in practice, though not in the worst case, where its growth factor\n",
    "can reach $2^{n-1}$ [@Higham.2002]; the SVD is backward stable\n",
    "outright; the naive variance formula is neither.\n",
    "\n",
    "What converts a small backward error into a possibly-large forward error is a\n",
    "property of the *problem*, and for linear systems it is exactly the condition\n",
    "number $\\kappa(\\mathbf{A}) = \\sigma_1/\\sigma_n$ of\n",
    "that section.\n",
    "\n",
    "**Proposition (forward error $\\le$ condition number $\\times$ backward\n",
    "error).** *Let $\\mathbf{A}$ be invertible, let\n",
    "$\\mathbf{A}\\mathbf{x} = \\mathbf{b}$, and suppose the computed\n",
    "$\\hat{\\mathbf{x}}$ exactly solves a nearby system,\n",
    "$(\\mathbf{A} + \\delta\\mathbf{A})\\,\\hat{\\mathbf{x}} = \\mathbf{b}$ with\n",
    "$\\|\\delta\\mathbf{A}\\| \\le \\varepsilon \\|\\mathbf{A}\\|$. Then*\n",
    "\n",
    "$$\n",
    "\\frac{\\|\\hat{\\mathbf{x}} - \\mathbf{x}\\|}{\\|\\hat{\\mathbf{x}}\\|}\n",
    "\\;\\le\\; \\kappa(\\mathbf{A})\\, \\varepsilon .\n",
    "$$\n",
    "\n",
    "**Proof.** Subtracting $\\mathbf{A}\\mathbf{x} = \\mathbf{b}$ from\n",
    "$(\\mathbf{A} + \\delta\\mathbf{A})\\hat{\\mathbf{x}} = \\mathbf{b}$ gives\n",
    "$\\mathbf{A}(\\hat{\\mathbf{x}} - \\mathbf{x}) = -\\delta\\mathbf{A}\\,\\hat{\\mathbf{x}}$,\n",
    "hence $\\hat{\\mathbf{x}} - \\mathbf{x} = -\\mathbf{A}^{-1}\\delta\\mathbf{A}\\,\\hat{\\mathbf{x}}$\n",
    "and\n",
    "\n",
    "$$\n",
    "\\|\\hat{\\mathbf{x}} - \\mathbf{x}\\|\n",
    "\\le \\|\\mathbf{A}^{-1}\\|\\, \\|\\delta\\mathbf{A}\\|\\, \\|\\hat{\\mathbf{x}}\\|\n",
    "\\le \\|\\mathbf{A}^{-1}\\|\\, \\|\\mathbf{A}\\|\\, \\varepsilon\\, \\|\\hat{\\mathbf{x}}\\|\n",
    "= \\kappa(\\mathbf{A})\\,\\varepsilon\\, \\|\\hat{\\mathbf{x}}\\|,\n",
    "$$\n",
    "\n",
    "using the operator-norm identities $\\|\\mathbf{A}\\| = \\sigma_1$ and\n",
    "$\\|\\mathbf{A}^{-1}\\| = 1/\\sigma_n$ from that section.\n",
    "$\\blacksquare$\n",
    "\n",
    "(The error here is measured relative to $\\hat{\\mathbf{x}}$; for small\n",
    "$\\varepsilon$ this matches the error relative to $\\mathbf{x}$ to first\n",
    "order.) This one inequality is the *division of labor* of numerical\n",
    "computing: the algorithm's job is to make $\\varepsilon$ small (backward\n",
    "stability gives $\\varepsilon \\approx u$), the problem's conditioning decides\n",
    "how much of that smallness survives, and the user's rule of thumb falls out\n",
    "by taking $\\log_{10}$ of both sides:\n",
    "\n",
    "$$\n",
    "\\textrm{correct digits in } \\hat{\\mathbf{x}}\n",
    "\\;\\approx\\; \\textrm{digits carried by the format} \\;-\\; \\log_{10} \\kappa(\\mathbf{A}).\n",
    "$$\n",
    "\n",
    "A backward-stable solve in float64 carries about $16$ digits, so\n",
    "$\\kappa = 10^{k}$ costs you $k$ of them, and at $\\kappa \\approx 10^{16}$\n",
    "the answer is pure noise even though the algorithm was flawless.\n",
    "\n",
    "### The Condition Number of a Linear System\n",
    "\n",
    "The same $\\kappa$ also governs sensitivity to errors in the right-hand side,\n",
    "the data in a least-squares problem: if $\\mathbf{A}\\mathbf{x} = \\mathbf{b}$\n",
    "and $\\mathbf{A}(\\mathbf{x} + \\delta\\mathbf{x}) = \\mathbf{b} + \\delta\\mathbf{b}$,\n",
    "then $\\|\\delta\\mathbf{x}\\|/\\|\\mathbf{x}\\| \\le\n",
    "\\kappa(\\mathbf{A})\\,\\|\\delta\\mathbf{b}\\|/\\|\\mathbf{b}\\|$, the perturbation\n",
    "bound the equation proved (together with the worst-case\n",
    "construction showing it is tight) in that section.\n",
    "\n",
    "Now we measure the digit loss on the most famously ill-conditioned matrix\n",
    "in the business: the **Hilbert matrix** $H_{ij} = 1/(i + j - 1)$, whose\n",
    "condition number grows exponentially with $n$. We solve\n",
    "$\\mathbf{H}\\mathbf{x} = \\mathbf{b}$ with the answer rigged to be\n",
    "$\\mathbf{x} = \\mathbf{1}$, and tabulate the forward error, the digits that\n",
    "survive, and the *backward* error, computed as the scaled residual\n",
    "$\\|\\mathbf{H}\\hat{\\mathbf{x}} - \\mathbf{b}\\| / (\\|\\mathbf{H}\\|\\,\\|\\hat{\\mathbf{x}}\\|)$.\n",
    "That this residual ratio equals the smallest relative perturbation of\n",
    "$\\mathbf{H}$ making $\\hat{\\mathbf{x}}$ exact is a classical theorem of\n",
    "Rigal--Gaches (see [@Higham.2002]), which is what lets a single\n",
    "computable number stand in for the definition's minimization:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-hilbert"
   },
   "source": [
    "print(' n      kappa   log10 kappa   forward error  correct digits  backward error')\n",
    "for n in [4, 6, 8, 10, 12]:\n",
    "    i = np.arange(n)\n",
    "    H = 1.0 / (1.0 + i[:, None] + i[None, :])   # Hilbert matrix\n",
    "    x_true = np.ones(n)\n",
    "    b = H @ x_true\n",
    "    x_hat = np.linalg.solve(H, b)\n",
    "    kappa = np.linalg.cond(H)\n",
    "    fwd = np.linalg.norm(x_hat - x_true) / np.linalg.norm(x_true)\n",
    "    bwd = (np.linalg.norm(H @ x_hat - b)\n",
    "           / (np.linalg.norm(H, 2) * np.linalg.norm(x_hat)))\n",
    "    print(f'{n:2d}  {kappa:9.1e}  {np.log10(kappa):8.1f}  {fwd:14.1e}  '\n",
    "          f'{-np.log10(fwd):11.1f}     {bwd:11.1e}')"
   ],
   "execution_count": 9,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      " n      kappa   log10 kappa   forward error  correct digits  backward error\n 4    1.6e+04       4.2         4.1e-14         13.4         0.0e+00\n 6    1.5e+07       7.2         1.4e-10          9.8         1.3e-16\n 8    1.5e+10      10.2         6.1e-08          7.2         5.7e-17\n10    1.6e+13      13.2         8.7e-05          4.1         1.0e-16\n12    1.8e+16      16.2         3.2e-01          0.5         9.8e-17\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-13 fw=all -->\n",
    "\n",
    "Read the table row by row against the rule of thumb. At $n = 4$,\n",
    "$\\log_{10}\\kappa \\approx 4.2$ and about $13$ digits survive of float64's\n",
    "$16$; by $n = 8$, $\\log_{10}\\kappa \\approx 10.2$ and about $7$ survive; at\n",
    "$n = 12$, $\\log_{10}\\kappa \\approx 16.2$ and barely one correct digit remains:\n",
    "the answer is essentially noise. (The trailing decimals of the error\n",
    "column vary with your LAPACK build; the staircase does not.) Meanwhile the backward-error column never\n",
    "leaves the $10^{-16}$ floor: *the solver is blameless at every row*. Each\n",
    "$\\hat{\\mathbf{x}}$ exactly solves a system one part in $10^{16}$ away from\n",
    "the one we posed; the Hilbert matrix simply maps that invisible perturbation\n",
    "to a visible one, exactly as the equation licenses it\n",
    "to. Geometrically, large $\\kappa$ means the level sets of\n",
    "$\\|\\mathbf{A}\\mathbf{x} - \\mathbf{b}\\|^2$ are extremely elongated ellipsoids:\n",
    "the same narrow valley of the figure, down\n",
    "which gradient descent zig-zags. Sensitivity of the solve and slowness of the\n",
    "descent are *one geometric fact* viewed from two sides.\n",
    "\n",
    "### Why Normal Equations Square the Pain\n",
    "\n",
    "Least squares offers a vivid demonstration that the *route* to a solution can\n",
    "ruin conditioning even when the destination is fine. The textbook route to\n",
    "$\\min_{\\mathbf{w}} \\|\\mathbf{A}\\mathbf{w} - \\mathbf{b}\\|^2$ is the normal\n",
    "equations $\\mathbf{A}^\\top\\mathbf{A}\\,\\mathbf{w} = \\mathbf{A}^\\top\\mathbf{b}$,\n",
    "which replace a solve governed by $\\kappa(\\mathbf{A})$ with one governed\n",
    "by $\\kappa(\\mathbf{A}^\\top\\mathbf{A})$. Recall from\n",
    "that section that this substitution is quadratically\n",
    "bad: for any matrix $\\mathbf{A}$ with full column rank,\n",
    "\n",
    "$$\n",
    "\\kappa(\\mathbf{A}^\\top\\mathbf{A}) = \\kappa(\\mathbf{A})^2 ,\n",
    "$$\n",
    "\n",
    "since $\\mathbf{A}^\\top\\mathbf{A} = \\mathbf{V}\\boldsymbol{\\Sigma}^2\\mathbf{V}^\\top$\n",
    "(that section) has singular values $\\sigma_i^2$. What is\n",
    "new here is the measurement of the cost in digits. By the rule of thumb, the\n",
    "normal equations lose $2\\log_{10}\\kappa$\n",
    "digits where an SVD- or QR-based solve, which works on $\\mathbf{A}$ directly,\n",
    "loses $\\log_{10}\\kappa$ (QR factors $\\mathbf{A} = \\mathbf{Q}\\mathbf{R}$ with\n",
    "$\\mathbf{Q}$ orthonormal and $\\mathbf{R}$ triangular, the factorization\n",
    "behind the `qr` demo of that section).\n",
    "With $\\kappa(\\mathbf{A}) = 10^5$ the predicted gap\n",
    "is five digits wide, large enough to measure:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-normal-equations"
   },
   "source": [
    "rng = np.random.default_rng(1)\n",
    "m, n = 100, 10\n",
    "U, _ = np.linalg.qr(rng.normal(size=(m, n)))    # random orthonormal columns\n",
    "V, _ = np.linalg.qr(rng.normal(size=(n, n)))\n",
    "sigma = np.logspace(0, -5, n)                   # kappa(A) = 10^5 by design\n",
    "A = U * sigma @ V.T                             # A = U diag(sigma) V^T\n",
    "w_true = rng.normal(size=n)\n",
    "b = A @ w_true\n",
    "print(f'kappa(A) = {np.linalg.cond(A):.1e}   '\n",
    "      f'kappa(A^T A) = {np.linalg.cond(A.T @ A):.1e}')\n",
    "w_ne = np.linalg.solve(A.T @ A, A.T @ b)        # normal equations\n",
    "w_svd = np.linalg.lstsq(A, b, rcond=None)[0]    # SVD-based solve\n",
    "for name, w in [('normal equations', w_ne), ('SVD (lstsq)     ', w_svd)]:\n",
    "    err = np.linalg.norm(w - w_true) / np.linalg.norm(w_true)\n",
    "    print(f'{name}: relative error {err:.1e}  '\n",
    "          f'({-np.log10(err):.1f} correct digits)')"
   ],
   "execution_count": 10,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "kappa(A) = 1.0e+05   kappa(A^T A) = 1.0e+10\nnormal equations: relative error 7.7e-08  (7.1 correct digits)\nSVD (lstsq)     : relative error 5.5e-14  (13.3 correct digits)\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-14 fw=all -->\n",
    "\n",
    "The printout confirms both the identity ($\\kappa$ of\n",
    "$\\mathbf{A}^\\top\\mathbf{A}$ is $10^{10}$, the square of $10^{5}$) and its\n",
    "consequence: the normal equations recover about seven correct digits, the SVD\n",
    "route about thirteen. Same problem, same data, same float64: five to six\n",
    "digits of accuracy, right at the predicted $\\log_{10}\\kappa = 5$, forfeited\n",
    "to a bad route. This is why `lstsq` exists, why numerical libraries solve\n",
    "least-squares subproblems by QR or SVD, and why\n",
    "that section built the pseudoinverse from the SVD\n",
    "rather than from $(\\mathbf{A}^\\top\\mathbf{A})^{-1}\\mathbf{A}^\\top$.\n",
    "\n",
    "### Ridge Regularization as Preconditioning\n",
    "\n",
    "When $\\kappa(\\mathbf{A})$ itself is the problem (nearly collinear\n",
    "features, a rank-deficient design), no choice of route saves the original\n",
    "problem. A **preconditioner** transforms a problem to reduce its condition\n",
    "number without changing its solution; the per-coordinate rescalings of\n",
    "that section apply the same idea inside an\n",
    "optimizer. Ridge regularization conditions the problem the way a\n",
    "preconditioner does, with one difference we return to below: it changes the\n",
    "problem, and it changes it in exactly the right direction. Minimizing\n",
    "$\\|\\mathbf{A}\\mathbf{w} - \\mathbf{b}\\|^2 + \\lambda \\|\\mathbf{w}\\|^2$ yields\n",
    "\n",
    "$$\n",
    "\\mathbf{w}_\\lambda = (\\mathbf{A}^\\top\\mathbf{A} + \\lambda \\mathbf{I})^{-1} \\mathbf{A}^\\top \\mathbf{b},\n",
    "$$\n",
    "\n",
    "and the added $\\lambda\\mathbf{I}$ acts directly on the spectrum.\n",
    "\n",
    "**Proposition (ridge improves conditioning monotonically).** *Let $\\mathbf{A}$\n",
    "have singular values $\\sigma_1 \\ge \\cdots \\ge \\sigma_n \\ge 0$. For every\n",
    "$\\lambda > 0$ the matrix $\\mathbf{A}^\\top\\mathbf{A} + \\lambda\\mathbf{I}$ is\n",
    "symmetric positive definite (hence invertible, even when $\\mathbf{A}$ is rank\n",
    "deficient), with*\n",
    "\n",
    "$$\n",
    "\\kappa(\\mathbf{A}^\\top\\mathbf{A} + \\lambda\\mathbf{I})\n",
    "= \\frac{\\sigma_1^2 + \\lambda}{\\sigma_n^2 + \\lambda},\n",
    "$$\n",
    "\n",
    "*which is strictly decreasing in $\\lambda$ whenever $\\sigma_1 > \\sigma_n$ and\n",
    "tends to $1$ as $\\lambda \\to \\infty$.*\n",
    "\n",
    "**Proof.** Writing $\\mathbf{A}^\\top\\mathbf{A} = \\mathbf{V}\\boldsymbol{\\Sigma}^2\\mathbf{V}^\\top$\n",
    "as above, $\\mathbf{A}^\\top\\mathbf{A} + \\lambda\\mathbf{I} =\n",
    "\\mathbf{V}(\\boldsymbol{\\Sigma}^2 + \\lambda\\mathbf{I})\\mathbf{V}^\\top$: the\n",
    "same eigenvectors, every eigenvalue shifted up to $\\sigma_i^2 + \\lambda \\ge\n",
    "\\lambda > 0$. Positive definiteness and the equation\n",
    "follow. For monotonicity, with $a = \\sigma_1^2 > b = \\sigma_n^2$,\n",
    "\n",
    "$$\n",
    "\\frac{d}{d\\lambda} \\frac{a + \\lambda}{b + \\lambda}\n",
    "= \\frac{b - a}{(b + \\lambda)^2} < 0,\n",
    "$$\n",
    "\n",
    "and as $\\lambda \\to \\infty$ the ratio tends to $1$. $\\blacksquare$\n",
    "\n",
    "Adding $\\lambda\\mathbf{I}$ lifts the floor of the spectrum while barely moving\n",
    "its ceiling: the narrow valley we met above rounds out\n",
    "into a bowl, as the figure shows, and the\n",
    "single shift pays twice.\n",
    "The solve in the equation becomes more accurate\n",
    "(fewer digits lost, by the rule of thumb) and gradient descent on the ridge\n",
    "objective, whose contraction factor $(\\kappa - 1)/(\\kappa + 1)$ we derived in\n",
    "that section, becomes faster. The cell below\n",
    "measures both at once: for each $\\lambda$ it computes\n",
    "the equation and *runs* gradient descent to a fixed\n",
    "relative tolerance of $10^{-6}$, counting iterations.\n",
    "\n",
    "![Level sets of the least-squares objective before and after adding the ridge term $\\lambda \\|\\mathbf{w}\\|^2$. The penalty lifts every eigenvalue of $\\mathbf{A}^\\top\\mathbf{A}$ by $\\lambda$, rounding the elongated valley into a bowl: the condition number drops from $\\sigma_1^2/\\sigma_n^2$ toward $1$, so linear solves lose fewer digits and gradient descent takes fewer steps.](../img/mdl-opt-conditioning-ellipse.svg)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "numerical-stability-conditioning-ridge"
   },
   "source": [
    "rng = np.random.default_rng(2)\n",
    "m, n = 200, 20\n",
    "U, _ = np.linalg.qr(rng.normal(size=(m, n)))\n",
    "V, _ = np.linalg.qr(rng.normal(size=(n, n)))\n",
    "sigma = np.logspace(0, -2, n)                   # kappa(A) = 100\n",
    "A = U * sigma @ V.T\n",
    "b = A @ rng.normal(size=n)\n",
    "\n",
    "def gd_iterations(lam, tol=1e-6):\n",
    "    \"\"\"Iterations for GD on the ridge objective to reach tol.\"\"\"\n",
    "    M, g = A.T @ A + lam * np.eye(n), A.T @ b\n",
    "    w_star = np.linalg.solve(M, g)\n",
    "    mu, L = sigma.min()**2 + lam, sigma.max()**2 + lam\n",
    "    eta, w = 2.0 / (L + mu), np.zeros(n)        # optimal fixed step size\n",
    "    for k in range(1, 200_000):\n",
    "        w -= eta * (M @ w - g)\n",
    "        if np.linalg.norm(w - w_star) <= tol * np.linalg.norm(w_star):\n",
    "            return k\n",
    "\n",
    "lams = np.logspace(-4, 0, 9)\n",
    "kappas = (sigma.max()**2 + lams) / (sigma.min()**2 + lams)\n",
    "iters = np.array([gd_iterations(lam) for lam in lams])\n",
    "print('iterations / kappa:', (iters / kappas).round(2))\n",
    "d2l.plot(lams, [kappas, iters.astype(float)], 'lambda', None,\n",
    "         legend=['kappa(A^T A + lambda I)', 'GD iterations to 1e-6'],\n",
    "         xscale='log', yscale='log')"
   ],
   "execution_count": 11,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "iterations / kappa: [5.89 5.9  5.94 5.99 6.06 6.15 6.19 6.25 6.  ]\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=\"234.3pt\" height=\"183.35625pt\" viewBox=\"0 0 234.3 183.35625\" 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:54:52.990919</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 183.35625 \nL 234.3 183.35625 \nL 234.3 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 31.8 145.8 \nL 227.1 145.8 \nL 227.1 7.2 \nL 31.8 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 40.677273 145.8 \nL 40.677273 7.2 \n\" clip-path=\"url(#pb34d60afc9)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_2\">\n      <defs>\n       <path id=\"mef58777f43\" d=\"M 0 0 \nL 0 3.5 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#mef58777f43\" x=\"40.677273\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_1\">\n      <!-- $\\mathdefault{10^{-4}}$ -->\n      <g transform=\"translate(28.927273 160.398438) 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        <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        <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-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-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-34\" transform=\"translate(186.855469 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_2\">\n     <g id=\"line2d_3\">\n      <path d=\"M 85.063636 145.8 \nL 85.063636 7.2 \n\" clip-path=\"url(#pb34d60afc9)\" 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=\"#mef58777f43\" x=\"85.063636\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_2\">\n      <!-- $\\mathdefault{10^{-3}}$ -->\n      <g transform=\"translate(73.313636 160.398438) 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-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-33\" transform=\"translate(186.855469 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_3\">\n     <g id=\"line2d_5\">\n      <path d=\"M 129.45 145.8 \nL 129.45 7.2 \n\" clip-path=\"url(#pb34d60afc9)\" 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=\"#mef58777f43\" x=\"129.45\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_3\">\n      <!-- $\\mathdefault{10^{-2}}$ -->\n      <g transform=\"translate(117.7 160.398438) 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-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(186.855469 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_4\">\n     <g id=\"line2d_7\">\n      <path d=\"M 173.836364 145.8 \nL 173.836364 7.2 \n\" clip-path=\"url(#pb34d60afc9)\" 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=\"#mef58777f43\" x=\"173.836364\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_4\">\n      <!-- $\\mathdefault{10^{-1}}$ -->\n      <g transform=\"translate(162.086364 160.398438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(186.855469 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_5\">\n     <g id=\"line2d_9\">\n      <path d=\"M 218.222727 145.8 \nL 218.222727 7.2 \n\" clip-path=\"url(#pb34d60afc9)\" 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=\"#mef58777f43\" x=\"218.222727\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_5\">\n      <!-- $\\mathdefault{10^{0}}$ -->\n      <g transform=\"translate(209.422727 160.398438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_6\">\n     <g id=\"line2d_11\">\n      <defs>\n       <path id=\"mfcac2261cc\" d=\"M 0 0 \nL 0 2 \n\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"33.801738\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_7\">\n     <g id=\"line2d_12\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"36.37579\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_8\">\n     <g id=\"line2d_13\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"38.646264\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_9\">\n     <g id=\"line2d_14\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"54.0389\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_10\">\n     <g id=\"line2d_15\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"61.85495\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_11\">\n     <g id=\"line2d_16\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"67.400526\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_12\">\n     <g id=\"line2d_17\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"71.70201\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_13\">\n     <g id=\"line2d_18\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"75.216577\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_14\">\n     <g id=\"line2d_19\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"78.188102\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_15\">\n     <g id=\"line2d_20\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"80.762153\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_16\">\n     <g id=\"line2d_21\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"83.032628\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_17\">\n     <g id=\"line2d_22\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"98.425263\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_18\">\n     <g id=\"line2d_23\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"106.241314\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_19\">\n     <g id=\"line2d_24\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"111.78689\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_20\">\n     <g id=\"line2d_25\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"116.088373\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_21\">\n     <g id=\"line2d_26\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"119.602941\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_22\">\n     <g id=\"line2d_27\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"122.574465\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_23\">\n     <g id=\"line2d_28\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"125.148517\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_24\">\n     <g id=\"line2d_29\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"127.418991\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_25\">\n     <g id=\"line2d_30\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"142.811627\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_26\">\n     <g id=\"line2d_31\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"150.627678\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_27\">\n     <g id=\"line2d_32\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"156.173254\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_28\">\n     <g id=\"line2d_33\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"160.474737\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_29\">\n     <g id=\"line2d_34\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"163.989304\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_30\">\n     <g id=\"line2d_35\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"166.960829\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_31\">\n     <g id=\"line2d_36\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"169.534881\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_32\">\n     <g id=\"line2d_37\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"171.805355\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_33\">\n     <g id=\"line2d_38\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"187.19799\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_34\">\n     <g id=\"line2d_39\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"195.014041\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_35\">\n     <g id=\"line2d_40\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"200.559617\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_36\">\n     <g id=\"line2d_41\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"204.8611\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_37\">\n     <g id=\"line2d_42\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"208.375668\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_38\">\n     <g id=\"line2d_43\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"211.347193\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_39\">\n     <g id=\"line2d_44\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"213.921244\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_40\">\n     <g id=\"line2d_45\">\n      <g>\n       <use xlink:href=\"#mfcac2261cc\" x=\"216.191719\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_6\">\n     <!-- lambda -->\n     <g transform=\"translate(110.714063 174.076563) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-6c\" d=\"M 603 4863 \nL 1178 4863 \nL 1178 0 \nL 603 0 \nL 603 4863 \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-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-62\" d=\"M 3116 1747 \nQ 3116 2381 2855 2742 \nQ 2594 3103 2138 3103 \nQ 1681 3103 1420 2742 \nQ 1159 2381 1159 1747 \nQ 1159 1113 1420 752 \nQ 1681 391 2138 391 \nQ 2594 391 2855 752 \nQ 3116 1113 3116 1747 \nz\nM 1159 2969 \nQ 1341 3281 1617 3432 \nQ 1894 3584 2278 3584 \nQ 2916 3584 3314 3078 \nQ 3713 2572 3713 1747 \nQ 3713 922 3314 415 \nQ 2916 -91 2278 -91 \nQ 1894 -91 1617 61 \nQ 1341 213 1159 525 \nL 1159 0 \nL 581 0 \nL 581 4863 \nL 1159 4863 \nL 1159 2969 \nz\n\" transform=\"scale(0.015625)\"/>\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      </defs>\n      <use xlink:href=\"#DejaVuSans-6c\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(27.783203 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6d\" transform=\"translate(89.0625 0)\"/>\n      <use xlink:href=\"#DejaVuSans-62\" transform=\"translate(186.474609 0)\"/>\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(249.951172 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(313.427734 0)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"matplotlib.axis_2\">\n    <g id=\"ytick_1\">\n     <g id=\"line2d_46\">\n      <path d=\"M 31.8 118.368124 \nL 227.1 118.368124 \n\" clip-path=\"url(#pb34d60afc9)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_47\">\n      <defs>\n       <path id=\"m8416a39281\" d=\"M 0 0 \nL -3.5 0 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m8416a39281\" x=\"31.8\" y=\"118.368124\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_7\">\n      <!-- $\\mathdefault{10^{1}}$ -->\n      <g transform=\"translate(7.2 122.167343) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_2\">\n     <g id=\"line2d_48\">\n      <path d=\"M 31.8 88.137123 \nL 227.1 88.137123 \n\" clip-path=\"url(#pb34d60afc9)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_49\">\n      <g>\n       <use xlink:href=\"#m8416a39281\" x=\"31.8\" y=\"88.137123\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_8\">\n      <!-- $\\mathdefault{10^{2}}$ -->\n      <g transform=\"translate(7.2 91.936342) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_3\">\n     <g id=\"line2d_50\">\n      <path d=\"M 31.8 57.906122 \nL 227.1 57.906122 \n\" clip-path=\"url(#pb34d60afc9)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_51\">\n      <g>\n       <use xlink:href=\"#m8416a39281\" x=\"31.8\" y=\"57.906122\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_9\">\n      <!-- $\\mathdefault{10^{3}}$ -->\n      <g transform=\"translate(7.2 61.705341) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-33\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_4\">\n     <g id=\"line2d_52\">\n      <path d=\"M 31.8 27.675121 \nL 227.1 27.675121 \n\" clip-path=\"url(#pb34d60afc9)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_53\">\n      <g>\n       <use xlink:href=\"#m8416a39281\" x=\"31.8\" y=\"27.675121\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_10\">\n      <!-- $\\mathdefault{10^{4}}$ -->\n      <g transform=\"translate(7.2 31.47434) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-34\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_5\">\n     <g id=\"line2d_54\">\n      <defs>\n       <path id=\"m1e5acd6927\" d=\"M 0 0 \nL -2 0 \n\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"139.498687\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_6\">\n     <g id=\"line2d_55\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"134.175272\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_7\">\n     <g id=\"line2d_56\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"130.398249\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_8\">\n     <g id=\"line2d_57\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"127.468562\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_9\">\n     <g id=\"line2d_58\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"125.074834\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_10\">\n     <g id=\"line2d_59\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"123.050966\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_11\">\n     <g id=\"line2d_60\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"121.297811\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_12\">\n     <g id=\"line2d_61\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"119.751419\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_13\">\n     <g id=\"line2d_62\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"109.267686\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_14\">\n     <g id=\"line2d_63\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"103.944271\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_15\">\n     <g id=\"line2d_64\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"100.167248\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_16\">\n     <g id=\"line2d_65\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"97.237561\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_17\">\n     <g id=\"line2d_66\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"94.843833\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_18\">\n     <g id=\"line2d_67\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"92.819965\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_19\">\n     <g id=\"line2d_68\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"91.06681\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_20\">\n     <g id=\"line2d_69\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"89.520418\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_21\">\n     <g id=\"line2d_70\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"79.036685\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_22\">\n     <g id=\"line2d_71\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"73.71327\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_23\">\n     <g id=\"line2d_72\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"69.936247\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_24\">\n     <g id=\"line2d_73\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"67.00656\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_25\">\n     <g id=\"line2d_74\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"64.612832\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_26\">\n     <g id=\"line2d_75\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"62.588964\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_27\">\n     <g id=\"line2d_76\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"60.835809\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_28\">\n     <g id=\"line2d_77\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"59.289417\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_29\">\n     <g id=\"line2d_78\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"48.805684\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_30\">\n     <g id=\"line2d_79\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"43.482269\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_31\">\n     <g id=\"line2d_80\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"39.705246\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_32\">\n     <g id=\"line2d_81\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"36.775559\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_33\">\n     <g id=\"line2d_82\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"34.381831\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_34\">\n     <g id=\"line2d_83\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"32.357963\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_35\">\n     <g id=\"line2d_84\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"30.604808\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_36\">\n     <g id=\"line2d_85\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"29.058416\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_37\">\n     <g id=\"line2d_86\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"18.574683\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_38\">\n     <g id=\"line2d_87\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"13.251268\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_39\">\n     <g id=\"line2d_88\">\n      <g>\n       <use xlink:href=\"#m1e5acd6927\" x=\"31.8\" y=\"9.474245\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n   </g>\n   <g id=\"line2d_89\">\n    <path d=\"M 40.677273 36.774247 \nL 62.870455 46.393968 \nL 85.063636 59.144342 \nL 107.256818 73.388921 \nL 129.45 88.137123 \nL 151.643182 102.885325 \nL 173.836364 117.129905 \nL 196.029545 129.880279 \nL 218.222727 139.5 \n\" clip-path=\"url(#pb34d60afc9)\" style=\"fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square\"/>\n   </g>\n   <g id=\"line2d_90\">\n    <path d=\"M 40.677273 13.5 \nL 62.870455 23.088871 \nL 85.063636 35.757834 \nL 107.256818 49.886151 \nL 129.45 64.482193 \nL 151.643182 79.036685 \nL 173.836364 93.200547 \nL 196.029545 105.823065 \nL 218.222727 115.974396 \n\" clip-path=\"url(#pb34d60afc9)\" style=\"fill: none; stroke-dasharray: 5.55,2.4; stroke-dashoffset: 0; stroke: #bf00bf; stroke-width: 1.5\"/>\n   </g>\n   <g id=\"patch_3\">\n    <path d=\"M 31.8 145.8 \nL 31.8 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 227.1 145.8 \nL 227.1 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 31.8 145.8 \nL 227.1 145.8 \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 31.8 7.2 \nL 227.1 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 38.8 140.8 \nL 198.854688 140.8 \nQ 200.854688 140.8 200.854688 138.8 \nL 200.854688 110.44375 \nQ 200.854688 108.44375 198.854688 108.44375 \nL 38.8 108.44375 \nQ 36.8 108.44375 36.8 110.44375 \nL 36.8 138.8 \nQ 36.8 140.8 38.8 140.8 \nz\n\" style=\"fill: #ffffff; opacity: 0.8; stroke: #cccccc; stroke-linejoin: miter\"/>\n    </g>\n    <g id=\"line2d_91\">\n     <path d=\"M 40.8 116.542188 \nL 50.8 116.542188 \nL 60.8 116.542188 \n\" style=\"fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square\"/>\n    </g>\n    <g id=\"text_11\">\n     <!-- kappa(A^T A + lambda I) -->\n     <g transform=\"translate(68.8 120.042188) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-6b\" d=\"M 581 4863 \nL 1159 4863 \nL 1159 1991 \nL 2875 3500 \nL 3609 3500 \nL 1753 1863 \nL 3688 0 \nL 2938 0 \nL 1159 1709 \nL 1159 0 \nL 581 0 \nL 581 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-70\" d=\"M 1159 525 \nL 1159 -1331 \nL 581 -1331 \nL 581 3500 \nL 1159 3500 \nL 1159 2969 \nQ 1341 3281 1617 3432 \nQ 1894 3584 2278 3584 \nQ 2916 3584 3314 3078 \nQ 3713 2572 3713 1747 \nQ 3713 922 3314 415 \nQ 2916 -91 2278 -91 \nQ 1894 -91 1617 61 \nQ 1341 213 1159 525 \nz\nM 3116 1747 \nQ 3116 2381 2855 2742 \nQ 2594 3103 2138 3103 \nQ 1681 3103 1420 2742 \nQ 1159 2381 1159 1747 \nQ 1159 1113 1420 752 \nQ 1681 391 2138 391 \nQ 2594 391 2855 752 \nQ 3116 1113 3116 1747 \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-41\" d=\"M 2188 4044 \nL 1331 1722 \nL 3047 1722 \nL 2188 4044 \nz\nM 1831 4666 \nL 2547 4666 \nL 4325 0 \nL 3669 0 \nL 3244 1197 \nL 1141 1197 \nL 716 0 \nL 50 0 \nL 1831 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-5e\" d=\"M 2988 4666 \nL 4684 2925 \nL 4056 2925 \nL 2681 4159 \nL 1306 2925 \nL 678 2925 \nL 2375 4666 \nL 2988 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-54\" d=\"M -19 4666 \nL 3928 4666 \nL 3928 4134 \nL 2272 4134 \nL 2272 0 \nL 1638 0 \nL 1638 4134 \nL -19 4134 \nL -19 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-20\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-2b\" d=\"M 2944 4013 \nL 2944 2272 \nL 4684 2272 \nL 4684 1741 \nL 2944 1741 \nL 2944 0 \nL 2419 0 \nL 2419 1741 \nL 678 1741 \nL 678 2272 \nL 2419 2272 \nL 2419 4013 \nL 2944 4013 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-49\" d=\"M 628 4666 \nL 1259 4666 \nL 1259 0 \nL 628 0 \nL 628 4666 \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-6b\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(56.160156 0)\"/>\n      <use xlink:href=\"#DejaVuSans-70\" transform=\"translate(117.439453 0)\"/>\n      <use xlink:href=\"#DejaVuSans-70\" transform=\"translate(180.916016 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(244.392578 0)\"/>\n      <use xlink:href=\"#DejaVuSans-28\" transform=\"translate(305.671875 0)\"/>\n      <use xlink:href=\"#DejaVuSans-41\" transform=\"translate(344.685547 0)\"/>\n      <use xlink:href=\"#DejaVuSans-5e\" transform=\"translate(413.09375 0)\"/>\n      <use xlink:href=\"#DejaVuSans-54\" transform=\"translate(496.882812 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(557.966797 0)\"/>\n      <use xlink:href=\"#DejaVuSans-41\" transform=\"translate(589.753906 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(658.162109 0)\"/>\n      <use xlink:href=\"#DejaVuSans-2b\" transform=\"translate(689.949219 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(773.738281 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6c\" transform=\"translate(805.525391 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(833.308594 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6d\" transform=\"translate(894.587891 0)\"/>\n      <use xlink:href=\"#DejaVuSans-62\" transform=\"translate(992 0)\"/>\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(1055.476562 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(1118.953125 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(1180.232422 0)\"/>\n      <use xlink:href=\"#DejaVuSans-49\" transform=\"translate(1212.019531 0)\"/>\n      <use xlink:href=\"#DejaVuSans-29\" transform=\"translate(1241.511719 0)\"/>\n     </g>\n    </g>\n    <g id=\"line2d_92\">\n     <path d=\"M 40.8 131.220313 \nL 50.8 131.220313 \nL 60.8 131.220313 \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_12\">\n     <!-- GD iterations to 1e-6 -->\n     <g transform=\"translate(68.8 134.720313) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-47\" d=\"M 3809 666 \nL 3809 1919 \nL 2778 1919 \nL 2778 2438 \nL 4434 2438 \nL 4434 434 \nQ 4069 175 3628 42 \nQ 3188 -91 2688 -91 \nQ 1594 -91 976 548 \nQ 359 1188 359 2328 \nQ 359 3472 976 4111 \nQ 1594 4750 2688 4750 \nQ 3144 4750 3555 4637 \nQ 3966 4525 4313 4306 \nL 4313 3634 \nQ 3963 3931 3569 4081 \nQ 3175 4231 2741 4231 \nQ 1884 4231 1454 3753 \nQ 1025 3275 1025 2328 \nQ 1025 1384 1454 906 \nQ 1884 428 2741 428 \nQ 3075 428 3337 486 \nQ 3600 544 3809 666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-44\" d=\"M 1259 4147 \nL 1259 519 \nL 2022 519 \nQ 2988 519 3436 956 \nQ 3884 1394 3884 2338 \nQ 3884 3275 3436 3711 \nQ 2988 4147 2022 4147 \nL 1259 4147 \nz\nM 628 4666 \nL 1925 4666 \nQ 3281 4666 3915 4102 \nQ 4550 3538 4550 2338 \nQ 4550 1131 3912 565 \nQ 3275 0 1925 0 \nL 628 0 \nL 628 4666 \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-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-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-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-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-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-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-47\"/>\n      <use xlink:href=\"#DejaVuSans-44\" transform=\"translate(77.490234 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(154.492188 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(186.279297 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(214.0625 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(253.271484 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(314.794922 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(355.908203 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(417.1875 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(456.396484 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(484.179688 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(545.361328 0)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(608.740234 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(660.839844 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(692.626953 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(731.835938 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(793.017578 0)\"/>\n      <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(824.804688 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(888.427734 0)\"/>\n      <use xlink:href=\"#DejaVuSans-2d\" transform=\"translate(949.951172 0)\"/>\n      <use xlink:href=\"#DejaVuSans-36\" transform=\"translate(986.035156 0)\"/>\n     </g>\n    </g>\n   </g>\n  </g>\n </g>\n <defs>\n  <clipPath id=\"pb34d60afc9\">\n   <rect x=\"31.8\" y=\"7.2\" width=\"195.3\" height=\"138.6\"/>\n  </clipPath>\n </defs>\n</svg>\n"
     },
     "metadata": {}
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-numerical-stability-conditioning-md-15 fw=all -->\n",
    "\n",
    "The two curves fall together across four orders of magnitude of $\\lambda$,\n",
    "from $\\kappa \\approx 5000$ and roughly $29{,}000$ iterations at\n",
    "$\\lambda = 10^{-4}$ down to $\\kappa = 2$ and a dozen iterations at\n",
    "$\\lambda = 1$, and the printed ratio sits between $5.89$ and $6.25$\n",
    "throughout: iteration count is a constant multiple of $\\kappa$, the constant\n",
    "being about $\\tfrac12 \\ln(1/\\textrm{tol}) \\approx 6.9$ predicted by the\n",
    "$(\\kappa - 1)/(\\kappa + 1)$ contraction (slightly less here because a random\n",
    "initial error is not perfectly aligned with the slowest eigendirection).\n",
    "Ridge thus conditions the problem exactly as a preconditioner would, with\n",
    "one difference: a true preconditioner leaves the minimizer unchanged, while\n",
    "ridge biases the solution, shrinking $\\mathbf{w}_\\lambda$ toward\n",
    "$\\mathbf{0}$; that section showed the\n",
    "precise sense in which the penalty $\\lambda\\|\\mathbf{w}\\|^2$ is the\n",
    "Lagrangian counterpart of a norm constraint. In practice, the $\\lambda$ you\n",
    "were going to add anyway, for statistical reasons, has been quietly\n",
    "stabilizing your arithmetic and accelerating your optimizer the whole time.\n",
    "\n",
    "## Summary\n",
    "\n",
    "* Floating point is scientific notation with a finite mantissa: relative\n",
    "  precision $\\varepsilon_{\\text{mach}}$ ($2^{-23}$ for fp32, $2^{-10}$ for\n",
    "  fp16, $2^{-7}$ for bfloat16), absolute gaps that double at every power of\n",
    "  two, and overflow/underflow thresholds ($e^x$ dies at $x \\approx 88.7$ in\n",
    "  fp32, $x \\approx 11.1$ in fp16). Mixed precision is engineering around\n",
    "  these limits: loss scaling fights fp16 underflow; bfloat16 trades mantissa\n",
    "  for fp32's exponent range; fp8 (E4M3/E5M2) tightens both budgets further\n",
    "  and makes per-tensor scaling mandatory.\n",
    "* Softmax is shift-invariant, so subtract the max before exponentiating;\n",
    "  shifted log-sum-exp is an exact identity that prevents exponential overflow\n",
    "  for finite logits (provided the final result itself is representable);\n",
    "  cross-entropy should be computed from logits as\n",
    "  $\\mathrm{lse}(\\mathbf{z}) - z_y$. Pass logits to your loss; the\n",
    "  from-probabilities route ends in `inf`, `NaN`, or (worse) a silently\n",
    "  clipped gradient.\n",
    "* Catastrophic cancellation: subtracting nearly equal numbers amplifies\n",
    "  existing relative error by $(|a|+|b|)/|a-b|$. The cure is reformulation,\n",
    "  not precision: `log1p`/`expm1` near zero, and Welford's recursion\n",
    "  $m_k = m_{k-1} + (x_k - m_{k-1})/k$,\n",
    "  $M_k = M_{k-1} + (x_k - m_{k-1})(x_k - m_k)$ in place of\n",
    "  $\\mathbb{E}[x^2] - \\mathbb{E}[x]^2$, which lost a factor of $400$ in\n",
    "  float64.\n",
    "* Backward error judges the algorithm, the condition number judges the\n",
    "  problem, and forward error $\\le \\kappa \\times$ backward error connects\n",
    "  them: a backward-stable solve loses about $\\log_{10}\\kappa$ digits, as the\n",
    "  Hilbert-matrix table showed digit for digit.\n",
    "* Normal equations square the condition number\n",
    "  ($\\kappa(\\mathbf{A}^\\top\\mathbf{A}) = \\kappa(\\mathbf{A})^2$; six digits\n",
    "  lost vs. the SVD route in our experiment); ridge regularization moves it\n",
    "  the other way, $\\kappa = (\\sigma_1^2 + \\lambda)/(\\sigma_n^2 + \\lambda)$,\n",
    "  monotonically improving both solve accuracy and the gradient-descent rate:\n",
    "  one number, two consequences, and one knob that helps both.\n",
    "\n",
    "## Exercises\n",
    "\n",
    "1. Compute $\\varepsilon_{\\text{mach}}$ for float32 by a halving loop: start\n",
    "   from $e = 1$ and halve until $1 + e == 1$ in float32. Why does the loop\n",
    "   exit at $e = 2^{-24}$ rather than at $2^{-23}$, and how does that relate\n",
    "   to the unit roundoff $u$ in the equation? Repeat the\n",
    "   loop with bfloat16 arithmetic (emulated or native) and confirm $2^{-7}$.\n",
    "2. Find all integer logits $x$ for which $e^x$ overflows in fp16 but not in\n",
    "   fp32. A network's final layer emits activations of size $\\approx 30$:\n",
    "   explain, with the numbers from this section, why training it in fp16\n",
    "   without loss scaling fails even though the softmax probabilities are\n",
    "   perfectly representable.\n",
    "3. Prove that $\\nabla\\, \\mathrm{lse}(\\mathbf{z}) = \\mathrm{softmax}(\\mathbf{z})$\n",
    "   and use that section to conclude that lse is convex.\n",
    "   Then show that the cross-entropy the equation has\n",
    "   gradient $\\mathrm{softmax}(\\mathbf{z}) - \\mathbf{e}_y$, another reason to\n",
    "   compute the loss from logits.\n",
    "4. Construct a small dataset (three numbers suffice) for which the naive\n",
    "   variance formula returns a strictly *negative* number in float64, and\n",
    "   verify that Welford's recursion the equation gets it right.\n",
    "   Then derive the *pairwise merge* rule: given $(m, M, k)$ for two disjoint\n",
    "   batches, express the combined statistics exactly; this is how the\n",
    "   computation parallelizes across devices.\n",
    "5. Show that if $a$ and $b$ carry relative errors of size $u$, the relative\n",
    "   error of the computed $a - b$ can be as large as\n",
    "   $(|a| + |b|)\\,u / |a - b|$, and that the subtraction itself adds no\n",
    "   rounding error when $a/2 \\le b \\le 2a$ (Sterbenz's lemma\n",
    "   [@Sterbenz.1974]; prove it for floats with the same exponent).\n",
    "6. Rewrite each of the following to avoid cancellation, and check one of them\n",
    "   numerically in float32: $\\sqrt{x + 1} - \\sqrt{x}$ for large $x$;\n",
    "   $1 - \\cos x$ for small $x$; the smaller root of $ax^2 + bx + c = 0$ when\n",
    "   $b^2 \\gg 4ac$.\n",
    "7. For the Hilbert experiment, compute the backward error of each solve with\n",
    "   respect to the right-hand side, $\\|\\mathbf{H}\\hat{\\mathbf{x}} - \\mathbf{b}\\|/\\|\\mathbf{b}\\|$,\n",
    "   and verify that the forward error is bounded by\n",
    "   the equation applied to that perturbation. Where\n",
    "   in the table is the bound tightest?\n",
    "8. Let $\\mathbf{A}$ have $\\sigma_1 = 1$ and $\\sigma_n = 10^{-3}$. Using\n",
    "   the equation and the GD contraction factor\n",
    "   $(\\kappa - 1)/(\\kappa + 1)$ from\n",
    "   that section, compute the condition\n",
    "   number and the predicted iteration count (to relative error $10^{-6}$)\n",
    "   for $\\lambda \\in \\{0, 10^{-4}, 10^{-2}\\}$. Then explain, via\n",
    "   that section, which constrained\n",
    "   problem each $\\lambda$ implicitly solves.\n",
    "\n",
    "## Discussions\n",
    "\n",
    "The condition number that set\n",
    "gradient descent's speed in that section is\n",
    "the same number that sets a linear solve's accuracy, and ridge regularization\n",
    "(the Lagrangian twin of a norm constraint, per\n",
    "that section) improves both at once.\n",
    "The stability toolkit travels well beyond this chapter: max-subtraction and\n",
    "log-sum-exp power every softmax and attention layer; log-space arithmetic is\n",
    "what makes naive Bayes (that section) and probabilistic\n",
    "inference generally feasible; the from-logits cross-entropy is the computation\n",
    "that section analyzes; and Welford-style running\n",
    "moments live inside batch normalization. One neighboring story is told\n",
    "elsewhere: the *depth* dimension of numerical stability (vanishing and\n",
    "exploding gradients as the conditioning of Jacobian products through many\n",
    "layers, and the initialization schemes that tame them) is the subject of\n",
    "that section in the main text, and\n",
    "that section takes up the continuous-time view of how such\n",
    "errors and noise propagate through training. For the full theory, error\n",
    "analysis of every algorithm in this book's linear-algebra substrate, the\n",
    "reference remains @Higham.2002 .\n",
    "\n",
    "[Discussions](https://d2l.discourse.group/t/numerical-stability-and-conditioning)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "d2l-jax",
   "display_name": "d2l (JAX)",
   "language": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
