fsum用于numpy数组的稳定求和

Question

我有许多包含小值的多维numpy.array需要相加，要保证数值误差很小。对于float类型，可以使用math.fsum进行求和（其实现方式在这里），一直以来都很好用。但是numpy.sum不够稳定。

那么如何对numpy.array进行稳定求和呢？

---

背景

这是针对{{link1：quadpy软件包}}的。小值数组是在（许多）间隔的特定点处对函数的评估，乘以它们的权重。它们的总和是该函数在间隔上的积分的近似值。

- Nico Schlömer

---

为什么numpy.sum或者your_array.sum()不能满足你的需求？ - ForceBru

谢谢回复。对我来说，numpy.sum不够稳定。我猜它只是一个接一个地添加（小）值，所以最终当大数和小数相加时，信息就会丢失。 - Nico Schlömer

我相信 numpy.sum 使用成对求和技术，因此优于单纯的求和方法。但是，如果您需要更多的功能，NumPy 不提供该功能。如果您的 NumPy 早于对成对求和技术的改进版本，则更新可能会有所帮助。 - user2357112

看起来GitHub上的这个线程可能会有用，尽管它说NumPy没有实现像math.fsum这样更精确的求和。 - ForceBru

Answer 1

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

这里是一个快速而简单的Kahan求和 numpy数组的实现。然而，请注意，对于病态求和，它并不是非常准确。

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

它能完成任务，但由于是使用Python循环遍历axis维度，所以速度较慢。

- Nico Schlömer