numpy.arrays的fsum，稳定求和

Question

我有一些numpy.array带有较小值的多维s，我需要将它们相加且数值误差很小。对于floats，这里有math.fsum（及其实现在此处），对我一直很好。numpy.sum不够稳定。

如何获得numpy.arrays 的稳定求和？

---

背景

这是用于Quadpy软件包的。小值数组是在（许多）时间间隔的特定点上乘以其权重的函数求值。这些总和是所述函数在间隔上的积分的近似值。

Answer 1

好吧，我已经实现了accupy，它提供了一些稳定的求和算法。

这是numpy 数组的Kahan 求和的快速而肮脏的实现。但是请注意，它对于病态和不是很准确。

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]):
        # /sf/answers/2997232731/
        y = a[(slice(None),) * axis + (i,)] - c
        t = s + y
        c = (t - s) - y
        s = t.copy()
    return s

它完成了这项工作，但速度很慢，因为它在axis第 -th 维度上进行 Python 循环。