performance - 如何向量化这个五点差代码以快速计算矩阵的导数？

Question

我正在尝试使用五点法计算关于二维矩阵中两个维度的偏一阶导数，即 dF/dx 和 dF/dy 。我通过遍历这些点成功地做到了这一点：

dF_dx = zeros(size(F));
dF_dy = zeros(size(F));

% Derivative with respect to y for each x value (Apply to all columns simultaneously)
dF_dy(2,1:(Nx-1)) = ( F(3,1:(Nx-1)) - F(1,1:(Nx-1)) )/(2*dy);
for m = 3:(Ny-2)
    dF_dy(m,1:(Nx-1)) = ( F(m-2,1:(Nx-1)) - F(m+2,1:(Nx-1)) + 8*F(m+1,1:(Nx-1)) - 8*F(m-1,1:(Nx-1)) )/(12*dy);
end
dF_dy(Ny-1,1:(Nx-1)) = ( F(Ny,1:(Nx-1)) - F(Ny-2,1:(Nx-1)) )/(2*dy);

% Derivative with respect to x for each y value (Apply to all rows simultaneously)
dF_dx(2:(Ny-1),2) = ( F(2:(Ny-1),3) - F(2:(Ny-1),1) )/(2*dx);
for n = 3:(Nx-2)
    dF_dx(2:(Ny-1),n) = ( F(2:(Ny-1),n-2) - F(2:(Ny-1),n+2) + 8*F(2:(Ny-1),n+1) - 8*F(2:(Ny-1),n-1) )/(12*dx);
end
dF_dx(2:(Ny-1),(Nx-1)) = ( F(2:(Ny-1),Nx) - F(2:(Ny-1),Nx-2) )/(2*dx);

是否有一种聪明的 Matlab 方法可以对其进行矢量化以使其执行速度更快？我已经阅读了几个似乎可以解决问题的函数（例如diff()、circshift()，甚至可能是kron ()），但不确定如何使用它们来解决这个问题。

（我的计划是稍后将其实现为 gpuArray，如果这与解决方案相关）。

谢谢！

---

第一次编辑

通过查看gradient()的源代码，我能够制作以下矢量化版本（即没有循环）：

F = rand(5,8);
Nx = size(F,2);
Ny = size(F,1);
dx = 2;
dy = 3;

dF_dx = zeros(size(F));
dF_dy = zeros(size(F));

dF_dx(2:(Ny-1),3:Nx-2) = (F(2:(Ny-1),1:Nx-4) - F(2:(Ny-1),5:Nx) + 8*F(2:(Ny-1),4:Nx-1) - 8*F(2:(Ny-1),2:Nx-3))/(12*dx);
dF_dx(2:(Ny-1),2) = ( F(2:(Ny-1),3) - F(2:(Ny-1),1) )/(2*dx);
dF_dx(2:(Ny-1),(Nx-1)) = ( F(2:(Ny-1),Nx) - F(2:(Ny-1),Nx-2) )/(2*dx);

dF_dy(3:Ny-2,1:(Nx-1)) = (F(1:Ny-4,1:(Nx-1)) - F(5:Ny,1:(Nx-1)) + 8*F(4:Ny-1,1:(Nx-1)) - 8*F(2:Ny-3,1:(Nx-1)))/(12*dy);
dF_dy(2,1:(Nx-1)) = ( F(3,1:(Nx-1)) - F(1,1:(Nx-1)) )/(2*dy);
dF_dy((Ny-1),1:(Nx-1)) = ( F(Ny,1:(Nx-1)) - F(Ny-2,1:(Nx-1)) )/(2*dy);

有没有办法在没有所有索引的情况下做到这一点（我怀疑这是最耗时的）？

---

第二次编辑

我现在已经实现了 Cris 和 chtz 的建议，并使用与 的卷积实现了这一点conv2()，如下所示：

F = rand(500,600);
Nx = size(F,2);
Ny = size(F,1);
dx = 2;
dy = 3;

 [dF_dx, dF_dy] = partial_derivatives(F, dx, dy, Nx, Ny)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [dF_dx, dF_dy] = partial_derivatives(F, dx, dy, Nx, Ny)

kernx = [-1 8 0 -8 1]/(12*dx); % Convolution kernel for x dimension
kerny = [-1;8;0;-8;1]/(12*dy); % Convolution kernel for y dimension

%%%%%%%% Partial derivative across x dimension %%%%%%%%
dF_dx = conv2(F, kernx, 'same') ; % Internal mesh points (five-point method)

% Second and penultimate mesh points (two-point method)
dF_dx(2:(Ny-1),2) = ( F(2:(Ny-1),3) - F(2:(Ny-1),1) )/(2*dx);
dF_dx(2:(Ny-1),(Nx-1)) = ( F(2:(Ny-1),Nx) - F(2:(Ny-1),Nx-2) )/(2*dx);
dF_dx(:,[1 Nx]) = 0; % Set boundary conditions
dF_dx([1 Ny],:) = 0; %


%%%%%%%% Partial derivative across x dimension %%%%%%%%
dF_dy = conv2(F, kerny, 'same') ; % Internal mesh points (five-point method)

% Second and penultimate mesh points (two-point method)
dF_dy(2,1:(Nx-1)) = ( F(3,1:(Nx-1)) - F(1,1:(Nx-1)) )/(2*dy);
dF_dy((Ny-1),1:(Nx-1)) = ( F(Ny,1:(Nx-1)) - F(Ny-2,1:(Nx-1)) )/(2*dy);
dF_dy(:,Nx) = 0;     % Set boundary conditions
dF_dy([1 Ny],:) = 0; %

end

作为 gpuArray，与第 1 版中的矢量化版本相比，这并没有提供任何明显的改进。有没有明显的改进方法？谢谢。

标签： performancematlabmatrixvectorizationdifference