Alec's Web Log

Often I like to get list of distances between all pairs of points in two sets. That is, if I have a set of points V and another set of points U, I’d like to build D, a #V by #U matrix of distance, where D(i,j) gives me the distance from point V(i) to point U(i).

For the sake of timings, construct some V and U values using:

V = rand(5000,3);
U = rand(5000,3);

Two for loops

You can achieve the all pairs distance computation between V and U in the most naive way using two for loops (like you would probably do it in C++):

  % Super-slow for loops method
  D = zeros(size(V,1),size(U,1));
  for i = 1:size(V,1)
    for j = 1:size(U,1)
      D(i,j) = sqrt(sum((V(i,:)-U(j,:)).^2,2));
    end
  end

Two repmats

You can also do this using two repmats in a one-liner, but this gets slow if you hit a memory wall:

  % Slow, especially if hitting memory paging
  D = ...
    squeeze(sqrt(sum((repmat(V,[1 1 size(U,1)])- ...
    permute(repmat(U,[1 1 size(V,1)]),[3 2 1])).^2,2)));

Single for loop, single repmat

You can try to alleviate this by using a single for loop and inside the loop a single repmat. This is faster, but will probably hit a memory wall later on, too:

  % slow single for loop, single repmat method, handles memory paging better
  D = zeros(size(V,1),size(U,1));
  for i = 1:size(V,1)
    D(i,:) = ...
      sqrt(sum((repmat(V(i,:),[size(U,1) 1]) - U).^2,2))';
  end

Single for loop, single bsxfun

The matlab function bsxfun will do the repeat virtually, so we can replace the minus operation inside the loop, giving us a pretty good speed up:

  % Fast single for loop, single bsxfun method, but matlab can handle both for
  % loops in bsxfun
  D = zeros(size(V,1),size(U,1));
  for i = 1:size(V,1)
    D(i,:) = ...
      sqrt(sum(bsxfun(@minus,V(i,:),U).^2,2));
  end

Single bsxfun

But there’s no reason for the outer for loop either, bsxfun is powerful enough to virtually repeat both list of points, resulting in a very fast, very concise, one-liner:

  % Fastest
  D = squeeze(sqrt(sum(bsxfun(@minus,V,permute(U,[3 2 1])).^2,2)));

Note: Contrary to what I thought would happen, the single bsxfun does hit a memory wall for very large input, so the single for loop, single bsxfun/repmat succeed here (see second column). It would be nice to have an automatic way of determining which code to run based on the machine’s memory.