# Find the nearest locations for multiple points

Assume that there are many holes with their locations fixed, and the same number of balls distributed randomly. What is the smallest total distance for the balls fitting into the holes on the precondition that each hole can only fit one ball? For instance, the holes(black dots) are regularly distributed, and the balls(red dots) are randomly distributed. The nearest holes of the individual balls are indicated by arrows.

holes = Tuples[Range[1, 2, 1], 2];
balls = RandomReal[{1, 2}, Dimensions[holes]];
Graphics[{PointSize[Large], Point[holes], Red,PointSize[Medium],Point[balls]}]


Thanks for all the helps and answers. The problem is called 'The Euclidean matching problem' or 'Euclidean minimum weight matching problem' 1. I recently found an approximate alogrithm which achieves nearly O(n) time complexity [2].

[2] A Near-Linear Constant-Factor Approximation for Euclidean Bipartite Matching

This is an instance of the assignment problem, which is a special case of the minimum-cost flow problem, which can be solved directly in Mathematica.

n = {5, 5};
SeedRandom[1234];
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;


Construct the bipartite graph between balls and holes with edge costs equal to the distances between them, and add two dummy "source" and "target" vertices. Strangely, this is the most time-consuming part.

graph = Graph[
Flatten@Table[
Property[ball[i] \[DirectedEdge] hole[j],
EdgeCost -> EuclideanDistance[balls[[i]], holes[[j]]]],
{i, Length@balls}, {j, Length@holes}]
~Join~
Table[Property[source \[DirectedEdge] ball[i], EdgeCost -> 0], {i, Length@balls}]
~Join~
Table[Property[hole[j] \[DirectedEdge] target, EdgeCost -> 0], {j, Length@holes}]];


Solve the minimum-cost flow problem.

assignments =
Cases[FindMinimumCostFlow[graph, source, target, "EdgeList"],
ball[_] \[DirectedEdge] hole[_]]
(*{ball[1] -> hole[18], ball[2] -> hole[15], ball[3] -> hole[1],
ball[4] -> hole[8], ball[5] -> hole[2], ball[6] -> hole[25],
ball[7] -> hole[16], ball[8] -> hole[11], ball[9] -> hole[10],
ball[10] -> hole[22], ball[11] -> hole[23], ball[12] -> hole[5],
ball[13] -> hole[6], ball[14] -> hole[24], ball[15] -> hole[12],
ball[16] -> hole[4], ball[17] -> hole[19], ball[18] -> hole[9],
ball[19] -> hole[21], ball[20] -> hole[13], ball[21] -> hole[3],
ball[22] -> hole[14], ball[23] -> hole[17], ball[24] -> hole[20],
ball[25] -> hole[7]} *)


Visualize the result.

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium], Point[balls],
Line[assignments /. ball[i_] \[DirectedEdge] hole[j_] :> {balls[[i]], holes[[j]]}]}]


• By the way, is there anything we can do about all the ugly instances of \[DirectedEdge] up there?
– user484
Oct 23, 2014 at 15:53
• -> instead of \[DirectedEdge]? Oct 23, 2014 at 18:29
• By the way, can Sqrt@Total[Outer[Plus,balls,-holes,1]]^2,{3}] increase the speed? Oct 23, 2014 at 19:40
• @ybeltukov: No, if I do edges = Flatten@Table[...]~Join~... it is instantaneous while graph = Graph[edges] is the slow part.
– user484
Oct 23, 2014 at 20:19
• Indeed... May be AdjacencyGraph? Because you almost construct the adjacency matrix. Your solution is very good but I believe it is possible to make it even better. Oct 23, 2014 at 20:35

A faster version of Rahul's Answer

This question and the answers herein have just helped me solve a very similar problem using a solution based on @Rahul's answer which I find very elegant.

However, as is discussed briefly in the comments on that answer the code given (though very easy to understand) starts to run slow for larger numbers of balls and holes. I needed to perform solve tens of assignment problems with upwards of 50 balls and holes in each. Rahul's code was taking around 12 seconds to construct the Graph for the 25 ball case, the solution below runs in 0.0012 seconds - I think faster than any of the others - I post it here for people looking for a fast solution in the future!

