加法时python/numpy中的浮点数精度分解

Question

我有一些问题，因为与numpy一起使用的数字非常低。我花了几个星期的时间来追溯我在数值积分中经常遇到的问题，当我在函数中添加浮点数时，float64精度就会丢失。使用乘积而不是总和执行数学上的恒等式计算，得到的值是正确的。

以下是代码示例和结果图：

from matplotlib.pyplot import *
from numpy import vectorize, arange
import math

def func_product(x):
    return math.exp(-x)/(1+math.exp(x))

def func_sum(x):
    return math.exp(-x)-1/(1+math.exp(x))

#mathematically, both functions are the same

vecfunc_sum = vectorize(func_sum)
vecfunc_product = vectorize(func_product)

x = arange(0.,300.,1.)
y_sum = vecfunc_sum(x)
y_product = vecfunc_product(x)

plot(x,y_sum,    'k.-', label='sum')
plot(x,y_product,'r--',label='product')

yscale('symlog', linthreshy=1E-256)
legend(loc='lower right')
show()

​

​

正如您所看到的，非常低的相加值分散在零附近，或者恰好是零，而相乘的值很好……

有没有人能帮忙/解释一下？非常感谢!

Answer 1

由于舍入误差，浮点精度对加/减非常敏感。最终，1+exp(x)变得如此之大，以至于将1与exp(x)相加得到与exp(x)相同的结果。在双精度下，这在exp(x) == 1e16附近

>>> (1e16 + 1) == (1e16)
True
>>> (1e15 + 1) == (1e15)
False

请注意，math.log(1e16)大约是37 --这大致就是您的图中出现问题的地方。

您可以有相同的问题，但规模不同：

>>> (1e-16 + 1.) == (1.)
True
>>> (1e-15 + 1.) == (1.)
False

对于你的制度中的绝大多数点，你的func_product实际上是在计算：

exp(-x)/exp(x) == exp(-2*x)

这就是为什么你的曲线图有一个很好的斜率-2。

把它推向另一个极端，你的其他版本正在计算(至少是近似的)：

exp(-x) - 1./exp(x) 

这大约是

exp(-x) - exp(-x)

Answer 2

这是catastrophic cancellation的一个例子。

让我们看一下计算出错的第一个点，当x = 36.0

In [42]: np.exp(-x)
Out[42]: 2.3195228302435691e-16

In [43]: - 1/(1+np.exp(x))
Out[43]: -2.3195228302435691e-16

In [44]: np.exp(-x) - 1/(1+np.exp(x))
Out[44]: 0.0

使用func_product的计算不会减去几乎相等的数字，因此它避免了灾难性的取消。

顺便说一句，如果你把math.exp改成np.exp，你可以去掉np.vectorize (它很慢)：

def func_product(x):
    return np.exp(-x)/(1+np.exp(x))

def func_sum(x):
    return np.exp(-x)-1/(1+np.exp(x))

y_sum = func_sum_sum(x)
y_product = func_product_product(x)

Answer 3

问题是您的func_sum是numerically unstable，因为它涉及两个非常接近的值之间的减法。

例如，在func_sum(200)的计算中，math.exp(-200)和1/(1+math.exp(200))的值是相同的，因为将1加到math.exp(200)是无效的，因为它超出了64位浮点的精度：

math.exp(200).hex()
0x1.73f60ea79f5b9p+288

(math.exp(200) + 1).hex()
0x1.73f60ea79f5b9p+288

(1/(math.exp(200) + 1)).hex()
0x1.6061812054cfap-289

math.exp(-200).hex()
0x1.6061812054cfap-289

这就解释了为什么func_sum(200)会给出零，但是位于x轴之外的点又如何呢？这也是由浮点不精确引起的；偶尔会发生math.exp(-x)不等于1/math.exp(x)的情况；理想情况下，math.exp(x)是最接近e^x的浮点值，而1/math.exp(x)是最接近math.exp(x)计算的浮点数的倒数的浮点值，不一定是e^-x。事实上，math.exp(-100)和1/(1+math.exp(100))非常接近，实际上只是在最后一个单元中有所不同：

math.exp(-100).hex()
0x1.a8c1f14e2af5dp-145

(1/math.exp(100)).hex()
0x1.a8c1f14e2af5cp-145

(1/(1+math.exp(100))).hex()
0x1.a8c1f14e2af5cp-145

func_sum(100).hex()
0x1.0000000000000p-197

因此，您实际计算的是math.exp(-x)和1/math.exp(x)之间的差异。您可以跟踪函数math.pow(2, -52) * math.exp(-x)的行，以查看它是否通过func_sum的正值传递(回想一下，52是64位浮点中有效位的大小)。