Is there a way to guarantee that calling np.sum(A), np.add.reduce(A), and A.sum() will generate precisely the same value (in an == sense) for a one-dimensional floating point numpy array A?
From the documentation of np.ufunc.reduce, it's explicitly stated that
add.reduce() is equivalent to sum()
and further specified that reduce is left-to-right associative; specifically, that 'op.reduce produces results equivalent to' the following code:
r = op.identity # op = ufunc
for i in range(len(A)):
r = op(r, A[i])
return r
This would be seem to be the expected behavior, as CPython's built-in sum() and functools.reduce() also explicitly guarantee left-to-right evaluation, which is relevant since floating point addition isn't a mathematically associative operation.
However, looking at the documentation of np.sum (which also helpfully mentions that add.reduce and ndarray.sum are equivalent methods) there is a seemingly contradictory note (emphasis mine):
For floating point numbers the numerical precision of sum (and
np.add.reduce) is in general limited by directly adding each number individually to the result causing rounding errors in every step. However, often numpy will use a numerically better approach (partial pairwise summation) leading to improved precision in many use-cases.
The description of when the 'pairwise summation' algorithm is preferred implies that a call to np.sum(A) as above would give a different evaluation order than np.add.reduce(A), and therefore a potentially different result. It's not at all clear to me whether this behavior can be disabled, if disabling the behavior would interfere with numpy's vectorization ability, and where in the source code this decision takes place (or even which functions are affected).
For some background knowledge, it was helpful to look at this question from 2013: What is the difference between np.sum and np.add.reduce?. However, that post (and all others I've seen) was asking about the performance differences, and the top answer by Warren Weckesser states that
when the argument is a numpy array, np.sum ultimately calls add.reduce to do the work.
While I'm not asking about the performance comparison, this (possibly outdated) answer added even more mystery, since it would also seem to imply that np.add.reduce should always give identical output to np.sum, and therefore that the documentation for ufunc.reduce is incorrect in stating that it produces results equivalent to left-to-right associativity since pairwise-summation could occur.
My goal isn't to sum floating point numbers in the most precise way, but rather to ensure that the value I get from using numpy to sum a one-dimensional, float64 numpy array A, evaluates as == to CPython's left-to-right sum(A). As such, I'm not exactly asking about the trade-offs of using either function, but about information regarding the actual evaluation order for these functions on the current version of numpy, given that the documentation mentions they are all 'equivalent' while also implying that they could give unequal results.