Fast hierarchical clustering algorithms in pure Crystal

This shard provides types and methods for fast hierarchical agglomerative clustering featuring efficient linkage algorithms.

The current implementation is heavily based on the work of Daniel Müllner [1] and derived from Müllner's own implementation found in the fastcluster C++ library [2], and parts of the implementation were also inspired by the kodama Rust crate (code for updating distances) and SciPy Python library (code for generating flat clusters).

The runtime performance of this library is on par with the reference implementations (see benchmark).

The most relevant types and methods for practical usage are the following:

• .cluster performs hierarchical clustering of a set of elements using the pairwise dissimilarities returned by the given block.
• .labels returns a mask of 0-based cluster indexes sorted by population.
• .centroids returns one representative element per cluster.
• .linkage performs hierarchical clustering based on a pairwise DistanceMatrix returning a Dendrogram as output.
• Rule is an enum that determines the linkage criterion or rule.
• Dendrogram represents a step-wise dendrogram that encodes the hierarchy of clusters obtained from hierarchical clustering.
• Dendrogram#flatten returns a list of flat clusters.
• Dendrogram#labels returns a mask of 0-based cluster indexes.

The available linkage rules are:

• Average or UPGMA method.
• Centroid or UPGMC method.
• Complete or farthest neighbor method.
• Median or WPGMC method.
• Single or nearest neighbor method.
• Ward's minimum variance method
• Weighted or WPGMA method.

This library is released under the MIT license.

Installation

1. Add the dependency to your shard.yml:

   dependencies:
     hclust:
       github: franciscoadasme/hclust

2. Run shards install

Usage

First define the data points to be clustered:

require "hclust"

coords = [
  [-0.30818828, 2.70462841, 1.84344886],
  [2.9666203, -1.39874721, 4.76223947],
  [3.21737027, 4.09489028, -4.60403434],
  [-3.51140292, -0.83953645, 2.31887739],
  [2.08457843, 4.24960773, -3.91378835],
  [2.88992367, -0.97659082, 0.75464131],
  [0.43808545, 3.70042294, 4.99126146],
  [-1.71676206, 4.93399583, 0.27392482],
  [1.12130963, -1.09646418, 1.45833231],
  [-3.45524705, 0.92812111, 0.15155981],
]

def euclidean(u, v)
  Math.sqrt (0...u.size).sum { |i| (u[i] - v[i])**2 }
end

The easiest way is to use .cluster. The following groups the coordinates by Euclidean distance with a cutoff of 4 using single linkage (the default):

clusters = HClust.cluster(coords, cutoff: 4) { |u, v| euclidean(u, v) }
clusters.size        # => 3
clusters.map(&.size) # => [5, 3, 2]

The block receives two elements and must return the dissimilarity between them. The result is an array of groups of the original elements.

Use a positional count to limit the number of clusters:

clusters = HClust.cluster(coords, 2) { |u, v| euclidean(u, v) }
clusters.size        # => 2
clusters.map(&.size) # => [8, 2]

.labels returns a mask of the same size as the input instead of grouped elements. Each position is the 0-based index of the corresponding cluster. Clusters are numbered by population (0 is the largest):

HClust.labels(coords, cutoff: 4) { |u, v| euclidean(u, v) }
# => [0, 1, 2, 0, 2, 1, 0, 0, 1, 0]

HClust.labels(coords, 2) { |u, v| euclidean(u, v) }
# => [0, 0, 1, 0, 1, 0, 0, 0, 0, 0]

.centroids returns one representative element per cluster (the member with the smallest average distance to the others):

centroids = HClust.centroids(coords, cutoff: 4) { |u, v| euclidean(u, v) }
centroids.size # => 3
centroids      # => [coords[0], coords[8], coords[2]]

The linkage rule can be passed as an extra argument. It changes how distances are updated when clusters are merged, so the grouping can change:

clusters = HClust.cluster(coords, cutoff: 4, rule: :centroid) { |u, v|
  euclidean(u, v)
}
clusters.size        # => 5
clusters.map(&.size) # => [2, 3, 2, 2, 1]

The same methods are available on any Indexable if you opt in to the core extensions:

require "hclust/core_ext"

coords.cluster(cutoff: 4) { |u, v| euclidean(u, v) }
coords.labels(2) { |u, v| euclidean(u, v) }
coords.centroids(cutoff: 4) { |u, v| euclidean(u, v) }

Alternatively, each step can be done manually. That is useful to inspect the dendrogram or to try different cutoffs without recomputing distances:

dism = HClust::DistanceMatrix.new(coords) { |u, v| euclidean(u, v) }
dendrogram = HClust.linkage(dism, :single)
clusters = dendrogram.flatten(4)
clusters.size # => 3
clusters      # => [[0, 3, 6, 7, 9], [1, 5, 8], [2, 4]]
dendrogram.labels(4)
# => [0, 1, 2, 0, 2, 1, 0, 0, 1, 0]

Dendrogram#flatten returns indexes into the original sequence. Dendrogram#labels returns the corresponding cluster mask. Refer to the API documentation for further details.

Benchmark

A Bash script is used to benchmark the code and be compared against reference implementations.

The script downloads the required libraries (except for SciPy) to run the corresponding code. The following programs are expected to be available: gcc for C++, cargo for Rust, and python for Python (SciPy must be installed in the current Python environment). Each benchmark will run for a number of times, and the best time will be printed out.

The following output was obtained on a machine with AMD® Ryzen 9 5950x under Pop!_OS 22.04 LTS using default benchmark values (see below):

$ bash bench/bench.sh
Testing Fastcluster (C++)...
Testing Kodama (Rust)...
Testing Scipy (Python)...
Testing HClust (Crystal)...
| name         | version | compiler     | time (ms) |
| ------------ | ------- | ------------ | --------- |
| fastcluster  | 1.2.6   | 11.2.0 (gcc) |     0.032 |
| kodama       | 0.2.3   | 1.61.0       |     0.041 |
| scipy        | 1.7.3   | 3.9.12       |     0.094 |
| hclust       | 1.0.0   | 1.4.1        |     0.067 |

The benchmark can be configured via the following environment variables:

• BENCH_SIZE sets the number of data points to bench (defaults to 100).
• BENCH_REPEATS sets the number of times to repeat the code during benchmark (defaults to 1000).
• BENCH_RULE sets the linkage rule to run: average, centroid, complete, median, single, ward (default), or weighted.
• BENCH_METHOD sets the clustering method to run: mst, chain, or generic (default).

Contributing

1. Fork it (https://github.com/franciscoadasme/hclust/fork)
2. Create your feature branch (git checkout -b my-new-feature)
3. Commit your changes (git commit -am 'Add some feature')
4. Push to the branch (git push origin my-new-feature)
5. Create a new Pull Request

Contributors

• Francisco Adasme - creator and maintainer