在 numpy 中对 ndarray 求和同时最小化浮点不准确的最有效方法是什么？

Question

我有一个大矩阵，其值的数量级变化很大。为了尽可能准确地计算总和，我的方法是将 ndarray 重塑为一维数组，对其进行排序，然后将其相加，从最小的条目开始。有没有更好/更有效的方法来做到这一点？

Answer 1

我认为，考虑到浮点精度问题，最适合您的任务的算法是Kahan summation。出于实际目的，卡汉求和的误差界限与被加数的数量无关，而朴素求和的误差界限随着被加数的数量线性增长。

NumPy 不使用 Kahan 求和，并且没有简单的方法可以在不牺牲性能的情况下实现它。但它使用了次优的方法，成对求和，在一些合理的假设下，例如被加数数量的对数的平方根，误差会增加。

因此，Numpy 本身很可能已经能够为您的问题提供足够好的精度。为了验证这一点，我实际上会通过 Kahan 求和运行一些示例案例（上面的维基百科链接中的伪代码可以简单地转换为 Python），并将其作为最佳的最佳结果，并将其与：

1. 按原样在矩阵上调用

   np.sum

   。
2. 在重塑为一维后在矩阵上调用

   np.sum

   ，如果您的矩阵在内存中不连续，这可能会给出更好的结果。
3. 在一维数组的排序版本上调用

   np.sum

   。

在大多数情况下，最后三个选项的行为应该类似，但了解的唯一方法是实际测试它。

Answer 2

我想对 numpy sum 添加一个简单的改进：使用 np.float128 （四倍精度）作为中间值，然后转换回标准 64 位浮点数：

np.float64(np.sum(array, dtype=np.float128))

这会将要求和的值转换为四倍精度，从而显着减少舍入误差。这比 64 位的 numpy.sum 慢，但比 fsum 或首先对数组进行排序要快得多。

作为示例，我正在评估标准柯西分布的 1000x1000 绝对值矩阵

matrix = np.abs(np.random.standard_cauchy([1000, 1000]))

。 128位结果与fsum相同，排序后的np.sum比没有排序要好，但仍然不太精确。

用

%timeit

测量的评估时间：

np.sum(matrix): 370 μs
np.float64(np.sum(matrix, dtype=np.float128)): 4.3 ms
math.fsum(matrix.flatten()): 60 ms
np.sum(np.sort(matrix.flatten())): 78 ms

备注：

1. float128 在 Windows 上不可用。
2. numba 目前还不支持 128 位浮点数。