# Bug or undocumented behaviour in Intersection

I have two lists consisting of element pairs:

a1 = Table[{i, i^2}, {i, 10}] // N;
a2 = Table[{2 i, Sqrt[i]}, {i, 5}] // N;


I want to pull out from both lists those pairs which have the same first element. Using Intersection with the SameTest option:

In:= d1 = Intersection[a1, a2, SameTest -> (First[#1] === First[#2] &)]

Out= {{2., 4.}, {4., 16.}, {6., 36.}, {8., 64.}, {10., 100.}}


gives a different result as compared to

In:= d2 = Intersection[a2, a1, SameTest -> (First[#1] === First[#2] &)]

Out= {{2., 1.}, {4., 1.41421}, {6., 1.73205}, {8., 2.}, {10., 2.23607}}


d1 contains the element pairs of the list a1 for which there exists also a similar first element list a2. d2 contains the element pairs of a2 for which there exists also a similar first element in list a1. Thus Intersection is not commutative.

I would like to get as "Intersection" all element pairs a1 and of a2 for which there is a match in the first element. I can of course call Intersection as shown above two times and then join the results d1 and d2. But is there not a more efficient way so that I need to call Intersection only once? I ask because in my real case the lists are very large and the SameTest is more complicated than the identity test used in this example.

• Technically cannot be done with 2 Intersections, proper approach is to gather up equivalence classes in a Association. In a meeting right now, cannot code right now – Manuel --Moe-- G Aug 20 '19 at 15:58
• Have you done benchmarking to see how slow the other operation actually is on a "real" dataset? These operations are fundamental enough where they should be pretty quick (outside of the fundamental algorithmic complexity) – user6014 Aug 20 '19 at 16:00
• @hippo3773 What, exactly, is the expected output in this example? something like {{2., 4., 1.}, {4., 16. 1.41}, ..., {10., 100., 2.23}}? – AccidentalFourierTransform Aug 20 '19 at 17:05
• @AccidentalFourierTransform: Ideally, I'd prefer the output to be just a list with two sublists, each containing the sorted element pairs which have corresponding first elements. In practice the original lists consist not just of element pairs, but are matrices, with possibly not the same row dimensions. Also the SameTest function is more complicated than Identity: e.g. Abs[#1-#2]<epsilon . – hippo3773 Aug 21 '19 at 17:33
• @hippo3773 Can you describe that giving the list explicitly instead of with words? What is the exact list you expect? As in, what numbers go inside {...}? – AccidentalFourierTransform Aug 21 '19 at 17:46

Maybe this produces does what you want (it is basically what Manuel --Moe-- G suggested).

f[a1_?MatrixQ, a2_?MatrixQ] := Thread /@ List @@@ Normal[
Select[
Merge[{
},
Identity
],
Length[#] >= 2 &
]
]

f[a1, a2]


{{{2., 4.}, {2., 1.}}, {{4., 16.}, {4., 1.41421}}, {{6., 36.}, {6., 1.73205}}, {{8., 64.}, {8., 2.}}, {{10., 100.}, {10., 2.23607}}}

# Alternate Method

One can obtain a more efficient implementation (without binning) with Nearest:

g[a1_?MatrixQ, a2_?MatrixQ] :=
Module[{idx, rowpointers, columnindices, A, i, j},
idx = Nearest[a2[[All, 1]] -> "Index", a1[[All, 1]], {\[Infinity], 0}];
rowpointers = Join[{0}, Accumulate[Length /@ idx]];
columnindices = Partition[Join @@ idx, 1];
A = SparseArray @@ {Automatic,
{Length[a1], Length[a2]}
,
0,
{1,
{rowpointers, columnindices},
ConstantArray[1, rowpointers[[-1]]]
}
};
{i, j} = Transpose[A["NonzeroPositions"]];
Union[a1[[i]], a2[[j]]]
]


Now

n = 1000000;
a1 = Transpose[{Range[1., 2 n], Range[1., 2 n]^2}];
a2 = Transpose[{Range[2., 2 n, 2.], Sqrt[Range[1., n]]}];

resultf = f[a1, a2]; // AbsoluteTiming // First
resultg = g[a1, a2]; // AbsoluteTiming // First
Union @@ resultf == resultg


14.7094

1.26785

True

You can obtain essentially the same binning (up to reordering) with GatherBy, but that takes a bit of extra time:

resultg2 = GatherBy[resultg, First]; // AbsoluteTiming // First
Sort[Sort /@ resultg2] == Sort[Sort /@ resultf]


1.26529

True

# Discussion

The second method is about 10 times faster for quite a range of values of n. Since Nearest employs a space partition tree that should have complexity $$O( n \, \log(n))$$ (instead of $$O(n^2)$$ for a naive approach), this tells me that also Merge uses some clever $$O( n \, \log(n))$$-method to perform the lookup. I've once read somewhere here on Mathematica.StackExchange that "Mathematica's" Association data type is based on a special kind of trie, so that is entirely possible. However, in this special case (where the keys are members of a suitable metric space), a space partition tree seems to be faster.

# Is Intersection's current behavior a bug?

I think it isn't. By assigning a SameTest, you effectively tell Mathematica that it should not care about which representative to choose within an equivalence class. So Mathematica chooses the first(?) representative it can find. Listing all representatives of the equivalence class is just not what you asked for.

• I would have embarrassed myself because no ability to code the association work so cleanly, very good indeed. – Manuel --Moe-- G Aug 20 '19 at 19:38
• Here's some implementation details on Association: youtu.be/xvOFEve2DMw?t=735 – Chip Hurst Aug 21 '19 at 18:58

After posting the question I found out that solutions to similar problems had already been posted on this website. For my problem I constructed the following solution:

[a1_?MatrixQ, a2_?MatrixQ] := Flatten[#, 1] & /@ Select[GatherBy[Flatten[{a1, a2}, 1], First], Length[#] == 2 &]

Faster than the f-function of @Henrik, but somewhat slower than the brilliant g-function of @Henrik. Works well with matrices that have different no of columns, and it is also straightforward to implement more complicated tests for "sameness" of the first elements.