Efficient computation of the sum of pairwise absolute differences

Hi,

I am working on a project that requires the computation of the sum of all pairwise absolute differences between elements at either end of a randomly placed vector with coordinates . For instance, starting with the following 3x3 matrix:

A = magic(3)

A = 3×3

8 1 6 3 5 7 4 9 2

A vector with coordinates will produce the absolute differences:

b = sum(abs([A(1,1)-A(1,2),A(1,2)-A(1,3),A(2,1)-A(2,2),A(2,2)-A(2,3),...
    A(3,1)-A(3,2),A(3,2)-A(3,3)]))

b = 28

Note that I took the differences between all "horizontal" matrix elements separated by a distance of .

The idea is to save this value in a matrix, let's say, C, in the position corresponding to the vector coordinate differences, i.e., with respect to the center of the matrix A. Thus, if we were to do the same with the vector defined by , we get:

C = zeros(size(A));
C(1,3) = sum(abs([A(2,1)-A(1,2),A(2,2)-A(1,3),A(3,1)-A(2,2),A(3,2)-A(2,3)]));

And adding in the previously calculated value:

C(2,3) = b

C = 3×3

0 0 6 0 0 28 0 0 0

As a proof of concept, I borrowed and adapted a function from another answer that does the computation. I have called it diffMap and is at the bottom of this question/script.

dMap = diffMap(A)

dMap = 5×5

6 2 16 6 2 2 18 20 6 6 8 28 0 28 8 6 6 20 18 2 2 6 16 2 6

As one can probably easily notice, this is just a modification of a rudimentary for-loopy cross-correlation.

I have been giving some thought to this problem in the hope of making this work with xcorr2 and discrete Fourier transforms. This is what I have got until now:

h = ones(size(A));
dMap_ = abs(xcorr2(A,h)-xcorr2(h,A))

dMap_ = 5×5

6 2 0 6 2 2 6 0 2 6 0 0 0 0 0 6 2 0 6 2 2 6 0 2 6

It would be great if this worked, because it could be easily adapted to work with fft2, ifft2, and fftshift thanks to the convolution theorem. However, this does not return what i was hoping for: it computes the absolute value of the sum of the differences and NOT the sum of the absolute differences.

The question: is there any way to make this work efficiently, ideally with conv2 or xcorr2? With big matrices four for-loops takes forever, obviously.

Thanks!!

Here's the diffMap function:

%% Function
function dMap = diffMap(A) 
%DIFFMAP Compute the absolute differences map

% Get the size of the matrix A
[r,c] = size(A);
% Initialize and index matrix
h = ones([r,c]);
% Padding
Rep = zeros(r + r*2-2, c + c*2-2);
Rep_ = zeros(r + r*2-2, c + c*2-2);
for x = r : r+r-1
    for y = c : c+r-1
        Rep(x,y) = A(x-r+1, y-c+1);
        Rep_(x,y) = h(x-r+1, y-c+1);
    end
end
dMap = zeros(r+r-1,c+c-1);
for x = 1 : r+r-1
    for y = 1 : c+c-1
        for i = 1 : r
            for j = 1 : c
                dMap(x, y) = dMap(x, y) + abs((Rep(x+i-1, y+j-1) * h(i, j)) -...
                    (Rep_(x+i-1, y+j-1) * A(i, j)));
            end
        end
    end
end
end

Answer 1

[Matt J]

I think parallelization might be the only way to accelerate things. Parfor seems to work well with the rewritten version of diffMap() below. With 12 cores, I can do a 250x250 matrix in about 2.7 sec., whereas the original diffMap code takes about 60 sec.

A = magic(250);

tic;
 dMap = diffMap(A);
toc

%% Function
function dMap = diffMap(A) 
%DIFFMAP Compute the absolute differences map

 [r,c] = size(A);

 [D1,D2]=deal(zeros(r,c));

  parfor k=1:r*c
      
       [i,j]=ind2sub([r,c],k);
        
        if i==1 && j==1, continue; end

        kernel=zeros(i,j);
        kernel([1,end])=[-1,1];
        
        D1(k)=sum(abs(conv2(A,kernel,'valid')),'all');
        D2(k)=sum(abs(conv2(A,flipud(kernel),'valid')),'all');
   end
  
   dMap=[flipud(D2);D1];
   dMap(end/2,:)=[];
   dMap=[rot90(dMap,2),dMap];
   dMap(:,end/2)=[];
  
   
end

[Santiago Benito]

Thank you very much for this version of the code! I tried replacing conv2 with conv2fft by Bruno Luong, but that actually takes longer. Probably because the matrices for the convolution computation in the initial iterations of the parfor loop are really small, or because that function is not really parallelizable.

In my current setup, a 480x680 matrix takes around 340 seconds to compute with conv2 and with conv2fft it took so long, that I had to stop it after 4000 seconds. For more info, I have four cores (AMD Ryzen 3 3200G).

Answer 2

[Matt J]

The question: is there any way to make this work efficiently, ideally with conv2 or xcorr2?

If you were taking the sum of squared differences, it could be done with xcorr2 in just a few lines:

W=ones(size(A));
tmp=xcorr2(A.^2,W);
dMap=tmp+rot90(tmp,2)-2*xcorr2(A);

However, because you are taking the sum of absolute differences, there is no simple connection with xcorr2.

With big matrices four for-loops takes forever, obviously.

Is there an upper bound on the size of the matrix that you will need to process?

[Santiago Benito]

Hello Matt, thank you very much for all of your answers, I really appreciate it!

I just read your comments and functions and tried some things on my own. I will post comments to your answers separately.

Firstly, I have at the moment matrices ranging from 480x640 to around 1000x800; they are relatively big for this kind of thing. 1200x1200 might be as big as it gets in this case.

Yeah, I wish I could stick to squared differences. I have that working relatively quickly already. I was trying to find a way to either:

• Compute the absolute differences independently but also using convolutions.
• Compute the absolute differences as a derivation of the squared differences map.

Arguably, these two options can be seen as the same thing, though.

Answer 3

[Matt J]

This GPU implementation may also be useful. I was able to process a 500x500 matrix in 20 seconds on the GTX 1080 Ti.

A = magic(500);
gd=gpuDevice;
tic;
dMap = diffMap(A);
wait(gd);
toc

%% Function
function dMap = diffMap(A) 
%DIFFMAP Compute the absolute differences map

% Get the size of the matrix A
[r,c] = size(A);

A=gpuArray(A);
T=gpuArray(toeplitz(1:r,[1,zeros(1,r-1)]));
t=nonzeros(T);
W=reshape(logical(T),r,1,[]);
R=double(W);

% Initialize and index matrix
% Padding
 
dMap = gpuArray.zeros(r,2*c-1);

 for j=1:c
 
     col=A(:,j);
     R(W)=col(t);
     tmp=permute(sum( abs(R-A).*W,1)  ,[3,2,1]);
     
     %ctr=c
     dMap(:,c-(j-1):2*c-j) = dMap(:,c-(j-1):2*c-j) +tmp;
    
 end

 
 dMap=[rot90(dMap,2);dMap];
 dMap(end/2,:)=[];
end

[Santiago Benito]

Again, thanks for the help! This function is, because of my setup, slower for me than the parallel one; I only have a GT 710. Just for reference, it takes 610 seconds to run, a little less than double the time of the parfor version.