How to use Compile to optimize the performance of a function calculating the distance between two points?

Question

I want to create a function that finds the real space distance between two three-dimensional points. Each point is a list of three real numbers.

I can create functions dist1, dist2, and distCompiled1 and test them on example points myPoint1 and myPoint2:

dist1[pt1_List, pt2_List] := Sqrt[(pt1[[1]] - pt2[[1]])^2 + (pt1[[2]] - pt2[[2]])^2
   + (pt1[[3]] - pt2[[3]])^2];
dist2[pt1_List, pt2_List] := Sqrt[(pt1 - pt2).(pt1 - pt2)];
distCompiled1 := Compile[{{pt1, _Real, 1}, {pt2, _Real, 1}},
   Sqrt[(pt1 - pt2).(pt1 - pt2)]];

myPoint1 = {1.0, 2.0, 3.0};
myPoint2 = {4.0, 5.0, 6.0};

dist1[myPoint1, myPoint2]
dist2[myPoint1, myPoint2]
distCompiled1[myPoint1, myPoint2]

5.19615

5.19615

5.19615

So, dist1, dist2, and distCompiled1 all give the same result.

But, I will be computing distances between many millions of points, so I want to optimize my function for speed. If I calculate the distance between 1 million pairs of points, for example, I get the following results:

dist1[pt1_List, pt2_List] := Sqrt[(pt1[[1]] - pt2[[1]])^2 + (pt1[[2]] - pt2[[2]])^2
   + (pt1[[3]] - pt2[[3]])^2];
dist2[pt1_List, pt2_List] := Sqrt[(pt1 - pt2).(pt1 - pt2)];
distCompiled1 := Compile[{{pt1, _Real, 1}, {pt2, _Real, 1}},
   Sqrt[(pt1 - pt2).(pt1 - pt2)]];
maxVal = 8.009469032;
numPts = 10^6;
SeedRandom[1234];
coords1 = Table[RandomReal[{0, maxVal}, 3], {numPts}];
SeedRandom[1235];
coords2 = Table[RandomReal[{0, maxVal}, 3], {numPts}];

(* Time the computations, and verify that the functions give the same results. *)
AbsoluteTiming[
dist1Result = Table[dist1[coords1[[i]], coords2[[i]]], {i, 1, numPts}];
]
AbsoluteTiming[
dist2Result = Table[dist2[coords1[[i]], coords2[[i]]], {i, 1, numPts}];
]
AbsoluteTiming[
distCompiled1Result = Table[distCompiled1[coords1[[i]], coords2[[i]]], {i, 1, numPts}];
]
(dist1Result == dist2Result) && (dist2Result == distCompiled1Result)

{6.725385, Null}

{5.900337, Null}

{30.191727, Null}

True

• Why is distCompiled1, the function that uses Compile, so much slower -- not faster -- than the others that don't use Compile?

• How can I rewrite distCompiled1 so that it ideally runs faster than dist1 and dist2?

Answer 1

distCompiled = 
  Compile[{{pt1, _Real, 1}, {pt2, _Real, 1}}, 
   Sqrt[(pt1[[1]] - pt2[[1]])^2 + 
   (pt1[[2]] - pt2[[2]])^2 + 
   (pt1[[3]] - pt2[[3]])^2],
   CompilationTarget -> "C", RuntimeAttributes -> {Listable},
     RuntimeOptions -> "Speed"];

dist1Result = Sqrt@Total[(coords1 - coords2)^2, {2}]; // AbsoluteTiming
dist2Result = distCompiled[coords1, coords2]; // AbsoluteTiming
dist1Result == dist2Result

{0.0334167, Null}

{0.0247759, Null}

True

Answer 2

Karl's answer is perfect as it is. If you know for sure that your input vectors have at least length 3, then it is safe to deactivate the bound checks that Mathematica writes into the C code. This can be done by the undocumented feature Compile`GetElement instead of Part. Note that you need RuntimeOptions -> "Speed" for this to work correctly.

