代码之家 › 专栏 › 技术社区 › Oskar Hofmann

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

integral optimization python

Oskar Hofmann · 技术社区 · 4 年前

我需要在二维网格上计算以下积分(x,y位置):

$\int_{0}^{r}\sqrt{1-e^{-s^2}} ds$

其中r=sqrt(x^2+y^2)并且2D网格以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找到),然后将它们保存在查找表中(hashmap?)或者我可以缓存函数值,这样只计算新值(通过@lru_cache)。但是,当我随后对函数进行矢量化时,这真的有效吗?
当积分从0到r时,积分通常是在已经计算过的区间上计算积分。例如,如果从0到1计算,然后从0到2计算,则只有从1到2的间隔是“新的”。但是,最好的利用方式是什么呢?使用scipy.integrate.quad,这真的会提高性能吗?

您是否有任何反馈或其他想法来优化此计算?

0 回复 | 直到 4 年前

Jérôme Richard 4 年前

您可以使用 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似乎没有提供向量化的方法 四边形 有效地处理您的案件。请注意,使用 nb.cfunc + scipy.LowLevelCallable 如@max9111所指出的,显著加快了执行速度。

由于网格以r=0为中心,因此网格上r的许多值都是相同的。由于对称性,所有值中只有约1/8是唯一的(对于方形网格)。一个想法是只计算唯一值的积分(通过np.unique找到),然后将它们保存在查找表中(hashmap?)或者我可以缓存函数值,这样只计算新值(通过@lru_cache)。但是,当我随后对函数进行矢量化时,这真的有效吗?

虽然不重新计算值确实是个好主意,但我不希望这种方法会更快。请注意,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 。