"Divide by zero encountered in log" 参数非常大

Question

我正在尝试评估此功能：

import scipy as sci
import numpy as np

np.pi*sci.special.gamma(2*i-1)/((2**(2*(i-1)))*(sci.special.gamma(i)**2))

这运行在“被零除”方面遇到了麻烦，所以我使用日志将其分成几部分以查看被零除出现的位置。

该方法在伪代码中看起来是：

log(f) = log(pi)+log(gamma(...))-log(2**...)-log(gamma(...))

问题出现在第三部分

我想在 i 的范围内求解这个项，所以有趣的是，在尝试循环时出现：

N = np.arange(2,21)
for i in N:
    f = np.log(2**(2*(i-1)))
    print(i, f)

Out:
    2 1.3862943611198906
    3 2.772588722239781
    4 4.1588830833596715
    5 5.545177444479562
    6 6.931471805599453
    7 8.317766166719343
    8 9.704060527839234
    9 11.090354888959125
    10 12.476649250079015
    11 13.862943611198906
    12 15.249237972318797
    13 16.635532333438686
    14 18.021826694558577
    15 19.408121055678468
    16 20.79441541679836
    17 -inf
    18 -inf
    19 -inf
    20 -inf
    __main__:2: RuntimeWarning: divide by zero encountered in log

看来，当参数变得太大 (i=17) 时，问题就出现了。

另一方面，当我只计算 i=17 的单个值时，它起作用了（对于更大的 i 也是如此）：

np.log(2**(2*(17-1)))
Out: 22.18070977791825

除以零的地方在哪里？当涉及到大参数时，我通常如何克服这个问题？

Answer 1

似乎 numpy 对大整数不满意。

如果您更改示例中的第一行

N = np.arange(2,21)

到

N = np.arange(2.0,21.0)

它应该工作得更好（浮点数而不是整数）

编辑：

之间的差异

import numpy as np
N = np.arange(2,100,dtype=int)
N = 2**(2*(N-1))
print(N)

和

import numpy as np
N = np.arange(2,100,dtype=float)
N = 2**(2*(N-1))
print(N)

Answer 2

由于 np.arange 的隐式类型，主要问题来自 整数溢出 。实际上，np.arange(2,21) 似乎会生成一个 int32 类型值的数组，因此 i 也是 int32 类型，因此 2**(2*(i-1)) 也是。 2**(2*(i-1)) 和 i>=17 大于 2**31（int32 类型）的限制。溢出导致值被截断为 0，因此在计算日志时会出现 Numpy 错误（因为没有为输入 0 定义日志函数）。

一个简单的解决方案是将类型更改为 int64 值，但这并不是很好，因为输出指数会快速增长。使用花车会好一点。但是，这远不是一个好主意，因为 i 值接近 500 或更高时会很快达到最大指数。问题的根源在于计算的数值稳定性。

一个好的解决方案是使用数学变换 来解决数值稳定性问题。事实上，log2(2**n) = n 和 log(n) = log2(n)*log(2)。因此，log(2**(2*(i-1))) = log2(2**(2*(i-1)))*log(2) = (2*(i-1))*log(2)。因此，以下代码应该给出等效的结果，但没有数字问题：

N = np.arange(2,21)
log2 = np.log(2)
for i in N:
    f = (2*(i-1))*log2
    print(i, f)

请注意，此代码实际上计算速度更快并且更简单。另请注意，使用此方法使用巨大的 i 值（例如大于 1,000,000 的值）没有问题！这是结果：

2 1.3862943611198906
3 2.772588722239781
4 4.1588830833596715
5 5.545177444479562
6 6.931471805599453
7 8.317766166719343
8 9.704060527839234
9 11.090354888959125
10 12.476649250079015
11 13.862943611198906
12 15.249237972318797
13 16.635532333438686
14 18.021826694558577
15 19.408121055678468
16 20.79441541679836
17 22.18070977791825
18 23.56700413903814
19 24.95329850015803
20 26.33959286127792

"Divide by zero encountered in log" 参数非常大

"Divide by zero encountered in log" with very large arguments

python

numpy

divide-by-zero