Would FindShortestTour be a better choice, (if I have read this correctly) as an example

AbsoluteTiming[SeedRandom[1]; pts = RandomReal[100, {100, 2}]; 
 pair = Partition[pts[[Last[FindShortestTour[pts]]]], 2]; 
 Print[Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}]]; 
 Total[Table[
   Sqrt[(pair[[q, 1, 1]] - pair[[q, 2, 1]])^2 + (pair[[q, 1, 2]] - 
      pair[[q, 2, 2]])^2], {q, 1, Length[pair]}]]]

or by rotating the points by 1 to get the other possible set of pairs

AbsoluteTiming[SeedRandom[1]; pts = RandomReal[100, {100, 2}]; 
 pair = Partition[pts[[Last[FindShortestTour[pts]]]], 2]; 
 pair = Partition[RotateLeft[Flatten[pair, 1]], 2]; 
 Print[Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}]]; 
 Total[Table[
   Sqrt[(pair[[q, 1, 1]] - pair[[q, 2, 1]])^2 + (pair[[q, 1, 2]] - 
      pair[[q, 2, 2]])^2], {q, 1, Length[pair]}]]]

in comparison to this method

AbsoluteTiming[SeedRandom[1]; pts = RandomReal[100, {100, 2}]; 
 g = CompleteGraph[Length[pts]]; 
 verGraph = VertexReplace[g, Thread[VertexList[g] -> pts]]; 
 weighGraph = 
  Graph[verGraph, 
   EdgeWeight -> EuclideanDistance @@@ EdgeList[verGraph]]; 
 pair = Partition[FindHamiltonianPath[weighGraph], 2]; 
 Print[Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}]]; 
 Total[Table[
   Sqrt[(pair[[q, 1, 1]] - pair[[q, 2, 1]])^2 + (pair[[q, 1, 2]] - 
      pair[[q, 2, 2]])^2], {q, 1, Length[pair]}]]]

The shortest tour option is also much faster, 3000 points is under a second on my machine.

Could I know how do you get that url(http://community.wolfram.com/groups/-/m/t/1075278) about a reply under a post?

Thanks for your endeavor for my this question.And I think a good message deserve to bring to you.I receive a good solution here based on FindShortestTour.

I have a method based on FindHamiltonianPath.It has a not very bad efficiency,which can find a pair(it is not always a global minimum pair as my test) from 100 points within 1s

SeedRandom[1]
pts = RandomReal[100, {100, 2}];
g = CompleteGraph[Length[pts]];
verGraph = VertexReplace[g, Thread[VertexList[g] -> pts]];
weighGraph = 
  Graph[verGraph, 
   EdgeWeight -> EuclideanDistance @@@ EdgeList[verGraph]];
pair = Partition[FindHamiltonianPath[weighGraph], 2];
Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}]

It is not a perfect solution,because a global minimum pair always is ecpected actually.

The blue ones always is the real minimum pair?

Brilliant try.Thanks very much.

But if I give some test pts for test,I found it not always give a right answer?Such as

SeedRandom[1]
pts = RandomReal[100, {8, 2}]

{{81.7389,11.142},{78.9526,18.7803},{24.1361,6.57388},{54.2247,23.1155},{39.6006,70.0474},{21.1826,74.8657},{42.2851,24.7495},{97.7172,82.5163}}

If we use a brute force method to find a definitely right answer:

First[pairs = 
  MinimalBy[Partition[#, 2] & /@ Permutations[pts], 
   Total[EuclideanDistance @@@ N[#]] &]]
First[Total[EuclideanDistance @@@ N[#]] & /@ pairs]

{{{81.7389,11.142},{78.9526,18.7803}},{{24.1361,6.57388},{42.2851,24.7495}},{{54.2247,23.1155},{97.7172,82.5163}},{{39.6006,70.0474},{21.1826,74.8657}}}

126.475(This is total distance)

But If we use your FindPairs

pts[[#]] & /@ FindPairs[pts]
Total[EuclideanDistance @@@ N[pts[[#]] & /@ FindPairs[pts]]]

{{{24.1361,6.57388},{54.2247,23.1155}},{{81.7389,11.142},{78.9526,18.7803}},{{39.6006,70.0474},{21.1826,74.8657}},{{42.2851,24.7495},{97.7172,82.5163}}}

141.565(This is total distance)