A square that evaluates negative

The two essays before this one end on verdicts that two builds disagree about, and in both cases there is a right answer: the determinant is 1, the matrix is definite or it is not, and an exact computation settles it. This one has no right answer to reach for, and that is the point of it.

Two traces of noise. The polynomial is (x − 1)⁶ expanded, so its value is a sixth power and cannot be negative; within 10⁻³ of the root it is under 10⁻¹⁸, which is below the rounding of the terms that produced it. One build returns a negative value at 179 of the 401 points and the other at 196, and the two disagree about the sign at 98 — a quarter of the window.

Why this is a different failure #

In one multiply the compiler removed the fused form is exactly right and the other is exactly wrong, because the algorithm was designed around the instruction. In a matrix that is definite on one machine neither form is reliably right, but one of them is right on each matrix and an exact computation says which.

Here neither trace is nearer the truth in any useful sense. The true value is 10⁻¹⁸ and both evaluations return quantities around 10⁻¹⁵, so both are wrong by three orders of magnitude, and asking which is closer to zero is asking which random number is smaller. The largest difference between the two evaluations across the window is 1.44·10⁻¹⁵ — six roundings, exactly what the bound allows — and both are entirely made of them.

What the picture shows is therefore not an accuracy comparison. It is the shape of the region in which a computed sign belongs to the build, and that region has a width, which is computable.

Where the window’s edge is #

The window is where |p(x)| falls below the rounding of the terms that make it up. Horner’s rule evaluating a degree-d polynomial commits about d roundings, each of size u times the largest partial value it holds, so the noise floor is roughly

u · max over the recurrence of |sₖ| ≈ u · Σ|cᵢ| |x|ⁱ,

which near x = 1 is u times the sum of the coefficients’ magnitudes: 1 + 6 + 15 + 20 + 15 + 6 + 1 = 64, giving about 1.4·10⁻¹⁴. The true value reaches that at |x − 1| = (1.4·10⁻¹⁴)^(1/6) ≈ 4.9·10⁻³.

Which is what the drag measures: at a window of 10⁻¹ the polynomial is 10⁻⁶ at the ends, both traces separate from zero and agree everywhere except in the middle; at 10⁻³ they are noise across the whole picture. The width is a property of the polynomial and the precision and is available before any evaluation is done.

That number — Σ|cᵢ| divided by |p(x)| — is exactly the summation condition number of the vector that hides it, applied to the terms of a polynomial. The essay’s census puts a residual at 10¹⁷ and this evaluation is the same object: an expression whose terms are enormous relative to their sum.

What a root finder does with it #

The consequence is not that a value is wrong. It is that a branch is taken differently, and the branch is inside every algorithm that locates something by sign.

A bisection asks for the sign at the midpoint and discards half the interval. Inside the window that sign is a property of the build, so two builds discard different halves, converge to different points, and return roots that differ by the width of the window. Their residuals are the same size — both are 10⁻¹⁵, which is what convergence looks like here — and neither is wrong by the standards a root finder applies to itself.

This is the same mechanism the roots are not the coefficients measures for a root’s conditioning: a double root moves like √ε, a root of multiplicity m like ε^(1/m), so half the digits are lost at a double root and five sixths at a sixth-order one, with nothing ill conditioned anywhere. The window here is that fact seen from the other side. The conditioning says the root cannot be located better than 10⁻³ from data at the rounding level; the window says that inside that 10⁻³ the arithmetic’s own decisions are what a bisection is following.

The two agree, and they are different statements. The conditioning is a property of the problem and it is the same for both builds. The window is a property of the arithmetic, and inside it the two builds take different paths through a computation that both complete successfully.

Counting the disagreements is the measurement #

The number worth quoting from the hero figure is not 179 or 196 — those are counts of a wrong sign, which is interesting once. It is 98, the number of points at which the two builds disagree with each other, because that is the quantity a person maintaining two machines actually meets.

It is worth seeing why it is neither 0 nor 401. It is not zero because the two evaluations differ by up to six roundings and the values are the size of one, so a disagreement is likely wherever the noise happens to straddle zero. It is not 401 because the two traces are strongly correlated: they share every input and differ only in how each of six multiply-adds was rounded, so most of the time they land on the same side of zero by accident of sharing most of their arithmetic.

A quarter is therefore the honest expectation, and it is stable across the window: at 10⁻⁴ and 10⁻⁵ the fraction is similar, and at 10⁻¹ it collapses to nearly nothing because the true value has come back above the noise. That stability is what makes the width of the window — rather than the count — the thing to compute in advance.

The general shape, since three essays have now produced it #

This is the third instance of one pattern and it is worth writing down before the field moves on to reductions again.

A computation produces a real number. The number is a difference of quantities much larger than itself, so its leading digits are rounding. A consumer compares it with something — zero, a tolerance, another number — and turns it into a discrete outcome. Two conforming builds compute the rounding differently, so the outcome differs, and every diagnostic anybody applies afterwards reports that both builds behaved correctly.

The three instances differ only in what is being compared and whether there is an adjudicator:

The third row is the one that closes the sequence, because it is where the site’s usual move stops working. This collection’s habit is to answer a disagreement by computing the quantity exactly and saying who was right. Here the exact value is 10⁻¹⁸, it is not zero, and it is positive — so the “right” answer is that every negative value in the figure is wrong, which is true and useless, because no build was ever going to return it and no algorithm that consumes the sign can act on it.

