# Parallelizing or optimizing a script for calculating the size of all pairwise intersections between collections of sets

I need to speed up or parallelize an operation that looks essentially like this:

comparisonSetsA = Round[Table[RandomReal[1., {100, 2}], {k, 1, 100}], 0.05];
comparisonSetsB = Round[Table[RandomReal[1., {100, 2}], {k, 1, 100}], 0.05];

intersectionSizes = Table[Length[Intersection[comparisonSetsA[[a]], comparisonSetsB[[b]]]], {a, 1, Length[comparisonSetsA]}, {b, 1, Length[comparisonSetsB]}] // AbsoluteTiming


Here, we're computing the length of intersections between all pairs of sets in comparisonSetsA and comparisonSetsB (containing real numbers), which, while less trivial in my actual application, are roughly about the same size as above. Simply writing ParallelTable makes things about 3.56 times worse in terms of timings on an 8-core machine.

Is there any way to efficiently parallelize or optimize the above process? Note that I cannot make assumptions about the size of any of the intersections, and the sets can be assumed to not contain duplicates.

Update - If it helps, I only need to know if the size of the intersection is greater than a smaller threshold value, e.g. $||A \cap B|| \geq 2$ or $3$. I don't need to know anything more than this.

Update 2 - I'd like to extend ybeltukov's to handle cases where the sets in A and B are not necessarily of equal size. Is there a simple way to do this that doesn't add much in the way of overhead? I'd asked a question about a possible "padding" method here: http://mathematica.stackexchange.com/questions/32770/equalizing-the-sizes-of-sets-in-an-array-using-non-intersecting-integer-elements

Update 3 - Let me be a bit more explicit about Update 2, since I did a poor job earlier communicating. I meant that I would like to extend ybeltukov's very nice result to the cases where internal to sets A and B, the sets are of different sizes. Say, for example, the case where list A has sets of size 5, 6, and 7, and set B has sets of size 8, 9, and 4.

I believe one calls these "ragged" lists.

-
(1) The sets you generate contain pairs rather than individual values, and they are reals. If this is representative of your actual problem then the intersections will all be empty. Unless you intend to use an equality test that allows for substantial granularity. – Daniel Lichtblau Sep 19 '13 at 19:10
(2) Do you know anything about the total size of intersection between A and B? If it is small compared to the total size of either collection then you can save time by processing those elements rather than sublist pairts (I can provide details, but the method will perform relatively poorly unless this size assumption holds). – Daniel Lichtblau Sep 19 '13 at 19:13
@DanielLichtblau Thanks for looking at my question - yes, I plan to coarse-grain the real-numbered values as per my previous question: mathematica.stackexchange.com/questions/32619/… – RM1618 Sep 19 '13 at 21:04
@DanielLichtblau Perhaps 50% of the intersections will be empty post coarse-graining (round to the nearest first decimal place with my actual data set) and the remainder will have something like a 20% overlap. – RM1618 Sep 19 '13 at 21:06

Mapping this problem to the matrix multiplication produces 250x speedup!

Let we have sets

$$a = \{2,3,5,7\},\\ b = \{3,5,5,9\}.$$

Intersection $a\cap b = \{3,5\}$ contain $2$ elements.

Let ${a_{\rm bin}}_n=1$ if $n\in a$ and $0$ otherwise ($b_{\rm bin}$ is similar)

$$a_{\rm bin} = (0,1,1,0,1,0,1,0,0),\\ b_{\rm bin} = (0,0,1,0,1,0,0,0,1).$$

Scalar product $a_\rm{bin}\cdot b_\rm{bin}=2$ which corresponds to $2$ elements in the intersection.

Let now we have sets of sets

$$A = \{\{2,3,5,7\},\{4,5,6,7\}\},\\ B = \{\{3,5,5,9\},\{4,7,8,9\}\}.$$

This corresponds to the matriсes

$$A_{\rm bin} = \left( \begin{array}{ccccccccc} 0 & 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 \end{array} \right)^T,\\ B_{\rm bin} = \left( \begin{array}{ccccccccc} 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 1 \end{array} \right)^T.$$ The lengths of pairwise intersections are the elements of the matrix product $$A_{\rm bin}^TB_{\rm bin} = \left( \begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right).$$

Original solution:

A = RandomReal[1., {500, 100, 2}];
B = RandomReal[1., {500, 100, 2}];
comparisonSetsA = Round[A, 0.05];
comparisonSetsB = Round[B, 0.05];

int = Table[
Length[Intersection[comparisonSetsA[[a]],
comparisonSetsB[[b]]]], {a, 1, Length[comparisonSetsA]}, {b, 1,
Length[comparisonSetsB]}]; // AbsoluteTiming


7.749266

Proposed solution:

1. Mapping to natural indexes

n = 20;
idA = 1 + Round[n A].{1, n + 1};
idB = 1 + Round[n B].{1, n + 1};

2. Searching nonempty bins

