Fastest way to differentiate only the sides of large matricies

I believe I have found one of the fastest ways to differentiate different portions of a large matrix using built in Matlab functions. However, there is a problem to my method, there is a remaining middle value between the side values. Allow me to demonstate with some simple code which makes use of gpuArray. For my case, I need the data on a gpu.

N = 300; % Number of points in x,y,z dimensions 
S = 10;  % Number of points from side of a dimension

for i = 1:20
B  = gpuArray(rand(N,N,N));      % Load some random data
An = diff(B(1:S+1,:,:),1,1);     % Differentation from left side
Ap = diff(B(end-S:end,:,:),1,1); % Differentation from right side

A1 = B([2:S+1,N-S+1:N],:,:)-B([1:S,N-S:N-1],:,:);
% The above is the finite difference which is applied to 
% both sides of the matrix at the same time. This method
% however, is slower than using built in diff

A2 = diff(B([1:S+1,N-S:N],:,:),1,1);
% The above method is the fastest method, however there is one
% additional data point between S+1 and N - S. Deleting this
% point would add extra computational time
end

However diff is coded on the backend of Matlab, and would need to be modified for my application. If this is true, how would I go about learning how to code up my own diff function for this special case?

Thank you

Answer 1

[Matt J]

This seems to be about as fast, at least on the GTX 1080 Ti.

T=reshape(B([1:S+1,N-S:N],:) ,S+1,2*N^2 );

A3 = reshape( diff(T,1,1), 2*S,N,N); 

[Matt J]

Here was my speed test:

N = 300; 
S = 10; 

B  = gpuArray.rand(N,N,N);     

gputimeit(@() baseLine(B,S,N))  %2.7684e-04 seconds

gputimeit(@() alternative(B,S,N)) %2.8324e-04 seconds
          
function baseLine(B,S,N)
 A2=diff(B([1:S+1,N-S:N],:,:),1,1);
end

function  alternative(B,S,N)
 T=reshape(B([1:S+1,N-S:N],:) ,S+1,2*N^2 );
 A3 = reshape( diff(T,1,1), 2*S,N,N); 
end

[Nathan Zechar]

@Matt J

Wow....This works, and it doesn't eat up a lot of memory to do this. It's kinda late and I wanted to respond before I pass out and say thanks. Tomorrow I'm going to have to sit down and wrap my head around what you just did here.

[Nathan Zechar]

@Matt J

Okay, so I believe I understand how to break down higher dimensional matracies using reshape, but I've come into a bit of a road bump when it comes to applying your technique along a different dimensional axis.

Say for instance I want to apply this along the 2nd or 3rd axis, this would require a permute, your technique, and then a permute back.

% along 1st axis (original
T  = reshape(B([1:S+1,N-S:N],:) ,S+1,2*N^2 );
A1 = reshape(diff(T,1,1), 2*S,N,N);

% along 2nd axis
B  = permute(B,[2 1 3])
T  = reshape(B([1:S+1,N-S:N],:) ,S+1,2*N^2 );
A  = reshape(diff(T,1,1), 2*S,N,N);
A2 = permute(A,[1 2 3])

% along 3rd axis
B  = permute(B,[3 2 1])
T  = reshape(B([1:S+1,N-S:N],:) ,S+1,2*N^2 );
A  = reshape(diff(T,1,1), 2*S,N,N);
A3 = permute(A,[1 2 3])  

These permutes, are computationally expensive. I wasn't sure if there was a better approach than permute. Additionally, I'm essentially applying something similiar to the built in grad function, however centeral differences are not being applied and there isn't an augmented algorithim for the ends of the matrix. This leads me to believe I might be able to achieve even more gains in efficiency.

[Matt J]

No, you don't need to permute.

% along 1st axis
T  = reshape(B([1:S+1,N-S:N],:) ,S+1,2,N,N );
A1 = reshape(diff(T,1,1), 2*S,N,N);

% along 2nd axis
T=reshape(B(:,[1:S+1,N-S:N],:),N,S+1,2,N);
A2=reshape( diff(T,1,2)  , N,2*S,N);

% along 3rd axis
T=reshape(B(:,:,[1:S+1,N-S:N]),N,N,S+1,2);
A3=reshape( diff(T,1,3)  , N,N,2*S);

[Nathan Zechar]

@Matt J, thank you again. I was hung up on what was happening on the inside of the first line in the original method.

B([1:S+1,N-S:N],:)

Reduces to 2 dimensions. I thought that was part of the magic and had to permute to that form to get everything to work. Thank you so much again.