Numpy 除法运行时警告: 在 double_scalars 中遇到无效值

我写了以下剧本:

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 把它们变成了零，但事实上除法有它的结果。

如何解决这类问题？

358006 次浏览
小开
最佳答案

You can't solve it. Simply answer1.sum()==0, and you can't perform a division by zero.

This happens because answer1 is the exponential of 2 very large, negative numbers, so that the result is rounded to zero.

nan is returned in this case because of the division by zero.

Now to solve your problem you could:

  • go for a library for high-precision mathematics, like mpmath. But that's less fun.
  • as an alternative to a bigger weapon, do some math manipulation, as detailed below.
  • go for a tailored scipy/numpy function that does exactly what you want! Check out @Warren Weckesser answer.

Here I explain how to do some math manipulation that helps on this problem. We have that for the numerator:

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)))

where above x=3* 1089 and y=3* 1093. Now, the argument of this exponential is

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

For the denominator you could proceed similarly but obtain that log(1+exp(-z+k)) is already rounded to 0, so that the argument of the exponential function at the denominator is simply rounded to -z=-3000. You then have that your result is

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

which is already extremely close to the result that you would get if you were to keep only the 2 leading terms (i.e. the first number 1089 in the numerator and the first number 1000 at the denominator):

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

For the sake of it, let's see how close we are from the solution of Wolfram alpha (link):

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

The difference between this number and the exponent above is +1.7053025658242404e-13, so the approximation we made at the denominator was fine.

The final result is 

'exp(-266.99999385580668) = 1.1050349147204485e-116

From wolfram alpha is (link)

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

and again, it is safe to use numpy here too.

小开

You can use np.logaddexp (which implements the idea in @gg349's answer):

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

Or you can use 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

提交答案