feat: Agglomerative clustering

[Mec-iS]

This is a draft PR to collect feedback on the preliminary research phase about implementing Agglomerative clustering via Muellner's "Modern hierarchical, agglomerative clustering algorithms", 2011.

The paper is interesting as presents an historical perspective for sequential, agglomerative, hierarchic, nonoverlapping methods (SAHN). As you may read in the documentation string in the module, different algorithms are presented in different versions; they are compared to what is called PRIMITIVE_CLUSTERING and formally proofed.

As the objective is to have MST-linkage (fastcluster) implemented, we need to go through some preliminary work with older algorithms as the objective implementation requires multiple procedures and data structures to be implemented:

1. LABEL (aka union-find), as in Fig.5 in Müllner, 2011
2. and a doubly-linked list to store the resulting dendrogram

These can be more easily implemented using Rust compared to C++; the latter requires custom implementations (see class doubly_linked_list and class auto_array_ptr that provide safer wrapper over C++ std library, those should be implemented in a more straightforward way using Rust stdlib); any suggestion or insight about how to proceed with these implementations are welcome.

There is a native-C++ implementation here that may provide a shortcut compared to the original Muellner code for fastcluster that implies R/Python bindings.

As a roadmap, I would suggest to start implementing PRIMITIVE_CLUSTERING as knowledge-building to tackle fastcluster that requires quite a lot of reasoning over the mentioned procedures of post-processing and data structures: (1.) is not quite clear to me yet and for (2.) we need to find the right Rust implementation to do it (I considered using petgraph as it looks like the standard way to handling graph-like structures). Once we have these latest things, we can move on to implement MST-linkage.

Please correct anything wrong found in the above construct.

[VolodymyrOrlov]

This is a very solid analysis!

I do not think that PRIMITIVE_CLUSTERING implementation will be of any use, since it will probably be very slow. I suggest to start right away from MST_linkage_core and NN_chain_core since in practice single, ward, complete and average are the most commonly used linkages. If N is not too big we can skip algorithm described in Figure 5 (A union-find data structure suited for the output conversion) (for now), since this is an optimization that helps to reduce space required to store node labels.

For doubly linked list can we use https://doc.rust-lang.org/std/collections/struct.LinkedList.html?

[VolodymyrOrlov]

one more suggestion. What if we remove parameters k and distance_threshold from method fit? I know Scikit-learn's AgglomerativeClustering.fit has both of these parameters but it is a bad design. Instead, we need a new method get_labels(k) and get_labels(distance_threshold) that will generate cluster labels. We can call this method anything, e.g. cluster_labels or split... The idea is that user fits algorithm to our data first with fit that returns Self. Then he/she can examine dendrogram and decide where to make a split or just take whatever hardcoded number of cluster he/she needs. What do you think?

[Mec-iS]

I will proceed implementing FastPair or MST-Linkage depending how difficult it is to implement the latest, or hopefully both. For the latest, I will use std::collections::LinkedList.

[Mec-iS]

@VolodymyrOrlov Please help about how to implement cmp for PairwiseDissimilarity<RealNumber>. cmp is implemented for RealNumber?

[Mec-iS]

FastPair implemented. Moving to integrating this first algorithm into agglomerative, adding the post-processing steps. Cleanup at the end.

[Mec-iS]

Please provide feedback on interface

[VolodymyrOrlov]

Please provide feedback on interface

Ideally, end user interface should be as simple as possible. It would be great to have a simple meta-class that represents agglomerative clustering algorithm, similar to AgglomerativeClustering. Particular implementation should be hidden behind a simple enum value and initialized within this meta-class, e.g. AgglomerativeClusteringAlgorithm with a value for FastPair. Linkage method should also be hidden behind enum (or interface with out-of-the-box implementations).

AgglomerativeClustering should have a static method fit that takes X, algorithm, a linkage method and returns an instance of itself. Then we also need a couple of non-static methods that will produce labels from number of classes (k) and from distance threshold (d). Of course, we must also provide a method that returns matrix with dendrogram or other way to select d

I think the rest of the classes should not be public or should be hidden behind pub (crate). This way we reserve a right to make any changes to these classes in the future

[Mec-iS]

Ideally, end user interface should be as simple as possible. It would be great to have a simple meta-class that represents agglomerative clustering algorithm, similar to AgglomerativeClustering. Particular implementation should be hidden behind a simple enum value and initialized within this meta-class, e.g. AgglomerativeClusteringAlgorithm with a value for FastPair. Linkage method should also be hidden behind enum (or interface with out-of-the-box implementations).

MST_Linkage works only with simple linkages (HCLUST_METHOD_SINGLE, HCLUST_METHOD_COMPLETE, HCLUST_METHOD_AVERAGE, HCLUST_METHOD_MEDIAN), so for now I would just avoid the enum for the linkage function.

I will proceed with a complete example so then we can iterate.

[Mec-iS]

dendrogram computation algorithm implemented.
Still missing the two post-processing steps (FindUnion and label) for fastcluster.

[VolodymyrOrlov]

@VolodymyrOrlov do we need to implement this step? if yes, what is the most straighforward way to do:

        // order = np.argsort(Z_arr[:, 2], kind='mergesort')
        // Z_arr = Z_arr[order]

?

Also, please provide some hints to implement threshold-based aggregation (currently it is implemented only with number of expected clusters).

Yes, I don't think dendrogram should be sorted

The best way to understand how threshold-based mechanism of splitting data into clusters works is to look at this example. The threshold is a parameter that is chosen by a user. It is basically a horizontal line in the dendrogram above which clusters are not merged. Look at the picture in the example I've share. If threshold is set by a user to 20 we should split data into 2 clusters. When the threshold is set to 10 data should be split into 4 clusters.

