3 ms·

Can't we just use Kahan summation [0] then divide by n? edit Someone on the reddit thread [1] had the same thought. [0] https://en.wikipedia.org/wiki/Kahan_su

Can't we just use Kahan summation [0] then divide by n? edit Someone on the reddit thread [1] had the same thought. [0] https://en.wikipedia.org/wiki/Kahan_summation_algorithm https://en.wikipedia.org/wiki/Kahan_summation_algorithm [1] https://old.reddit.com/r/perl/comments/2zneat/how_you_average_numbers_matters/ https://old.reddit.com/r/perl/comments/2zneat/how_you_averag...

m4r35n357 5y ago
Or can we sum more accurately by taking the stable_mean and multiplying by n?
SuchAnonMuchWow 5y ago
Kahan summation is better than naive summation, but not perfect. The stable_mean algorithm in the article is not perfect either. It just has more precision than the naive one, but it more or less has the same precision as twice-as-precise floats. >>> 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
- dataflow 5y ago
  Specifically https://code.activestate.com/recipes/393090/ https://code.activestate.com/recipes/393090/
- Majromax 5y ago
  In this case, it seems as if stable_mean is less accurate than the naive sum?
  - SuchAnonMuchWow 5y ago
    The two behaves differently, so you can not strictly compare the two. Most of the time stable_mean will give more accurate results, but this is not always the case. For example you can create test cases where the naive approach will return an exact result but stable_mean will have some rounding error. Kahan summation however is always more precise than naive summation.
- vitus 5y ago
  > Note that in python there is a math.fsum function which return a perfectly rounded summation of floats. 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 Namely, the "naive" mean yielded the smallest error (due to cancellation lining up just right such that it's capped after two rounds), fsum yielded a nominally larger one, and stable_mean yielded the biggest error. 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 edit: Naturally, if you want accurate decimal calculations, you should consider using decimal.Decimal or fractions.Fraction (albeit at the cost of significant performance). Using numpy.float128 will help, but it's still subject to rounding errors due to imprecise representation.
  - MauranKilom 5y ago
    Some of these results are indeed very surprising. The "correct" result of (double(0.3) - double(0.2) - double(0.1)) is -1 * 2^-54 = -5.5511151231257827e-17, and dividing by three comes out to -1.8503717077085942e-17. 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.
    - vitus 5y ago
      I think you might be off by one for the "correct" result (i.e. 2^-55), which is what fsum arrives at. >>> 0.3 - 0.2 - 0.1 -2.7755575615628914e-17 >>> (0.3 - 0.2 - 0.1) / 3 -9.25185853854297e-18 But either way, the fundamental problem in this case (as you noted) is actually due to loss of precision from numeric representation.
gpderetta 5y ago
Interestingly enough I was given this problem (well, numerically stable online variance, but close enough) a while ago as a take-home job interview and I ended up researching this quite a bit. 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.
- pixel_fcker 5y ago
  How did you force using the x87 fpu?
  - gpderetta 5y ago
    with GCC long double is always an 80 bit precision type.

m4r35n357 5y ago

Or can we sum more accurately by taking the stable_mean and multiplying by n?

SuchAnonMuchWow 5y ago

Kahan summation is better than naive summation, but not perfect. The stable_mean algorithm in the article is not perfect either. It just has more precision than the naive one, but it more or less has the same precision as twice-as-precise floats. >>> 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

dataflow 5y ago
Specifically https://code.activestate.com/recipes/393090/ https://code.activestate.com/recipes/393090/
Majromax 5y ago
In this case, it seems as if stable_mean is less accurate than the naive sum?
- SuchAnonMuchWow 5y ago
  The two behaves differently, so you can not strictly compare the two. Most of the time stable_mean will give more accurate results, but this is not always the case. For example you can create test cases where the naive approach will return an exact result but stable_mean will have some rounding error. Kahan summation however is always more precise than naive summation.
vitus 5y ago
> Note that in python there is a math.fsum function which return a perfectly rounded summation of floats. 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 Namely, the "naive" mean yielded the smallest error (due to cancellation lining up just right such that it's capped after two rounds), fsum yielded a nominally larger one, and stable_mean yielded the biggest error. 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 edit: Naturally, if you want accurate decimal calculations, you should consider using decimal.Decimal or fractions.Fraction (albeit at the cost of significant performance). Using numpy.float128 will help, but it's still subject to rounding errors due to imprecise representation.
- MauranKilom 5y ago
  Some of these results are indeed very surprising. The "correct" result of (double(0.3) - double(0.2) - double(0.1)) is -1 * 2^-54 = -5.5511151231257827e-17, and dividing by three comes out to -1.8503717077085942e-17. 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.
  - vitus 5y ago
    I think you might be off by one for the "correct" result (i.e. 2^-55), which is what fsum arrives at. >>> 0.3 - 0.2 - 0.1 -2.7755575615628914e-17 >>> (0.3 - 0.2 - 0.1) / 3 -9.25185853854297e-18 But either way, the fundamental problem in this case (as you noted) is actually due to loss of precision from numeric representation.

gpderetta 5y ago

Interestingly enough I was given this problem (well, numerically stable online variance, but close enough) a while ago as a take-home job interview and I ended up researching this quite a bit. 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.

pixel_fcker 5y ago
How did you force using the x87 fpu?
- gpderetta 5y ago
  with GCC long double is always an 80 bit precision type.