The well on the far side of the band

The numbers below the smallest one drew the subnormal band as a staircase that stops halving, counted the pairs that flush-to-zero makes subtract to nothing, and then drew a careful line around what the band is worth. The identity $x - y = 0 \Rightarrow x = y$ is real. The claim that subnormals save computations mostly is not: “a computation whose intermediate values have decayed into the subnormal range is a computation in trouble for reasons that a few extra bits at the bottom will not repair.” Subnormals turn a cliff into a slope, and a slope down to zero is still a slope down to zero.

That sentence is right about a quantity that ends in the band. It says nothing about a quantity that passes through the band on its way somewhere else — that is small for a few steps and then large again, and whose large values are computed from the small ones. Such quantities are not exotic. They are what a recurrence produces whenever the thing it tracks has to cross a barrier, and the simplest one in linear algebra is a stationary distribution.

A chain with two wells #

A birth–death chain moves one state left or right at each step, with probability $p_i$ of going up from state $i$ and $q_{i+1}$ of coming back down from $i+1$. Its stationary distribution satisfies detailed balance, $q_{i+1}\,\pi_{i+1} = p_i\,\pi_i$, so fixing $\pi_0 = 1$ gives every other state by one multiplication and one division:

$$\pi_{i+1} = \pi_i \, \frac{p_i}{q_{i+1}}.$$

That is forward substitution on a lower bidiagonal system — the balance equations written as $L\pi = e_1$ — and it is how every textbook computes the distribution of a queue, a population or a molecule hopping between conformations. The chain here has 48 states and a potential shaped like two wells, $\sin^2(\pi s)$ plus a slight tilt towards the near side, with Metropolis rates: the walker climbs out of a well with probability proportional to the Boltzmann factor of the climb. The temperature is set so that the top of the barrier has probability exactly $2^{-d}$ of the near well’s bottom, and $d$ — the depth, in octaves — is the only knob.

The quantity asked for is ordinary: the share of the probability in the far well, the right half of the chain. At a depth of ten octaves it is 0.337, and it falls slowly, to 0.245 at sixteen and 0.101 at thirty, because the tilt matters more as the temperature drops. Nothing about it is small.

The rates are stored in half precision, and the reference is the same recurrence in double precision on those stored rates, so every difference measured below is the arithmetic’s and none of it is the data’s. Five arithmetics are compared: half precision with gradual underflow; half precision flushing any result below the smallest normal to zero, as a fast-math mode does in hardware; a format with no subnormals that rounds correctly to what it does have; gradual half precision started at $\pi_0 = 1024$ rather than at one; and the same eleven-bit significand with no limit on the exponent at all.

What crossing the band looks like #

At twelve octaves the whole distribution stays above the smallest normal, $2^{-14} \approx 6.1\cdot 10^{-5}$, and the arithmetics agree to half precision. Turn the dial to sixteen and the barrier’s top sits two octaves into the band. Gradual underflow carries the chain through it a few bits short and brings it back up the other side: its far well holds 0.2454, the exact figure to four digits. Flush to zero sets the first state below $2^{-14}$ to zero, and every state after it is zero too, because every one of them is a multiple of that state — thirty-one states of forty-eight are exactly zero, the far well is empty, and a normalised answer puts the whole of the probability on the near side of the barrier.

That answer is not a NaN, not an infinity and not a warning. It is a probability vector that sums to one, satisfies the balance equations to the precision the norm can see (as measured below), and says with certainty that a walker never visits a well it actually spends a quarter of its time in.

The deepest barrier each arithmetic survives #

The hero figure at the top of the page is the sweep. With flush to zero the far well is right to half precision up to a depth of 13 octaves and exactly zero from 14 — the depth at which the barrier’s top first falls below the smallest normal. There is no transition: one octave the answer is fine and the next it is gone.

Gradual underflow stays within ten units of roundoff to sixteen octaves and then degrades a little under a bit an octave: 2.7·10⁻² at 18, 5.3·10⁻² at 19, 0.27 at 22, 0.46 at 23. At 24 and 25 it overshoots — the far well comes out 2.2 and 3.5 times too large, which the section on rounding at the floor explains — and from 26 its far well is zero as well. With its far well right to within half, it survives to 23 octaves against flush-to-zero’s 13.

Ten octaves. That is not a coincidence of this chain. Half precision’s significand has eleven bits, so its subnormal band spans $p - 1 = 10$ octaves, from $2^{-14}$ down to $2^{-24}$, and the recurrence survives exactly as long as the barrier’s top has somewhere representable to go.

Holding the exponent field at five bits and varying the significand from four bits to eleven makes the same point as a law. Flush to zero survives 12 or 13 octaves at every width, because what ends it is the smallest normal, which the significand does not move. Gradual underflow survives 16, 17, 18, 19, 20, 21, 22 and 23. What the subnormals add over flushing is 4, 5, 6, 6, 8, 9, 9 and 10 octaves, within one of $p - 1$ at every width and equal to it at eleven bits; the octave of slack is the rates themselves, rounded to $p$ bits, moving the barrier’s top a little from $2^{-d}$. So the band’s width, which the earlier essay derived as a fact about the format, is also a fact about what a computation can get through — and it grows with the precision, as that essay said the band did.

