How to compute the average of doubles, so that the total error is minimal?

Question

Assume we have a long array of doubles, say, N == 1000000.

array<double, N> arr;

There are two naive approaches to compute the average. First

double result = 0;
for (double x : arr) {
    result += x;
}
result /= arr.size();

This may be inaccurate when the sum of values is very big. Floating point numbers lose precision then.

Another approach is:

double result = 0;
for (double x : arr) {
    result += x / arr.size();
}

This may lose precision when the numbers are small.

Is there any fail-safe way to calculate a simple average of floating point numbers? Solutions, which use only the standard library are appreciated.

Answer 1

If you want to squeeze more accuracy out of doubles, you can use Kahan summation and finally division by number of elements. There is however no standard library implementation of Kahan summation I know of.

An easy, standard way (almost like cheating) would of course be calculation using long doubles, basically using your first implementation and only converting the result back to double precision.

Answer 2

The so-called naive ways are not naive. What do the data mean, and how accurately can you measure those values? Unless the answer is something very unusual, the simple method with doubles is fine. However floats are a bit under-powered for general use.

If you add the small absolute values first you might get an extra bit or so of precision. That requires a sort. If the data are all above a certain threshold, subtracting the minimum may also give you another bit.

You can also store a partial total, and a partial mean, and check at each stage that partial mean * number processed is within a certain tolerance of the partial total. That won't give you any extra accuracy, but it will tell you if the fpu is too inaccurate for your purposes.

You can also use long double, or even code your own extended-precision floating point library (or use someone else's). However the solutions get increasingly heroic.

Answer 3

One way to reduce loss of precision would be to sort the doubles and then add them together in sorted order, starting with the smallest values and then at the end divide the final sum by the number of doubles.

So the tools you need would be std::sort and std::accumulate and plain old division /.