有没有numpy数值稳定性的文件？

Question

我环顾四周寻找一些关于numpy/scipy函数在数值稳定性方面的表现的文档,例如,是否有任何改善数值稳定性的方法,或者是否存在其他稳定实现.

我特别感兴趣的+是浮点数组的加法(运算符)numpy.sum(),numpy.cumsum()和numpy.dot().在所有情况下,我基本上总结了大量的浮点数,我担心这种计算的准确性.

有没有人知道在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指出有许多不同版本的标准例程实现可用 - 它们之间的速度优化和精度会有所不同。

您的示例没有引入太多舍入，“+”没有不必要的舍入，当较小的数字具有无法在答案中表示的低位位时，精度完全取决于浮点数据类型中隐式的舍入。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