The principle is as @ybeltukov suggests that rather than generate the graph we simply write down the weighted adjacency matrix for the situation and pass it to FindMinimumCostFlow as a cost matrix:

    SourceTargetCostMatrix[pointsA_, pointsB_] :=
Module[{lA = Length[pointsA], lB = Length[pointsB]},
ArrayFlatten@{
{0, ConstantArray[1, {1, lA}], ConstantArray[0, {1, lB}], 0},
{ConstantArray[0, {lA, 1}], ConstantArray[0, {lA, lA}],
Outer[EuclideanDistance, pointsA, pointsB, 1],
ConstantArray[0, {lA, 1}]},
{ConstantArray[0, {lB, 1}], ConstantArray[0, {lB, lA}],
ConstantArray[0, {lB, lB}], ConstantArray[1, {lB, 1}]},
{0, ConstantArray[0, {1, lA}], ConstantArray[0, {1, lB}], 0}
}
]

costMatrix = SourceTargetCostMatrix[balls, holes];

assignments = Cases[
FindMinimumCostFlow[costMatrix, 1, Length[costMatrix], "EdgeList"],
x_ \[DirectedEdge] y_ /; x != 1 && y != Length[costMatrix]
];

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],  Point[balls],
Line[assignments /. i_ \[DirectedEdge] j_ :> {balls[[i - 1]], holes[[j - Length[balls] - 1]]}]
}]


It can solve a 25x25 grid in just over a second, roughly a quarter of the time is for SourceTargetCostMatrix, the remaining is FindMinimumCostFlow:

Different Cost Functions

My problem required more priority placed on assigning those balls nearest to holes to that corresponding hole at the cost of having a few balls a very long way from a hole. I therefore used the Log of the EuclideanDistance which worked very well. To do so obviously simply replace EuclideanDistance in SourceTargetCostMatrix by whatever cost function you would like.

• A small bug detected~ When a point p exists in both balls and holes then the EuclideanDistance between this two ps would be 0, indicating there's no connection between them! It's better to add a small distance offset in the distance matrix: Outer[EuclideanDistance, pointsA, pointsB, 1]+1.*^-10
– Wjx
Oct 11, 2017 at 8:14
• There is a bug in MMA 11.3 and earlier that for certain ball/hole matrices, the FindMinimumCostFlow never finishes its calculations [CASE:4156292]. Sep 11, 2018 at 2:24
• The problem occurs when there are identical balls or holes. The work-around is before generating the costMatrix to do If[Length[balls]>Length[Union[balls]],balls=balls+RandomReal[{-0.001,0.001},Dimensions[balls]]]. And likewise for the holes. Sep 13, 2018 at 15:24

I will crib shamelessly from example and code for illustrating by @ybeltukov.

The example:

n = {5, 5};
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;


We can solve this as a linear programming problem. It looks like an integer linear program, but these are known to be solvable as relaxations thereof, that is, solutions to the relaxed LP will be integer valued (provided the solution is unique).

We set up the problem as below, to use FindMinimum. That will invoke linear programming. I find it easier to formulate in terms of variables rather than explicit matrix and vector constraints.

len = Length[holes];
vars = Array[x, {len, len}];
fvars = Flatten[vars];
c2 = Thread[Total[vars, {2}] == 1];
c3 = Map[0 <= # <= 1 &, fvars];
dists = Table[
vars[[j, k]]*EuclideanDistance[balls[[j]], holes[[k]]], {j,
len}, {k, len}];
obj = Total[dists, 2];


Now we solve it. The option setting is for speed. When the problem size is a bit larger than this, it will use interior point anyway, but for this size the automatic mode makes a slower choice.

SetOptions[LinearProgramming, Method -> "InteriorPoint"];

{min, vals} = FindMinimum[{obj, Flatten[{c1, c2, c3}]}, fvars];

res = Position[Round[vars /. vals], 1, 2]

(* {{1, 10}, {2, 8}, {3, 7}, {4, 23}, {5, 20}, {6, 2}, {7, 6}, {8,
11}, {9, 18}, {10, 24}, {11, 1}, {12, 19}, {13, 25}, {14, 14}, {15,
17}, {16, 22}, {17, 4}, {18, 15}, {19, 5}, {20, 12}, {21, 21}, {22,
13}, {23, 16}, {24, 3}, {25, 9}} *)


The picture:

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ res]}]


