How to remove duplicates from set of machine precision 2D points?

Question

I have a set of 2D points with machine precision coordinates. I need to remove all duplicates. Performance is important.

This is the most obvious fast solution:

Union[points]

Unfortunately it turns out that unlike Equal and SameQ, Union has no tolerance for any difference with machine precision numbers, which makes it unusable for this purpose.

In[1]:= 
a = 1.;
b = a + $MachineEpsilon;

In[3]:= a == b
Out[3]= True

In[4]:= a === b
Out[4]= True

In[5]:= Union[{a, b}]
Out[5]= {1., 1.}

Using the SameTest option of Union works, but it is not a good option because it makes the complexity $n^2$.

So what is the fastest way to get rid of duplicates while allowing for a tolerance of comparisons, and preferably being able to control it (Internal`$EqualTolerance, Internal`$SameQTolerance)?

Answer 1

Here's another compiled implementation. It's only very slightly faster than the code of s0rce, which is due to this version making only one array access (rather than two) per inner loop iteration. It also uses the Internal`Bag, which may be advantageous for memory consumption in case there are many duplicates.

It must be said that this is still not as fast as Union, but at least it acknowledges Internal`$EqualTolerance while being faster than Split. The value of Internal`$EqualTolerance is actually hard-coded into the bytecode on compilation, so it will be necessary to recompile if a different tolerance is required.

compiledUnion = Compile[{{r, _Real, 2}},
   Block[{
     sorted = Sort[r], output,
     seen, current
    },
    output = Internal`Bag[seen = First[sorted], 1];
    Do[
     If[i != seen, Internal`StuffBag[output, seen = i, 1]],
     {i, sorted}
    ];
    Partition[Internal`BagPart[output, All], Length[seen]]
   ],
   RuntimeOptions -> {"Speed", "CompareWithTolerance" -> True},
   CompilationTarget -> "C"
  ];

Test case, with the default value of Internal`$EqualTolerance == 7 Log[2]/Log[10]:

compiledUnion[{
  {1., 1., 1.}, {1. + $MachineEpsilon, 1., 1.}, {1., 1. + $MachineEpsilon, 1.},
  {2., 1., 1.}
 }]
(* -> {{1., 1., 1.}, {2., 1., 1.}} *)

Performance comparison, in ascending order of run-time:

r = RandomReal[1, {10000000, 2}];
r = RandomChoice[r, 10000000];

DeleteDuplicates[r]; (* 1.813 seconds; incorrect result (no tolerance) *)
DeleteDuplicates[r, Equal]; (* same timing and result--incorrect special-casing *)
Sort[r]; (* 4.297 seconds; base case for Union-like approaches *)
Union[r]; (* 4.609 seconds; incorrect result (no tolerance) *)

compiledUnion[r]; (* 5.875 seconds *)
With[{sorted = Sort[r]}, Pick[sorted, deleteDuplicatesC[sorted], 1]]; (* 6.109 seconds *)

Split[Sort[r]][[All, 1]]; (* 11.953 seconds; unpacks + additional memory overhead *)
DeleteDuplicates[r, Equal[##] &]; (* 294.8 days; unpacks *)

Answer 2

One possible solution is using Split, which obeys Internal`$SameQTolerance

r = RandomReal[1, {10000000, 2}];
r = RandomChoice[r, 10000000];

Split[Sort[r]][[All, 1]]; // AbsoluteTiming

(* ==> {12.756712, Null} *)

This is about 3-4 times slower than Union:

Union[r]; // AbsoluteTiming

(* ==> {4.306363, Null} *)

Answer 3

How about a compiled solution. I think this solves the transitivity issues with Split.

deleteDuplicatesC = 
  Compile[{{v, _Real, 2}}, 
   Block[{i, len = Length[v], output = Table[1, {i, len}]},
    Do[If[Compile`GetElement[v, i] == Compile`GetElement[v, i - 1], 
      output[[i]] = 0], {i, 2, len}];
    output], 
   RuntimeOptions -> {"Speed", "CompareWithTolerance" -> True}, 
   CompilationTarget -> "C"];