What the exact value does establish is where the window is, which is the thing to compute.

Three things that do not fix it #

Higher precision. It moves the edge: quadruple precision would put the noise floor at 10⁻³² and the window at |x − 1| ≈ 5·10⁻⁶. The window is narrower and it is still there, and inside it the same thing happens. This is the shape the whole field keeps producing — precision moves a boundary and does not remove one.

A tolerance. Declaring |p| < 10⁻¹⁴ to be zero converts a sign into three outcomes, positive, negative and cannot tell, which is genuinely better — but a bisection cannot use cannot tell without an answer for what to do next, and the honest answer is to stop and report an interval. Most implementations do not have that outcome and the ones that do return an interval of width 10⁻³, which is the conditioning again.

Turning off the contraction. It makes the trace reproducible between builds that share the flag, and it does not make it correct, and it makes it the less accurate of the two forms elsewhere. It is the option that most looks like a fix and buys the least.

What does work #

Two things, and both are changes of formulation rather than of arithmetic.

Do not expand the polynomial. (x − 1)⁶ evaluated as a sixth power of (x − 1) is exact to a rounding at every point, has no cancellation anywhere, and is never negative. The whole failure is a consequence of writing a polynomial in a basis in which its value is a difference of large terms — which is what a monomial expansion does to any polynomial with a root away from zero, and what the roots are not the coefficients is about when it measures a Chebyshev basis against a monomial one.

Or evaluate with an error-free transformation. Compensated Horner — the same twoProduct and twoSum machinery the reproducible summation policies use — evaluates a polynomial as though in twice the precision, at about four times the cost, and pushes the noise floor to u² instead of u. That does not remove the window either; it narrows it to |x − 1| ≈ 4·10⁻⁶, which for most purposes is narrow enough that the conditioning is the binding constraint rather than the arithmetic.

Both are the same recommendation this site keeps arriving at from different directions: when a computed quantity is a difference of large terms, the repair is upstream of the arithmetic.

The window is not a defect of Horner’s rule #

Worth heading off, because the obvious response to the figure is to blame the algorithm.

Horner’s rule is the best of the standard evaluations by every measure the subject applies. It uses d multiply-adds where the naive sum of powers uses about 2d operations, its backward error is optimal — the computed value is the exact value of a polynomial whose coefficients differ from the given ones by a relative γ₂d — and it is what every library uses. The failure here is not a failure of Horner’s rule; it is the failure the backward error statement predicts.

A polynomial whose coefficients are perturbed by 10⁻¹⁵ relative has roots that move by (10⁻¹⁵)^(1/6) ≈ 3·10⁻³ for a sixth-order root, and inside that neighbourhood the perturbed polynomial takes both signs. Horner returns the exact value of one such polynomial, and the fused build returns the exact value of a different one. Both are correct, both are answers to slightly different questions, and the questions differ by less than the data can distinguish.

That is the site’s founding identity applied to an evaluation rather than to a solve, and it puts this essay’s window exactly where three errors and one number would predict: the conditioning of the root is enormous, the backward error of the evaluation is at the rounding level, and their product is the width of the region in which nothing can be said.

Where else the sign of a nearly-zero quantity is read #

Naming them, because a reader will recognise more of these than they expect.

A discriminant. b² − 4ac decides whether a quadratic’s roots are real, and near a double root it is the difference of two nearly equal products — exactly the determinant of one multiply the compiler removed with different letters.

A pivot. When symmetry is not enough chooses between a 1 × 1 and a 2 × 2 block by comparing computed quantities whose difference decides the branch.

A convergence test on a difference of iterates. |xₖ₊₁ − xₖ| is a difference of nearly equal things by construction, and its comparison with a tolerance is a tolerance that reads its own residual.

A determinant used to decide orientation, a curvature used to decide a step, a sign used to decide a reflection. Every one is a real number computed from cancelling terms and consumed as a bit.

The field’s general statement covers all of them: a stopping test is a race prices the third, and the pattern is that the arithmetic’s disagreement is small, the consumer is discrete, and the composition is a branch two machines take differently.

What to report instead of a sign #

The recommendation the three FMA essays share, stated once here because this is where it is unavoidable.

A routine that computes a quantity and consumes its sign should return the quantity and a scale, not the sign. For the polynomial: return p(x) and Σ|cᵢ||x|ⁱ, from which a caller computes the noise floor and decides whether the sign means anything. For a Cholesky: return the pivot and the matrix norm. For a determinant: return the value and the sum of the magnitudes of its two products.

In every case the second number costs one accumulation the loop was nearly doing anyway, and it converts an unanswerable question — is this sign real — into an arithmetic one. That is the same move rank is a decision makes for a rank: the honest output is the spectrum with the threshold drawn on it, and the count above the threshold is the caller’s.

The reason it is not standard is that a sign has an obvious return type and a scale does not have an obvious name. sgn is one bit and (value, scale) is a struct, and a routine that returns a struct where the mathematics says sign looks like it is hedging. It is not hedging: the sign genuinely is not determined by the inputs at this precision, and returning one is a claim the arithmetic cannot support.

At other settings #