[Mec-iS]

After passing code reviewer's screening, I would prefer to merge this as the code is providing fastpair+fastcluster and returnining the dendrogram. We can move the implementation of the rest of the interface to the next PR.

PS. There are some clippy warnings from other modules that we should figure out how to fix or silence

[Mec-iS]

@VolodymyrOrlov

When the threshold is set to 10 data should be split into 4 clusters.

shouldn't this to be 3 clusters? the parent cluster is split into two children, the parent one still count?

[Mec-iS]

Ready for merging the algorithms. Clustering is performed with FastPair + fastcluster.

I will add the remaining high-level methods (fit_k and fit_threshold) soon later this month.

[VolodymyrOrlov]

Great work, Lorenzo! I think this PR is very close to finishing line but there are some changes in the interface that I would like to see before we can merge it.
I've also attached a notebook that demonstrates capabilities of Scikit's AgglomerativeClustering that we have to match if we want our algorithm to be useful to end users. Let me know if you need help refactoring your code.

https://gist.github.com/VolodymyrOrlov/84548cef1a0fefd83ff9e57f0de5e03e

[VolodymyrOrlov]

Can you put fastpair behind pub(crate) for now?

Suggested change

	pub mod fastpair;
	pub(crate) mod fastpair;

[Mec-iS]

I have tried but then the code in agglomerative.rs cannot reach it.
If you can provide the code for this so I can see how it works.

[VolodymyrOrlov]

can you move method fit and labels directly into AggregativeFastPair and remove this train for now? We do need a broader system of interface for all our algorithms but I think it would be best to have it in a separate PR.

[Mec-iS]

I can do this

[Mec-iS]

Just a note to clear: FastPair do not perform clustering, it is only a faster implementation to find the closest neighbour

[VolodymyrOrlov]

we need 2 methods that calculate cluster labels. One method will calculate labels for fixed number of clusters, k, while the other method produce labels from a threshold parameter. Take a look at the Jupyter notebook with demonstration that I've attached to this PR.

I am not sure what would be the best way to call these methods, labels_by_k, labels_by_threshold?

[Mec-iS]

I can do it but I need the pseudo-code for the methods mentioned, as I have no idea how to compute the outcome from the dendrogram. The current implementation compute the dendrogram for the dense matrix in the shape (for the Iris dataset):

       let expected = [
            //  index, index, distance, label
            [0.0, 4.0, 0.01999999999999995, 2.0],
            [7.0, 19.0, 0.029999999999999964, 3.0],
            [1.0, 9.0, 0.030000000000000037, 2.0],
            [14.0, 18.0, 0.05999999999999993, 2.0],
            [2.0, 3.0, 0.0600000000000001, 2.0],
            [6.0, 23.0, 0.06999999999999998, 3.0],
            [10.0, 12.0, 0.07000000000000003, 2.0],
            [11.0, 16.0, 0.07000000000000013, 2.0],
            [8.0, 24.0, 0.08999999999999982, 4.0],
            [21.0, 27.0, 0.09000000000000011, 6.0],
            [20.0, 28.0, 0.10999999999999982, 9.0],
            [22.0, 26.0, 0.1799999999999998, 4.0],
            [25.0, 30.0, 0.2600000000000009, 6.0],
            [13.0, 29.0, 0.27000000000000013, 10.0],
            [5.0, 32.0, 0.3800000000000002, 11.0],
            [15.0, 33.0, 0.54, 12.0],
            [31.0, 34.0, 0.6899999999999995, 18.0],
            [17.0, 35.0, 0.7000000000000001, 19.0],
            [5.0, 17.0, 4.460000000000001, 0.0],
        ];

From this the methods should compute the "by cluster number" and "by threshold", can you provide the procedures to do so or provide a link to code?

[VolodymyrOrlov]

@Mec-iS Take a look at how it is done in Scikit's _agglomerative.py. Start from _hc_cut method. If I understand this method correctly, the basic idea is to take k records from the dendrogram with the smallest index values - these are your top nodes in k clusters. All you need to do next is to iterate through descendants of each top node and assign all nodes at the lower level label of its parent.
Graphically it will look as if you are making a horizontal cut at some level that covers exactly k trees. For example, on this picture

the horizontal cut will produce k=3 clusters. Now you need to follow descendants of these three nodes and assign them labels of their parent.

Once you are done with "by cluster number" you can easily implement method "by threshold" using this logic

[VolodymyrOrlov]

Do you plan to add different types of linkage (“ward”, “complete”, “average”, “single”) later? Even if you plan to support only "ward" at this moment we still should have an enum parameter linkage with a single item "ward"

[Mec-iS]

FastPair works only and only with single. To add different linkages we need also to add an alternative algorithm like NearestNeighbour that can work with all the linkages.

PS. on another side: what about using different distances functions to compute the pairwise distances?

[VolodymyrOrlov]

I feel that it should be possible. Take a look at this paper:

"We can perform median, centroid, Ward’s, or other
bottom-up clustering methods in which clusters are represented by objects, in time
O(n^2 log^2n) and space O(n)."

[Mec-iS]

I don't know if I understood wrong but I think that paper refers to fastcluster that is the clustering algorithm; FastPair is a closest neighbour algorithm. The correspective of FastPair in that paper is the Nearest Neighbour algorithm that can be used with other linkages. FastPair instead can be used only with single as, if I remember correctly, stated in the FastPair paper.

[Mec-iS]

moved on with #142