12
$\begingroup$

Suppose I have two matrices Nx2, Mx2 representing N, M 2d vectors respectively. Is there a simple and good way to calculate distances between each vector pair (n, m)?

The easy but inefficient way is of course:

d = zeros(N, M);
for i = 1:N,
  for j = 1:M,
    d(i,j) = norm(n(i,:) - m(j,:));
  endfor;
endfor;

The closest answer I've found is bsxfun, used like so:

bsxfun(inline("x-y"),[1,2,3,4],[3;4;5;6])

ans =
  -2 -1  0  1
  -3 -2 -1  0
  -4 -3 -2 -1
  -5 -4 -3 -2
$\endgroup$
  • $\begingroup$ I took a look at this and I couldn't do much better than vectorizing the computation. I think this computation is a pretty good candidate for writing an external C/Fortran function. $\endgroup$ – Aron Ahmadia Feb 26 '12 at 11:05
  • 1
    $\begingroup$ I bet you could make a 2xNxM matrix which you populate with an outer product, then square each of the entries and sum along the zeroth axis and square root. In Python this would look like: distance_matrix = (n[:,:,nexaxis] * m[:,newaxis,:]); distance_matrix = distance_matrix**2; distance_matrix = sqrt(distance_matrix.sum(axis=1)); If you only need to know the closest n-vectors there are much better ways to do this! $\endgroup$ – meawoppl Feb 27 '12 at 20:46
  • 3
    $\begingroup$ @meawoppl (New to Octave) I found out how to use the linear-algebra package in Octave, which provides cartprod, so now I can write: (1) x = cartprod(n(:,1), m(:,1)); (2) y = cartprod(n(:,2), m(:,2)); (3) d = sqrt((x(:,1)-x(:,2)).^2+(y(:,1)-y(:,2)).^2) ..which runs much faster! $\endgroup$ – Kelley van Evert Feb 27 '12 at 23:36
  • $\begingroup$ How about octave.sourceforge.net/statistics/function/pdist.html $\endgroup$ – Nemo Mar 7 '15 at 10:37
6
$\begingroup$

Vectorizing is straightforward in these situations using a strategy like this:

eN = ones(N,1);
eM = ones(M,1);
d  = sqrt(eM*n.^2' - 2*m*n' + m.^2*eN');

Here's an example that vectorizes the for loop with a 15x speedup for M=1000 and N=2000.

n = rand(N,2);
m = rand(M,2);
eN = ones(N,2);
eM = ones(2,M);

tic;
d_vect  = sqrt(eN*m.^2' - 2*n*m' + n.^2*eM);
vect_time = toc;

tic;
for i=1:N
  for j=1:M
     d_for(i,j) = norm(n(i,:)-m(j,:));
  end
end
for_time = toc; 

assert(norm(d_vect-d_for) < 1e-10*norm(d_for)) 
$\endgroup$
  • $\begingroup$ David, nice to see you on scicomp! I shamelessly edited your code fragment and expanded it some, please revert if my edits went the wrong direction from clarifying what you intended. $\endgroup$ – Aron Ahmadia Apr 9 '12 at 21:56
2
$\begingroup$

From Octave 3.4.3 and later the operator - does automatic broadcasting (uses bsxfun internally). So you can proceed in this way.

Dx = N(:,1) - M(:,1)';
Dy = N(:,2) - M(:,2)';
D = sqrt (Dx.^2 + Dy.^2);

You can do the same using a 3d matrix but I guess like this is more clear. D is a NxM matrix of distances, every vector in N against every vector in M.

Hope this helps

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.