"It fails for x-y when x and y are close"

I don't think that's quite right. Each individual operation is accurate within epsilon, the trouble with cancellation is that the error from operations prior to the subtraction is much larger than the result of the subtraction.

Example, if I do 1 + 2^30 - 2^30 in single precision I get zero. But each individual operation is correct. 1 + 2^30 = 2^30, so the error is 1 < 2^30 * epsilon. Then 2^30 - 2^30 = 0, which is exactly correct.

Yes. The big errors appear when you remember x is actually (1+ε)x, as explained towards the end of the article. If you assume the inputs are exact then you get an accurate result, but if you are doing error analysis that does not help you.

This is the right answer. Substraction is accurate, as much as it can be given the its inputs. It’s everything (combined) leading up to it that isn’t.

Every year I struggle to teach this simple concept in class :-) :-(

This is due to the non-associative nature of FP. (a+b) + c != a + (b+c).

I disagree that the "axiom" as stated is fundamental. My main argument with it is that I don't see an easy way to go from the axiom to usable theorems about FP.

Happy to be wrong on this, but I am missing this at this time.

I'd say not only is fundamental, it's pretty much the only tool you have to start the analysis of an FP algorithm. If you want to analyze a naive sum algorithm for instance you do that by recursively applying the a ⊕ b = (a + b)(1 + ε) rule and figure out how the epsilons bubble up - the result being that the error is linear with the sum length.

Yep the term for this is catastrophic cancellation, which doesn’t invalidate the axiom.

> Kahan’s summation algorithm

Radford Neal (2015) has an approach to computing sums exactly (to the correct rounded floating point number) in less than twice the time to do it naively. In some situations it can do work in parallel and saturate memory bandwidth, so takes roughly the same time: http://www.cs.toronto.edu/~radford/xsum.abstract.html

"Kahan’s (1965) method for reducing summation error (but without producing the exact result) was also tested on all systems. On modern 64-bit processors, computing the exact sum with the large superaccumulator method was faster than Kahan’s method for summations of more than about 1000 terms. Kahan’s method was significantly faster than the small superaccumulator method only for summations of less than about 100 terms."

I think you're misinterpretating the fundamental axiom because it's still true for `x-y`. What you're thinking is that the fundamental axiom applies to the composition of operations `closest_float(x) + closest_float(y)`, but the axiom only applies to each individual operation.

The axiom can be stated more rigorously "If f is a (possibly multivariate) function on the set of floats to the reals, the IEEE-compliant approximation of f must give the closest float for all float values of its arguments."

> "If f is a (possibly multivariate) function on the set of floats to the reals, the IEEE-compliant approximation of f must give the closest float for all float values of its arguments."

It should be noted that this is not exactly true for some IEEE functions, due to the table-maker's dilemma: computing the "closest float approximation" to an arbitrary function can be computationally infeasible if the true result is close to being halfway in-between two consecutive floats. (This applies to the default rounding mode, but can be extended by analogy to other rounding modes.)

They're correct that the "In mathematical terms" version of it is wrong (even for addition, if one is negative). It's missing fl()'s around the other a's and b's. Your statement is greatly superior than the poor man's version presented in the article.

Kahan summation is unfortunately slower than naive summation, which can be problematic. Pairwise summation [1], which roundoff error is proportional to log n, can be used as a drop-in replacement of naive summation as it performs more or less equally.

[1] https://en.wikipedia.org/wiki/Pairwise_summation

That's not true. The axiom is correct. Say you're using double (i.e., you have 16 decimal digits of accuracy, roughly) and you are doing ((1.0 + 1e-20) - 1.0).

That becomes (1.0 - 1.0) = 0.0 because 1e-20 was dropped during the addition. The substraction was perfectly accurate given its inputs (1.0 - 1.0 is exactly 0.0).

The issue is that float(1.0 + 1e-20) = 1.0. But that's also very accurate ! (All the 16 decimal digits of 1.0 are correct). But yeah the 20th (decimal) digit got dropped.

You see how each individual computation is maximally accurate. So the axiom is correct. What it boils down to is that you are manipulating (adding/subtracting) numbers of very different magnitudes. So cancellations, which would be negligible if all numbers where of similar magnitude, becomes catastrophic. Hence the name catastrophic cancellations.

The key is the framing; as stated, x and y are _exact_ floating-point values, which makes the statement true[1]. In fact, it's embarrassingly trivial, since it is just a weakening of the definition of rounding used by IEEE 754. In particular, if cancellation occurs in x - y, then the result is exact (Sterbenz' lemma).

[1] Under the assumption that no underflow or overflow occurs (which is the usual framing for this sort of naive error analysis). It's easy enough to correct for underflow by adding an absolute error term as well.

In other words, while we all know that addition and subtraction over real numbers is associative, addition and subtraction over floating point numbers is not.

> Subtraction is the easiest way to get big floating point errors.

That's why most modern floating point representations allow subnormal numbers: without them, subtracting two near, normal numbers will result in zero even if the numbers aren't equal.

> It fails for x-y when x and y are close. Subtraction is the easiest way to get big floating point errors.

And note: addition might be subtraction, while subtraction might be addition.

Consider x=1.0 and y=-0.99999.

x + y = 0.00001

You did addition, but you got cancellation error. The only solution is to sort the numbers from largest to smallest dynamically during runtime.