5
$\begingroup$

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?

$\endgroup$
3
  • 2
    $\begingroup$ Meta question: does you end use actually require computing all point-point-pair distances? Or might you only really require distances between "close" pairs? $\endgroup$ Commented Nov 2, 2023 at 19:32
  • $\begingroup$ @DanielLichtblau You're right. I only need distances between "close" pairs. But I'm not experienced enough to know the best way to do that; I need to do some research. A simple way might be to divide space into segments, and examine the current and neighboring segments first. Do you have any hints? $\endgroup$
    – Andrew
    Commented Nov 2, 2023 at 20:07
  • 2
    $\begingroup$ Nearest can deliver all points within a specified radius of a given one. I’ll post an example when I’m back at my desk. $\endgroup$ Commented Nov 3, 2023 at 0:28

5 Answers 5

8
$\begingroup$
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

$\endgroup$
2
  • 1
    $\begingroup$ Wow, this is an amazing demonstration. It shows the capabilities, in terms of performance, of listability in Mathematica. $\endgroup$
    – Andrew
    Commented Nov 2, 2023 at 20:27
  • $\begingroup$ @Karl Did you really use numPts = 10^6? Even with the performance improvement with listability, you must have a fast computer. $\endgroup$
    – Andrew
    Commented Nov 2, 2023 at 21:04
6
$\begingroup$

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 insteadl 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.

$\endgroup$
4
$\begingroup$

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}

$\endgroup$
4
$\begingroup$

[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}}} *)
$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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