r = RandomReal[1, {10000000, 2}];
r = RandomChoice[r, 10000000];

selected = 
   With[{sorted = Sort[r]}, 
    Pick[sorted, deleteDuplicatesC[sorted], 1]]; // AbsoluteTiming

(* {4.033231, Null} *)

Looks like it works:

a = 1.;
b = a + $MachineEpsilon;
test = {{a, a}, {b, b}, {b, a}, {a, b}};
Pick[test, deleteDuplicatesC[test], 1]

(* {{1., 1.}} *)

Answer 4

Tally also obeys Internal`$SameQTolerance:

a = RandomReal[{-1, 1}, 100000];

Tally[a][[All, 1]] // Length

100000

Internal`$SameQTolerance = 12.0;
Tally[a][[All, 1]] // Length

33818

I have used this to handle such problems in the past, but be aware that there remains a kind of instability.

Answer 5

Here's a different approach that I would find acceptable in many situations in which numbers within a certain radius or tolerance are to be considered equal. Divide the number line up into bins of a width relative to the power of 2 of the location on the number line, and pick a representative from each bin. The boundary cases at powers of 2 can be a problem if they are important. I doubt they are in many applications. Similarly if the data contain 0., which is not a power of 2, problems might arise. This approach is basically what floating-point numbers do anyway, each number being a representative of a bin of relative radius $MachineEpsilon/2. To semi-emulate Equal, the idea would be to round off the lowest seven bits.

One way to do it is to use Round, and another way to round off an integral number of bits is to add and subtract from a number x the value 10.^Internal`$EqualTolerance * $MachineEpsilon/2 * x.

SeedRandom[0];
r = RandomReal[{1.9999999997, 2.0000000003}, {10000000, 2}]

Block[{basetol = 10.^Internal`$EqualTolerance*$MachineEpsilon/2}, 
 DeleteDuplicates@Round[r,
     Clip[   (* $MinMachineNumber is a dodge to handle extremely small numbers and 0. *)
      basetol*(2.^Floor@Log2[$MinMachineNumber + Abs@r]),
      {$MinMachineNumber, $MaxMachineNumber}
      ]
     ] // Dimensions // AbsoluteTiming
 ]
(*  {3.31892, {9936697, 2}}  *)

(* assumes tolerance is an integral number of bits *)
Block[{tol = Internal`$EqualTolerance/Log10[2] - 1 // Round, shift},
  shift = 2.^(tol + Floor@Log2[$MinMachineNumber + Abs@r]);
  DeleteDuplicates@ ((r + shift) - shift)
  ] // Dimensions // AbsoluteTiming
(*  {3.1248, {9936776, 2}}  *)

Compared to a straight DeleteDuplicates:

DeleteDuplicates[r] // Dimensions // AbsoluteTiming
(*  {2.54819, {9999992, 2}}  *)

Well, the two methods are not too bad compared to the presumptive optimum of DeleteDuplicates.

The difference in outcomes of the two methods is explained by rounding error in second argument to Round:

(* assumes tolerance is an integral number of bits *)
Block[{basetol = Log2[10.^Internal`$EqualTolerance*$MachineEpsilon] - 1 // Round}, 
 DeleteDuplicates@Round[r,
     Clip[
      2.^(basetol + Floor@Log2[$MinMachineNumber + Abs@r]),
      {$MinMachineNumber, $MaxMachineNumber}
      ]
     ] // Dimensions // AbsoluteTiming
 ]
(*  {3.33764, {9936776, 2}}  *)

It is slightly faster not to worry about the power of 2 precisely in figuring out the shift, with somewhat different rounding.

