Numpy floating point rounding errors

Question

While searching for some numpy stuff, I came across a question discussing the rounding accuracy of numpy.dot():

Numpy: Difference between dot(a,b) and (a*b).sum()

Since I happen to have two (different) Computers with Haswell-CPUs sitting on my desk, that should provide FMAand everything, I thought I'd test the example given by Ophion in the first answer, and I got a result that somewhat surprised me:

After updating/installing/fixing lapack/blas/atlas/numpy, I get the following on both machines:

>>> a = np.ones(1000, dtype=np.float128)+1e-14
>>> (a*a).sum()
1000.0000000000199999
>>> np.dot(a,a)
1000.0000000000199948

>>> a = np.ones(1000, dtype=np.float64)+1e-14
>>> (a*a).sum()
1000.0000000000198
>>> np.dot(a,a)
1000.0000000000176

So the standard multiplication + sum() is more precise than np.dot(). timeit however confirmed that the .dot() version is faster (but not much) for both float64 and float128.

Can anyone provide an explanation for this?

edit: I accidentally deleted the info on numpy versions: same results for 1.9.0 and 1.9.3 with python 3.4.0 and 3.4.1.

Answer 1

It looks like they recently added a special Pairwise Summation to ndarray.sum for improved numerical stability.

From PR 3685, this affects:

all add.reduce calls that go over float_add with IS_BINARY_REDUCE true
so this also improves mean/std/var and anything else that uses sum.

See here for code changes.