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,:));

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


ans =
  -2 -1  0  1
  -3 -2 -1  0
  -4 -3 -2 -1
  -5 -4 -3 -2
  • $\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$ Feb 26, 2012 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, 2012 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$ Feb 27, 2012 at 23:36
  • $\begingroup$ How about octave.sourceforge.net/statistics/function/pdist.html $\endgroup$
    – Nemo
    Mar 7, 2015 at 10:37

2 Answers 2


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);

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

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

assert(norm(d_vect-d_for) < 1e-10*norm(d_for)) 
  • $\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$ Apr 9, 2012 at 21:56

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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