Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
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
related :… – 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.) – J. M. Mar 23 '13 at 16:06
up vote 32 down vote accepted

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.

Here is the code I will adopt from this answer:


(*  {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];


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[HierarchicalClustering`DistanceMatrix[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.

share|improve this answer
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 :) – R. M. 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
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? – Ray Koopman Nov 10 '13 at 7:42

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

distmatrix = Compile[{{pts, _Real, 2}},
   Block[{x, y, list = Internal`Bag[Most[{0.}]]},
    For[x = 1, x <= Length[pts], x++,
     For[y = x + 1, y <= Length[pts], y++,
         Abs[Compile`GetElement[pts, x, 1] - 
            Compile`GetElement[pts, y, 1]] + 
          Abs[Compile`GetElement[pts, x, 2] - 
            Compile`GetElement[pts, y, 2]]];
    Internal`BagPart[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.

share|improve this answer
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:… &… – s0rce Mar 24 '13 at 16:58

The DistanceMatrix function, newly introduced in version 10.3, is very fast for Euclidean distances. It is still marked as [[Experimental]] in the version 10.3 documentation.

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

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


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 = HierarchicalClustering`DistanceMatrix[pts, DistanceFunction -> EuclideanDistance]; // AbsoluteTiming
(* {0.582991, Null} *)

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

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


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 = 
    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" -> CCompilerDriver`GenericCCompiler`GenericCCompiler, 
   "CompilerInstallation" -> "/opt/local/bin", 
   "CompilerName" -> "g++-mp-5", 
   "SystemCompileOptions" -> 
    "-O3 -m64 -fPIC -framework Foundation -framework mathlink"};

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[
    {LFun["distMat", {{Real, 2, "Constant"}, {Real, 2, "Constant"}}, {Real, 2}]}

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


obj = Make["DistMatrix"];

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

dm1 == dm4
(* True *)
share|improve this answer

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];
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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