A norm that overflows before it is a norm

This site has a thread for pairs of expressions that a textbook derivation cannot tell apart and a computer can. Classical against modified Gram–Schmidt. The normal equations against QR. Elimination with and without a row swap. In every case the algebra is identical, the arithmetic is not, and the difference is a loss of digits — a number that comes back with four correct figures instead of fifteen.

This is the same thread with a different failure at the end of it. The two expressions here differ in whether they return a number at all, and they differ over half the format’s range.

The two expressions #

The Euclidean norm of a vector is √(Σxᵢ²), and written that way it is one line:

s = 0
for each x:  s = s + x*x
return sqrt(s)

The alternative divides by the largest entry first:

m = max |x|;  if m = 0 return 0
s = 0
for each x:  t = x/m;  s = s + t*t
return m * sqrt(s)

They compute the same thing. Both are backward stable, both make the same number of rounding errors to within a small constant, and the second costs one extra pass over the data and one division per entry.

The second is what hypot does, what every BLAS nrm2 does, and what almost nobody writes.

Why squaring is expensive in a way that has nothing to do with cost #

Multiplication adds exponents. So xᵢ² has twice the exponent of xᵢ, and a value comfortably inside the top of the range has a square that is not.

fp16’s largest finite number is 65,504, so anything above 256 has a square that overflows. The Euclidean norm of a vector whose entries are all 1,000 is 4,000 — representable, with three significant figures to spare — and the sum of squares is 1.6·10⁷, which is not a number fp16 has.

Half the range is gone, in the logarithm, and it is gone from the top. The measurement says how much:

Every one of the format’s normal exponents is tried: the same unit vector is rescaled by every power of two the format holds, both expressions are evaluated, and the answer is compared against the true norm. The scaled expression is right at every normal exponent in all three formats. The naive one is right over between 52% and 60% of them.

The 60% for fp16 rather than 50% is the vector length doing its work at the small end: with n = 8 the sum of squares of tiny entries is eight times a single square, which lifts the bottom of the naive window slightly above where the squares alone would put it.

The pairing the site could not draw before #

The comparison this whole extension exists for takes one line and it is a vector of sixteen entries, each 1,000.

fp16: the naive expression returns infinity. tf32: the naive expression returns 4,000, which is right. bfloat16: the naive expression returns 3,968, which is right to within its precision.

fp16 and tf32 have the same eleven significand bits. Every measurement this site made before the maturity phase gave them identical answers, because the simulation had no exponent. Here one of them returns the answer and the other returns an infinity, on identical data at identical precision.

And the third line is the one worth carrying away. bfloat16 is the least precise of the three and the one that succeeds where the most precise of them fails. It returns 3,968 rather than 4,000 — an error of 0.8%, which is exactly what eight significand bits are worth — while fp16 returns nothing usable at all. That is the entire argument for bfloat16 in three numbers, and it cannot be made on an axis that measures only precision.

The site asserts all three: infinite in fp16, finite and correct in tf32, and finite in bfloat16 with a larger error than tf32’s — the last one in that direction, because “the less precise format did better” is a sentence that would be false if read without the qualifier.

Where this bites #

The Euclidean norm is not an unusual computation. It is what every iterative method computes at every step, to decide whether to stop.

Conjugate gradients compute ‖r‖ each iteration. GMRES computes it as part of the Arnoldi orthogonalisation. Every stopping test in the rate the condition number predicts and the spectrum that predicts nothing is a norm ratio. Householder QR computes a norm to build every reflector — which is a reflection cannot stop being one’s subject, and where a spurious infinity would produce a reflector that is not a reflector.

At double precision none of this matters and that is why it is invisible. binary64’s range is 616 decades and half of it is 308, which no ordinary computation approaches. The problem appears exactly when the format narrows, which is exactly what the hardware has spent the last decade doing.

It is also why the standard library functions exist. hypot(x, y) is in every language’s maths library and it is there for one reason: √(x² + y²) written directly overflows when the answer does not. The scaled norm is that function generalised to n entries, and the reason nobody writes it by hand is that at double precision the naive version works and the habit never forms.

What the scaled version costs, and where it is not free #

The extra pass over the data is not nothing on a large vector, and the division per entry is worse than the multiply it accompanies. On a memory-bound computation the second pass roughly doubles the cost of the norm.

Which is why real implementations do not do exactly what the code above does. The standard approach is a single pass with a running scale: carry both a scale factor and a running sum, and rescale only when an entry arrives that is larger than the current scale. On ordinary data the rescale happens once or twice and the cost is close to the naive loop’s; on adversarial data it happens more often and the routine is still correct. That is a genuinely better algorithm than either of the two above and it is not drawn here, because the argument this essay is making is about the two-expression comparison and a third expression that is better than both would obscure it.

There is a second approach worth naming because it is what the fast paths actually do. Compute the naive sum, and check the exception flag afterwards: if nothing overflowed, the answer is right and was obtained at full speed; if something did, fall back to the scaled routine. IEEE 754’s sticky flags exist for exactly this pattern. It costs one branch per call rather than one branch per entry, which is the difference between a technique and a tax.

