pdist2 equivalent in MATLAB version 7

Question

I need to calculate the euclidean distance between 2 matrices in matlab. Currently I am using bsxfun and calculating the distance as below( i am attaching a snippet of the code ):

for i=1:4754
test_data=fea_test(i,:);
d=sqrt(sum(bsxfun(@minus, test_data, fea_train).^2, 2));
end

Size of fea_test is 4754x1024 and fea_train is 6800x1024 , using his for loop is causing the execution of the for to take approximately 12 minutes which I think is too high. Is there a way to calculate the euclidean distance between both the matrices faster?

I was told that by removing unnecessary for loops I can reduce the execution time. I also know that pdist2 can help reduce the time for calculation but since I am using version 7. of matlab I do not have the pdist2 function. Upgrade is not an option.

Any help.

Regards,

Bhavya

score 12 · Answer 1 · answered Oct 14 '11 at 22:41

12

Here is vectorized implementation for computing the euclidean distance that is much faster than what you have (even significantly faster than PDIST2 on my machine):

D = sqrt( bsxfun(@plus,sum(A.^2,2),sum(B.^2,2)') - 2*(A*B') );

It is based on the fact that: ||u-v||^2 = ||u||^2 + ||v||^2 - 2*u.v

Consider below a crude comparison between the two methods:

A = rand(4754,1024);
B = rand(6800,1024);

tic
D = pdist2(A,B,'euclidean');
toc

tic
DD = sqrt( bsxfun(@plus,sum(A.^2,2),sum(B.^2,2)') - 2*(A*B') );
toc

On my WinXP laptop running R2011b, we can see a 10x times improvement in time:

Elapsed time is 70.939146 seconds.        %# PDIST2
Elapsed time is 7.879438 seconds.         %# vectorized solution

You should be aware that it does not give exactly the same results as PDIST2 down to the smallest precision.. By comparing the results, you will see small differences (usually close to eps the floating-point relative accuracy):

>> max( abs(D(:)-DD(:)) )
ans =
  1.0658e-013

On a side note, I've collected around 10 different implementations (some are just small variations of each other) for this distance computation, and have been comparing them. You would be surprised how fast simple loops can be (thanks to the JIT), compared to other vectorized solutions...

answered Oct 14 '11 at 22:41

Amro

123,847
25
243
454

1

@Alex: thanks. I see your solution is basically the same, only I avoid creating temp matrices in memory by using BSXFUN (instead of using REPMAT and the like) – Amro Jun 29 '12 at 17:27
Agree, bsxfun is one of the underrated and unknown functions in MATLAB. It doesn't waste memory like repmat. More users need to learn how to use it when 'vectorizing' their code. – Alexey Jun 29 '12 at 19:31
Just a quick note, I compared both PDIST2 and your method using bsxfun and serially bsxfun is much faster. However, inside of a parfor loop performance actually got worse using your bsxfun method for DD. – brown.2179 Mar 02 '14 at 17:14
but vectorized code can be used to have even better performance if GPUs are around and for loops aren't...right? – Charlie Parker Mar 24 '16 at 05:08
@CharlieParker I guess, you'd have to benchmark and see.. My point was that for-loops in MATLAB had gotten a bad reputation in the past, but they're not as bad nowadays. – Amro Mar 24 '16 at 05:13
@amro it also surprises me that there is such a big difference (even if its only on one machine). I'm surprised mathworks wouldn't have changed that function already to ur implementation. – Charlie Parker Mar 24 '16 at 15:10
@CharlieParker please note that comparison was done in 2011, worth repeating it with the latest version... Also the thing is with vectorized code is that tends to use more memory than an equivalent for-loop... That's one downside – Amro Mar 24 '16 at 23:22
so I ran test on my laptop (quad core i7 running Win10 and R2015b), this time I got `8.5183` vs. `1.1109` seconds (still in favor of the vectorized solution). – Amro Mar 25 '16 at 02:50
Another update, in 2016b, Win 10, i7-3770K, I get `6.686229` and `0.950623` still favoring the vectorized solution. – jodag Dec 13 '17 at 03:49

Jonas · Accepted Answer · 2011-10-08T14:34:55.470

2

You could fully vectorize the calculation by repeating the rows of fea_test 6800 times, and of fea_train 4754 times, like this:

rA = size(fea_test,1);
rB = size(fea_train,1);

[I,J]=ndgrid(1:rA,1:rB);

d = zeros(rA,rB);

d(:) = sqrt(sum(fea_test(J(:),:)-fea_train(I(:),:)).^2,2));

However, this would lead to intermediary arrays of size 6800x4754x1024 (*8 bytes for doubles), which will take up ~250GB of RAM. Thus, the full vectorization won't work.

You can, however, reduce the time of the distance calculation by preallocation, and by not calculating the square root before it's necessary:

rA = size(fea_test,1);
rB = size(fea_train,1);
d = zeros(rA,rB);

for i = 1:rA
    test_data=fea_test(i,:);
    d(i,:)=sum( (test_data(ones(nB,1),:) -  fea_train).^2, 2))';
end

d = sqrt(d);

edited Oct 08 '11 at 14:34

answered Oct 08 '11 at 13:02

Jonas

74,690
10
137
177

3

`repmat` is never good for performance. You're better off inlining in this case: replace `repmat(test_data,nB,1)` with `test_data(ones(1,n8), :)`. – Nzbuu Oct 08 '11 at 13:31
@Jones thank you for this reply, but this does not reduce the number of for loops i have in my code. I already have a for loop within which this loop is found and the execution time is still the same. – bhavs Oct 08 '11 at 15:01
2

@BhavyaPH: The solution that reduces the number of `for` loops will most likely not run on your computer due to not enough RAM. The other solution should reduce the execution time. Have you used the profiler to compare? – Jonas Oct 08 '11 at 15:49
1

Removing for loops isn't always fastest due to the JIT. See http://stackoverflow.com/questions/7569368/matlab-for-loop-vs-vectorization, and many other questions on SO. – Nzbuu Oct 08 '11 at 17:32

Alexey · Answer 3 · 2012-06-29T20:23:12.817

Try this vectorized version, it should be pretty efficient. Edit: just noticed that my answer is similar to @Amro's.

function K = calculateEuclideanDist(P,Q)
% Vectorized method to compute pairwise Euclidean distance
% Returns K(i,j) = sqrt((P(i,:) - Q(j,:))'*(P(i,:) - Q(j,:)))

[nP, d] = size(P);
[nQ, d] = size(Q);

pmag = sum(P .* P, 2);
qmag = sum(Q .* Q, 2);

K = sqrt(ones(nP,1)*qmag' + pmag*ones(1,nQ) - 2*P*Q');

end

pdist2 equivalent in MATLAB version 7

3 Answers3

Linked