如何加快 scipy.integrate.quad？答案

【问题标题】：How can I speed up scipy.integrate.quad?如何加快 scipy.integrate.quad？
【发布时间】：2021-03-23 22:37:42
【问题描述】：

考虑以下代码：

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.85s）。如何加快速度？

【问题讨论】：

不要使用np.vectorize。它并不比列表理解快。至于quad 本身，除了确保disFunc 快速之外，您无能为力。我不认为quad 有转向参数，但我可能错了。

标签： python performance scipy numerical-integration

【解决方案1】：

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

请注意，您希望集成的功能实际上等同于Lower Incomplete Gamma Function。 scipy 包括 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 编译。这在我的机器上运行速度快了约 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]），因此您可能需要单独计算并将其添加到数组的前面。

【讨论】：

你的递归封闭形式很棒，但我应该提到 n 可以带小数。是否可以修改您的 integration_by_parts 函数来处理 n 可以有小数的情况（例如，5.5）？
@user1363251 啊，不幸的是，对于非自然数的 n 值，基本函数中不存在积分的封闭形式，因此无法修改函数以适应这一点，需要数值积分.我建议使用scipy.special.gamma 来执行近似，而不是使用集成模块，因为它可能专门针对这个积分进行了更好的优化。
更好的是：lower incomplete gamma function 应该正是您要查找的内容，您只需将其乘以完整的 gamma 函数即可。
@user1363251 我已经更新了我的答案，以展示如何使用下不完全伽马函数计算积分。

【解决方案2】：

首先，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！

【讨论】：