How to average floating point numbers (2015)
nu42.com
nu42.com
>>> from interval import interval
>>> data = (interval(1_000_000_000.1), interval(1.1)) * 50_000
>>> sum(data) / len(data)
interval([500000000.5971994, 500000000.601347])
The result is definitely between the 2 values. If you try another method to do the calculation, (e.g. sorting first) you will know which one was better without the answer having to be obvious as it is in this case.
You might find one method is better at the lower bound, the other is better at the higher bound, which means you can narrow the value further than using either individually. It's a really powerful system that I think should be much better known.
edit Someone on the reddit thread [1] had the same thought.
[0] https://en.wikipedia.org/wiki/Kahan_summation_algorithm
[1] https://old.reddit.com/r/perl/comments/2zneat/how_you_averag...
>>> data = (1_000_000.1, 1.1, -1_000_000.1, -1.1) * 25_000
>>> sum(data) / len(data) == -2.3283153183228934e-16
>>> sum(sorted(data)) / len(data) == 2.7122441679239272e-11
>>> stable_mean(data) == -2.571225228716577e-13
>>> float(sum(numpy.array(data, dtype=numpy.float128)) / len(data)) == 2.2648549702353195e-19
Note that in python there is a math.fsum function which return a perfectly rounded summation of floats. There are multiple different techniques for that, but they are more complex than Kahan summation.
>>> math.fsum(data) / len(data) == 0.0
Kahan summation however is always more precise than naive summation.
Quirks of the underlying data can yield surprising results.
I was playing around with the "well-known" floating point error result of 0.3 - 0.2 - 0.1 = some small number (instead of 0), and found:
>>> sum((0.3, -0.2, -0.1) * 1_000_000) / 3e6
-2.775557561562891e-23
>>> math.fsum((0.3, -0.2, -0.1) * 1_000_000) / 3e6
-9.25185853854297e-18
>>> stable_mean((0.3, -0.2, -0.1) * 1_000_000)
-1.387894586067203e-17
Of course, if you sort the values (such that we're accumulating -0.2 + ... -0.2 + -0.1 + ... + -0.1 + 0.3 + ... + 0.3, then sum performs much worse, stable_mean is a little worse, and fsum performs exactly the same.
>>> sum(sorted((0.3, -0.2, -0.1) * 1_000_000)) / 3e6
-1.042215244884126e-12
>>> math.fsum(sorted((0.3, -0.2, -0.1) * 1_000_000)) / 3e6
-9.25185853854297e-18
>>> stable_mean(sorted((0.3, -0.2, -0.1) * 1_000_000))
-3.0046468865706775e-16
Hence, as far as those functions are concerned, the naive mean yielded the worst error in both tests, but stable_mean coincidentally happens to be the best in the first test. I'm slightly surprised that none of them got the correct result.
>>> 0.3 - 0.2 - 0.1
-2.7755575615628914e-17
>>> (0.3 - 0.2 - 0.1) / 3
-9.25185853854297e-18
It has been a while, but if I remember correctly pairwise sum was in practice as good or better as Kahan summation (and both more stable than Walford's) but faster.
But by far the fastest and still very stable solution was to use the legacy 80 bit x87 fpu as accumulator. Which is not surprising as Kahan explicitly added 80 bit precision to x87 for intermediate computations.
>>> data = (1_000_000_000.1, 1.1) * 50_000
>>> sum(data) / len(data)
500000000.60091573
>>> sum(sorted(data)) / len(data)
500000000.6004578
>>> sum(sorted(data, reverse=True)) / len(data)
500000000.6012391
>>> import statistics
>>> statistics.mean(data)
500000000.6
>>> X = [10_000_000_000 + k for k in (4, 7, 13, 16)]
>>> statistics.variance(X)
30 >> RUBY_VERSION
⇒ "2.7.2"
>> data = [1_000_000_000.1, 1.1] * 50_000 ; nil
⇒ nil
>> data.sum
⇒ 50000000060000.0
>> data.sum / data.length
⇒ 500000000.6
>> data.sort.sum / data.length
⇒ 500000000.6
>> data.sort.reverse.sum / data.length
⇒ 500000000.6
We don’t do that in TruffleRuby and so don’t produce quite such an accurate result.
Here's the internal _sum function for example, which you can see maps _exact_ratio over the data before summing: https://github.com/python/cpython/blob/3.9/Lib/statistics.py...
>>> data = (1_000_000_000.1, 1.1) * 50_000
>>> statistics.fmean(data)
500000000.6
>>> import boltons.statsutils
>>> data = (1_000_000_000.1, 1.1) * 50_000
>>> boltons.statsutils.mean(data)
500000000.60091573
>>> sum(data) / len(data)
500000000.60091573The algorithm uses 2n bits for the sum. https://radfordneal.wordpress.com/2015/05/21/exact-computati...
My issue is that I've come across errors due to floating point in weird places. One example, searching a list of strings for another string. If the sort function decides those strings look like numbers it may coerce them into doubles (ints being too small to hold them). Then sort wont match correctly due to the least significant bits being rounded off.
The fact is: such kind of SNAFU does happen.
You sometimes even see that with APIs inconsistently wrapping behind quotes or not (seemingly at random) some values. And then it gets trickier because you can have the parser correctly parse into, say, "big decimals" some of the values, while others may be automagically, imprecisely, parsed into floating-point numbers.
The problem is of course made worse by people overall overusing decimal numbers where integers would be totally fine.
One of the most beautiful floating-point bugs ever is this one (from about ten years ago as I recall it): language preference in browsers could be set as a decimal number between 0 and 1. Java (and other languages) had a broken floating-point implementation where parsing a certain floating-point number would lead to an infinite loop. So you could set one core of any webserver with that bug by simply requesting a page, with the language preference set to that floating-point number. Rinse and repeat a few times and all cores were spinning at 100%, effectively DDoS'ing the webserver.
After all of this, I still think that IEEE 754 is pretty darn good. The fact that it is available and fast in pretty much all hardware is a major selling point.
That said, I absolutely hate -0 (minus zero). It's a stupid design wart that I've never seen any good use case for. But it is observable and therefore cannot be ignored, leading to even more icky special cases.
E.g. I've seen it argued that +0 and -0 represent limits of functions that approach zero from the negative side or the positive side. But that makes no sense, because no other floating point number represents a limit going from above or below the number, so why should 0? Also, all this walking on eggshells around a possible -0 (because that represents the limit of 0 coming from below 0) is just lost when someone writes 0. Admittedly, I've not seen the guts of most numerical code, but the only time I have ever seen a minus 0 written into code was to test the minus-zero handling of something.
/shrug
Just use a programming language that doesn't coerce strings into numbers. Or use a comparison function/operator that explicitely compares arguments as strings and not as numbers.
Like what? Packing the real number line in 32 (or 64) bits is an incredible trick. The roughly scale-invariant distribution of represented values means you can do math without regard to choice of units, which would be impossible with any fixed point scheme.
Floating point is overall pretty good, sometimes I wish there was a version without some many ways to store NAN(for 16 bit floats it is so wasteful), and maybe versions of sqrt/div/rcp/rsqrt that would clamp invalid inputs to 0, but other than that it is hard to see how you make any major improvements given the fixed # of bits.
Actually one other possible improvement would be a way to request the hardware to dither when rounding down to the fixed bit count, although I'm not sure how it could do anything besides pure white noise, but better than nothing..
[Ed. smaller meaning closer to zero]
Except this doesn't work if you must handle negative numbers. So what matters is not the smallest numbers but the number nearer from zero.
Summing [4, 5, 6, 7] would yield the series (4, )9, 15, 22 with "sort then add", but what the previous poster described would first sum 4+5 (leaving the list [9, 6, 7]), then 6+7 (leaving the list [9, 13]) and finally arriving at 22.
That algorithm could be improved (for accuracy, not speed) by replacing sums of the smallest number with multiplications as I suppose that multiplication of more than 3 identical numbers or more is loosing less precision than multiple additions.
For example:
[4, 3, 7, 7, 8] could be calculated as:
(3+4) => [7, 7, 7, 8]
(7*3) => [21, 8]
(8+21) => 29If you do, in fact, need to worry about floating point error, then there is not some magical "it will never happen if you have this much precision and use this neat algorithm". Heck, even fixed/decimal/rational numbers aren't perfect solutions, they just make it easier to tell if you're staying within your precision.
Without infinite storage, you can't retain infinite precision, so the only real question is how much error you are willing to tolerate. And that depends very much on what you're doing - I can tolerate a lot more error in rolls on a MMO than in my bank account balance...
So, you are, in fact, gonna need an algorithm for summation that controls that error, if you're not willing to tolerate very much.
YAGNI is a prediction, not a rule. Just like premature optimization, alternatives to floating point - or using longer floating point types - have real costs, which must be balanced against the needed precision.
As the OP says: "One saving grace of the real world is the fact that a given variable is unlikely to contain values with such an extreme range"
My point is that 99.999% algorithm is going to have a cost over the 99.9% algorithm (more complex code, more resource usage), and depending on what you're doing, you may have to consider that cost.
double fsum(const double *p, size_t n) {
size_t i;
double s;
for (s = i = 0; i < n; ++i) s += p[i];
return s;
}
double fsum(const double *p, size_t n) {
size_t i;
double s;
if (n > 8) return fsum(p, n / 2) + fsum(p + n / 2, n - n / 2);
for (s = i = 0; i < n; ++i) s += p[i];
return s;
}
What we need to consider is that sometimes there is no tradeoff.
If you are averaging bazillions of numbers and they range from weeny to huge and you thing 50_000_000.000006 != 50_000_000.06 for your purposes then it is hard for me to see that you are doing anything sensible.
This is, except in very specialised cases that I cannot imagine, a non issue.
Just add them up and divide, (and check for overflow, that could ruin your day!)
Even the naive floating implementation can't possibly get it this wrong for just 4 numbers.
>>> import numpy as np
>>> arr = np.array( (1_000_000_000.1, 1.1) * 50_000)
>>> arr.mean()
500000000.6000003Better: recursive bucket accumulators using the exponent
Clarification: I do not say that treating floats careful is not important, and I know doing stuff like subtracting large numbers can lead to inaccurate leading digits. This is clear and part and parcel of any numerical work. But this is not what the article is about. I am claiming instead the article is looking for a problem where there is none (or let's say I fail to see) .
You wouldn't want those kind of rounding errors in financial applications, in rocket science applications, etc ...
Many of us have had their attitudes readjusted by a bug where these "insignificant" errors became anything but.
We're often wrong in ways we can't imagine, so it is wise to remain humble and try to strive for correctness over our potentially misguided ideas how things work.
For example, you'd better still have multiple layers of security even though we've already made some other layer "impenetrable". And still check for those "impossible" error cases before they create nearly undebuggable issues.
From TFA:
"Now, in the real world, you have programs that ingest untold amounts of data. They sum numbers, divide them, multiply them, do unspeakable things to them in the name of “big data”. Very few of the people who consider themselves C++ wizards, or F# philosophers, or C# ninjas actually know that one needs to pay attention to how you torture the data. Otherwise, by the time you add, divide, multiply, subtract, and raise to the nth power you might be reporting mush and not data.
One saving grace of the real world is the fact that a given variable is unlikely to contain values with such an extreme range. On the other hand, in the real world, one hardly ever works with just a single variable, and one can hardly every verify the results of individual summations independently.
Anyway, the point of this post was not to make a grand statement, but to illustrate with a simple example that when using floating point, the way numbers are added up, and divided, matters."
Yes, real data is noisy but testing needs do be precise and repeatable. For example, if we need to test that a value is <=100, then it must pass at 100 and fail at 100.00000001. And yes, we use margins and fixed point numbers too but sometimes, we need to be precise with floating point numbers.
It also matters in some calculations, for example when doing GCD/LCM with periods and frequencies. For that, we eventually switched to rational numbers because precision of source data was all over the place.
We all know that parts per billion rarely make sense but where do we draw the line? Sometimes, it really matters (ex: GPS clocks), sometimes, PI=3 is fine. So in doubt, use the highest precision possible, you can relax later if you can characterize your error.
And even if you use these techniques butterfly effect will kick in eventually.
Imagine three redundant systems on a plane that want to compute the same from the same input array (and then check whether they all agree). You want to implement this, and figure that (since there's a lot of data) it makes sense to parallelize the sum. Each thread sums some part and then in the end you accumulate the per-thread results. No race conditions, clean parallelization. And suddenly, warning lights go off, because the redundant systems computed different results. Why? Different thread team assignments may lead to different summation orders and different results, because floating point addition is not associative.
More generally, whenever you want fully reproducible results (like "the bits hash to the same value"), you need to take care of these kinds of "irrelevant" problems.
Real image data, roughly 1 megapixel per tile at that time, and 32-bit floats.