13
$\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$
9
  • 4
    $\begingroup$ Try EuclideanDistance[x1,#]&/@y1! $\endgroup$ Feb 2, 2013 at 17:42
  • 2
    $\begingroup$ You need to read up on Map. $\endgroup$
    – Szabolcs
    Feb 2, 2013 at 17:51
  • 1
    $\begingroup$ Related: Calculating the distance between the coordinates in a list $\endgroup$
    – Yves Klett
    Feb 2, 2013 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$ Feb 2, 2013 at 18:27
  • 2
    $\begingroup$ Also related $\endgroup$ Feb 2, 2013 at 18:27

5 Answers 5

11
$\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$
5
  • $\begingroup$ Just a suggestion: formatting code inline with text makes copy&pasting less convenient compared to separate code lines. $\endgroup$
    – Yves Klett
    Feb 2, 2013 at 20:38
  • $\begingroup$ @Yves, sorry I'm new to StackExchange, I'll try to edit it. $\endgroup$
    – RunnyKine
    Feb 2, 2013 at 20:58
  • $\begingroup$ You don't even need the Table. Just put the entire list of x1 into Outer. $\endgroup$
    – Jens
    Feb 2, 2013 at 21:01
  • $\begingroup$ No worries whatsoever, welcome to the party! $\endgroup$
    – Yves Klett
    Feb 2, 2013 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, 2013 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$
4
  • $\begingroup$ +1 Perhaps add link to your LTemplate (to save time searching for a package that saves time...)? $\endgroup$ Nov 8, 2017 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, 2017 at 13:01
  • 2
    $\begingroup$ @bobthechemist Just added on RPi benchmark for you. $\endgroup$
    – Szabolcs
    Nov 8, 2017 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$ Nov 8, 2017 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$
11
  • $\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, 2017 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, 2017 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, 2017 at 15:53
  • $\begingroup$ Yes, that's what I meant, I just didn't phrase it very well. $\endgroup$
    – Szabolcs
    Nov 8, 2017 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, 2017 at 16:31

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.