python - mpmath 中的元素运算比 numpy 慢及其解决方案

Question

我有一些涉及阶乘的计算，它们的爆炸速度非常快，因此我决定使用任意精度库 mpmath。

我的代码如下所示:

import numpy as np
import mpmath as mp
import time

a    = np.linspace( 0, 100e-2, 100 )
b    = np.linspace( 0, np.pi )
c    = np.arange( 30 )

t    = time.time()
M    = np.ones( [ len(a), len(b), len(c) ] )
A, B = np.meshgrid( a, b, indexing = 'ij' )
temp = A**2 + B
temp = np.reshape( temp, [ len(a), len(b), 1 ] )
temp = np.repeat( temp, len(c), axis = 2 )
M   *= temp
print 'part1:      ', time.time() - t
t    = time.time()

temp = np.array( [ mp.fac(x) for x in c ] )
temp = np.reshape( temp, [ 1, 1, len(c) ] )
temp = np.repeat(  temp, len(a), axis = 0 )
temp = np.repeat(  temp, len(b), axis = 1 )
print 'part2 so far:', time.time() - t
M   *= temp
print 'part2 finally', time.time() - t
t    = time.time()

似乎花费最多时间的是最后一行，我怀疑这是因为 M 有一堆 float 和 temp 有一堆 mp.mpf。我尝试使用 mp.mpfs 初始化 M，但随后一切都变慢了。

这是我得到的输出:

part1:        0.00429606437683
part2 so far: 0.00184297561646
part2 finally 1.9477159977

有什么想法可以加快速度吗？

Answer 1

对于此类计算，

gmpy2 明显快于 mpmath。以下代码在我的机器上运行速度大约是原来的 12 倍。

import numpy as np
import gmpy2 as mp
import time

a = np.linspace(0, 100e-2, 100)
b = np.linspace(0, np.pi)
c = np.arange(30)

t = time.time()
M = np.ones([len(a), len(b), len(c)])
A, B = np.meshgrid( a, b, indexing = 'ij' )
temp = A**2+B
temp = np.reshape(temp, [len(a), len(b), 1])
temp = np.repeat(temp, len(c), axis=2)
M *= temp
print 'part1:', time.time() - t
t = time.time()

temp = np.array([mp.factorial(x) for x in c])
temp = np.reshape(temp, [1, 1, len(c)])
temp = np.repeat(temp, len(a), axis=0)
temp = np.repeat(temp, len(b), axis=1)
print 'part2 so far:', time.time() - t
M *= temp
print 'part2:', time.time() - t
t = time.time()

mpmath 是用 Python 编写的，通常使用 Python 的 native 整数进行计算。如果gmpy2可用，它将使用gmpy2提供的更快的整数类型。如果您只需要 gmpy2 直接提供的功能之一，那么直接使用 gmpy2 通常会更快。

更新

我进行了一些实验。实际发生的事情可能不是你所期望的。计算 temp 时，值可以是整数(math.factorial、gmpy.fac 或 gmpy2.fac)或浮点值(gmpy2.factorial、mpmath.fac)。当 numpy 计算 M *= temp 时，temp 中的所有值都会转换为 64 位 float 。如果该值是整数，则转换会引发 OverflowError。如果该值为 float ，则转换返回无穷大。您可以通过将 c 更改为 np.arange(300) 并在末尾打印 M 来看到这一点。如果您使用gmpy.fac或math.factorial，您将得到OverflowError。如果您使用 mpmath.factorial 或 gmpy2.factorial，您不会得到 OverflowError，而是得到 M将包含无穷大。

如果您试图避免 OverflowError，则需要使用浮点值计算 temp，以便转换为 64 位 float 将导致无穷大。

如果您没有遇到 OverflowError，那么 math.factorial 是最快的选择。

如果您想避免 OverflowError 和无穷大，那么您需要使用 mpmath.mpf 或 gmpy2.mpfr 浮点类型。 (不要尝试使用gmpy.mpf。)

更新#2

这里是一个使用精度为 200 位的 gmpy2.mpfr 的示例。使用c=np.arange(30)，它比原始示例快约 5 倍。我使用 c = np.arange(300) 来显示它，因为这会生成 OverflowError 或无穷大。较大范围的总运行时间与原始代码大致相同。

import numpy as np
import gmpy2
import time

from gmpy2 import mpfr

gmpy2.get_context().precision = 200

a = np.linspace(mpfr(0), mpfr(1), 100)
b = np.linspace(mpfr(0), gmpy2.const_pi())
c = np.arange(300)

t = time.time()
M = np.ones([len(a), len(b), len(c)], dtype=object)
A, B = np.meshgrid( a, b, indexing = 'ij' )
temp = A**2+B
temp = np.reshape(temp, [len(a), len(b), 1])
temp = np.repeat(temp, len(c), axis=2)
M *= temp
print 'part1:', time.time() - t
t = time.time()

temp = np.array([gmpy2.factorial(x) for x in c], dtype=object)
temp = np.reshape(temp, [1, 1, len(c)])
temp = np.repeat(temp, len(a), axis=0)
temp = np.repeat(temp, len(b), axis=1)
print 'part2 so far:', time.time() - t
M *= temp
print 'part2:', time.time() - t
t = time.time()

免责声明:我维护 gmpy2。