bin = Transpose@{#, ConstantArray[#2, Length[#]]} &[Tally[#][[All, 1]], #2[[1]]] &;
binA = Join @@ MapIndexed[bin, idA];
binB = Join @@ MapIndexed[bin, idB];

3. Converting to sparse arrays

spA = SparseArray[binA -> ConstantArray[1, Length[binA]], {(n + 1)^2, Length[A]}];
spB = SparseArray[binB -> ConstantArray[1, Length[binB]], {(n + 1)^2, Length[B]}];

4. Matrix multiplication

int2 = Normal[Transpose[spA].spB]; // AbsoluteTiming


0.028975

5. Result checking

Max@Abs[int2 - int]


0

P.P.S You can use any mapping to natural integers. Rounding is the simplest possibility but not unique. The maximum value of indexes (max in Daniel Lichtblau answer) is irrelevant because I use sparse arrays.

-
Beautiful method, and quite a surprising speedup. – RM1618 Sep 20 '13 at 22:22
How might I handle the case where the sets are not of equal size? I'm getting nonrectangular tensor errors, which makes sense. – RM1618 Sep 21 '13 at 6:48
I mean the sets in A and B, respectively. Can I "pad" them to be of equal size while avoiding intersections? – RM1618 Sep 21 '13 at 6:56
@user9603 I thought that my solution is general. I will test it in several hours. Currently I'm not near the computer. – ybeltukov Sep 21 '13 at 10:47
I think the trouble is with the dot product during the mapping? – RM1618 Sep 21 '13 at 10:53

Here is a method based on iterating over all values that appear in both the A and B sets. I avoided the rounding stuff by working with integers, but you might be able to apply something similar to your case. I'll illustrate with sets of single ints (not pairs), with 1000 A and B sets each containing 1000 elements (possibly with repeats though one can easily remove that). I chose to have all lie between 0 and 10^5 so as to emulate the case where intersection sizes will be small.

Set it up:

n1 = 1000;
n2 = 1000;
max = 100000;
SeedRandom[11111]

comparisonSetsA = RandomInteger[max, {n1, n2}];
comparisonSetsB = RandomInteger[max, {n1, n2}];


Find intersection of all elements in the A sets with all in the B sets.

t1 = Timing[
setA = Union[Flatten[comparisonSetsA]];
setB = Union[Flatten[comparisonSetsB]];
setAB = Intersection[setA, setB];
Scan[(keep[#] = True) &, setAB];
]

(* Out[313]= {0.610000, Null} *)


Create "hash" sparse sets that record, for each value in setAB, what sets in comparisonSetsA and comparisonSetsB contain it. By "hash" I mean that we use Mathematica DownValues to record these (don't fret if this description makes no sense-- the code can still be understood without recourse to the verbiage).

Clear[h1, h2];
h1[_] := {}
h2[_] := {}

t2 = Timing[
Do[Scan[If[TrueQ[keep[#]], h1[#] = {j, h1[#]}] &,
comparisonSetsA[[j]]], {j, Length[comparisonSetsA]}];
Do[Scan[If[TrueQ[keep[#]], h2[#] = {j, h2[#]}] &,
comparisonSetsB[[j]]], {j, Length[comparisonSetsB]}];
]

(* Out[317]= {10.830000, Null} *)

t3 = Timing[
Scan[(h1[#] = Union[Flatten[h1[#]]]) &, setAB];
Scan[(h2[#] = Union[Flatten[h2[#]]]) &, setAB];
]

(* Out[318]= {1.680000, Null} *)

Now set up the intersection table. Initially everything is an empty set. We iterate over the values in setAB, find all A and B subsets each lies in, and use those to fill in the intersection table.

t4 = Timing[
itable =
ConstantArray[{}, {Length[comparisonSetsA],
Length[comparisonSetsB]}];
Do[aindices = h1[val]; bindices = h2[val];
Do[itable[[i, j]] = {itable[[i, j]], val}, {i, aindices}, {j,
bindices}],
{val, setAB}];
itable = Map[Flatten, itable, {2}];
]

(* {17.610000, Null} *)


Compare to the pairwise subset intersection approach.

Timing[
itable2 =
Table[Intersection[comparisonSetsA[[a]], comparisonSetsB[[b]]], {a,
1, Length[comparisonSetsA]}, {b, 1, Length[comparisonSetsB]}];]

(* Out[333]= {155.630000, Null} *)

itable === itable2

(* Out[334]= True *)


I do not know what the method by @ybeltukov might do with this example, as I did not want to sort through the various differences in input in order to check (read: I do not understand that code and don't wish to admit it). Looks like it could be significantly faster still.

Also what I show is tailored for smaallish pairwise intersections and not likely to be useful in cases where the pairwise intersection sets tend to be large. For example, if I change max to to^6 then the main loop drops from 17.6 to 4.5 seconds (and with max at 10^7 that loop is no longer the bottleneck), whereas if I instead drop max to 10^4, that time goes up to 139.2 seconds, which is no gain at all over the brute force method.

-
I hope that my approach is now clearer, see update. – ybeltukov Sep 20 '13 at 16:47
@ybeltukov Thanks for the update. Very nice method. It is in a way similar to what I did (I had already realized that much but didn't quite see how the matrix product fit in). The surprise, for me, is that since only a length is wanted, a (sparse) matrix product will give exactly that, and with considerable speed. This is in part because I lost track early on of the fact that only a length is needed... – Daniel Lichtblau Sep 20 '13 at 18:29
...which raises the following possibility. Sonce we really do not need ragged arrays in the intersection table, it should possible to gain speed by using Compile. I think your matrix product would still be faster though. – Daniel Lichtblau Sep 20 '13 at 18:31