The error made at the bottom is carried #

What gradual underflow guarantees is absolute rather than relative: a value in the band is rounded to a multiple of $2^{-24}$, so its error is at most half of that however small the value is. At the bottom of a twenty-octave barrier, where the exact probability is $9.5\cdot 10^{-7}$, a spacing of $6\cdot 10^{-8}$ is a relative error of several per cent, and eighteen states sit in the band.

The profile shows what happens to that error next. On the near side of the band every state is right to $8\cdot 10^{-4}$ in every arithmetic, the ordinary rounding of eleven bits. On the far side, the state-by-state error of the gradual route is flat at 2.1 to 2.2·10⁻² — every state in the far well carries the error the bottom made, no more and no less. The climb back up multiplies each state by a ratio above one, and multiplication by an exact ratio carries a relative error through unchanged; it neither amplifies it nor repairs it. The unbounded route and the one started at 1024 never enter the band and never acquire the error.

This is the structural reason the earlier essay’s reading fails here. A slope down to zero is a slope down to zero only if the answer lives at the bottom. When the answer lives on the far side, whatever the bottom costs in relative terms is billed to every state beyond it, and the far well is billed in full.

The norm calls every one of these answers exact #

The standard certificate for a computed solution is its backward error: how small a change to the matrix makes the computed vector exact. The exact answer to a nearby problem is the collection’s account of why that is the right question, and on almost every problem measured here it is.

On this one the normwise backward error, $\|r\|_\infty / (\|L\|_\infty\|\hat\pi\|_\infty + 1)$, is below $6.1\cdot 10^{-5}$ at every depth in every arithmetic — beneath half precision’s unit roundoff of $4.9\cdot 10^{-4}$, including every answer whose far well is exactly zero. By the norm, flush to zero’s empty far well is the exact stationary distribution of a chain whose rates differ from the true ones in the fifth digit.

The reason is arithmetic, not luck. The distribution is scaled so its largest entry is one; the residual row where a state was flushed is the product of a rate and a value below the smallest normal, so it is at most about $6\cdot 10^{-5}$; and $6\cdot 10^{-5}$ divided by a norm of order one is below the unit roundoff. Anything that happens below $u$ times the largest entry is invisible to a normwise measure, and the whole subnormal band of half precision is below that line. A normwise backward error cannot see underflow in a vector whose largest entry is near one, by construction.

The componentwise backward error, which divides each row’s residual by that row’s own $(|L||\hat\pi|)_i$, sees it at once: it is one — as large as it can be — from the first depth at which flush to zero empties the far well, and for gradual underflow it rises through the band and reaches one at 26 octaves, where its far well empties too. Small compared to what is the essay about this choice of denominator; here it is the whole difference between a certificate and a false one.

And the condition number blames the problem #

A second explanation offers itself for the empty far well, and it is the textbook’s default. The balance matrix is badly conditioned: its normwise condition number $\kappa_\infty(L)$ is $1.1\cdot 10^4$ at a depth of ten octaves, $3.7\cdot 10^6$ at twenty and $1.45\cdot 10^9$ at thirty, growing by a little under an octave per octave of barrier. The standard bound multiplies it by the unit roundoff, which passes one at seven octaves, so the bound says that past seven octaves half precision cannot be expected to get anything right.

It gets everything right. With the exponent unbounded, the same eleven-bit significand returns the far well to between $7\cdot 10^{-6}$ and $3.3\cdot 10^{-3}$ at every depth to forty, where the normwise bound has passed $10^8$. The componentwise condition number of every state’s probability — Skeel’s number read one component at a time — is exactly 95 at every depth: state $i$’s is $2i+1$, one contribution for each rate on the path to it, and the last state is state 47. A barrier forty octaves deep does not make any state’s probability harder to determine from the rates. It makes the matrix’s norm a bad summary of the rates, because the norm compares entries of $L^{-1}$ that differ by $2^{40}$.

Two condition numbers of one matrix found the same split on a scaled system: the normwise number is a worst case over perturbations that put their whole budget on the smallest entry, and a problem whose data are all positive and multiply rather than cancel is never that worst case. What loses the far well is not the conditioning and not the significand. It is the exponent field, and only the exponent field — which is why the unbounded route, whose one difference is the exponent, has no trouble. The condition number is an amplifier is right that κ multiplies a perturbation; the perturbation here is not one κ describes.

Rounding at the floor goes up as well as down #

“No subnormals” names two different arithmetics, and they fail differently. A flush-to-zero mode replaces any result below the smallest normal with zero. A format that simply lacks subnormals, and rounds every result to the nearest number it can represent, sends a value between half the smallest normal and the smallest normal up to the smallest normal.

On this recurrence the second one does something neither of the others does at that depth: it inflates the far well. At fourteen octaves it returns 0.431 against an exact 0.2745; at sixteen 0.722 against 0.245; at eighteen 0.899 against 0.218, 4.1 times too large; and from nineteen, zero. The mechanism is a sticking point. A state that should be 0.7 of the smallest normal rounds up to it; the next state is that value times a ratio a little under one, which is again above half the smallest normal and rounds up to it again; and the walk sits on the floor for several steps while the exact values keep falling, then climbs out of the barrier from a height it should never have had.

