# 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[], pts[]], {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. Mar 22, 2013 at 19:42
• Mar 22, 2013 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? Mar 23, 2013 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. Mar 23, 2013 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, 2013 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
];//AbsoluteTiming

(*  {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):

Sqrt[HierarchicalClusteringDistanceMatrix[pts, DistanceFunction -> EuclideanDistance]];// 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. In recent versions of Mathematica it is optimized for several possible DistanceFunction settings, including ManhattanDistance.

• 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 Mar 22, 2013 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. Mar 22, 2013 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, 2013 at 20:38
• Pre-complementing saves some time: With[{tr = -Transpose@pts}, Total@Abs[# + tr]& /@ pts]. Also, for "small" problems, say < 1700 or so points (on my system -- YMMV), the built-in DistanceMatrix function in the HierarchicalClustering package seems to be faster. I'm puzzled why it hasn't been mentioned. Are people not aware of it? Nov 10, 2013 at 7:42
• Thanks for the comment, @Ray. I did not see any particular speedup from pre-complementing on my system, but DistanceMatrix with the default DistanceFunction is indeed faster in the range you specified, for Euclidean distances. I wasn't aware of that. But also, this is only true when we can avoid explicitly specifying DistanceFunction - please see my edit. Nov 10, 2013 at 9:50

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 listability 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, let's 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} *)


It seems that the simple procedural implementation is a bit less than twice as fast, and 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 ...) Mar 22, 2013 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. Mar 22, 2013 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. Mar 24, 2013 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/… Mar 24, 2013 at 16:58

The DistanceMatrix function, newly introduced in version 10.3, is very fast for Euclidean distances.

Here's a speed comparison with Leonid's fast solution.

pts = RandomReal[1, {5000, 2}];


Euclidean

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

(* {0.952627, Null} *)

dm2 = DistanceMatrix[pts, pts]; // AbsoluteTiming
(* {0.212496, Null} *)

dm3 = HierarchicalClusteringDistanceMatrix[pts, DistanceFunction -> EuclideanDistance]; // AbsoluteTiming
(* {0.582991, Null} *)

dm1 == dm2 == dm3
(* True *)


Note that HierarchicalClusteringDistanceMatrix is a built-in function, not one provided by the package. Most performance-critical functions of this "package" are in reality highly optimized built-ins. Also note that the default distance for this function is not EuclideanDistance, but that square of that. So we needed to specify EuclideanDistance explicitly.

Manhattan

Let's test if these functions are special-cased for the Manhattan distance.

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

dm2 =
DistanceMatrix[pts, pts,
DistanceFunction -> ManhattanDistance]; // AbsoluteTiming
(* {0.211771, Null} *)

dm3 =
HierarchicalClusteringDistanceMatrix[pts,
DistanceFunction -> ManhattanDistance]; // AbsoluteTiming
(* {0.621123, Null} *)

dm1 == dm2 == dm3
(* True *)


A final note

For accurate benchmarking is it very important to use AbsoluteTiming and not Timing here. In recent versions of Mathematica all of these operations are internally parallelized and Timing would measure the total CPU time spent by each core, added up, instead of the wall time.

Just for fun, here's a C++ version using LTemplate. This is specialized for 2D points!

<< LTemplate


I'm on a Mac where the system compiler doesn't support OpenMP. I'll use gcc from MacPorts to be able to use OpenMP.

$CCompiler = {"Compiler" -> CCompilerDriverGenericCCompilerGenericCCompiler, "CompilerInstallation" -> "/opt/local/bin", "CompilerName" -> "g++-mp-5", "SystemCompileOptions" -> "-O3 -m64 -fPIC -framework Foundation -framework mathlink"}; SetDirectory[$TemporaryDirectory];
code = "
#include <cmath>

inline double sqr(double x) { return x*x; }

struct DistMatrix {
mma::RealTensorRef distMat(mma::RealMatrixRef a, mma::RealMatrixRef b) {
mma::RealMatrixRef mat = mma::makeMatrix<double>(a.rows(), b.rows());
#pragma omp parallel for
for (mint i=0; i < a.rows(); ++i)
for (mint j=0; j < b.rows(); ++j)
mat(i, j) = std::hypot(a(i,0) - b(j,0), a(i,1) - b(j,1));
return mat;
}
};
";
Export["DistMatrix.h", code, "String"];

tem = LClass[
"DistMatrix",
{LFun["distMat", {{Real, 2, "Constant"}, {Real, 2, "Constant"}}, {Real, 2}]}
];

CompileTemplate[tem,
"CompileOptions" -> {"-std=c++14", "-fopenmp"}]

obj = Make["DistMatrix"];

dm4 = obj@"distMat"[pts, pts]; // AbsoluteTiming
(* {0.062397, Null} *)

dm1 == dm4
(* True *)

• Hello, Szabolcs. I always wanted to learn Wolfram Library link technique. However, for the beginer(like me), I discovered that it is not easy(hard) for me to learn that via the Wolfram LibraryLink User Guide. Could you give some suggestions to beginer to let them understand that technique more easily? Thanks a bunch:)
– xyz
May 29, 2016 at 2:33

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


This is a good case that can sped up with FunctionCompile. My starting point is Leonid's code from his accepted answer:

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

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


Gives: {1.7097, Null}

Write the code as a function, with typed arguments:

f =
Function[{Typed[pts, TypeSpecifier["PackedArray"]["Real64", 2]]},
With[{tr = Transpose[pts]}, Map[Sqrt[Total[(# - tr)^2]] &, pts]]];

AbsoluteTiming[d2 = f[pts];]


Gives a similar timing: {1.70066, Null}

Use FunctionCompile to compile it:

cf = FunctionCompile[f];

AbsoluteTiming[d3 = cf[pts];]


Gives a 2-3x speedup: {0.51376, Null}

Check that all the results are the same:

d1 == d2 == d3
`

Gives: True