# Fastest way to calculate matrix of pairwise distances

It is a very common problem that given a distance function $d(p_1,p_2)$ and a set of points pts, we need to construct a matrix mat so that mat[[i,j]] == d[ pts[[i]], pts[[j]] ].

What is the most efficient way to do this in Mathematica?

Let's assume that the points are in $\mathbb{R}^n$ for simplicity, and because that's the case I'm dealing with now, but theoretically the points could be any type of object, e.g. strings with $d$ being an edit distance.

For the specific problem I have right now I need to calculate the EuclideanDistance and ManhattanDistance of 2D points.

The simplest way to do this is

pts = RandomVariate[NormalDistribution[], {1000, 2}];

mat = Outer[ManhattanDistance, pts, pts, 1]; // AbsoluteTiming

(* ==> {0.595327, Null} *)


This obviously calculates all distances twice, which is wasteful. So I was hoping for an easy $2\times$ speedup, but it isn't as easy as one would hope. Doing the same operation the same number of times in a Do loop takes considerably longer (probably because of indexing):

Do[ManhattanDistance[pts[[10]], pts[[20]]], {Length[pts]^2}]; // AbsoluteTiming

(* ==> {1.902417, Null} *)


So what programming pattern do you typically use when calculating such a distance matrix and which one would you recommend for this specific problem?

-
There are of course lots of different ways to do the calculation. Unfortunately it is not at all obvious which is likely to be the fastest. –  Szabolcs Mar 22 '13 at 19:42
–  andre Mar 22 '13 at 19:55
Good question, just what I was looking for at the moment! I also wonder, how much speed can be gained if only the distance functions is reimplemented & compiled? –  István Zachar Mar 23 '13 at 6:16
I just want to point out that most of the answers (using the vectorized speedup) won't work for the more general case (e.g. Hamming distance of vectors, edit distance of strings, etc.). If one does not want to desing a specific, highly optimized method for each such case, I found that using Outer is generally the fastest way. –  István Zachar Mar 23 '13 at 9:11
Why not generate all the entries in, say, the lower triangle (with something like Table[(* stuff *), {j, n}, {k, j}] // Flatten) and feed this to SymmetrizedArray[]? (I do not have version 9, so I can't test this.) –  Ｊ. Ｍ. Mar 23 '13 at 16:06

Using Outer is here one of the worst methods, and not just because it computes the distance twice, but because you can't leverage vectorization in this approach. This is actually a common issue and an important point to stress: Outer works pairwise and is unable to utilize the possible vectorized nature of the operation it is performing on an element-by-element basis.

distances=
With[{tr=Transpose[pts]},
Function[point,Sqrt[Total[(point-tr)^2]]]/@pts
];//Timing

(*  {0.046875,Null} *)


which is an order of magnitude faster. You can Compile it with a C target which may improve the performance further. Also, essentially the same approach I used in this recent answer, with good performance.

For Manhattan distance, use

distances =
With[{tr = Transpose[pts]},
Function[point, Total@Abs[(point - tr)]] /@ pts];


EDIT

As noted by Ray Koopman in comments, the function DistanceMatrix from the package HierarchicalClustering may be faster for Euclidean distance, for small and medium data size (up to a couple of thousands):

<< HierarchicalClustering
Sqrt[DistanceMatrix[pts]];// AbsoluteTiming

(* {0.019351, Null} *)


Note, however, that this is only true for the particular case of Euclidean distance, or perhaps other distances which don't require to set the DistanceFunction option explicitly on the top-level. In other cases (for example, for Manhattan distance), it will be quite slow, because when DistanceFunction is set explicitly, one can not leverage vectorization any more, once again.

-
I had not taken apart the previous code. This is impressive. Even when you consider that it is performing the calculation twice as often as required. +1 –  rcollyer Mar 22 '13 at 20:29
@rcollyer Thanks. What "previous code" do you mean though? I actually did not want to get the absolute fastest (otherwise would go with Compile), but wanted to stress this issue with Outer and explain why it is slow here. –  Leonid Shifrin Mar 22 '13 at 20:31
Seeing as I goofed up the dimensions, I pulled out an old MATLAB code of mine for a vectorized solution to aim for speed (I had forgotten what I had done) and it was more or less identical to this :) –  rm -rf Mar 22 '13 at 20:38
@rm-rf Well, that's a cool compliment, since it comes from Matlab guru :-). –  Leonid Shifrin Mar 22 '13 at 20:42
@LeonidShifrin I really think Mathematica is too hard to learn right and that hinders adoption. –  Rolf Mertig Mar 22 '13 at 22:28

