numpy.arrays的fsum，稳定求和 [英] fsum for numpy.arrays, stable summation

问题描述

我有许多多维numpy.array，它们的值很小 我需要加起来几乎没有数值误差.对于float，有 math.fsum (带有它的实现此处)，这对我一直很有益. numpy.sum 不够稳定.

如何获得numpy.array的稳定总和?

---

背景

这是针对 quadpy软件包的.小值数组是在(许多)时间间隔的特定点上乘以其权重的函数求值.这些总和是该函数在整个区间上的积分的近似值.

解决方案

那么，我已经实现了 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]):
        # 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

它可以完成工作，但是速度很慢，因为它在第axis个维度上进行了Python循环.

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.

解决方案

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.

这篇关于numpy.arrays的fsum，稳定求和的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！