15

我四处寻找一些关于 numpy/scipy 函数如何在数值稳定性方面表现的文档，例如，是否采取了任何措施来提高数值稳定性，或者是否有替代的稳定实现。

我+对浮点数组numpy.sum()、numpy.cumsum()和的加法（运算符）特别感兴趣numpy.dot()。在所有情况下，我本质上都是对大量浮点数求和，我担心这种计算的准确性。

有谁知道 numpy/scipy 文档或其他来源中对此类问题的任何引用？

4

1 回答 1

2

短语“稳定性”是指算法。如果您的算法一开始就不稳定，那么提高精度或减少组件步骤的舍入误差不会获得太多收益。

像“solve”这样更复杂的 numpy 例程是 ATLAS/BLAS/LAPACK 例程的包装器。您可以参考那里的文档，例如“dgesv”使用具有部分旋转和行交换的 LU 分解求解实值线性方程组：LAPACK 的基础 Fortran 代码文档可以在此处查看http://www.netlib.org /lapack/explore-html/但http://docs.scipy.org/doc/numpy/user/install.html指出有许多不同版本的标准例程实现可用 - 速度优化和精度将因它们而异。

您的示例没有引入太多舍入，“+”没有不必要的舍入，当较小的数字具有无法在答案中表示的低位时，精度完全取决于浮点数据类型中隐含的舍入。Sum 和 dot 仅取决于评估的顺序。Cumsum 不能轻易地重新排序，因为它输出一个数组。

对于“cumsum”或“dot”函数期间的累积舍入，您可以选择：

在 Linux 64 位上，numpy 提供了对高精度“long double”类型 float128 的访问，您可以使用它来减少中间计算中的精度损失，但会以性能和内存为代价。

但是在我的 Win64 上安装“numpy.longdouble”映射到“numpy.float64”一个普通的 C 双精度类型，所以你的代码不是跨平台的，检查“finfo”。（Canopy Express Win64 上不存在真正更高精度的 float96 或 float128）

log2(finfo(float64).resolution)
> -49.828921423310433
actually 53-bits of mantissa internally ~ 16 significant decimal figures

log2(finfo(float32).resolution)
> -19.931568                          # ~ only 7 meaningful digits

由于sum()并将dot()数组减少为单个值，因此使用内置函数可以轻松最大化精度：

x = arange(1, 1000000, dtype = float32)
y = map(lambda z : float32(1.0/z), arange(1, 1000000))
sum(x)                     #  4.9994036e+11
sum(x, dtype = float64)    #  499999500000.0
sum(y)                     #  14.357357
sum(y, dtype = float64)    #  14.392725788474309
dot(x,y)                             # 999999.0
einsum('i,i', x, y)                  # * dot product = 999999.0
einsum('i,i', x, y, dtype = float64) # 999999.00003965141
  • 请注意，在这种情况下，“点”内的单精度舍入取消，因为每个近似整数都被舍入为精确整数

优化舍入取决于您要相加的事物类型 - 首先添加许多小数字可以帮助延迟舍入，但不会避免存在大数字但相互抵消的问题，因为中间计算仍然会导致精度损失

显示评估顺序依赖性的示例...

x = array([ 1., 2e-15, 8e-15 , -0.7, -0.3], dtype=float32)
# evaluates to
# array([  1.00000000e+00,   2.00000001e-15,   8.00000003e-15,
#   -6.99999988e-01,  -3.00000012e-01], dtype=float32)
sum(x) # 0
sum(x,dtype=float64) # 9.9920072216264089e-15
sum(random.permutation(x)) # gives 9.9999998e-15 / 2e-15 / 0.0
于 2013-09-10T00:38:25.503 回答