Gradual underflow does exactly the same thing one band lower. At 24 and 25 octaves the bottom of the barrier is near the smallest subnormal, a value rounds up to one unit of $2^{-24}$ and sticks, and the far well comes out 2.2 and 3.5 times too large before it empties. Rounding to nearest sticks at whatever floor a format has, as long as each step’s ratio is above one half. Only flush to zero never overshoots — it is biased strictly towards the answer being zero — and the figure shows its ratio pinned at one and then on the floor.

That changes what a comparison of two runs can say. The earlier essay measured a two-run detector — compare a run with subnormals against one without — and found it fires only while values are inside the band. On this chain the two runs differ at every depth from fourteen octaves to twenty-five, and the difference does not say which run is nearer the truth: at 25 octaves the run with subnormals is 3.5 times too large and the run without is zero, and the exact answer lies between them. From 26 the two agree again, both at zero, which is the earlier essay’s silent case.

Starting high buys more than the subnormals do #

The fourth arithmetic changes nothing about the format. It starts the recurrence at $\pi_0 = 1024$ instead of at one, using ten of the sixteen octaves half precision has above one and that the normalised recurrence leaves idle. It survives to 33 octaves, ten more than gradual underflow alone, and it stays within ten units of roundoff to 26.

The figure at the top of the page shows the three remedies side by side. Subnormals bought ten octaves; starting high bought ten more; an unbounded exponent — which in practice means holding the values near one by exact power-of-two rescaling, or working with logarithms of the probabilities — bought all forty and would buy any depth. The last one costs a test and a shift per step, or a logarithm, and requires knowing in advance that the computation needs it.

The other half of a format treated the exponent field as a budget to be spent, and this is the same budget from the other end: a recurrence normalised at its largest value spends none of the range above it, and half precision’s range above one is sixteen octaves out of forty. The caveat is on the same page as the gain — 48 states at up to 1024 each sum to at most 49,152, inside the format’s largest number of 65,504 but not by much, so the headroom above is spent by the sum rather than by the values, and a longer chain would have to start lower. A norm that overflows before it is a norm is that trap at the top, and the scaling that buys ten orders is the general form of the repair.

What the earlier reading got right #

None of this contradicts the earlier essay’s account of the band — it confirms it at the level of single numbers, where precision falls a bit an octave and the identity $x - y = 0 \Rightarrow x = y$ is bought with a slow case in the hardware. What it corrects is the inference from that account to computations. The inference assumed the small values were the answer. When they are a waypoint, the band’s slope is not a slope down to zero. It is a slope down and back up, and the far side inherits the bottom’s relative error in full — which gradual underflow keeps at a few per cent for ten octaves longer than flush to zero keeps it at all.

The detail worth carrying is the asymmetry of the certificates. The one quantity that would tell a user of the half-precision answer that something went wrong, the componentwise backward error, is the one that is almost never computed; the normwise one that is computed reports every answer on this page, including the empty far well, as exact. Eight bits and a format that breaks the rules priced a format whose smallest normal is $2^{-6}$; by the law above, the same chain would lose its far well there past a barrier of six octaves without subnormals and past about nine with them.

What one chain does not show #

One potential, one length, one set of rates. A chain whose path to the far well crosses several barriers pays the bottom’s error at each one; whether the costs add or compound is a measurement this sweep does not make. The recurrence is the simplest solver for a birth–death chain and not the only one; a general chain’s stationary vector is usually computed by an elimination, and a solver that renormalises as it goes would never enter the band. Half precision is the format measured, because it is the common format whose smallest normal a real computation reaches. Nothing here measures what subnormals cost in time, which is the reason flush to zero exists and which the earlier essay describes. And the survival depths are read at an error of one half in the far well; at a tighter tolerance every route survives less far, gradual underflow by more than the others, since its losses start sooner.

Still open: an elimination that keeps its sign, and formats with less band #

The elimination every Markov-chain code uses. The Grassmann–Taksar–Heyman variant of Gaussian elimination computes a general chain’s stationary vector without a single subtraction, which is why it keeps every component to high relative accuracy where ordinary elimination does not. That guarantee is stated for an arithmetic without an exponent limit. The prediction with a sign is that it fails exactly as the recurrence does — the band decides, not the elimination — and that its componentwise backward error, which it can report cheaply, says so.

A barrier crossed twice. Three wells in a row put two barriers on the path to the last one. If each barrier’s bottom costs its own relative error and the climb carries both, the errors should multiply, so a chain with two sixteen-octave barriers would lose more than twice what one costs. That is a two-parameter sweep of the same recurrence.

Eight-bit formats. E4M3 carries four significant bits and so a subnormal band of three octaves; E5M2 carries three and two. The law above predicts that on those formats the subnormals are worth almost nothing and starting high is worth almost everything, and both formats are used for exactly the probability-like quantities in which a recurrence of this kind appears.