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

Question

我有兴趣进行二维数值积分。现在我正在使用 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) 以获取 q、z 和 t 向量和 return 使用 np.tensordot 的 f 值的 3D 数组，然后乘以您的面积元素 (dtdz)使用函数值并使用 np.sum 对它们求和。如果您的面积元素不是常量，则必须制作面积元素数组并使用 np.einsum 要考虑积分限制，您可以使用屏蔽数组在汇总之前屏蔽积分限制之外的函数值。请注意，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 对 151s），但唯一的用途是有一个简单的例子来优化。

使用Numba简单编译函数

要达到运行，您还需要 Numba and 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 的出色展示了一种简单的方法来执行此操作，包括附加参数。

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)
#3.2645838260650635

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

What would be the computationally faster way to implement this 2D numerical integration?

python

integration

numpy

scipy

numerical-integration

您可以使用 Numba 或低级调用