添加数字时在 python/numpy 中浮点精度细分

Float precision breakdown in python/numpy when adding numbers

提问人:Dändän 提问时间:1/11/2013 最后编辑:Dändän 更新时间:1/11/2013 访问量:5390

问:

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

enter image description here

如您所见,相当低的总和值分散在零附近或正好为零,而乘法值则很好......

拜托,有人可以帮忙/解释吗?多谢!

Python NumPy 精度

评论

2赞 eumiro 1/11/2013
如果还不够,还有.它有同样的问题吗?np.float64np.float128
0赞 Dändän 1/11/2013
是的,float96 也有同样的问题......
0赞 jeff7 1/12/2013
为什么不在对数域中做这些求和以避免精度问题呢?lingpipe-blog.com/2012/02/16/......

答:

6赞 mgilson 1/11/2013 #1

由于舍入误差,浮点精度对加法/减法非常敏感。最终,它变得如此之大,以至于将 1 添加到 exp(x) 会得到与 exp(x) 相同的结果。在双精度中,大约是:1+exp(x)exp(x) == 1e16

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

请注意,这大约是 37 -- 这大致是你的情节中事情变得疯狂的地方。math.log(1e16)

您可能会遇到相同的问题,但规模不同:

>>> (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)
3赞 ecatmur 1/11/2013 #2

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

例如,在计算 时,和具有相同的值,因为加到 没有效果,因为它超出了 64 位浮点的精度:func_sum(200)math.exp(-200)1/(1+math.exp(200))1math.exp(200)

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

这就解释了为什么给出零,但是位于 x 轴之外的点呢?这些也是由浮点不精确引起的;偶尔会发生不等于的情况;理想情况下,是最接近 的浮点值,并且是最接近 的浮点数的浮点值,计算公式为 ,不一定是 。事实上,并且非常接近,实际上仅在最后一个单元中有所不同:func_sum(200)math.exp(-x)1/math.exp(x)math.exp(x)e^x1/math.exp(x)math.exp(x)e^-xmath.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

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

4赞 unutbu 1/11/2013 #3

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

让我们看一下计算出错的第一点,当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.expnp.expnp.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)