与numpy及其解决方案相比，mpmath中的元素操作较慢 [英] Elementwise operations in mpmath slow compared to numpy and its solution

问题描述

我进行了一些涉及阶乘分解的计算，因此我决定使用任意精度库mpmath.

I have some calculations that involve factorials that explode pretty fast so I resolved to use the arbitrary precision library 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 s，而temp具有一堆mp.mpf s.我尝试用mp.mpf s初始化M，但是随后一切变慢了.

The thing that seems to take the most time is the very last line and I suspect it is because M has a bunch of floats and temp has a bunch of mp.mpfs. I tried initializing M with mp.mpfs but then everything slowed down.

这是我得到的输出:

part1:        0.00429606437683
part2 so far: 0.00184297561646
part2 finally 1.9477159977

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

Any ideas how I can speed this up?

推荐答案

gmpy2明显快于mpmath.以下代码在我的计算机上的运行速度提高了约12倍.

gmpy2 is significantly faster that mpmath for this type of calculation. The following code runs about 12x faster on my machine.

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的本机整数进行计算.如果gmpy2可用，它将使用gmpy2提供的更快的整数类型.如果只需要gmpy2直接提供的功能之一，那么直接使用gmpy2通常会更快.

mpmath is written in Python and normally uses Python's native integers for its computations. If gmpy2 is available, it will use the faster integer type provided by gmpy2. If you just need one of the functions that is provided directly by gmpy2, then using gmpy2 directly is usually faster.

更新

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

I ran a few experiments. What's actually happening may not be what you expect. When you calculate temp, the values can either be an integer (math.factorial, gmpy.fac, or gmpy2.fac) or a floating-point value (gmpy2.factorial, mpmath.fac). When numpy computes M *= temp, all the values in temp are converted to a 64-bit float. If the value is an integer, the conversion raises an OverflowError. If the value is a floating point number, the conversion returns infinity. You can see this by changing c to np.arange(300) and print M at the end. If you use gmpy.fac or math.factorial, you will get and OverflowError. If you use mpmath.factorial or gmpy2.factorial, you won't get an OverflowError but the resulting M will contain infinities.

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

If you are trying to avoid the OverflowError, you will need to calculate temp with floating point values so the conversion to a 64-bit float will result in infinity.

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

If you aren't encountering OverflowError, then math.factorial is the fastest option.

如果您尝试同时避免使用OverflowError和无穷大，则需要在整个过程中使用mpmath.mpf或gmpy2.mpfr浮点类型. (不要尝试使用gmpy.mpf.)

If you are trying to avoid both OverflowError and infinities, then you'll need to use either the mpmath.mpf or the gmpy2.mpfr floating point types throughout. (Don't try to use gmpy.mpf.)

更新#2

以下是使用gmpy2.mpfr且精度为200位的示例.使用c=np.arange(30)，它比原始示例快5倍.我使用c = np.arange(300)展示它，因为那样会生成OverflowError或无穷大.较大范围的总运行时间与原始代码大致相同.

Here is an example that uses gmpy2.mpfr with a precision of 200 bits. With c=np.arange(30), it is ~5x faster than your original example. I show it using c = np.arange(300) since that would either generate an OverflowError or infinities. The total running time for the larger range is about the same as your original code.

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.

这篇关于与numpy及其解决方案相比，mpmath中的元素操作较慢的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！