fsum for numpy.arrays, stable summation

Question

I have a number of multidimensional numpy.arrays with small values that I need to add up with little numerical error. For floats, there is math.fsum (with its implementation here), which has always served me well. numpy.sum isn't stable enough.

How can I get a stable summation for numpy.arrays?

---

Background

This is for the quadpy package. The arrays of small values are the evaluations of a function at specific points of (many) intervals, times their weights. The sum of these is an approximation of the integral of said function over the intervals.

Answer 1

Alright then, I've implemented accupy which gives a few stable summation algorithms.

Here's a quick and dirty implementation of Kahan summation for numpy arrays. Notice, however, that it is not not very accurate for ill-conditioned sums.

def kahan_sum(a, axis=0):
    '''Kahan summation of the numpy array along an axis.
    '''
    s = numpy.zeros(a.shape[:axis] + a.shape[axis+1:])
    c = numpy.zeros(s.shape)
    for i in range(a.shape[axis]):
        # https://stackoverflow.com/a/42817610/353337
        y = a[(slice(None),) * axis + (i,)] - c
        t = s + y
        c = (t - s) - y
        s = t.copy()
    return s

It does the job, but it's slow because it's Python-looping over the axis-th dimension.