(Disclosure: Had this looked incorrect, I would have thrown transposes into the formulation of the objective until I got it right.)

• It seems to fail on big problems (above 20 by 20 seems to be enough): "Could not process the constraints..." and "Lists {min,vals} and FindMinimum[{obj,Flatten[{c1,c2,c3}]},fvars] are not the same shape. >>" Jan 5, 2015 at 12:34
• @P.Fonseca Some experimentation indicates you might do better using NMinimize in place of FindMinimum. I do not know if the result is guaranteed to be the same and I am looking into this. Also I would use Dispatch[vals] in the replacement because otherwise it is quite slow. Jan 5, 2015 at 19:55
• @P.Fonseca What version of Mathematica are you using? Jan 5, 2015 at 19:56
• 10.0.2 win64 sp1. I don't think i copied something wrong (since for less than 20, it worked fine), but if you can't repeat my problem (I first started with 50, and then went down to a rouded working limit), i can give it another try tomorrow, for a better insight. Jan 5, 2015 at 21:08
• @P.Fonseca I'm looking into some issues that arose in 10.0.2, most notably a slowdown. I am not sure if the failure to converge that I can make happen (with some obscure settings at least) is indicative of anything wrong or just a difficult optimization. At 20x20 we have 400^2 variables, which is not necessarily overwhelming for IP, but certainly not a small problem. Jan 6, 2015 at 15:31

Minimization of the total Euclidean distance from balls to holes is quite difficult optimization problem. As a start point I propose a simple greedy algorithm:

• Find the nearest hole for every ball.
• Put balls to holes starting from the closest ball-hole pair while the corresponding holes are empty
• Repeat

It is not the best algorithm and may be I will update my post in the future.

n = {5, 5};
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium], Point[balls]}]


res = Flatten[#, 2] &@ Last@Reap@
Module[{h = holes, b = balls, nh, nb, nf, num, put, fill,
nondup},
nb = nh = Range@Length@h;
While[h != {},
nf = Nearest[h -> Automatic];
num = nf /@ b // Flatten;
put = Ordering@Total[(h[[num]] - b)^2, {2}];
fill = num[[put]];
nondup = Floor@BinarySearch[Range@Length@fill, 1/2,
1 - Boole@DuplicateFreeQ@fill[[;; #]] &];
put = put[[;; nondup]];
fill = fill[[;; nondup]];
Sow@Transpose@{nb[[put]], nh[[fill]]};
b = Delete[b, Transpose@{put}];
nb = Delete[nb, Transpose@{put}];
h = Delete[h, Transpose@{fill}];
nh = Delete[nh, Transpose@{fill}];
]
]
(* {{8, 24}, {13, 5}, {19, 16}, {23, 10}, {1, 23}, {3, 17}, {5,
14}, {16, 21}, {17, 3}, {9, 8}, {7, 4}, {18, 11}, {4, 22}, {21,
15}, {11, 20}, {2, 6}, {14, 2}, {10, 7}, {15, 9}, {25, 1}, {12,
25}, {6, 18}, {22, 19}, {20, 13}, {24, 12}} *)

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ res]}]


• @ ybeltukov, your answer is helpful and enlightening.The first step is ‘Closest pair of points problem’. But a simple counterexample is for 2 hole-ball pairs with the given coordinates: holes = {{0, 0}, {0, 2}}; balls = {{0, 1.01}, {2, 2}}; I think the farest balls should be considered first. Oct 23, 2014 at 8:32
• @novice Anyway, greedy algorithm is bad. It can't get the global minimum. I found that Hungarian algorithm is exactly what you want. It is a polynomial algorithm. Unfortunately I have no time to implement it now. This Wolfram demonstration can also be helpful. Oct 23, 2014 at 11:33
• The Hungarian algorithm solves what is known as the assignment problem. Apparently one can solve such problems directly in Mathematica since version 9.
– user484
Oct 23, 2014 at 11:55
• @RahulNarain I knew that it should correspond to a graph problem! Very nice :) Oct 23, 2014 at 18:31
• Because of triangle inequality, no two edges of the optimal matching can cross each other. /*A Survey on Algorithms for Euclidean Matching*/ Jan 6, 2015 at 0:53