Here is a little procedural implementation using Bag compiled to C:

distmatrix = Compile[{{pts, _Real, 2}},
Block[{x, y, list = InternalBag[Most[{0.}]]},
For[x = 1, x <= Length[pts], x++,
For[y = x + 1, y <= Length[pts], y++,
InternalStuffBag[list,
Abs[CompileGetElement[pts, x, 1] -
CompileGetElement[pts, y, 1]] +
Abs[CompileGetElement[pts, x, 2] -
CompileGetElement[pts, y, 2]]];
];
];
InternalBagPart[list, All]
], CompilationTarget -> "C", RuntimeOptions -> "Speed"];

distmatrix[pts]; // AbsoluteTiming

(*

{0.009000, Null}

*)


edit: an even better performing solution is based on Mr. Wizard's vectorization approach and relying on the lisability and parallelizability of compiled functions and as a nice touch it doesn't rely on undocumented functions.

distmatrix2 =
Compile[{{point, _Real, 1}, {tr, _Real, 2}},
Total@Abs[point - tr], CompilationTarget -> "C",
RuntimeOptions -> "Speed", RuntimeAttributes -> {Listable},
Parallelization -> True];


For comparison against Leonid's method, lets use more points.

pts = RandomVariate[NormalDistribution[], {10000, 2}];

distmatrix[pts]; // AbsoluteTiming

distances =
With[{tr = Transpose[pts]},
Function[point, Total@Abs[(point - tr)]] /@ pts]; // AbsoluteTiming

distmatrix2[pts, Transpose[pts]]; // AbsoluteTiming

(* {0.865050, Null}, {1.562089, Null},  {0.319018, Null} *)


Seems that the simple procedural implementation is a bit less then twice as fast, not really worth the extra work/complexity. The listable/parallelized compiled solution is simpler and about 5x faster.

-
Can you run a comparison against Leonid's code, as my timings differ from yours? (Yours is still faster ...) –  rcollyer Mar 22 '13 at 21:12
+1 still. The reason why it is hard to make the code much more efficient particularly for large lists is that the larger the list, the closer is the performance of my code to the low-level C performance. This is because what matters always is to make fast inner loops - and in my case the inner loop is replaced by a vectorized operation. The longer the list, the more the bottleneck is in the vectorized operation, and you can't get much faster here since those are highly optimized. –  Leonid Shifrin Mar 22 '13 at 21:56
Where can I find information about Internal*, Compile* etc.? It doesn't appear to be in any standard documentation. I find it unfair that WRI hides so much these useful internal features. They will always be able to produce more efficient code if they have such tools available. –  Federico Mar 24 '13 at 16:32
@Federico, unfortunately they are undocumented. I figured out how to use them from this site, see: mathematica.stackexchange.com/questions/1934/… & mathematica.stackexchange.com/questions/8650/… –  s0rce Mar 24 '13 at 16:58

This is another vectorized approach which is an order of magnitude faster than using Outer, but about 1.5 times slower than Leonid's answer:

dist = With[{c = ConstantArray[Dot[#, #] & /@ pts, {Length@pts}]},
c + Transpose@c - 2 pts . Transpose@pts // Sqrt];
`
-