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

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)? • A single argument DeleteDuplicates seems faster than single argument Union here – Rojo Mar 26, 2013 at 2:03 • @Rojo For a list of what size with how many duplicates? Mar 26, 2013 at 5:08 • For the same test data of your answer. Union takes about 7 seconds and DeleteDuplicates slightly less than 3. So, 10000000 length and 3677844 duplicates this time – Rojo Mar 26, 2013 at 5:21 • Would you mind adding the test case in the question? For me, DeleteDuplicates is usually the fastest solution. There are some threads on mathgroup why this should be so, too... Mar 26, 2013 at 7:50 • @Szabolcs a random comment - congratulations on getting the 1337th nice question badge! :) Mar 26, 2013 at 15:28 ## 6 Answers 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 InternalBag, 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 = InternalBag[seen = First[sorted], 1]; Do[ If[i != seen, InternalStuffBag[output, seen = i, 1]], {i, sorted} ]; Partition[InternalBagPart[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 *)

DeleteDuplicates[r, Equal[##] &]; (* 294.8 days; unpacks *)

• Also including the DeleteDuplicates timing would be great. Mar 26, 2013 at 16:31
• @YvesKlett okay, I will add it later. However, it should be noted that DeleteDuplicates, while faster, does not actually produce the correct result. Mar 26, 2013 at 16:33
• Oleksandr, any info on that would also be welcome. I did a run that matched the result of Union (after being sorted). Mar 26, 2013 at 16:50
• @YvesKlett yes--contrast e.g. DeleteDuplicates[{1., 1. + $MachineEpsilon}, SameQ[##] &]. From this we can see that SameQ is special-cased internally by DeleteDuplicates, but not correctly so. Mar 27, 2013 at 7:51 • @YvesKlett I added the DeleteDuplicates timings. I couldn't be bothered to wait for 9 months for the last one to finish, so that's an extrapolation. Mar 27, 2013 at 9:12 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} *)

• Since Split only compares adjacent elements, it is not obvious to me that one iteration will always be enough. You may need something like FixedPoint[Split[#][[All,1]]&,Sort[r]]. Mar 25, 2013 at 22:00
• @LeonidShifrin If I have {a,b,c,d,e} and a == band b==c, then Split will give {{a,b,c}, ...}, so it should be fine. The problem could be that it removes a bit too much because comparison with tolerance is not transitive, so it's possible that a==b and b==c but a != c. Mar 25, 2013 at 22:02
• Actually, yes, I meant this non-transitivity, but did not make the right conclusion. Mar 25, 2013 at 22:04
• I might not get it, but for this one DeleteDuplicates is about 2.5 times faster than Union... Mar 26, 2013 at 7:55
• @Szabolcs, Split has the disadvantage that it unpacks and thus needs more memory. Why exactly can you not use DeleteDuplicates? Seems faster.
– user21
Mar 26, 2013 at 11:09

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[CompileGetElement[v, i] == CompileGetElement[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.}} *)  • Compiled Equal is a call to InternalCompareNumeric. The tolerance is $InternalEqualTolerance (7 bits) if RuntimeOptions -> "CompareWithTolerance" -> True is set and zero otherwise. It seems like it would be possible to set the tolerance explicitly either by changing Internal$EqualTolerance before compiling the code or by editing the compiled bytecode directly. Note that you must manually set "CompareWithTolerance" -> True here because the default for RuntimeOptions -> "Speed" is "CompareWithTolerance" -> False! Mar 25, 2013 at 23:28 • Thanks @OleksandrR. Once you posted your first comment I found in the docs that I couldn't specify "Speed" think this works for the example given now. Mar 25, 2013 at 23:39 • You actually can specify RuntimeOptions -> {"Speed", "CompareWithTolerance" -> True}, FWIW. Mar 25, 2013 at 23:42 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. 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}}  *)


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.

• I seem to get that deleteDuplicatePoints[r1, $MachineEpsilon] is shorter than r1, and that doesn't seem right to me. Dec 16, 2017 at 18:06 • @MichaelE2 I think you meant that Union removes more duplicates than deleteDuplicatePoints does for the largish dataset with lots of duplicates? I think this is because Union uses chessboard distance instead of Euclidean distance. By that I mean, with Union both the x and y coordinates can be off by $MachineEpsilon, while with deleteDuplicatePoints only x or y but not both can be off by $MachineEpsilon. Note that the difference between relative tolerance and absolute tolerance doesn't play much of a role because of the range I chose. Dec 16, 2017 at 18:20 • I meant iterating deleteDuplicatePoints: i.stack.imgur.com/uHlTS.png -- I guess in part it has to do with nontransitivity, but the first run seems to miss a lot? Dec 16, 2017 at 18:27 • @MichaelE2 I tweaked the Nearest call, so that more points are pruned initially. Further work is needed, though, in determining which points can be deleted. For example, in one dimension we could have neighbors sets with indices$(1, 2)$,$(1, 2, 3)$and$(2, 3)\$. We shouldn't delete points 2 and 3, since the final set {1} would not be a duplicate of 3. I'll update my answer with some thoughts on this when I get a chance. Dec 16, 2017 at 19:12