我怎样才能加快 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 Function. scipy includes the scipy.special.gammainc 函数，用于逼近下不完全 Gamma 函数，由（完整）Gamma 函数正则化（即：其输出除以 Gamma( 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 编译。这运行s 在我的机器上快了 ~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 运行ge [a, b] 上的积分等价于 [a, c] 上的积分加上 [c, b] 上的积分，对于a和b之间的一些c。因此，我们不是每次都计算 (0, b * t) 之间的积分，而是只计算 (b * prev_t, b * t) 之间的积分（其中 prev_t 是我们使用的最后一个 t )，然后运行一个累加和。这只会让它更快，因为在较小的运行ge.

上执行积分所需的近似迭代次数要少得多

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

我怎样才能加快 scipy.integrate.quad？

How can I speed up scipy.integrate.quad?

python

performance

scipy

numerical-integration

如果你只是想集成你在这里给出的功能：

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