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

Question

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[1]
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[20];
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.

Answer 1

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[1]
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[#[[1]] +  #2[[1]] == 0,
       -bound,
       If[
        #[[1]] #2[[1]] < 0 
        ,
        bound
        ,
        EuclideanDistance[#[[2 ;;]], #2[[2 ;;]]]
        ]] &), Method -> "IntegerLinearProgramming"];

resultLength = nn*bound + N@First@tour

30.8036

pairs = pts[[#]] & /@ Partition[tour[[2]], 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

Answer 2

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[1]
pts = RandomInteger[20, {10, 2}];
g = CompleteGraph[10]

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

Answer 3

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

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

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

o = FindShortestTour[pts][[2]];

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

pairs = MinimalBy[%, tdist, 1][[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

Answer 4

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[1]
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

Answer 5

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[1]
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 ;; -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];
c1 = Thread[
   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[{#[[1]], #[[2]] + #[[1]]} &, 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.