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

Question

我需要在二维网格（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 分钟。现在我想减少那些计算中的一些冗余以加快计算速度。即：

由于网格以 r = 0 为中心，因此网格上的许多 r 值都是相同的。由于对称性，所有值中只有约 1/8 是唯一的（对于方形网格）。一个想法是只计算唯一值的积分（通过 np.unique 找到），然后将它们保存在查找 table（哈希表？）中，或者我可以缓存函数值，以便只有新的计算值（通过@lru_cache）。但是，当我之后对函数进行矢量化时，这真的有效吗？
随着积分从 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 所指出的那样。

As the grid is centered around r = 0, many values for r on the grid are the same. Due to symmetry only ~1/8 of all values are unique (for a square grid). One idea was to calculate the integral only for the unique values (found via np.unique) and then save them in a look-up table (hashmap?) Or I could cache the function values so that only new values are calculated (via @lru_cache). But does that actually work when I vectorize the function afterwards?

虽然不重新计算值确实是个好主意，但我不希望这种方法明显更快。请注意，hashmap 和 np.unique 一样非常慢。我建议您仅 select 输入数组的四分之一 R。类似于 R[0:R.shape[0]//2, 0:R.shape[1]//2]。形状不对要小心

As the integral goes from 0 to r, the integral is often calculating integrals over intervals it has already calculated. E.g. if you calculate from 0 to 1 and afterwards from 0 to 2, only the interval from 1 to 2 is "new". But what would be the best way to utilize that? And would that even be a real performance boost using scipy.integrate.quad?

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

的可选参数降低计算 sub-intervals 的精度

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

Reducing redundancy for calculating large number of integrals numerically

python

optimization

integral