10

由于 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()

在此处输入图像描述

如您所见，非常低的总和值分散在零附近或恰好为零，而相乘的值很好......

请，有人可以帮助/解释吗？非常感谢！

4

3 回答 3

6

由于舍入误差，浮点精度对加法/减法非常敏感。最终，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)
于 2013-01-11T14:10:17.163 回答
4

这是灾难性取消的一个例子。

让我们看看计算出错的第一点，当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)
于 2013-01-11T14:14:55.203 回答
3

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

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 位浮点中有效数的大小）。

于 2013-01-11T14:14:12.197 回答