Faster "Closest Pair of Points Problem" implementation?

Question

The closest pair of points problem is a common computational geometry problem: given n points, find a pair of points with the smallest distance between them. A naive algorithm of finding distances between all pairs of points and selecting the minimum requires O(n²) time. The problem may be solved in O(n log n) time in a Euclidean space. Rosetta Code has only the brute force solution for Mathematica. I coded up a O(n log n) version, but my implementation is slow, and only beats the O(n²) for n > 10,000 on my machine. How can this code be made faster? Is there a built-in Mathematica function for this?

Here is my O(n²) code:

nearestPairN[
  data_] := (*This is 10x faster than \
https://rosettacode.org/wiki/Closest-pair_problem#Mathematica_.2F_\
Wolfram_Language.*)
 Last @@ Module[{npt}, 
   MinimalBy[
    Table[npt = 
      N[First@Nearest[Drop[data, i], data[[i]], 
         1]]; {EuclideanDistance[data[[i]], npt], {data[[i]], 
       npt}}, {i, Length[data] - 1}], First]]

Here is my O(n log n) code:

closestPair[ptsIn_] := Module[{xP, yP, pts},
  (*Top level function. 
  Sorts the pts by x and by y and then calls closestPairR[] *)
  
  pts = N[ptsIn];
  xP = Sort[pts, #1[[1]] < #2[[1]] &];
  yP = Sort[pts, #1[[2]] < #2[[2]] &];
  closestPairR[xP, yP]
  ]

closestPairR[xP_, yP_] := 
 Module[{n, mid, xL, xR, xm, yL, yR, dL, pairL, dmin, pairMin, yS, nS,
    closest, closestP, k, cDist},
  (*where xP is P(1).. P(n) sorted by x coordinate, and
            
  yP is P(1).. P(n) sorted by y coordinate (ascending order) *)
  
  n = Length[xP];
  If[ n <= 3, (*Brute Force*)
   Piecewise[{   
     {{\[Infinity], {}}, n < 2},
     {{EuclideanDistance[xP[[1]], xP[[2]]], {xP[[1]], xP[[2]]}}, 
      n == 2},
     {Last@MinimalBy[{
         {EuclideanDistance[xP[[1]], xP[[2]]], {xP[[1]], xP[[2]]}},
         {EuclideanDistance[xP[[1]], xP[[3]]], {xP[[1]], xP[[3]]}},
         {EuclideanDistance[xP[[3]], xP[[2]]], {xP[[3]], xP[[2]]}}
         }, First], n == 3}
     }],
   mid =  Ceiling[n/2];
   xL = xP[[1 ;; mid]];
   xR = xP[[mid + 1 ;; n ]];
   xm = xP[[mid]];
   yL = Select[yP, #[[1]] <= xm[[1]] &]; 
   yR = Select[yP, #[[1]] > xm[[1]] &]; 
   {dL, pairL} = closestPairR[xL, yL];
   {dmin, pairMin} = closestPairR[xR, yR];
   If[dL < dmin, {dmin, pairMin} = {dL, pairL}];
   yS = Select[yP, Abs[#[[1]] - xm[[1]]] <= dmin &];
   nS = Length[yS];
   {closest, closestP} = {dmin, pairMin} ;
   Table[
    k = i + 1;
    While[(k <= nS) && (yS[[k, 2]] - yS[[i, 2]] < dmin), 
     cDist = EuclideanDistance[ yS[[k]], yS[[i]]];
     If[cDist < closest, 
      {closest, closestP} = {cDist, { yS[[k]], yS[[i]]}}
      ];
     k = k + 1]
    , {i, 1, nS - 1}];
   {closest, closestP}
   ](*end if*)]

And test code:

pts = RandomReal[1, {10000, 2}];
Timing[nearestPairN[pts]]
Timing[closestPair[pts]]

One sample run gave me:

{1.58398, {{0.696976, 0.770051}, {0.697001, 0.770091}}}

{1.54305, {0.000047247, {{0.696976, 0.770051}, {0.697001, 0.770091}}}}

[Edit] The three submitted solutions are compared below, using RepeatedTiming for n from 10k to 500k. @Carl's solution is the fastest, it uses built-in functions, and is succinct.

Answer 1

Here's a refinement of @Henrik's approach. The key difference is that using Nearest[pts->"Distance"] is over an order of magnitude faster than using Nearest[pts->{"Index", "Distance"}]:

SeedRandom[1];
pts = RandomReal[1, {10000,2}];

AbsoluteTiming[
    data=Nearest[pts->"Distance"][pts,2][[All,2]];
    d=Min[data];
    i=First @ Ordering[data,1];
    j=Nearest[pts->"Index"][pts[[i]],2][[2]];
    {d, pts[[i]],pts[[j]]}
]

{0.008782, {0.000108351, {0.644732, 0.0760254}, {0.64463, 0.0759882}}}

Compared to @Henrik's answer:

AbsoluteTiming[
    data=Nearest[pts->{"Index","Distance"}][pts,2][[All,2]];
    i=OrderingBy[data,Last,1][[1]];
    j=data[[i,1]];
    dist=data[[i,2]];
    result={dist,pts[[i]],pts[[j]]}
]

{0.149502, {0.000108351, {0.644732, 0.0760254}, {0.64463, 0.0759882}}}

Answer 2

Maybe this works for you:

AbsoluteTiming[
 data = Nearest[pts -> {"Index", "Distance"}][pts, 2][[All, 2]];
 i = OrderingBy[data, Last, 1][[1]];
 j = data[[i, 1]];
 dist = data[[i, 2]];
 result = {dist, pts[[i]], pts[[j]]}
 ]

i and j are the indices of the pair of points that are closest to each other; dist is their distance.

Nearest[pts -> {"Index", "Distance"}][pts, k] finds the for every point its k nearest neighbors. The first "neighbor" will always be the point itself. That's why I take Nearest[pts -> {"Index", "Distance"}][pts, 1][[All, 2]] to get the nearest "true" neighbor.

I am pretty sure that this does quite the same like C.E.'s NearestNeighborGraph implementation, but it is about 3 times faster on my machine and alsmost twice as fast as b3m2a1's compiled implementation (on my machine and for 10000 points). This is only faster than NearestNeighborGraph because it does not rely on Mathematica's Graph implementation and all its overhead. Yet another example where avoiding Graph boosts performance. But I have to say that NearestNeighborGraph catches up a bit when one increases the number of points.

The performance degradations start with Graph using a fancy interface that enforces conversion of packed lists of edge indices to unpacked lists of UndirectedEdges and DirectedEdges. And it goes on with all the fancy vertex and edge attributes. Once I had to dig through the half-hidden implementation details of some Graph feature for debugging. I got cold shivers seeing that one had to pass a dozen of parser layers before one arrives at the actual algorithms...

Answer 3

Straight-forward application of Compile to your code. No effort has been taken towards optimization. You'll need to evaluate this twice for Compile to pick up on the recursion. I flattened the return structure and made it explicitly 2D, but you can easily make it higher-dimensional if you want to by turning the 2;;5 into some calculated parameter.

simpleMin =
  With[
   {
    EuclideanDistance = (Sqrt@Total[(# - #2)^2] &)
    },
   Compile[
    {
     {n, _Integer}, 
     {xP, _Real, 2}
     },
    Module[{d1, d2, d3},
     Which[
      n == 2,
      Join[{EuclideanDistance[xP[[1]], xP[[2]]]}, xP[[1]], xP[[2]]],
      n == 3,
      d1 = EuclideanDistance[xP[[1]], xP[[2]]];
      d2 = EuclideanDistance[xP[[1]], xP[[3]]];
      d3 = EuclideanDistance[xP[[3]], xP[[2]]];
      If[d1 <= d2 && d1 <= d3,
       Join[{EuclideanDistance[xP[[1]], xP[[2]]]}, xP[[1]], xP[[2]]],
       If[d2 <= d1 && d2 <= d3,
        Join[{EuclideanDistance[xP[[1]], xP[[3]]]}, xP[[1]], xP[[3]]],
        Join[{EuclideanDistance[xP[[2]], xP[[3]]]}, xP[[2]], xP[[3]]]
        ]
       ],
      True,
      {-1., -1., -1., -1., -1.}
      ]
     ]
    ]
   ];
closestPairRC =
  With[
   {
    EuclideanDistance = (Sqrt@Total[(# - #2)^2] &),
    simpleMin = simpleMin
    },
   Compile[
    {
     {xP, _Real, 2},
     {yP, _Real, 2}
     },
    Module[{
      n, mid,
      xL, xR, xm, yL, yR,
      dL, pairL, dmin, pairMin,
      yS, nS, closest, closestP,
      k, cDist
      },
     (*where xP is P(1).. P(n) sorted by x coordinate,
     and yP is P(1).. P(n) sorted by y coordinate (ascending order)*)
     
     n = Length[xP];
     dL = -1.;
     pairL = {-1., -1., -1., -1., -1.};
     dmin = -1.;
     pairMin = {-1., -1., -1., -1., -1.};
     If[n <= 3,
      simpleMin[n, xP],
      (* hard to make the recursion work well inside `Compile`. 
      Might be a case for `FunctionCompile` *)
      mid = n;
      mid = Ceiling[n/2];
      xL = xP[[1 ;; mid]];
      xR = xP[[mid + 1 ;; n]];
      xm = xP[[mid]];
      yL = Select[yP, #[[1]] <= xm[[1]] &];
      yR = Select[yP, #[[1]] > xm[[1]] &];
      pairL = closestPairRC[xL, yL];
      dL = pairL[[1]];
      pairL = pairL[[2 ;; 5]];
      pairMin = closestPairRC[xR, yR];
      dmin = pairMin[[1]];
      pairMin = pairMin[[2 ;; 5]];
      If[dL < 0 || dL < dmin, dmin = dL; pairMin = pairL;];
      yS = Select[yP, Abs[#[[1]] - xm[[1]]] <= dmin &];
      nS = Length[yS];
      closest = dmin;
      closestP = pairMin;
      Do[
       k = i + 1;
       While[
        (k <= nS) && (yS[[k, 2]] - yS[[i, 2]] < dmin), 
        cDist = EuclideanDistance[yS[[k]], yS[[i]]];
        If[cDist < closest,
         closest = cDist; closestP = Join[yS[[k]], yS[[i]]]
         ];
        k = k + 1
        ],
       {i, 1, nS - 1}
       ];
      Join[{closest}, closestP]
      ](*end if*)
     ]
    ]
   ];

Corresponding uncompiled version

closestPairR =
  With[
   {
    EuclideanDistance = (Sqrt@Total[(# - #2)^2] &),
    simpleMin = simpleMin
    },
   Function[
    {xP, yP},
    Module[{
      n, mid,
      xL, xR, xm, yL, yR,
      dL, pairL, dmin, pairMin,
      yS, nS, closest, closestP,
      k, cDist
      },
     (*where xP is P(1).. P(n) sorted by x coordinate,
     and yP is P(1).. P(n) sorted by y coordinate (ascending order)*)
     
     n = Length[xP];
     dL = -1.;
     pairL = {-1., -1., -1., -1., -1.};
     dmin = -1.;
     pairMin = {-1., -1., -1., -1., -1.};
     If[n <= 3,
      simpleMin[n, xP],
      (* hard to make the recursion work well inside `Compile`. 
      Might be a case for `FunctionCompile` *)
      mid = n;
      mid = Ceiling[n/2];
      xL = xP[[1 ;; mid]];
      xR = xP[[mid + 1 ;; n]];
      xm = xP[[mid]];
      yL = Select[yP, #[[1]] <= xm[[1]] &];
      yR = Select[yP, #[[1]] > xm[[1]] &];
      pairL = closestPairR[xL, yL];
      dL = pairL[[1]];
      pairL = pairL[[2 ;; 5]];
      pairMin = closestPairR[xR, yR];
      dmin = pairMin[[1]];
      pairMin = pairMin[[2 ;; 5]];
      If[dL < 0 || dL < dmin, dmin = dL; pairMin = pairL;];
      yS = Select[yP, Abs[#[[1]] - xm[[1]]] <= dmin &];
      nS = Length[yS];
      closest = dmin;
      closestP = pairMin;
      Do[
       k = i + 1;
       While[
        (k <= nS) && (yS[[k, 2]] - yS[[i, 2]] < dmin), 
        cDist = EuclideanDistance[yS[[k]], yS[[i]]];
        If[cDist < closest,
         closest = cDist; closestP = Join[yS[[k]], yS[[i]]]
         ];
        k = k + 1
        ],
       {i, 1, nS - 1}
       ];
      Join[{closest}, closestP]
      ](*end if*)
     ]
    ]
   ];

I'm dropping the wrapper so we can directly compare the two

pts = BlockRandom[N@RandomReal[1, {10000, 2}]];
sortPts = {Sort[pts], SortBy[pts, Last]};
RepeatedTiming[closestPairRC @@ sortPts]
RepeatedTiming[closestPairR @@ sortPts]

{0.0898, {0.0000934546, 0.41575, 0.168734, 0.415831, 0.168781}}

{0.830, {0.0000934546, 0.41575, 0.168734, 0.415831, 0.168781}}

So it's like a factor of 10 faster. You might be able to get better performance with FunctionCompile since I know it's able to deal with recursion cleanly. You can also play with compilation options or any number of other things.