(* assumes tolerance is an integral number of bits *)
Block[{basetol = 10.^Internal`$EqualTolerance/2 // Round, shift},
  shift = basetol*r;
  s3 = DeleteDuplicates@((r + shift) - shift);
  s3 // Dimensions
  ] // AbsoluteTiming
(*  {2.77393, {9937043, 2}}  *)

Answer 6

It seems to me that all of the other answers produce incorrect results. Consider:

data = {{1., 1.}, {1.+$MachineEpsilon,-1.}, {1.+$MachineEpsilon,1.}};

Clearly the first and last points should be considered duplicates, yet none of the other answers remove this duplicate. For example:

SetOptions[EvaluationNotebook[], PrintPrecision->20];

compiledUnion[data]
Split[Sort[data]][[All, 1]]
Pick[data, deleteDuplicatesC[data], 1]
Tally[data][[All, 1]]

{{1., 1.}, {1.0000000000000002, -1.}, {1.0000000000000002, 1.}}

{{1., 1.}, {1.0000000000000002, -1.}, {1.0000000000000002, 1.}}

{{1., 1.}, {1.0000000000000002, -1.}, {1.0000000000000002, 1.}}

{{1., 1.}, {1.0000000000000002, -1.}, {1.0000000000000002, 1.}}

The other answers don't work because they basically sort and then remove adjacent duplicates. In the example above, the duplicates are not adjacent after sorting. The only suggested answer that works for this example is:

Union[data, SameTest->SameQ]

{{1., 1.}, {1.0000000000000002, -1.}}

However, as mentioned by the OP, this method has $O(n^2)$ complexity, and moreover it uses a relative tolerance instead of an absolute tolerance. For example:

data2 = {{1., 1.}, {1.*^-10, 1.*^-10}, {1.*^-10+$MachineEpsilon, 1.*^-10}};

One might wish that the second and third points are considered duplicates, yet:

Union[data2, SameTest->SameQ]

{{1.*10^-10, 1.*10^-10}, {1.0000022204460493*10^-10, 1.*10^-10}, {1., 1.}}

This is because the absolute difference is $MachineEpsilon, while the relative difference is only $\approx 10^{-6}$, and SameQ is based on relative difference.

Here is an approach that uses absolute differences, and is orders of magnitude faster than Union (with the SameTest option) for large datasets:

deleteDuplicatePoints[data_, tol_] := With[{n = Nearest[data->"Index", data, {All, tol}]},
    data[[DeleteDuplicates[Min/@n]]]
]

There are two important differences between deleteDuplicatePoints and Union. First, deleteDuplicatePoints uses an absolute tolerance, while Union uses a relative tolerance. Second, deleteDuplicatePoints uses euclidean distance, while Union uses chessboard distance (see ChessboardDistance).

Let's check deleteDuplicatePoints on our two previous examples:

deleteDuplicatePoints[data, $MachineEpsilon]
deleteDuplicatePoints[data2, $MachineEpsilon]

{{1., 1.}, {1.0000000000000002, -1.}}

{{1., 1.}, {1.*10^-10, 1.*10^-10}}

In both cases, the duplicate is removed. What about larger datasets? Here is a comparison between deleteDuplicatePoints and Union with a SameTest option for largish dataset with lots of duplicates:

data = RandomReal[{1, 1 + 100 $MachineEpsilon}, {10^5, 2}];
r1 = deleteDuplicatePoints[data, $MachineEpsilon];//AbsoluteTiming
r2 = Union[data, SameTest->SameQ]; //AbsoluteTiming

{0.053882, Null}

{1.065719, Null}

And a largish datasets with few duplicates:

data = RandomReal[{1, 1 + 10^5 $MachineEpsilon}, {10^4, 2}];
r1 = deleteDuplicatePoints[data, $MachineEpsilon];//AbsoluteTiming
r2 = Union[data, SameTest->SameQ]; //AbsoluteTiming

{0.005078, Null}

{10.17658, Null}

When there are lots of duplicates, then Union with a SameTest option becomes competitive on timing because many results get pruned, and so less comparisons are needed. When there aren't a lot of duplicates, then the $O\left(n^2\right)$ behavior becomes clear. My tests indicate that the Nearest approach appears to be $O(n \ln(n))$.

Finally, a speed comparison with the other answers. I will just compare deleteDuplicatePoints with Union (without the SameTest option):

data = RandomReal[{1, 1 + 10^7 $MachineEpsilon}, {10^7, 2}];
r1 = deleteDuplicatePoints[data, $MachineEpsilon];//AbsoluteTiming
r2 = Union[data]; //AbsoluteTiming

{10.612259, Null}

{3.048912, Null}

So, not too much slower than Union.