实现此2D数值积分的计算上更快的方法是什么？

Question

我有兴趣进行2D数值积分。现在我正在使用scipy.integrate.dblquad，但是它非常慢。请参见下面的代码。我需要用完全不同的参数来评估这100次积分。因此，我想使处理尽可能快和高效。代码是：

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) * (1 / (np.sqrt(2 * np.pi) * 2)) * 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)

花费时间是

212.96751403808594

太多了。请提出一种更好的方法来实现我的目标。在来到这里之前，我曾尝试进行搜索，但没有找到任何解决方案。我读过quadpy可以更好，更快地完成这项工作，但是我不知道如何实现同样的目的。请帮忙。

Answer 1

通常，通过矩阵运算进行求和要比使用scipy.integrate.quad（或dblquad）快得多。您可以重写f（q，z，t）以获取aq，z和t向量，并使用np.tensordot返回f值的3D数组，然后将area元素（dtdz）乘以函数值并求和他们使用np.sum。如果您的area元素不是恒定的，则必须制作一个area-elements数组并使用np.einsum。要考虑积分限制，可以在汇总之前使用masked数组来遮盖积分限制之外的函数值。请注意，np.einsum会忽略掩码，因此，如果使用einsum，则可以使用np.where将积分限制之外的函数值设置为零。示例（具有恒定面积的元素和简单的积分限制）：

import numpy as np
import scipy.special as ss
import time

def f(q, t, z):

    # Making 3D arrays before computation for readability. You can save some time by
    # Using tensordot directly when computing the output
    Mq = np.tensordot(q, np.ones((len(t), len(z))), axes=0)
    Mt = np.tensordot(np.ones(len(q)), np.tensordot(t, np.ones(len(z)), axes = 0), axes = 0)
    Mz = np.tensordot(np.ones((len(q), len(t))), z, axes = 0)

    return Mt * 0.5 * (ss.erf((Mt - Mz) / 3) - 1) * (Mq * Mt) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
     -0.5 * ((Mz - 40) / 2) ** 2)

q = np.linspace(0.03, 1, 1000)
t = np.linspace(0, 50, 250)
z = np.linspace(10, 60, 250)

#if you have constand dA you can shave some time by computing dA without using np.diff
#if dA is variable, you have to make an array of dA values and np.einsum instead of np.sum
t0 = time.process_time()
dA = np.diff(t)[0] * np.diff(z)[0]
func_vals = f(q, t, z)
I = np.sum(func_vals * dA, axis=(1, 2))
t1 = time.process_time()

这在我的2012年macbook pro（2.5GHz i5）上耗时18.5秒，dA = 0.04。通过这种方式进行操作，还可以使您轻松地在精度和效率之间进行选择，并将dA设置为知道函数功能时有意义的值。

但是，值得注意的是，如果您想要更多的点，则必须拆分积分，否则就有可能使内存最大化（1000 x 1000 x 1000）加倍需要8GB内存。因此，如果您要进行非常高精度的大型集成，则值得在运行之前快速检查所需的内存。

Answer 2

您可以使用Numba或低级可拨打电话

几乎是您的示例

我只是将函数直接传递给scipy.integrate.dblquad，而不是使用lambda生成函数的方法。

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(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)
#143.73969149589539

这已经快了一点点（143 vs. 151s），但是唯一的用途是有一个简单的示例进行优化。

仅使用Numba编译功能

要使其运行，您还需要Numba和numba-scipy。 numba-scipy的目的是提供scipy.special中的包装函数。

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()

#error_model="numpy" -> Don't check for division by zero
@nb.njit(error_model="numpy",fastmath=True)
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)
#8.636585235595703

使用低级可调用项

scipy.integrate函数还提供了传递C回调函数而不是Python函数的可能性。这些函数可以用例如C，Cython或Numba编写，在本示例中将使用它们。主要优点是，在函数调用时无需Python解释程序交互。

@ Jacques Gaudin的出色answer显示了包括附加参数在内的简便方法。

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 = jit(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)
#3.2645838260650635