python - 是否有任何关于 numpy 数值稳定性的文档？

Question

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

我对 float 组的加法(+ 运算符)、numpy.sum()、numpy.cumsum() 特别感兴趣和 numpy.dot()。在所有情况下，我基本上都是对大量 float 求和，我担心此类计算的准确性。

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

Answer 1

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

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

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

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

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

但是在我的 Win64 安装中，“numpy.longdouble”映射到“numpy.float64”一个普通的 C double 类型，所以你的代码不是跨平台的，检查“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