Top k closest pairs in a set of million 128-dimensional points I have a set of 1 million points in 128-dimensional space. Among all trillion pairs of the points in the set, I need to get a subset of 100 million pairs whose cosine distances are less than that of every pair outside the subset.

I tried measuring cosine distances of all trillion pairs, sort them and get first 100 million of sorted pairs. But this program is estimated to take few years in single thread.

from scipy.spatial.distance import cosine

points = load_points_list()         # List of points, where each point is tuple of 128 floats
M = len(points)                     # 1 million points

cosine_distances, pairs = [], []
for i in range(M):
    for j in range(M):
        cosine_distances.append(cosine(points[i], points[j]))
        pairs.append((i, j))

# Sort pairs based on cosine distances
cosine_distances, pairs = (list(s) for s in zip(*sorted(zip(cosine_similarities, pairs_indices))))

top_100_million_closest_pairs = pairs[:100000000]


Is there a more efficient algorithm or GPU-enabled acceleration for a problem at this scale?
First, cut your work load in half.  You only need to consider pairs where i &lt; j, so:

for i in range(M-1):
    for j in range(i+1, M):
        # ...


If you actually need both (i, j) and (j, i), then you only need the 50 million closest pairs, and can duplicate &amp; swap every pair in the list as the last step.   More savings!

Record &amp; sort the first 100 million (50 million?) pairs.  Then, as new pairs are generated, if closer than the worst stored pair, use bisect.bisect() to find insertion spot, add the new pair, and kick out the worst pair.   

Time: \$O(n^2 log k)\$, 
Space: \$O(k)\$, n is number of points, k is desired number of pairs.