Numpy float mean calculation precision

Question

I happen to have a numpy array of floats:

a.dtype, a.shape
#(dtype('float64'), (32769,))

The values are:

a[0]
#3.699822718929953
all(a == a[0])
True

However:

a.mean()
3.6998227189299517

The mean is off by 15th and 16th figure.

Can anybody show how this difference is accumulated over 30K mean and if there is a way to avoid it?

In case it matters my OS is 64 bit.

Answer 1

This is due to incapability of float64 type to store the sum of your float numbers with correct precision. In order to get around this problem you need to use a larger data type of course^*. Numpy has a longdouble dtype that you can use in such cases:

In [23]: np.mean(a, dtype=np.longdouble)                                                                                                                                                                    
Out[23]: 3.6998227189299530693

Also, note:

In [25]: print(np.longdouble.__doc__)                                                                                                                                                                       
Extended-precision floating-point number type, compatible with C
    ``long double`` but not necessarily with IEEE 754 quadruple-precision.
    Character code: ``'g'``.
    Canonical name: ``np.longdouble``.
    Alias: ``np.longfloat``.
    Alias *on this platform*: ``np.float128``: 128-bit extended-precision floating-point number type.

_{* read the comments for more details.}

Answer 2

The mean is (by definition):

a.sum()/a.size

Unfortunately, adding all those values up and dividing accumulates floating point errors. They are usually around the magnitude of:

np.finfo(np.float).eps
Out[]: 2.220446049250313e-16

Yeah, e-16, about where you get them. You can make the error smaller by using higher-accuracy floats like float128 (if your system supports it) but they'll always accumulate whenever you're summing a large number of float together. If you truly want the identity, you'll have to hardcode it:

def mean_(arr):
    if np.all(arr == arr[0]):
        return arr[0]
    else:
        return arr.mean()

In practice, you never really want to use == between floats. Generally in numpy we use np.isclose or np.allclose to compare floats for exactly this reason. There are ways around it using other packages and leveraging arcane machine-level methods of calculating numbers to get (closer to) exact equality, but it's rarely worth the performance and clarity hit.

Answer 3

Here is a rough approximation of a bound on the maximum error. This will not be representative of average error, and it could be improved with more analysis.

Consider calculating a sum using floating-point arithmetic with round-to-nearest ties-to-even:

sum = 0;
for (i = 0; i < n; ++n)
    sum += a[i];

where each a[i] is in [0, m).

Let ULP(x) denote the unit of least precision in the floating-point number x. (For example, in the IEEE-754 binary64 format with 53-bit significands, if the largest power of 2 not greater than |x| is 2^p, then ULP(x) = 2^p−52. With round-to-nearest, the maximum error in any operation with result x is ½ULP(x).

If we neglect rounding errors, the maximum value of sum after i iterations is i•m. Therefore, a bound on the error in the addition in iteration i is ½ULP(i•m). (Actually zero for i=1, since that case adds to zero, which has no error, but we neglect that for this approximation.) Then the total of the bounds on all the additions is the sum of ½ULP(i•m) for i from 1 to n. This is approximately ½•n•(n+1)/2•ULP(m) = ¼•n•(n+1)•ULP(m). (This is an approximation because it moves i outside the ULP function, but ULP is a discontinuous function. It is “approximately linear,“ but there are jumps. Since the jumps are by factors of two, the approximation can be off by at most a factor of two.)

So, with 32,769 elements, we can say the total rounding error will be at most about ¼•32,769•32,770•ULP(m), about 2.7•10⁸ times the ULP of the maximum element value. The ULP is 2⁻⁵² times the greatest power of two not less than m, so that is about 2.7•10⁸•2⁻⁵² = 6•10⁻⁸ times m.

Of course, the likelihood that 32,768 sums (not 32,769 because the first necessarily has no error) all round in the same direction by chance is vanishingly small but I conjecture one might engineer a sequence of values that gets close to that.

An Experiment

Here is a chart of (in blue) the mean error over 10,000 samples of summing arrays with sizes 100 to 32,800 by 100s and elements drawn randomly from a uniform distribution over [0, 1). The error was calculated by comparing the sum calculated with float (IEEE-754 binary32) to that calculated with double (IEEE-754 binary64). (The samples were all multiples of 2⁻²⁴, and double has enough precision so that the sum for up to 2²⁹ such values is exact.)

The green line is c n √n with c set to match the last point of the blue line. We see it tracks the blue line over the long term. At points where the average sum crosses a power of two, the mean error increases faster for a time. At these points, the sum has entered a new binade, and further additions have higher average errors due to the increased ULP. Over the course of the binade, this fixed ULP decreases relative to n, bringing the blue line back to the green line.