Low memory version of single linkage clustering

[lmcinnes]

Currently in scipy single linkage clustering on vector data starts by calling pdist to generate the full distance matrix, and then runs an MST calculation on that. While this is efficient it is very memory intensive for large datasets. It is relatively straightforward to to calculate the MST generating distances only as they are needed through the algorithm, allowing efficient computation without filling memory. It would be beneficial if such an option was available in scipy.

Reproducing code example:

import numpy as np
import scipy.cluster.hierarchy

data = np.random.normal(size=(64000,2))
tree = scipy.cluster.hierarchy.linkage(data, method='single')

Error message:

On my mac laptop this simply freezes the whole machine. On Linux a MemoryError with no traceback results. If the (lack of) error is not reproducible on your machine simply make data a larger array.

Scipy/Numpy/Python version information:

1.0.0 1.11.3 sys.version_info(major=3, minor=5, micro=4, releaselevel='final', serial=0)

[rgommers]

This is both a defect (no MemoryError on macOS) and an enhancement request. Discussion with @lmcinnes at the SciPy'18 sprint - changing to his single linkage implementation seems feasible. Does require a Cython interface to the distance functions in spatial though (gh-9032), so cannot easily be integrated right now.

[Mortal]

Actually, it turned out that it's rather easy to combine scipy's qhull bindings with scipy's MST code to implement single-linkage Euclidean clustering in N log N time and linear memory. I have a proof of concept in https://github.com/Mortal/singlelinkage - it easily handles 64000 points in 2d per the reproduction of @lmcinnes. It uses the edges of the Delaunay triangulation to limit the size of the input to the MST algorithm.

I would be happy to implement this into the scipy library as a replacement for the current Euclidean single-linkage implementation. Any pointers where to start?

EDIT TO ADD: Based on the comment from @lmcinnes I have added to the linked repo a quadratic-time, linear-memory algorithm that doesn't use scipy's Delaunay triangulation or MST algorithm, but instead implements Prim's algorithm directly.

[j-bowhay]

@Mortal it is first best to put your proposal to the mailing list to see if there is support for it

[lmcinnes]

This looks like a good implementation. The biggest catches are that it only works for Euclidean distance, and since it relies on Delaunay triangulation it only works well for low dimensional data (Delaunay in high dimensions becomes expensive very quickly).

Since we ultimately just want the MST, and assuming we use Prims algorithm for that, it is easy enough to only generate the distances required for each step of Prims (i.e. distances from the last added point to all remaining points) which is linear in memory and quadratic in time complexity.

If we are only interested in low(ish) dimensional data then Dual-tree Boruvka will scale better than Delaunay triangulation with dimension, and assuming we use Ball-Trees from scipy, it should support a wide range of distances beyond just Euclidean. Unfortunately Dual-tree Boruvka is rather more complicated to implement than the solution using Delaunay triangulation. Let's see what the mailing list thinks of things ...