减少数值计算大量积分的冗余

Question

我需要在 2D 网格（x，y 位置）上计算以下积分：

与 r = sqrt(x^2 + y^2) 且二维网格以 x=y=0 为中心。 实现很简单：

import numpy as np
from scipy import integrate

def integralFunction(x):
    def squareSaturation(y):
        return np.sqrt(1-np.exp(-y**2))
    return integrate.quad(squareSaturation,0,x)[0]

#vectorize function to apply function with integrals on np-array
integralFunctionVec = np.vectorize(integralFunction)

xmax = ymax = 5
Nx = Ny = 1024

X, Y = np.linspace(-xmax, xmax, Nx), np.linspace(-ymax, ymax, Ny)
X, Y = np.meshgrid(X, Y)

R = np.sqrt(X**2+Y**2)
Z = integralFunctionVec(R)

但是，我目前正在处理 1024x1024 网格，计算大约需要 1.5 分钟。现在这些计算中有一些冗余，我想减少这些冗余以加快计算速度。即：

1. 由于网格以 r = 0 为中心，因此网格上的许多 r 值是相同的。由于对称性，所有值中只有约 1/8 是唯一的（对于方形网格）。一种想法是仅计算唯一值的积分（通过 np.unique 找到），然后将它们保存在查找表中（哈希图？）或者我可以缓存函数值，以便只计算新值（通过 @ lru_缓存）。但是当我之后对函数进行矢量化时，这实际上有效吗？

2. 当积分从 0 到 r 时，积分通常会在已计算的区间内计算积分。例如，如果您从 0 到 1 进行计算，然后再从 0 到 2 进行计算，则只有从 1 到 2 的区间是“新的”。但利用它的最佳方法是什么？使用 scipy.integrate.quad 是否能真正提升性能？

您有任何反馈或其他想法来优化此计算吗？

Answer 1

您可以使用Numba来加速quad的计算。下面是一个例子：

import numpy as np
import numba as nb
from scipy import integrate

@nb.cfunc('float64(float64)')
def numbaSquareSaturation(y):
    return np.sqrt(1-np.exp(-y**2))

squareSaturation = scipy.LowLevelCallable(numbaSquareSaturation.ctypes)

def integralFunction(x):
    return integrate.quad(squareSaturation,0,x)[0]

integralFunctionVec = np.vectorize(integralFunction)

xmax = ymax = 5
Nx = Ny = 1024

X, Y = np.linspace(-xmax, xmax, Nx), np.linspace(-ymax, ymax, Ny)
X, Y = np.meshgrid(X, Y)

R = np.sqrt(X**2+Y**2)
Z = integralFunctionVec(R)

在我的机器上这大约快了 25 倍。该代码仍然不是最佳的，因为 squareSaturation 调用会带来很大的开销，但 SciPy 似乎没有提供一种针对您的情况有效矢量化 quad 的方法。请注意，使用 nb.cfunc+scipy.LowLevelCallable 可以显着加快执行速度，如 @max9111 所指出的。

由于网格以 r = 0 为中心，因此网格上的许多 r 值是相同的。由于对称性，所有值中只有约 1/8 是唯一的（对于方形网格）。一种想法是仅计算唯一值的积分（通过 np.unique 找到），然后将它们保存在查找表中（哈希图？）或者我可以缓存函数值，以便只计算新值（通过 @ lru_缓存）。但是当我之后对函数进行矢量化时，这实际上有效吗？

尽管不重新计算值确实是一个好主意，但我并不认为这种方法会明显更快。请注意，hashmap 和 np.unique 一样非常慢。我建议您只选择输入数组 R 的四分之一。类似于 R[0:R.shape[0]//2, 0:R.shape[1]//2]。如果形状奇怪，请小心。

当积分从 0 到 r 时，积分通常会在已计算的区间内计算积分。例如，如果您从 0 到 1 进行计算，然后再从 0 到 2 进行计算，则只有从 1 到 2 的区间是“新的”。但利用它的最佳方法是什么？使用 scipy.integrate.quad 是否能真正提升性能？

这可能会有所帮助，因为积分的域更小并且函数应该更平滑。这意味着 Scipy 的计算速度应该更快。即使它不会自动执行此操作，您也可以使用 quad。