Ok this may end up being a greedy algorithm, though it was based on an old Dijkstra algorithm that I modified, but since balls and holes are not connected in any way is just greedy.

Since I don't have V10 installed here, I don't have the same BinarySearch[] functions (Must be different from Combinatorica), I couldn't compare against @ybeltukov. But I bet my is slower due to the indexing of distances, as a N^2 matrix, and finding the nearest using Min on it. However it might be educative to see another one.

Using his sample data,

n = {5, 5};
SeedRandom[1234];
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;


If you find another way to index the distance, it may be faster.

dist = Table[
Norm[holes[[i]] - balls[[j]]], {i, 1, Length[holes]}, {j, 1,
Length[balls]}];
filledholes = ConstantArray[False, Length[holes]];
balldroped = ConstantArray[False, Length[balls]];
parm = Max[dist];

a = Reap[
While[And @@ filledholes != True,
posmin = First@Position[dist, Min[dist]];
If[filledholes[[posmin[[1]]]] != True &&
balldroped[[posmin[[2]]]] != True,
filledholes[[posmin[[1]]]] = True;
balldroped[[posmin[[2]]]] = True;
Sow[{posmin[[2]], posmin[[1]]}];
dist[[First@posmin, Last@posmin]] = parm + 1.0;
,
dist[[First@posmin, Last@posmin]] = parm + 1.0;
Continue[];
];
];]

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ (Last@Last@a)]}]


The total distance being

Plus @@ (Norm[balls[[#]] - holes[[#2]]] & @@@ (Last@Last@a))


22.9346

With Timing it is

0.039872

In case you want maximize the distance, then easy to replace Min for Max, vis versa, and the adjustment to the index. In ybeltukov answer you would weighted the results with 1/nf.

dist = Table[
Norm[holes[[i]] - balls[[j]]], {i, 1, Length[holes]}, {j, 1,
Length[balls]}];
filledholes = ConstantArray[False, Length[holes]];
balldroped = ConstantArray[False, Length[balls]];
parm = Min[dist];
b = Reap[
While[And @@ filledholes != True,
posmax = First@Position[dist, Max[dist]];
If[filledholes[[posmax[[1]]]] != True &&
balldroped[[posmax[[2]]]] != True,
filledholes[[posmax[[1]]]] = True;
balldroped[[posmax[[2]]]] = True;
Sow[{posmax[[2]], posmax[[1]]}];
dist[[First@posmax, Last@posmax]] = 0.99*parm;
,
dist[[First@posmax, Last@posmax]] = 0.99*parm;
Continue[];
];
];]

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ (Last@Last@b)]}]


The total distance being

Plus @@ (Norm[balls[[#]] - holes[[#2]]] & @@@ (Last@Last@b))


101.006

This may help:

locations=RandomInteger[{-10,10},{4,2}];
holes=RandomInteger[{-10,10},{4,2}];
ListPlot[{locations, holes}, PlotLegends -> {"locations", "holes"}]


distanceFunc = Nearest[holes, DistanceFunction -> EuclideanDistance];
N@Total[EuclideanDistance[#, First@distanceFunc[#]] & /@ locations]

• thanks for your answer, what if each hole can only take one ball? Oct 22, 2014 at 7:49
• @novice ahh, I miss that Oct 22, 2014 at 8:06