如何加快应用 scipy 两个二维数组的积分？

Question

我正在使用 scipy.integrate 将积分函数应用于两个二维数组。这是示例：

from scipy import integrate
import numpy as np

def integrate_lno2(top, bottom, peak_height, peak_width):
    return integrate.quad(lambda x: np.exp(-np.power(x - peak_height, 2)/(2*np.power(peak_width, 2))), top, bottom)[0]

# change row and col to test speed
row = 100; col = 100; peak_height=300; peak_width=60

top = np.linspace(100, 200, row*col).reshape(row, col)
bottom = np.linspace(800, 900, row*col).reshape(row, col)

res = np.zeros((row, col))

for i in range(row):
    for j in range(col):
        res[i, j] = integrate_lno2(top[i, j], bottom[i, j], peak_height, peak_width)

如果二维数组的形状增加，for 循环可能会变慢。我找到了numba integrand example，但是它不接受上下限

Answer 1

就像在之前的回答中一样，您可以使用 Numba 来加速由于 Numpy 开销很大而非常慢的 lambda 调用（Numpy 未针对标量进行优化，因此执行此操作非常慢).更好的是：你可以告诉 Numba 生成一个 C 函数，它可以直接从 Scipy 调用，开销很小（因为它几乎完全消除了慢速 CPython 解释器的开销）。您还可以 pre-compute 除以变量并将其转换为乘法（更快）。

这是结果代码：

import numba as nb
import numpy as np
import scipy as sp

factor = -1.0 / (2 * np.power(peak_width, 2))

# change row and col to test speed
row = 100; col = 100; peak_height=300; peak_width=60

@nb.cfunc('float64(float64)')
def compute_numba(x):
    return np.exp(np.power(x - peak_height, 2) * factor)

compute_c = sp.LowLevelCallable(compute_numba.ctypes)

def integrate_lno2(top, bottom):
    return sp.integrate.quad(compute_c, top, bottom)[0]

top = np.linspace(100, 200, row*col).reshape(row, col)
bottom = np.linspace(800, 900, row*col).reshape(row, col)

res = np.zeros((row, col))

for i in range(row):
    for j in range(col):
        res[i, j] = integrate_lno2(top[i, j], bottom[i, j])

计算循环在我的机器上快了大约 100 倍。

如何加快应用 scipy 两个二维数组的积分？

How to speed up applying scipy integrate of two 2D arrays?

python

numpy

scipy

numba