distCompiled2 = Compile[{{pt1, _Real, 1}, {pt2, _Real, 1}},
   Block[{dx, dy, dz},
    dx = Compile`GetElement[pt1, 1] - Compile`GetElement[pt2, 1];
    dy = Compile`GetElement[pt1, 2] - Compile`GetElement[pt2, 2];
    dz = Compile`GetElement[pt1, 3] - Compile`GetElement[pt2, 3];
    Sqrt[dx dx + dy dy + dz dz]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

On my machine, this makes the function up to 20% faster.

Answer 3

[Answer to a different question...]

The compiled methods are fine if, for each point in one set, you only want to check distances to a small number of points in the other set (the examples have "small" set to 1). If instead you want to know which points in the second are within a specified distance to the first, Nearest is the function to use.

maxVal = 8.009469032;
numPts = 10^6;
SeedRandom[1234];
coords1 = Table[RandomReal[{0, maxVal}, 3], {numPts}];
SeedRandom[1235];
coords2 = Table[RandomReal[{0, maxVal}, 3], {numPts}];

Find all coords2 neighbors to coords1 points that are within distance of 0.1.

AbsoluteTiming[
 nbrs = Nearest[coords1, coords2, {Infinity, .1}];
 ]

(* Out[81]= {8.00687, Null} *)

Check the first few neighbor sets.

In[82]:= nbrs[[1 ;; 4]]

(* Out[82]= {{{4.22894, 7.34831, 3.22588}, {4.2806, 7.37094, 
   3.30835}, {4.22592, 7.2422, 3.29461}, {4.25438, 7.32368, 
   3.34095}}, {{0.93075, 5.86791, 4.62192}, {0.877728, 5.88219, 
   4.70008}, {0.928627, 5.96515, 4.66545}, {0.984484, 5.91222, 
   4.69613}, {0.901294, 5.97538, 4.64348}, {0.860855, 5.95889, 
   4.67252}, {0.944867, 5.83116, 4.61225}, {0.924464, 5.80529, 
   4.69967}}, {{4.34926, 3.20961, 5.60807}, {4.36852, 3.25073, 
   5.60718}, {4.34147, 3.28848, 5.63214}, {4.35976, 3.29842, 
   5.56703}, {4.37305, 3.30205, 5.58017}, {4.35996, 3.29815, 
   5.53587}, {4.34908, 3.24295, 5.67379}, {4.30285, 3.31839, 
   5.53733}, {4.30903, 3.25791, 5.49811}, {4.2832, 3.27267, 
   5.66872}, {4.37751, 3.17006, 5.59264}, {4.40635, 3.28206, 
   5.61345}, {4.3599, 3.16009, 5.608}, {4.39588, 3.18371, 
   5.57433}}, {{1.48945, 3.21571, 4.98408}, {1.51413, 3.21257, 
   5.02518}, {1.44867, 3.23852, 5.0191}, {1.50369, 3.21718, 
   5.04696}, {1.47354, 3.15735, 4.98367}, {1.48652, 3.21397, 
   5.06329}, {1.50156, 3.30154, 5.01192}, {1.43768, 3.2166, 
   4.91252}, {1.48613, 3.26883, 4.90126}}} *)

Answer 4

What follows is a more detailed demonstration of Karl's answer on this page. This demonstrates the performance improvement (amazing, in my opinion) offered by Listable functions in Mathematica:

(* Initialization *)
maxVal = 8.009469032;
numPts = 10^6;
SeedRandom[1234]; coords1 = Table[RandomReal[{0, maxVal}, 3], {numPts}];
SeedRandom[1235]; coords2 = Table[RandomReal[{0, maxVal}, 3], {numPts}];

(* Functions *)
distAndrew[pt1_List, pt2_List] := Sqrt[(pt1[[1]] - pt2[[1]])^2
    + (pt1[[2]] - pt2[[2]])^2
    + (pt1[[3]] - pt2[[3]])^2];

distKarl[listOfPts1_List, listOfPts2_List] := 
  Sqrt@Total[(listOfPts1 - listOfPts2)^2, {2}];

distKarlCompiled = Compile[{{pt1, _Real, 1}, {pt2, _Real, 1}},
   Sqrt[(pt1[[1]] - pt2[[1]])^2
     + (pt1[[2]] - pt2[[2]])^2
     + (pt1[[3]] - pt2[[3]])^2],
   CompilationTarget -> "C", RuntimeAttributes -> {Listable},
   RuntimeOptions -> "Speed"];

(* Performance testing *)
distAndrewResult = Table[distAndrew[coords1[[i]], coords2[[i]]], {i, 1, 
     numPts}]; // AbsoluteTiming
distAndrewWithMapResult = 
   Map[Apply[distAndrew, #] &, Transpose[{coords1, coords2}]]; // AbsoluteTiming
distKarlResult = distKarl[coords1, coords2]; // AbsoluteTiming
distKarlCompiledResult = distKarlCompiled[coords1, coords2]; // AbsoluteTiming

(distAndrewResult == distAndrewWithMapResult) && 
(distAndrewWithMapResult == distKarlResult) &&
(distKarlResult == distKarlCompiledResult)

{7.979455, Null}

{8.709511, Null}

{0.116064, Null}

{0.062322, Null}

Answer 5

If you don't want to bother with Compile, you could use either of the built-in functions EuclideanDistance or Norm:

AbsoluteTiming[
 edResult = Table[EuclideanDistance[coords1[[i]], coords2[[i]]], {i, 1, numPts}];]
AbsoluteTiming[n2Result = Norm /@ (coords1 - coords2);]

Even the second one is slower than distKarlCompiled by a factor of 5 though.