编辑3：

最终(我认为)版本,更清晰,更快速地融入max9111’s answer的创意.

import numpy as np

from numba import as nb

@nb.njit()

def func1_jit(a, b, c, d):

# Precompute

exp_min = 5 - (a + b + c + d)

exp_max = b

exp = 2. ** np.arange(exp_min, exp_max + 1)

fact_e = np.empty((a + b - 2))

fact_e[0] = 1

for ei in range(1, len(fact_e)):

fact_e[ei] = ei * fact_e[ei - 1]

# Loops

B = 0

for ai in range(0, a):

for bi in range(0, b):

for ci in range(0, c):

for di in range(0, d):

for ei in range(0, ai + bi):

for fi in range(0, ci + di):

B += exp[ei - fi - ai - ci - di + 1 - exp_min] * (ei * ei - 2 * (ei * fi) - 7 * di) * fact_e[ei]

return B

这已经比以前的任何选项都快得多,但我们仍然没有利用多个CPU.一种方法是在函数本身内,例如并行化外循环.这会在每次调用时增加一些开销来创建线程,因此对于小输入实际上有点慢,但对于更大的值应该明显更快：

import numpy as np

from numba import as nb

@nb.njit(parallel=True)

def func1_par(a, b, c, d):

# Precompute

exp_min = 5 - (a + b + c + d)

exp_max = b

exp = 2. ** np.arange(exp_min, exp_max + 1)

fact_e = np.empty((a + b - 2))

fact_e[0] = 1

for ei in range(1, len(fact_e)):

fact_e[ei] = ei * fact_e[ei - 1]

# Loops

B = np.empty((a,))

for ai in nb.prange(0, a):

Bi = 0

for bi in range(0, b):

for ci in range(0, c):

for di in range(0, d):

for ei in range(0, ai + bi):

for fi in range(0, ci + di):

Bi += exp[ei - fi - ai - ci - di + 1 - exp_min] * (ei * ei - 2 * (ei * fi) - 7 * di) * fact_e[ei]

B[ai] = Bi

return np.sum(B)

或者,如果您有许多要评估函数的点,也可以在该级别进行并行化.这里a_arr,b_arr,c_arr和d_arr是要评估函数的值的向量：

from numba import as nb

@nb.njit(parallel=True)

def func1_arr(a_arr, b_arr, c_arr, d_arr):

B_arr = np.empty((len(a_arr),))

for i in nb.prange(len(B_arr)):

B_arr[i] = func1_jit(a_arr[i], b_arr[i], c_arr[i], d_arr[i])

return B_arr

最佳配置取决于您的输入,使用模式,硬件等,因此您可以根据您的情况组合不同的想法.

编辑2：

实际上,忘记我之前说过的话.最好的是JIT编译算法,但是以更有效的方式.首先计算昂贵的部分(我采用指数和阶乘),然后将其传递给编译的循环函数：

import numpy as np

from numba import njit

def func1(a, b, c, d):

exp_min = 5 - (a + b + c + d)

exp_max = b

exp = 2. ** np.arange(exp_min, exp_max + 1)

ee = np.arange(a + b - 2)

fact_e = scipy.special.factorial(ee)

return func1_inner(a, b, c, d, exp_min, exp, fact_e)

@njit()

def func1_inner(a, b, c, d, exp_min, exp, fact_e):

B = 0

for ai in range(0, a):

for bi in range(0, b):

for ci in range(0, c):

for di in range(0, d):

for ei in range(0, ai + bi):

for fi in range(0, ci + di):

B += exp[ei - fi - ai - ci - di + 1 - exp_min] * (ei * ei - 2 * (ei * fi) - 7 * di) * fact_e[ei]

return B

在我的实验中,这是迄今为止最快的选项,并且只占用很少的额外内存(只有预先计算的值,输入上的大小为线性).

a, b, c, d = 4, 6, 3, 4

# The original function

%timeit func1_orig(a, b, c, d)

# 2.07 ms ± 33.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# The grid-evaluated function

%timeit func1_grid(a, b, c, d)

# 256 µs ± 25 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

# The precompuation + JIT-compiled function

%timeit func1_jit(a, b, c, d)

# 19.6 µs ± 3.25 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)

那么总是可以选择网格评估整个事情：

import numpy as np

import scipy.special

def func1(a, b, c, d):

ai, bi, ci, di, ei, fi = np.ogrid[:a, :b, :c, :d, :a + b - 2, :c + d - 2]

# Compute

B = (2.) ** (ei - fi - ai - ci - di + 1) * (ei ** 2 - 2 * (ei * fi) - 7 * di) * scipy.special.factorial(ei)

# Mask out of range elements for last two inner loops

m = (ei < ai + bi) & (fi < ci + di)

return np.sum(B * m)

print(func1(4, 6, 3, 4))

# 21769947.844726562

显而易见,随着参数的增加,其内存成本将快速增长.代码实际上执行的计算比必要的多,因为两个内部循环具有不同的迭代次数,因此(在此方法中)您必须使用最大的,然后删除您不需要的.希望是矢量化将弥补这一点.一个小的IPython基准：

a, b, c, d = 4, 6, 3, 4

# func1_orig is the original loop-based version

%timeit func1_orig(a, b, c, d)

# 2.9 ms ± 110 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# func1 here is the vectorized version

%timeit func1(a, b, c, d)

# 210 µs ± 6.34 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

编辑：

请注意,以前的方法也不是一个全有或全无的事情.您可以选择仅对某些循环进行网格评估.例如,两个最里面的循环可以像这样矢量化：

def func1(a, b, c, d):

B = 0

e = np.arange(a + b - 2).reshape((-1, 1))

f = np.arange(c + d - 2)

for ai in range(0, a):

for bi in range(0, b):

ei = e[:ai + bi]

for ci in range(0, c):

for di in range(0, d):

fi = f[:ci + di]

B += np.sum((2.) ** (ei - fi - ai - ci - di + 1) * (ei ** 2 - 2 * (ei * fi) - 7 * di) * scipy.special.factorial(ei))

return B

这仍然有循环,但它确实避免了额外的计算,并且内存要求低得多.哪一个最好取决于我猜的输入大小.在我的测试中,使用原始值(4,6,3,4),这甚至比原始函数慢;此外,对于这种情况,似乎在每个循环上为ei和fi创建新数组比在预先创建的循环上操作更快.但是,如果将输入乘以4(14,24,12,16),那么这比原始(约x5)快得多,尽管仍然比完全矢量化的(约x3)慢.另一方面,我可以计算输入的值,用十(40,60,30,40)来缩放这个(在~5分钟内)而不是前一个因为内存(我没有测试如何)它需要与原始功能一起使用).使用@ numba.jit有点帮助,虽然不是很大(由于阶乘函数不能使用nopython).您可以尝试使用或多或少的循环向量化,具体取决于输入的大小.