卷积中心差分

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阵列时，卷积IMHO将是一种过度杀伤力，因为切片阵列实际上在内存和性能上都是免费的。由于计算量大，可以研究numexpr以利用多核