如何加速scipy.integrate.quad？

Question

考虑下面的代码：

def disFunc(x,n):
        y = np.power(x, n-1)*np.exp(-x)
        return y
    

def disFunc2(t,n,beta):
    y = integrate.quad(disFunc, a = 0, b = beta*t, args = (n))
    return y[0]   

def disFunc3(func):
    y = np.vectorize(func)
    return(y)

maxCounter = 6000
stepSize = .001
n = 5
beta = 25
t = np.cumsum(np.ones(maxCounter)*stepSize)-stepSize
x = disFunc3(disFunc2)
start_time = time.time()
y = x(t,n,beta)
time_taken = (time.time() - start_time)
print (time_taken)

工作起来像个符咒，但速度太慢了（1.85秒）。我怎样才能加快速度

Answer 1

首先，np.vectorize通常非常慢。也许您可以尝试更改代码，使其不必依赖于此

否则，您可以尝试Numba/Numba-Scipy的jit（即时编译）。这大大加快了使用NumPy和Scipy的代码的速度： http://numba.pydata.org/

然后可以导入jit装饰器，并在函数上使用它：

from numba import jit

@jit  # or jit(nopython=False) in other scenarios
def disFunc(x,n):
        y = np.power(x, n-1)*np.exp(-x)
        return y
    
> ...

确保同时安装Numba和Numba Scipy

Answer 2

如果您只想集成此处提供的功能：

请注意，您希望集成的函数实际上相当于Lower Incomplete Gamma Functionscipy包括^{}函数，用于近似较低的不完全伽马函数，由（完全）伽马函数正则化（即：其输出除以伽马（n））。因此，我们可以使用这些特殊函数更有效地近似积分：

from scipy import special
import time

# ...

start_time = time.time()
y = special.gammainc(n, beta * t)
y *= special.gamma(n)  # Reverse the regularisation
time_taken = (time.time() - start_time)
print(time_taken)
print(y)

对于不同的n值，所需的时间会有所不同，但对于n=5.5，在我的机器上运行的时间约为0.0005秒，而下面的方法为0.3秒

如果要集成不存在闭合形式的任意函数：

这里有一个想法，它只是改进了您在数学上使用的方法，而不是使用JIT编译。这在我的机器上运行速度快了约10倍（您的代码为我运行大约需要2.2秒，这需要约0.3秒）：

import numpy as np
from scipy import integrate
import time


def disFunc(x, n):
    y = np.power(x, n - 1) * np.exp(-x)
    return y


def disFunc2(t_prev, t_cur, n, beta):
    y = integrate.quad(disFunc, a=beta * t_prev, b=beta * t_cur, args=(n))
    return y[0]


def disFunc3(func):
    y = np.vectorize(func)
    return y


maxCounter = 6000
stepSize = 0.001
n = 5
beta = 25
t = np.cumsum(np.ones(maxCounter) * stepSize) - stepSize
x = disFunc3(disFunc2)
start_time = time.time()
y = x(t[:-1], t[1:], n, beta)
time_taken = time.time() - start_time
print(time_taken)
print(np.cumsum(y))

这个想法利用了积分的性质：对于a和b之间的一些c，范围[a，b]上的积分等于[a，c]上的积分加上[c，b]上的积分。因此，我们不是每次计算（0，b*t）之间的积分，而是计算（b*prev\t，b*t）之间的积分（其中prev\t是我们使用的最后一个t），然后运行一个累积和。这只会使它更快，因为在较小范围内执行积分所需的近似迭代次数要少得多

应该注意的是，这确实跳过了第一个积分（从0到beta*t[0]），因此您可能希望单独计算该积分并将其添加到数组的前面