使用RuntimeWarning进行numpy除法:double_scalars中遇到无效值

Question

我写了以下脚本:

import numpy

d = numpy.array([[1089, 1093]])
e = numpy.array([[1000, 4443]])
answer = numpy.exp(-3 * d)
answer1 = numpy.exp(-3 * e)
res = answer.sum()/answer1.sum()
print res

但是我得到了这个结果并且发生了错误:

nan
C:\Users\Desktop\test.py:16: RuntimeWarning: invalid value encountered in double_scalars
  res = answer.sum()/answer1.sum()

似乎输入元素太小,以至于python将它们变成零,但实际上除法有其结果.

如何解决这类问题？

Answer 1

你无法解决它.简单地说answer1.sum()==0,你不能执行除零.

这answer1是因为是2个非常大的负数的指数,所以结果四舍五入为零.

nan 在这种情况下返回是因为除以零.

现在要解决您的问题,您可以:

• 去找一个高精度数学的图书馆,比如mpmath.但那不那么有趣.
• 作为更大武器的替代品,做一些数学操作,如下所述.
• 寻求量身定制的scipy/numpy功能,完全符合您的要求!看看@Warren Weckesser的回答.

在这里,我解释如何做一些有助于解决这个问题的数学运算.我们有这个分子:

exp(-x)+exp(-y) = exp(log(exp(-x)+exp(-y)))
                = exp(log(exp(-x)*[1+exp(-y+x)]))
                = exp(log(exp(-x) + log(1+exp(-y+x)))
                = exp(-x + log(1+exp(-y+x)))

以上x=3* 1089和y=3* 1093.现在,这个指数的论证是

-x + log(1+exp(-y+x)) = -x + 6.1441934777474324e-06

对于分母,您可以进行类似的操作,但获得log(1+exp(-z+k))已经舍入到0的分母,以便分母处的指数函数的参数简单地舍入到-z=-3000.然后你就得到了你的结果

exp(-x + log(1+exp(-y+x)))/exp(-z) = exp(-x+z+log(1+exp(-y+x)) 
                                   = exp(-266.99999385580668)

如果你只保留2个主要术语(即1089分子中的第一个数字和分母的第一个数字1000),那么你已经非常接近结果:

exp(3*(1089-1000))=exp(-267)

为了它,让我们看看我们与Wolfram alpha(链接)的解决方案有多接近:

Log[(exp[-3*1089]+exp[-3*1093])/([exp[-3*1000]+exp[-3*4443])] -> -266.999993855806522267194565420933791813296828742310997510523

这个数字和上面的指数之间的差异是+1.7053025658242404e-13,所以我们在分母处的近似值很好.

最终的结果是

'exp(-266.99999385580668) = 1.1050349147204485e-116

来自wolfram alpha是(链接)

1.105034914720621496.. × 10^-116 # Wolfram alpha.

而且,在这里使用numpy也是安全的.

Answer 2

你可以使用np.logaddexp(在@ gg349的答案中实现这个想法):

In [33]: d = np.array([[1089, 1093]])

In [34]: e = np.array([[1000, 4443]])

In [35]: log_res = np.logaddexp(-3*d[0,0], -3*d[0,1]) - np.logaddexp(-3*e[0,0], -3*e[0,1])

In [36]: log_res
Out[36]: -266.99999385580668

In [37]: res = exp(log_res)

In [38]: res
Out[38]: 1.1050349147204485e-116

或者您可以使用scipy.special.logsumexp:

In [52]: from scipy.special import logsumexp

In [53]: res = np.exp(logsumexp(-3*d) - logsumexp(-3*e))

In [54]: res
Out[54]: 1.1050349147204485e-116