12
$\begingroup$

I want to find the Euclidean distance between one point (x1) and a list of points (y1), which contains a lot of coordinates

x1 = killer[[2]]

{6.05102, 5.87667}

y1 = victim[[2 ;;]]

{{1.40687, 4.92494}, {0.419206, 1.70406}, {6.29657,0.577941}, {4.12022, 4.94952},
{2.04784, 5.94545}, {1.29192,1.43152}, {3.26737, 1.90134}, {4.27274, 0.528028},
{2.79659,1.37788}, {5.43955, 1.81355}}

Is it possible for me to find the EuclideanDistance between x1 and y1, where it will show all results between x1 and each elements in y1.

$\endgroup$
  • 4
    $\begingroup$ Try EuclideanDistance[x1,#]&/@y1! $\endgroup$ – PlatoManiac Feb 2 '13 at 17:42
  • 2
    $\begingroup$ You need to read up on Map. $\endgroup$ – Szabolcs Feb 2 '13 at 17:51
  • 1
    $\begingroup$ Related: Calculating the distance between the coordinates in a list $\endgroup$ – Yves Klett Feb 2 '13 at 18:12
  • 5
    $\begingroup$ @PlatoManiac If efficiency is important, one can gain an order of magnitude speedup by using Sqrt[Total[(x1 - #)^2]] & /@ y1 instead, due to auto-compilation (EuclideanDistance is not compilable). Slightly faster still can be a vectorized solution like Sqrt[Total[(Transpose[y1] - x1)^2]]. $\endgroup$ – Leonid Shifrin Feb 2 '13 at 18:27
  • 2
    $\begingroup$ Also related $\endgroup$ – Leonid Shifrin Feb 2 '13 at 18:27
10
$\begingroup$

Here is what you want to do. Let the first point's coordinate be stored as a list of list as follows:

x1 = {{6.05102, 5.87667}}    

and the second set of coordinates

y1 = {{1.40687, 4.92494}, {0.419206, 1.70406}, {6.29657,0.577941}, {4.12022, 4.94952},
      {2.04784, 5.94545}, {1.29192,1.43152}, {3.26737, 1.90134}, {4.27274, 0.528028},
      {2.79659,1.37788}, {5.43955, 1.81355}} 

Now to compute the Euclidean distances between x1 and every element in y1 use Outer, your best friend from now on.

Outer[EuclideanDistance, x1, y1, 1]//Flatten

This then gives you

{4.74067, 7.00914, 5.30442, 2.14187, 4.00377, 6.51217, 4.85304, 5.63651, 5.55252, 4.10887}. 

Hope this helps. In fact, you can loop this through various x1's as follows.

Table[Outer[EuclideanDistance, {killer[[k]]}, y1, 1], {k, 1, n}]

Where n is the Length of the killer list. This is fast and compact code.

$\endgroup$
  • $\begingroup$ Just a suggestion: formatting code inline with text makes copy&pasting less convenient compared to separate code lines. $\endgroup$ – Yves Klett Feb 2 '13 at 20:38
  • $\begingroup$ @Yves, sorry I'm new to StackExchange, I'll try to edit it. $\endgroup$ – RunnyKine Feb 2 '13 at 20:58
  • $\begingroup$ You don't even need the Table. Just put the entire list of x1 into Outer. $\endgroup$ – Jens Feb 2 '13 at 21:01
  • $\begingroup$ No worries whatsoever, welcome to the party! $\endgroup$ – Yves Klett Feb 2 '13 at 21:09
  • $\begingroup$ @Jens, yes you're right. In this case you don't. But suppose you want to do further analysis of a large data set, like calculate the Max or Min of the distances from each of the x1's, then the table definitely comes in handy for Memory conservation. Otherwise Mathematica will store all of those unnecessary data in memory, because Outer will have to finish before you can Map "Min" or "Max" to the entire List. $\endgroup$ – RunnyKine Feb 2 '13 at 21:10
13
$\begingroup$

I would use DistanceMatrix for this. Here is a comparison of DistanceMatrix with the other answers:

SeedRandom[1];
data = RandomReal[1, {10^7, 2}];

r1 = First @ DistanceMatrix[{{1., 1.}}, data]; //AbsoluteTiming
r2 = distance[{1,1}, data]; //AbsoluteTiming (* Mahdi *)
r3 = Outer[EuclideanDistance, {{1., 1.}}, data, 1] //Flatten; //AbsoluteTiming (* RunnyKine *)

r1 == r2 == r3

{0.384095, Null}

{0.88146, Null}

{9.64241, Null}

True

$\endgroup$
9
$\begingroup$

I understand that this question already has an accepted answer, but couldn't resist posting my answer. This answer is mainly useful if you have a very large dataset. I define the following function using Compile function in Mathematica.

distance = Compile[{{n, _Real, 1}, {z, _Real, 2}},
Sqrt[Total[(# - n)^2]] & /@ z, RuntimeOptions -> "Speed", 
Parallelization -> True, CompilationTarget -> "C", 
RuntimeAttributes -> {Listable}
];

Then we map the function over a two dimensional set with 10000000 (=$10^7$) points.

data = RandomReal[{0, 1}, {10000000, 2}];

The result is:

new[{0, 0}, data] // AbsoluteTiming // First
(* 0.609684 s *)

which is 10 times faster compared to the accepted answer:

Outer[EuclideanDistance, {{0, 0}}, data, 1] // AbsoluteTiming // First
(* 5.927118 s *)
$\endgroup$
6
$\begingroup$

Here's a vectorized solution:

SeedRandom[1];
data = RandomReal[1, {10^7, 2}];
pt = {1., 1.};

res2 = Sqrt@Total[(Transpose[data] - pt)^2]; // AbsoluteTiming
(* {0.472308, Null} *)

Compare to Carl's DistanceMatrix, Mahdi's compiled version and jkuczm's compiled version:

res1 = First@DistanceMatrix[{pt}, data]; // AbsoluteTiming
(* {0.399328, Null} *)

res3 = distance[{1, 1}, data]; // AbsoluteTiming
(* {0.958243, Null} *)

res5 = distanceLib[{1,1 }, data]; // AbsoluteTiming
(* {0.065559, Null} *)

res1 == res2 == res3
(* True *)

The lesson is that while built-in functions are usually the best (or at least have the highest potential for running fast), vectorization is usually the next best thing, and often beats an auto-parallelized compiled implementation.

I suspect that this is because vectorization can make use of not only parallelization, but also SIMD instructions.

Finally, vector arithmetic probably relies on the MKL, a highly optimized library that auto-generated C code just can't compete with.


As an experiment, we can also compare the performance of hand-written C++ code. While naive C++ code won't be parallelized, we have the advantage that we know in advance that we are going to be working with 2D points, so we can specialize the code for that case.

Using LTemplate (to save time),

Needs["LTemplate`"]
SetDirectory[$TemporaryDirectory];

tem = LClass["Dist",
   {LFun["dist", {{Real, 1, "Constant"} (* point *), {Real, 2, "Constant"} (* point list *)}, {Real, 1}(* distance list *)]}
   ];

code = "
  #include <cmath>

  using namespace mma;

  struct Dist {
    RealTensorRef dist(RealTensorRef pt, RealMatrixRef ptlist) {
        massert(pt.size() == 2 && ptlist.cols() == 2);
        auto res = makeVector<double>(ptlist.rows());
        for (int i=0; i < res.size(); ++i)
            res[i] = std::hypot(ptlist(i,0) - pt[0], ptlist(i,1) - pt[1]);
        return res;
    }   
  };
  ";
Export["Dist.h", code, "String"];

CompileTemplate[tem]
LoadTemplate[tem]
obj = Make[Dist]

res4 = obj@"dist"[pt, data]; // AbsoluteTiming
(* {0.05077, Null} *)

res1 == res4
(* True *)

That's a 7.5x speedup over the DistanceMatrix function with only a trivial amount of naïve C++ code. I think this illustrates well the advantages of LTemplate. It makes LibraryLink development easy enough that I can write libraries like this in as little as 5 minutes. Using plain LibraryLink, it would take several times longer, which is enough of a deterrent that I would not do this nearly as frequently as LTemplate allows me to.

Update: It's only marginally faster than jkuczm's distanceLib.


Since people using a slow system like the Raspberry Pi may be particularly interested in performance, here's a benchmark with $10^6$ points on a first-gen RPi / Raspbian Jessie / gcc 4.9.2 / Mathematica 11.0.1:

DistanceMatrix: 1.12 s; Vectorization: 1.44 s; Library function: 0.67 s.

The speedup is not as dramatic as it was on OS X / x86_64. I am not sure if gcc 6.3 in Raspbian Stretch would bring any improvements.

$\endgroup$
  • $\begingroup$ +1 Perhaps add link to your LTemplate (to save time searching for a package that saves time...)? $\endgroup$ – bobthechemist Nov 8 '17 at 12:53
  • $\begingroup$ @bobthechemist I thought I had already added the link! Thanks for pointing out that it was missing. $\endgroup$ – Szabolcs Nov 8 '17 at 13:01
  • 2
    $\begingroup$ @bobthechemist Just added on RPi benchmark for you. $\endgroup$ – Szabolcs Nov 8 '17 at 13:09
  • $\begingroup$ Thanks - I suspect there is only so much that can be done with the caps on RAM and a 32 bit OS. Even small improvements are worth the effort, IMO. $\endgroup$ – bobthechemist Nov 8 '17 at 14:12
4
$\begingroup$

Mahdi's compiled function is not really using runtime Listable attribute, since it accepts rank 2 tensor as second argument.

To use runtime listability and Parallelization we need to define function accepting lower rank tensor:

distanceListable = Compile[{{x0, _Real, 1}, {x, _Real, 1}},
    Sqrt[
        Subtract[Compile`GetElement[x, 1], Compile`GetElement[x0, 1]]^2+
        Subtract[Compile`GetElement[x, 2], Compile`GetElement[x0, 2]]^2
    ],
    CompilationTarget -> "C", RuntimeOptions -> "Speed",
    Parallelization -> True, RuntimeAttributes -> {Listable}
];

If we accept rank 2 tensor and manually loop over it, we don't need runtime attributes, and can completely get rid of CompiledFunction overhead by directly using underlying LibraryFunction (Last element of CompiledFunction):

distanceLib = Last@Compile[{{x0, _Real, 1}, {x, _Real, 2}},
    Table[
        Sqrt[
            Subtract[Compile`GetElement[x, i, 1], Compile`GetElement[x0, 1]]^2+
            Subtract[Compile`GetElement[x, i, 2], Compile`GetElement[x0, 2]]^2
        ],
        {i, Length@x}
    ],
    CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

On my machine vectorized solution has same speed as DistanceMatrix solution but uses less memory, listable compiled function is about two times faster and uses least amount of memory, library function, extracted from compiled function with manual looping, is almost seven times faster.

SeedRandom@1;
x0 = {1., 1.} // Developer`ToPackedArray;
data = RandomReal[1, {10^7, 2}];

r1 = distance[x0, data]; // MaxMemoryUsed // RepeatedTiming
(* {0.946,  240003488} *)
r2 = First@DistanceMatrix[{x0}, data]; // MaxMemoryUsed // RepeatedTiming
(* {0.56,   490785120} *)
r3 = Sqrt@Total[Subtract[Transpose@data, x0]^2]; // MaxMemoryUsed // RepeatedTiming
(* {0.56,   320001120} *)
r4 = distanceListable[x0, data]; // MaxMemoryUsed // RepeatedTiming
(* {0.351,   80004264} *)
r5 = distanceLib[x0, data]; // MaxMemoryUsed // RepeatedTiming
(* {0.0815, 160000440} *)

r1 === r2 === r3 === r4 === r5
(* True *)
$\endgroup$
  • $\begingroup$ I just realized that I didn't use constant passing for my library function. Doing so gave a much larger speedup than I expected. Maybe that could have something to do with why the library function extracted form the compiled function is so fast? (But it doesn't make sense.) $\endgroup$ – Szabolcs Nov 8 '17 at 15:35
  • $\begingroup$ I kind of wish you had benchmarked the LibraryLink version too, especially if you're on Windows :-) $\endgroup$ – Szabolcs Nov 8 '17 at 15:48
  • $\begingroup$ @Szabolcs distanceLib does use "Constant" argument passing. I thought that difference between CompiledFunction and extracted from it LibraryFunction comes from the fact that CompiledFunction is doing something with function result which requires making copy of it. $\endgroup$ – jkuczm Nov 8 '17 at 15:53
  • $\begingroup$ Yes, that's what I meant, I just didn't phrase it very well. $\endgroup$ – Szabolcs Nov 8 '17 at 15:59
  • $\begingroup$ @Szabolcs I'm on Linux, LTemplate version is surprisingly slow on my computer: 0.230 s, only 2.4 times faster than DistanceMatrix. Maybe I have "weak" default compiler optimizations. $\endgroup$ – jkuczm Nov 8 '17 at 16:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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