Octave: calculate distance between two matrices of vectors

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

• 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. – Aron Ahmadia Feb 26 '12 at 11:05
• 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! – meawoppl Feb 27 '12 at 20:46
• @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! – Kelley van Evert Feb 27 '12 at 23:36
• – Nemo Mar 7 '15 at 10:37

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

• 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. – Aron Ahmadia Apr 9 '12 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