# How to group 10k points into closest pairs?

I was first faced with a problem of grouping 1000 pairs into unique pairs to minimise the total distance - since this was a small number I solved it using linear optimisation:

Part 1 - Basic/Naive approach:

A bit more detail on it here but the gist is:

Create an objective function which wants to minimise all the entries in the lower triangle of a distance matrix (multiplied by the respective dummy). For example v1_10 means that point 1 and point 10 are connected. Such that each store is connected to only one (other store).

The result of the optimisation takes us from the first picture to the second. The structure of the code is thus:

``````shops = [[id, latitude, longitude],...]
``````

The helper functions to create the optimisation:

``````def store_connect(shops):
"""
Cross stores and return connection e.g. "v1-2" by ID
"""
for i in range(len(shops)-1):
for j in range(i+1,len(shops)):
yield 'v'+str(shops[i][0])+'-'+str(shops[j][0])

def distances(shops):
"""
Return distance in miles for crosses
"""
for i in range(len(shops)-1):
for j in range(i+1,len(shops)):
yield haversine([shops[i][1],shops[i][2]],
[shops[j][1],shops[j][2]])

def all_constraints(shops):
"""
Return constraint row
"""
for a in shops:
cons = []
for b in shops:
if a[0]<b[0]:
cons.append("v"+str(a[0])+"-"+str(b[0]))
elif a[0]>b[0]:
cons.append("v"+str(b[0])+"-"+str(a[0]))
yield cons
``````

And then the optimisation itself:

``````prob = LpProblem("Minimise distance",LpMinimize)

# Variables
x = []
xnames = []
for one in store_connect(shops):
xnames.append(one)
x.append(LpVariable(one,0,1,LpInteger))
print("Created %d variables" % len(x))

# Objective
prob += pulp.lpSum([i[1]*x[i[0]] for i in enumerate(distances(shops))])
print("Created objective")

# Constraints
counter = 0
for constraint_rw in all_constraints(shops):
prob += pulp.lpSum(x[xnames.index(one_cn)] for one_cn in constraint_rw) == 1
counter += 1
if counter % 10 == 0:
print("Added %d contraints out of %d" % (counter,len(shops)))
``````

Which is straight-forward to solve directly or upload the .lp file to CPLEX:

``````#Using PULP - if not too complicated ... otherwise CPLEX:
prob.solve()
# The status of the solution is printed to the screen
print("Status:", LpStatus[prob.status])

# Optimum variable values
sol = []
for v in prob.variables():
if v.varValue:
sol.append(v.name)
print("Total minimised distance: ", value(prob.objective))

shp_names = [row[0] for row in shops]
edges = []
for x in sol:
v_list = (x)[1:].split("_") #remove the first letter 'v' and split
conn = [it for sub in [shops[shp_names.index(float(v_list[0]))][1:3],
shops[shp_names.index(float(v_list[1]))][1:3]] for it in sub]
edges.append(conn)

# Print diagram
plot_points(shops,edges)
``````

Part 2 - Better approach (for lots of points):

If I expand the problem to 10,000 points I get roughly 50,000,000 variables. I was thinking of the following ways of simplifying the problem:

1) Use clustering to group the 10,000 points into something like 100 groups and then solve 100 times within those groups (however I'm not sure how to check whether the solution would be optimal).

2) Project the points; put them into a quad-tree and then instead of allowing connections between one point and any other point -> only allow connections to the closest 100 points (e.g.): however, I am not sure how to properly set this (the cut-off). I thought perhaps if a point is chosen that is the cut-off point then the algorithm will restart with a larger cut-off.

3) Run something like a minimum spanning tree or concorde tsp algorithm and then back-track the solution from that?

• Are you looking for a Code Review or is there a specific something that is not working? The latter would suit the focussed Q&A that we do here but I think the former is too broad.
– PolyGeo
Commented Feb 5, 2016 at 10:24
• @PolyGeo I don't think code review is appropriate because I have something that 'works' only for a small number of points. However I would like an approach that works for 10k points (or maybe even 100k points); which I think turns this into a different beast entirely. Commented Feb 5, 2016 at 10:26
• It still seems to me that you are describing something that works for a small number of points but has performance that degrades severely as larger numbers of points are thrown at it. Code Review says that "If you are looking for feedback on a specific working piece of code from your project in the following areas… •Performance" then your question would be on-topic there.
– PolyGeo
Commented Feb 5, 2016 at 10:32
• @PolyGeo I think the problem is that this works in n^2 time. So after 1000 points it doesn't just 'degrade' but stops working at all. So I was looking for a completely different approach which would scale to 10,000 points (it take a long time but it would work ... not like 100 years). However, if you still think that this should be Code Review I can repost - or can this be moved? Commented Feb 5, 2016 at 10:49

My suggestion is based on a method I apply to group smaller subcatchments into larger.

There are no streams, thus 1st step to compute Euclidean minimum spanning tree:

I guess it is some sort of optimisation already.

Pick your sink and create directed graph. Picture below shows "Flow Accumulation" in links, i.e. count of nodes discharging into it:

Pick links with flow greater or equal to desired "subcatchment" size, 10 in this case. Prune graph (remove nodes) at "source" branches. Picture show remaining nodes:

Nodes of each pruned subgraph form new "subcatchments". As one can see at 1st iteration 4 larger subcatchments were formed.

Repeat calculation of "Flow Accumulation" and pruning until no nodes left. Nodes at the picture below labelled by their parent subcatchment size. As one can see it ranges from 5 to 15, although target was set to 10.

Start making your pairs, going downstream. If count of nodes is even you are lucky, if not just add downstream subcatchment to the group.

Anyway it seem reasonable to use triangulation inside groups to compute distance matrix:

Hopefully it will reduce number of parameters to optimise

UPDATE ON ORIGINAL SOLUTION

It took me 15 minutes to regroup original split of 500 points:

Into series each containing even count of points:

• This seems to rely on luck, guesswork, and intuition. It also seems it could leave an awful lot of orphan (unpaired) points behind. Could you supply any reasons why it should work, if only approximately? Commented Feb 6, 2016 at 16:30
• @whuber yes it will most likely produce what they call sub-optimal solution. Some manual intervention required to eliminate orphans using simple rule (reason you want?) that is total (or average) of 2 odd numbers result in even number. See updated answer Commented Feb 7, 2016 at 4:45