# Building a large distance matrix

I am trying to build a distance matrix for around 600,000 locations for which I have the latitudes and longitudes. I want to use this distance matrix for agglomerative clustering. Since this is a large set of locations, calculating the distance matrix is an extremely heavy operation. Is there any way to opimize this process while keeping in mind that I am going to use this matrix for clustering later. Below is the code I am using.

from scipy.spatial.distance import pdist
import time
start = time.time()
# dist is a custom distance function that I wrote
y = pdist(locations[['Latitude', 'Longitude']].values, metric=dist)
end = time.time()
print(end - start)


The simplest solution for such a task coming to my mind is to do a simple kmeans clustering (or batch variants) using the exact same metrics as planned for the later hierarchical clustering step (in your case eucledian / minkowski with p=2). For the initial kmeans step you chose the number of clusters k such that a distance computations on those cluster centers will be feasable. With this you basically initialize the hierarchical clustering one level down the hierarchy.

import numpy as np
from sklearn.cluster import KMeans
from scipy.spatial.distance import pdist, squareform
from sklearn.neighbors import NearestNeighbors
x = np.random.rand(600000*2).reshape((600000, -1))
kmeans = KMeans(k=int(len(x)/100), n_jobs=16)
kmeans.fit(x)
x_ = kmeans.cluster_centers_
D = squareform(pdist(x_))


Here x is your data (simulated here via random coordinates in a square) and the reduction factor for the first step is 100. Still this operation is quite expensive in terms of computing time. A faster solution follows:

Note that in case your data is quite uniform consider a random preselection of a subset of your data (maybe with a distance criteria similar to what happens in Poisson disc sampling) instead of the kmeans. That would be super fast:

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import pairwise_distances
# input data
x = np.random.rand(600000*2).reshape((600000, -1))

# reduction of the size of the set of samples using uniformity
# chose a couple of times randomly from the input data and compute clusters in that subset
# this allows to avoid expensive clustering on the huge input data set
repetitions = 20
reduction = 1000
x_ = [x[np.random.choice(len(x), size=int(len(x)/reduction), replace=False)] for i in range(repetitions)]
x_ = np.reshape(x_, (repetitions*int(len(x)/reduction),2))
kmeans = KMeans(int(len(x)/reduction), n_jobs = 16)
kmeans.fit(x_)
# Compute the distance matrix to start with for the hierarchical clustering
D = pairwise_distances(kmeans.cluster_centers_, n_jobs=16)
# Compute the assignments from each data point to any of the level-0 clusters
D_data2clusters = pairwise_distances(x, kmeans.cluster_centers_, n_jobs=16)
cluster_labels = np.argmin(D_data2clusters,1)

from matplotlib import pyplot as plt
for cluster in np.unique(cluster_labels):
plt.scatter(*x[cluster_labels == cluster].T)
plt.scatter(*kmeans.cluster_centers_.T)
plt.savefig("clustertest")


Here kmeans.cluster_centers_ (green dots) or distance matrix D could be the input for your hierarchical clustering.

from scipy.cluster.hierarchy import linkage, dendrogram
plt.figure()
dendrogram(z)
plt.savefig("dendrogram")


Another possible workaround in such scenarios can be to compute an incomplete distance matrix using neighbor tree based approaches. This would basically be your approximation of the distance matrix.

To this end you first fit the sklearn.neighbors.NearestNeighbors tree to your data and then compute the graph with the mode "distances" (which is a sparse distance matrix). You will need to push the non-diagonal zero values to a high distance (or infinity). In that sparse matrix basically only the information about the closer neighborhood of each data is stored and larger distances are not even computed and put into that matrix. However, for your scenario memory for a float matrix of size 600000^2 would have to be allocated - that is 2.62 TiB which is unrealistic.

Aussuming you have such a distance matrix you can try and play around whether any hierarchical clustering approach handles that sort of incomplete distance matrix appropriately but as was pointed out by the earlier answer it will be very expensive. I therefore suggest in such cases to use the very efficient kmeans (on the data itself as shown above - also batch versions might be useful in that case) or kmedoid (on the sparse distance matrix) which you can construct and apply in a hierarchical fashion as well.

Generally if you can reformulate an algorithm such that it does not require access to the full distance matrix at once you may use sklearn's pairwise_distances_chunked.

Have you considered that the following steps will be even worse?

The standard algorithm for hierarchical clustering scales O(n³). You just don't want to use it on large data. It does not scale.

Why don't you do a simple experiment yourself: measure the time to compute the distances (and do the clustering) for n=1000,2000,4000,8000,16000,32000 and then estimate how long it will take you to process the entire data set assuming that you had enough memory... You will see that it is not feasible to use this algorithm on such big data.

You should rather reconsider your approach...

• You just don't want to use it on large data. But how do I do it on large data sets. There must some way. Can you suggest any other approach? I am a bit new to this. Commented Mar 21, 2019 at 7:24
• Karthik: compute how much memory you would need. You can try to implement this yourself by storing the distances on disk, as this will be several TB. But trust me, do the experiment I suggested first and estimate how long it will take you. Are you willing to wait weeks for the result? Commented Mar 21, 2019 at 16:14
• If the experiment shows your runtime increases by 4 with each doubling the size, going from 32k to 600k means you'll need about 350x as long. With the expected O(n³) increase, it will take 6600x as long. Then you can estimate if it's worth trying. Maybe add a factor of 10x additionally for working on disk instead of in-memory. You'll need about 1.341 TB disk space to store the matrix, and as much working space. That is doable. You'll need to read this matrix many many times though, so even with a SSD this will take several days just for the IO. Commented Mar 21, 2019 at 16:31