如何提高积分速度？

Question

我有一个想要集成的函数：

def f(z, t, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * np.exp(-0.5 * ((z - 40) / 2) ** 2)

我以这个函数为例来理解和展示不同集成方法之间的区别。根据通过该平台上的各种答案以及文档中获得的早期知识，我尝试了不同的方法来提高该功能的集成速度。我将在这里简要解释和比较它们：

1. 仅使用 python 和 scipy.quad：（需要 202.76 秒）

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def f(q, z, t):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t)  * np.exp(-0.5 * ((z - 40) / 2) ** 2)


y = np.empty([len(q)])
for n in range(len(q)):
    y[n] = integrate.dblquad(lambda t, z: f(q[n], z, t), 0, 50, lambda z: 10, lambda z: 60)[0]

end = time.time()
print(end - start)

#202.76 s

正如预期的那样，这不是很快。

2. 使用低级可调用并使用 Numba 编译函数：（需要 6.46 秒）

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb
from numba import cfunc
from numba.types import intc, CPointer, float64
from scipy import LowLevelCallable

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def jit_integrand_function(integrand_function):
    jitted_function = nb.njit(integrand_function, nopython=True)

    #error_model="numpy" -> Don't check for division by zero
    @cfunc(float64(intc, CPointer(float64)),error_model="numpy",fastmath=True)
    def wrapped(n, xx):
        ar = nb.carray(xx, n)
        return jitted_function(ar[0], ar[1], ar[2])
    return LowLevelCallable(wrapped.ctypes)

@jit_integrand_function
def f(t, z, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):
    return 10.

def upper_inner(z):
    return 60.


y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)

#6.46 s

这比以前的方法更快。

3. 使用scipy.quad_vec：（需要 2.87 秒）

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def f(z, t, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * np.exp(-0.5 * ((z - 40) / 2) ** 2)

ans2 = integrate.quad_vec(lambda t: integrate.quad_vec(lambda z: f(z,t,q),10, 60)[0],0,50)[0]

end = time.time()
print(end - start)
#2.87s

涉及的方法scipy.quad_vec是最快的。我想知道如何让它更快？scipy.quad_vec不幸的是不接受低级可调用函数；有没有办法像 一样有效地矢量化scipy.quad或函数？或者我可以使用其他方法吗？任何形式的帮助将非常感激。scipy.dblquadscipy.quad_vec

Answer 1

这是速度提高 10 倍的版本

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def fa(z, t, q):
    return (erf((t - z) / 3) - 1) * np.exp(-0.5 * ((z - 40) / 2) ** 2)

ans3 = 0.5 * integrate.quad_vec(lambda t: t * j0(q * t) * 
         integrate.quad_vec(lambda z: fa(z,t,q),10, 60)[0],0,50)[0]

end = time.time()
print(end - start)

这里的区别在于去掉j0(q*t)了内部积分。