# How to get a pair with smallest total distance from one list in high-efficiency method

This is similar to this post, but actually more complicated. I have two lists in that question, and I want to get a pair. But in here, I want to get a pair from one list (even points). Given these 20 points:

SeedRandom
pts = RandomInteger[20, {20, 2}]


{{5,0},{7,0},{2,3},{0,0},{16,14},{3,8},{19,5},{18,16},{12,0},{19,4},{7,3},{0,4},{20,3},{5,12},{19,8},{11,2},{3,10},{4,2},{17,11},{15,6}}

I can use FindIndependentEdgeSet to get one pair like this:

g = CompleteGraph;
List @@@ FindIndependentEdgeSet[VertexReplace[g, Thread[VertexList[g] -> pts]]]


{{{5,0},{7,0}},{{2,3},{3,8}},{{0,0},{16,14}},{{19,5},{4,2}},{{18,16},{3,10}},{{12,0},{7,3}},{{19,4},{11,2}},{{0,4},{20,3}},{{5,12},{19,8}},{{17,11},{15,6}}}

The total distance is

Total[EuclideanDistance @@@ pair] // N


113.859

But I'm sure it is not the smallest pair. Actually, I think I need a FindIndependentEdgeSet of edge weight version, but it seems the FindIndependentEdgeSet regards weighted graph as unweighted directly. Can anyone give me advice? Of course, I'm happy to know other methods that can do this which aren't based on Graph Theory.

This post related this question.

• Are you interested only in a perfect solution, i.e. the minimal weight pairing, or would an optimization approach that produces a "good" pairing be useful too? – Mr.Wizard Apr 9 '17 at 7:36
• @Mr.Wizard I just expected the optimal paring and hope to speed up it. Thanks. :) – yode Apr 9 '17 at 9:21
• @Mr.Wizard A brute force look shows the min distance with SeedRandom to be 30.8036. This takes ~100s to run for the 20 point case and scales badly (Binomial[2 n, n])... – Quantum_Oli Apr 9 '17 at 11:43

There is a solution based on IGLargestIndependentVertexSets,which from Szabolcs's package IGraphM,when the point is less 15 or much less.I mean,the effeciency it very poor.I will produce just 10 points to test it.

### Built a weight graph with EuclideanDistance

SeedRandom
pts = RandomInteger[20, {10, 2}];
g = CompleteGraph ### Find a smallest independent vertex sets

