2

There are currently a bunch of pages about efficiently calculating pairwise distance in MATLAB (some even by me!). I'm looking to do something a little different.

Rather than compute the overall distance between pairs of rows in two matrices, say X and Y:

X = [1 1 1; 2 2 2];
Y = [1 2 3; 4 5 6];

I want to calculate a 3 dimensional matrix, storing the raw column-wise difference between each pair of rows. In the above example, this matrix would have two rows (for the 2 observations in X), 3 columns, and 2 slices in the 3rd dimension (for the 2 observations in Y):

diffs(:,:,1) =

     0    -1    -2
     1     0    -1


diffs(:,:,2) =

    -3    -4    -5
    -2    -3    -4

I have so far come up with two ways accomplishing this as a general case, but I want to find something that is elegant, transparent, and efficient.

Repmat+Permute Approach

% Set up data
X = rand(100,10);
Y = rand(200,10);

timer = tic;
X_tiled = repmat(X,[1 1 size(Y,1)]);
Y_tiled = repmat(permute(Y,[3,2,1]),[size(X,1),1,1]);
diffs = X_tiled - Y_tiled;
toc(timer)
% Elapsed time is 0.001883 seconds.

For-Loop Approach

timer = tic;
diffs = zeros(size(X,1),size(X,2),size(Y,1));
for i = 1:size(X,1)
    for  j =1:size(Y,1)
        diffs(i,:,j) = X(i,:) - Y(j,:);
    end
end
toc(timer)
% Elapsed time is 0.028620 seconds.

Does anyone else have something better than what I've got?

Community
  • 1
  • 1
Nolan Conaway
  • 2,639
  • 1
  • 26
  • 42

1 Answers1

3

You could use bsxfun to replace that messy repmat for under-the-hood broadcasting after using permute on Y to send the first dimension to the third position keeping the second dimension in its place to match up with the second dimension of X. This could be achieved with permute(Y,[3 2 1]. Thus, the solution would be -

diffs = bsxfun(@minus,X,permute(Y,[3 2 1]))

Benchmarking

Benchmarking Code -

% Set up data
X = rand(100,10);
Y = rand(200,10);

% Setup number of iterations
num_iter = 500;

%// Warm up tic/toc.
for iter = 1:50000
    tic(); elapsed = toc();
end

disp('---------------- With messy REPMAT')
timer = tic;
for itr = 1:num_iter
    X_tiled = repmat(X,[1 1 size(Y,1)]);
    Y_tiled = repmat(permute(Y,[3,2,1]),[size(X,1),1,1]);
    diffs = X_tiled - Y_tiled;
end
toc(timer)

disp('---------------- With sassy BSXFUN')
timer = tic;
for itr = 1:num_iter
    diffs1 = bsxfun(@minus,X,permute(Y,[3 2 1]));
end
toc(timer)

Output -

---------------- With messy REPMAT
Elapsed time is 3.347060 seconds.
---------------- With sassy BSXFUN
Elapsed time is 0.966760 seconds.
Divakar
  • 218,885
  • 19
  • 262
  • 358