与卷积的中心差异

Question

所以基本上我试图在二维数组上做有限差分而不做太多 for 循环。我想要数组的 Hessian 矩阵和梯度。所以我需要数组的一阶和二阶导数。

这可以通过计算数组上的以下等式来实现。

为了处理边界，我们只为内部点计算边界，因此此导数的代码可能类似于以下内容

arr = np.random.rand(16).reshape(4,4)
result = np.zeros_like(arr)
w, h = arr.shape
for i in range(1, w-1):
    for j in range(1, h-1):
        result[i,j] = (arr[i+1, j] - arr[i-1, j]) / (2*dx)

这给出了正确答案，但与 nu numpy 操作相比可能非常慢，所以我心想。这基本上只是与内核的卷积，看起来像这样

kernel = [1, 0 , -1]

所以我们执行下面的代码

from scipy.sigmal import convolve
result = np.pad((convolve(arr,kernel,mode='same', 
             method = 'direct')/(2*dx))[1:-1, 1:-1], 1).T

因为我们只处理内部点，所以我们将它们切掉并在之后用零填充，以模仿在之前的简单情况下会发生什么。

这有效！但是对于某些数组，朴素案例和卷积案例之间的均方误差飙升。因此，在某些情况下，数值误差似乎会增加很多。

我想要通过卷积获得的速度与原始案例的稳定性。有什么帮助吗？

Answer 1

我们可以简单的切片操作。因此，在输出初始化之后，执行 -

result[1:-1,1:-1] = (arr[2:,1:-1] - arr[:-2,1:-1])/(2*dx)

卷积恕我直言，在使用 NumPy 数组时会有点矫枉过正，因为切片数组实际上不占用内存和性能。由于计算量大，可以研究 numexpr 以利用多核。

与卷积的中心差异

Central difference with Convolution

python

numpy

vectorization

python-3.x