I'm going to recommend that you use cosine similarity to solve this. The idea is that each element is turned into a cluster of one element. Then you use cosines to keep on finding the two most similar clusters, and merge them. Continue until you're happy with your groups.
Here is the basic math behind that.
Suppose that C is a cluster of elements. We can turn it into a ~10000 dimension vector V = [v0, v1, ... vn] by making each list a dimension, and just putting down the number of elements of the cluster in that list.
The "dot product" C o D of two clusters is just the sum of vi * wi. The "length" of a cluster || C || is sqrt(C o C). The "cosine" cos(C, D)between two clusters is (C o D) / || C || / || D ||. See, for example, this explanation of why. When cosine is close to 1, two clusters are very similar. When it is close to 0, they are very different.
Now suppose that we have clusters C, D and E, and decide to merge C and D into a bigger cluster. Then || CuD ||^2 = (CuD o CuD) = (C o C) + (D o D) + 2 * || C || * || D || * cos(B, C) from the law of cosines. And, finally, CuD o E = (C o E) + (D o E). This allows us to combine clusters without having to recalculate everything from scratch.
OK, enough math, now let's talk about programming.
You have lists of ~1000 elements with approximately ~500,000 pairs each. You have 10,000 of these. For approximately 5,000,000,000 data points of interest. That isn't a huge data set, but it is enough we need to not be naive about how to do it.
The technique you need is called MapReduce. A map-reduce has the following phases:
- For each input record, emit some
(key, value) pairs through a map function.
- The framework groups by
key.
- Then your reduce function gets
(key, [value1, value2, value3, ...]) and does something with it.
This works well for distributed programming because you can run as many mappers as you want in parallel, the framework can do the heavy lifting for step 2, and then you can run your reducers in parallel. And this paradigm can handle a wide variety of problems.
But for a problem of your size, you can ALSO just fake it. You "emit" pairs by writing key value lines to a file. The heavy lifting for step 2 is done by a Unix command:
LC_ALL=C sort file_in > file_out
And now your keys have been grouped together, letting you easily gather them together and reduce them.
So given a list L we simply:
def mapper (L):
for i in range(len(L)):
for j in range(i, len(L)):
emit((L[i], L[j]), 1)
And then the reducer will get:
def reducer (key, values):
elt1, elt2 = key
dot_product = len(values)
... record this somewhere ...
And now we have calculated all non-zero dot products between elements. Any missing pairs are zeros. If you have ~5000 elements, you'll get ~12,500,000 dot products. From which we can calculate ~5000 lengths and ~12,500,000 cosines.
Next programming concept, a priority queue which can be implemented with a heap. (We'll need to use some properties of a heap.) So the basic idea is as follows:
put all pairs of distinct elements in a priority queue, ordered by max cosine.
put all elements in a lookup saying they are valid clusters
for some number of times:
take the highest cosine pair
if both are still valid:
create a new cluster that is the union of the two
calculate its lengths and cosines with all other valid clusters
(The necessary math is above.)
stick the cosine pairs into the queue
take the original clusters out of the valid lookup
add the new cluster as valid
Now this has a problem. Its problem is that every time we merge two groups we correctly record what is valid and what is not. But we have thousands of pairs that are now garbage, and thousands that are added to the priority queue. The queue is rapidly growing.
We solve this by doing garbage collection every so often, such as every time it doubles in size. For that we simply go through the array underlying the heap, and make a new array of the valid ones. That new array is not a heap, but we can call the heapify operation on it to turn it back into a heap, and therefore a priority queue. And then go back to work.
And now you simply run this until you hit some sort of stopping condition. Such as having the desired number of groups, or the cosine getting too small.
And then your final list of valid groups is your answer.