iden = IGLargestIndependentVertexSets[LineGraph[g]];
pair = List @@@ First[MinimalBy[EdgeList[g][[##]] & /@ iden,
Total[EuclideanDistance @@@ N[#]] &]]


{{{5,0},{0,0}},{{7,0},{12,0}},{{2,3},{3,8}},{{16,14},{18,16}},{{19,5},{19,4}}}

N[Total[EuclideanDistance @@@ pair]]


18.9274

I hope this is a good start..

• @Szabolcs Could I expected IGLargestIndependentEdgeSets in future? – yode Apr 9 '17 at 7:07

FindShortestTour appears to be a sub-optimal but fast method.

SeedRandom
pts = RandomInteger[20, {20, 2}];

tdist = N@Total[EuclideanDistance @@@ #] &;

o = FindShortestTour[pts][];

Partition[# @ pts[[o]], 2] & /@ {Most, Rest};

pairs = MinimalBy[%, tdist, 1][]

tdist @ pairs

ListLinePlot @ pairs

{{{5, 0}, {4, 2}}, {{2, 3}, {0, 0}}, {{0, 4}, {3, 8}}, {{3, 10}, {5, 12}},
{{16, 14}, {18, 16}}, {{17, 11}, {19, 8}}, {{19, 5}, {20, 3}},
{{19, 4}, {15, 6}}, {{12, 0}, {11, 2}}, {{7, 3}, {7, 0}}}

32.0483 The same code with pts = RandomInteger[20, {10, 2}]; gives:

{{{7, 0}, {12, 0}}, {{19, 4}, {19, 5}}, {{18, 16}, {16, 14}},
{{3, 8}, {2, 3}}, {{0, 0}, {5, 0}}}

18.9274

• @yode I did not look at your self-answer closely enough! Why does your output only have ten points in it? Shouldn't it have all 20? I think I misunderstand your problem. :-( – Mr.Wizard Apr 9 '17 at 8:57
• Sorry,my solution is very slow,so I use RandomInteger[20, {10, 2}] to instead of 20. – yode Apr 9 '17 at 9:00
• @yode I see. :-) I guess you are interested in only the optimal paring and not one that is "close" -- is that correct? If so I should delete my answers as neither method is optimal. – Mr.Wizard Apr 9 '17 at 9:04
• Yes,I just expected the optimal paring and hope to speed up it. :) – yode Apr 9 '17 at 9:06
• @yode No, I do not think this will always be the shortest pairing, which is why I wrote sub-optimal. I suspect that to find the true shortest pairing every time will require exhaustive testing and be very slow, as this is true for the Traveling Salesman Problem to the best of my knowledge. – Mr.Wizard Apr 9 '17 at 9:36

I don't know if this is of any value at all but I was curious to see how an optimization approach would work, and this is what I came up with:

SeedRandom
pts = RandomInteger[20, {20, 2}];

tdist = N @ Total[EuclideanDistance @@@ #] &;

fn[pt_List, x_ /; VectorQ[x, NumberQ]] := tdist @ Partition[pt[[Ordering @ x]], 2]

vars = Unique @ Table["p", {Length @ pts}];

{min, sol} =
NMinimize[fn[pts, vars], vars, Method -> "SimulatedAnnealing",
MaxIterations -> 1000, PrecisionGoal -> 1];

pts[[Ordering[vars /. sol]]] ~Partition~ 2

tdist @ %

{{{12, 0}, {11, 2}}, {{0, 0}, {0, 4}}, {{2, 3}, {4, 2}}, {{3, 10}, {5, 12}},
{{18, 16}, {16, 14}}, {{17, 11}, {19, 8}}, {{3, 8}, {7, 3}}, {{15, 6},
{19, 5}}, {{5, 0}, {7, 0}}, {{19, 4}, {20, 3}}}

31.675


The same code starting with pts = RandomInteger[20, {10, 2}]; gives:

{{{19, 5}, {19, 4}}, {{2, 3}, {3, 8}}, {{12, 0}, {7, 0}},
{{0, 0}, {5, 0}}, {{18, 16}, {16, 14}}}

18.9274


Might not be terribly efficient but this can be done with integer linear programming. It is quite similar to a method shown here but unfortunately the relaxed problem in this case need not deliver an integer solution.

Here is the example.

SeedRandom
len = 20;
pts = RandomInteger[20, {len, 2}];


For each 1<=i<j<=20 we'll create a variable to denote whether or not to create an edge between the ith and jth points. So these are intended as 0-1 variables. We impose constraints that every point be hit by exactly one edge. The objective function will be the sum of all distances between points linked by an edge.

ovars = Array[x, {len, len}];
vars = MapIndexed[#1[[#2[] + 1 ;; -1]] &, ovars];
fvars = Flatten[vars];
zeroRule = Thread[Complement[Flatten[ovars], fvars] -> 0];

dists = Table[
vars[[j, k]]*EuclideanDistance[pts[[j]], pts[[j + k]]], {j,
len - 1}, {k, 1, Length[vars[[j]]]}];
obj = Total[dists, 2];
Map[Total, vars] + (Map[Total, Transpose[ovars]] /. zeroRule) == 1];
c2 = Map[0 <= # <= 1 &, fvars];


We solve it with FindMinimum.

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


Now deduce which edges get used.

ores = Position[Round[vars /. vals], 1, 2];
res = Map[{#[], #[] + #[]} &, ores]


Again cribbing from someone else's graphics, we show the result.

Graphics[{PointSize[Large], Point[pts], Red, PointSize[Medium],
Arrow[{pts[[#]], pts[[#2]]} & @@@ res]}] This method handles 150 points in a half minute or so. One bottleneck is the computation of the objective function.

• I have seen your answer many times,give me some time to check the FindMinimum always give a global minimal pair.I will response you after some days. – yode May 4 '17 at 0:18

This gives the same (optimal) result as Daniel Lichtblau's answer for the given example. Integer linear programming is also used, but in this answer this is done by FindShortestTour.

max = 20;
nn = 20;
SeedRandom
pts = RandomInteger[max, {nn, 2}]
boolPts =
Join[
Transpose@Join[{-Range[nn]}, Transpose@ pts],
Table[{i, 0, 0}, {i,  nn}]];
bound = 2 max^2;
tour = FindShortestTour[
boolPts,
DistanceFunction -> (If[#[] +  #2[] == 0,
-bound,
If[
#[] #2[] < 0
,
bound
,
EuclideanDistance[#[[2 ;;]], #2[[2 ;;]]]
]] &), Method -> "IntegerLinearProgramming"];

resultLength = nn*bound + N@First@tour

30.8036

pairs = pts[[#]] & /@ Partition[tour[], 2][[1 ;; -1 ;; 2]]
resultLength == Total[EuclideanDistance @@@ pairs]

{{{5, 0}, {7, 0}}, {{5, 12}, {3, 10}}, {{0, 0}, {2, 3}}, {{11,
2}, {12, 0}}, {{19, 4}, {20, 3}}, {{17, 11}, {19, 8}}, {{15,
6}, {19, 5}}, {{16, 14}, {18, 16}}, {{7, 3}, {4, 2}}, {{3, 8}, {0,
4}}}
True

• Some question should be listed.Firstly,the boolPts can be Catenate[MapIndexed[{PadRight[#2, 3], Prepend[#1, -First[#2]]} &, pts]]?Secondly,why we set the bound is 2 max^2?.Thirdly,Do you sure the FindShortestTour always will give Global minimal pairs and have you see this?fourthly,could you explain your DistanceFunction? – yode May 4 '17 at 0:39
• The result is not same with this comment.And how to get the pair as your tour? – yode May 4 '17 at 0:55
• @yode 1) yes, you can use that Catenate code if you prefer. 2) bound is set heuristically to 2 max^2. It can be seen as a reward/penalty that has is so big that if you can get it as a reward, you will always choose to do so and similarly you always avoid it as a penalty. If you do not include a reward, my procedure may find triplets instead of pairs. It is also used as a penalty to "remove edges" between vertices. 3) I think it will always find the optimum, I think the reduction to a FindshortestTour (traveling salesman) problem is correct and... – Jacob Akkerboom May 4 '17 at 7:02
• FindShortestTour with the option Method -> "IntegerLinearProgramming" will give the optimal solution for the problem it is considering. 4) See my explanation under "2)". Re: "The result is not same...", well I get the same decimals right (30.8036)? Re: "how to get the pair(s)..", it is stored in the variable pairs. – Jacob Akkerboom May 4 '17 at 7:05
• @yode it was probably a bit silly of me to not use a Graph. I felt it might be more efficient to not have to calculate all distances between all points, but perhaps the time this takes is negligible anyway. If I make a graph the intended connections between vertices may be more apparent. – Jacob Akkerboom May 4 '17 at 7:10