# Faster "Closest Pair of Points Problem" implementation?

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(n2) 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(n2) 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(n2) 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}}}}


 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.

• Are you aware of DistanceMatrix[]? Jan 8, 2021 at 23:09
• I tried that, but my implementation using DistanceMatrix[] is much slower (slower even than my nearestPairN), and seems to build the O(n_^2) array of distances. I particularly want a O(_n log (n) ) implementation. Jan 8, 2021 at 23:52
• I do notice that you are building a whole Table, then throwing it away. Does using Do help any? Jan 9, 2021 at 0:57
• @b3m2a1 Uses the Do to speed it up. Thanks! Jan 9, 2021 at 15:48

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}}}

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}}}

• Hah! (+1) And I was under the impression that I had found a new "good" feature for myself! ;) But honestly: Why the heck is constructing the Kd-tree twice still faster?! I realize that Nearest[pts->{"Index","Distance"}] does not allow packing. But still, this seems to be a design flaw. (Indeed, I was originally quite surprised that the Kd-tree was so slow...) Jan 9, 2021 at 17:26
• @HenrikSchumacher You might want to file a suggestion report. Unless the slowness is in a single unpacking, it can very likely be improved. Jan 10, 2021 at 0:11
• @DanielLichtblau Done. Jan 10, 2021 at 0:27
• @HenrikSchumacher Thanks. On a more general note, I am on the fence about how to approach the speed gap between your response and this one by Carl Woll. If one ignores underlying realities of packed array efficiencies, yours "should" be optimal. So one could view this behavior as a deficiency. But maybe we should instead take Carl's improvements, which are significant, as a paradigm for best practice when speed is paramount. In any case, having the suggestion on record, especially from an external user, is a good thing. Jan 10, 2021 at 15:59

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

• "seeing that one had to pass a dozen of parser layers" - I imagine it's for flexibility purposes, but is otherwise not great for squeezing out performance. Jan 9, 2021 at 13:40
• I upvoted but I'm not sure why. It might have been because this saved me the time of writing what would have been substantially similar code. But then I think, maybe it was the fortitude shown in piercing the Graph implementation armor to get, eventually, to its core (which, tellingly, I first mistyped as "code"). Jan 9, 2021 at 14:32
• As @J.M. said, flexibility. General-purpose, exciting a user when ˢhe finds hᵉᵢʳₘself allowed to use a Graph as a node in another Graph. Encapsulation, freeing you from (or otherwise interpreted, hindering you from) looking at the [commercial,] internal core of the implementation. Jan 29, 2021 at 7:10

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.