Neither of those is available in this site’s simulation, since a JavaScript number has no flags to check. What is available is the measurement of which expression works where, and that is the claim the figures make.

Two failures that are not the same failure #

It is worth being precise about why overflow gets its own essay when cancellation takes the answer already covers “an expression that is bad in a way the algebra hides”.

Cancellation is a relative error and it is bounded. Subtracting two nearly equal numbers returns a result whose relative error is large and whose absolute error is not: the answer is still within an ulp of the difference of the stored values, and the trouble is that those stored values were already rounded. Every backward error analysis in this site’s error field rests on that boundedness.

Overflow is not an error of any size. The result is an infinity, so the relative error is infinite, and infinity is absorbing: a subsequent subtraction gives a NaN and the NaN propagates through everything downstream. A single overflow does not degrade the answer, it removes it.

The practical consequence is the one worth having. There is no way to compute more accurately than the format permits — the response to cancellation is a better-conditioned expression or a longer significand. But overflow can always be engineered around, by scaling, at a cost of one extra pass or one exception check. A failure that has a fix and a failure that does not are different kinds, and the fix is three lines.

The scaling that is not a scaling #

There is a variant of this problem which the two expressions above do not cover and which is worth naming, because it is where the fix stops being three lines.

The scaled norm works because every entry can be divided by the largest one without loss: the ratios are at most one, so nothing overflows, and they are at least the smallest ratio in the data, so nothing underflows unless the data already spans the whole range. That argument needs the entries to be comparable. If a vector holds one entry at 10³⁵ and one at 10⁻³⁵, dividing by the first sends the second to 10⁻⁷⁰, which is not a binary32 number, and the scaled expression returns the larger entry alone.

That is the right answer to fifteen digits — a Euclidean norm is dominated by its largest entry when the spread is that wide — so nothing is lost here. It stops being the right answer in the computations where the small entries are the point, which is most of what this site is about: the residual of a converged iteration, the trailing singular values, the entries below a rank cutoff.

The general fix is to carry more than one accumulator — one for the entries near the top of the range, one for the middle, one for the bottom — and combine them at the end. That is what the carefully written library routines do, and it is why a well-implemented nrm2 is a hundred lines rather than five. The measurement here does not reach it, because a unit vector rescaled uniformly has no spread at all, and building a test with the spread would be building a test of a routine this site does not implement.

What the site does record is where its own scaled expression stops being right: four exponents at the bottom of each format, all of them subnormal, where the entries have already lost significand bits before any expression is applied. Reporting those rather than restricting the measurement to the range where the claim holds is the difference between a check and an illustration.

What is asserted here #

The scaled expression is right at every normal exponent of fp16, bfloat16 and binary32 — 30, 254 and 254 of them, with no exceptions.

The naive one covers between 45% and 62% of them, asserted in both directions. The upper bound is the claim; the lower one exists because an expression that covered a quarter of the range would be a different and much worse expression, and asserting only the upper bound would not notice.

And the scaled expression loses accuracy only among the subnormals, at most four exponents in each format, where the entries themselves have already lost significand bits before any norm is taken. That is reported rather than excluded — a routine described as covering “the whole range” and measured only over the normal part is a claim measured where it holds.

fp16 returns infinity where tf32 returns 4,000, at the same precision, and bfloat16 returns a finite answer less accurate than tf32’s.

The refusals: that fp16 and tf32 compute the same thing, which is the claim this site was making by omission for two phases; and that overflow is a rounding error near the answer, fed 300² in fp16.

The habit this suggests, and the one it does not #

The advice that follows is narrower than “always use the scaled version”, and worth stating carefully because the broad version is wrong.

At double precision, write whichever is clearer. Six hundred decades of range is more than any physical quantity spans, the naive expression covers three hundred of them, and a computation that leaves that window has a scaling problem the norm is not going to fix. Rewriting every norm in a double-precision code for this reason would be work spent on a failure that is not going to happen.

At sixteen bits, do not write either — call the library. Twelve decades of range with six usable by the naive expression is a window narrow enough that ordinary data leaves it, and the correct routine is more subtle than the scaled loop above. This is the same conclusion the road that squares the problem reaches about the normal equations, and for the same reason: an expression that is fine in one regime and disastrous in another is exactly the kind of thing a library exists to get right once.

And in between, know which regime the data is in. That is not an evasion. The question “what is the largest magnitude this computation will produce” has an answer for most real codes, it is usually knowable in advance, and knowing it settles the question completely.

What the essay does not support is the stronger claim that the naive expression is a mistake. It is the right expression at double precision, it is faster, and it is clearer. What makes it a mistake is writing it into a routine that will later be compiled at a narrower precision by somebody who was not thinking about its range — which is a description of most of what has happened to numerical software in the last five years, and the reason this essay is in the maturity phase rather than the foundation.