rustyml 0.14.0

A high-performance machine learning & deep learning library in pure Rust, offering ML algorithms and neural network support
Documentation
# 5.3. Clustering Metrics

Clustering has no residual to square and no confusion matrix, unlike [regression](./5.1._Regression_Metrics.md) and [classification](./5.2._Classification_Metrics.md). A clustering never names its groups. It only splits the samples into parts. This fact shapes every metric on this page. External metrics compare your partition against a ground-truth partition. Internal metrics score the geometry of the partition alone. Both kinds must ignore the names of the clusters. RustyML's clustering metrics live in `rustyml::metrics::clustering`. They are re-exported flat from `rustyml::metrics`, and also through the [prelude](../Chapter-01/1.5._The_Prelude_and_Imports.md). One line, `use rustyml::metrics::*;`, brings every function below into scope.

Like the rest of the `metrics` leaf module, these functions **panic** on precondition violations. Examples are mismatched lengths, empty input, and an out-of-range cluster count. They do not return the crate's `Error` type. This is a deliberate choice. The `metrics` module depends only on `ndarray` and `ahash`. It mirrors `ndarray`'s own rule of panicking on a dimension mismatch, instead of adding the error machinery from [1.6. Error Handling](../Chapter-01/1.6._Error_Handling.md). Validate your labels before you pass them to a metric, if a panic would be fatal in your program.

## 5.3.1. Two families of metric

The external metrics take two label arrays, `labels_true` and `labels_pred`, both of type `isize`. They measure how well the predicted partition reproduces the reference partition. Every metric on this page takes `isize` labels, the same type scikit-learn's `labels_` uses. This lets any clustering estimator in the crate feed any metric on this page. [KMeans](../Chapter-02/2.7._KMeans_Clustering.md), [Mean Shift](../Chapter-02/2.9._Mean_Shift.md), and [DBSCAN](../Chapter-02/2.8._DBSCAN.md) all return `Array1<isize>`, so no cast is ever needed. DBSCAN's and Mean Shift's `-1` noise label still needs attention (see [5.3.7](#537-handling-dbscans-noise-sentinel)), but it needs no type conversion. The internal metrics take the feature matrix `x` plus one label array. They score the partition against the geometry of the data, with no ground truth at all. Use internal metrics when you cluster unlabeled data, the usual case. Use external metrics when you have a reference partition and want to benchmark an algorithm against it.

| Function | Family | Inputs | Range | Perfect score | Chance-corrected |
|---|---|---|---|---|---|
| `adjusted_rand_index` | external | two `isize` label arrays | `[-0.5, 1.0]` | `1.0` | yes |
| `adjusted_mutual_info` | external | two `isize` label arrays | `[-1.0, 1.0]` typ. | `1.0` | yes |
| `normalized_mutual_info` | external | two `isize` label arrays | `[0.0, 1.0]` | `1.0` | **no** |
| `v_measure_score` | external | two `isize` label arrays | `[0.0, 1.0]` | `1.0` | no |
| `homogeneity_score` | external | two `isize` label arrays | `[0.0, 1.0]` | `1.0` | no |
| `completeness_score` | external | two `isize` label arrays | `[0.0, 1.0]` | `1.0` | no |
| `fowlkes_mallows_score` | external | two `isize` label arrays | `[0.0, 1.0]` | `1.0` | no |
| `silhouette_score` | internal | `x`, labels, metric | `[-1.0, 1.0]` | `1.0` (higher better) | n/a |
| `davies_bouldin_score` | internal | `x`, labels | `>= 0.0` | `0.0` (lower better) | n/a |
| `calinski_harabasz_score` | internal | `x`, labels | `>= 0.0` | higher better | n/a |

All the external metrics are symmetric except 2. Swapping `labels_true` and `labels_pred` leaves the score unchanged for the rest. The 2 exceptions are `homogeneity_score` and `completeness_score`. They are duals of each other, so swapping the arguments swaps the two scores. Their harmonic mean, `v_measure_score`, is symmetric again.

## 5.3.2. Why label matching is the wrong tool

Classification metrics give you the instinct to line up two label vectors element by element and count the matches. This instinct fails for clustering, because cluster identifiers are arbitrary. Suppose one run labels a group `0` and another run labels the same group `7`. Both describe the same partition. An element-wise "accuracy" still reports total disagreement. This is not a rare edge case. It happens every time you rerun [KMeans](../Chapter-02/2.7._KMeans_Clustering.md) with a different seed, because KMeans numbers its clusters by initialization order.

The metrics on this page avoid the problem. They work from the **contingency table** of the two labelings: how many samples fall into each (true cluster, predicted cluster) cell. They also work from pair-agreement counts derived from that table. Both stay unchanged under any renaming of the clusters. The example below swaps the names of the two clusters in a partition. A naive accuracy score collapses to 0. ARI and NMI stay pinned at `1.0`.

```rust
use ndarray::array;
use rustyml::metrics::{adjusted_rand_index, normalized_mutual_info};

fn main() {
    let truth = array![0isize, 0, 1, 1];
    // Same partition, cluster names swapped: {0,1} -> label 1, {2,3} -> label 0.
    let relabeled = array![1isize, 1, 0, 0];

    // Element-wise "accuracy" collapses to 0 even though the grouping is identical.
    let naive_accuracy = truth
        .iter()
        .zip(relabeled.iter())
        .filter(|&(a, b)| a == b)
        .count() as f64
        / truth.len() as f64;

    let ari = adjusted_rand_index(&truth, &relabeled);
    let nmi = normalized_mutual_info(&truth, &relabeled);

    println!("naive accuracy = {naive_accuracy}"); // 0.0
    println!("ARI = {ari}, NMI = {nmi}"); // 1.0, 1.0

    assert_eq!(naive_accuracy, 0.0);
    assert!((ari - 1.0).abs() < 1e-12);
    assert!((nmi - 1.0).abs() < 1e-12);
}
```

This permutation invariance means you never need to solve an assignment problem before scoring, such as Hungarian matching of predicted clusters to true classes. The metric already handles that step.

## 5.3.3. External metrics and chance correction

A subtler trap exists. A metric can be permutation-invariant and still mislead you, because random labelings do not score 0. The Rand index counts the fraction of the `C(n, 2)` sample pairs where two partitions agree. Agreement means both partitions group the pair together, or both partitions separate the pair. The expected value of the Rand index under random labeling is not 0. It increases as you add clusters. A bare Rand index of 0.7 means nothing on its own. The **adjusted** Rand index subtracts that expected value and rescales the result. Independent labelings then score about `0.0`. A perfect match scores `1.0`. Agreement that is systematically worse than random can go negative, down to about `-0.5`.

Mutual information has the same problem, in a sharper form. MI keeps rising as you cut the data into more clusters. It reaches the full entropy when every sample forms its own cluster. Comparing raw MI across candidate values of K is therefore meaningless. `adjusted_mutual_info` subtracts the **expected mutual information** (EMI). RustyML computes EMI exactly, under a hypergeometric model of random partitions with the same cluster sizes. It evaluates every binomial coefficient in log space, from a shared log-factorial table. AMI is therefore the mutual-information analogue of ARI. It scores about `0.0` for independent labelings, `1.0` for identical ones, and occasionally a small negative value.

`normalized_mutual_info` does **not** correct for chance. It only rescales MI into `[0.0, 1.0]`. It divides MI by the arithmetic mean of the two clusterings' entropies, `(H_true + H_pred) / 2`. This normalizer is a fixed convention in RustyML. There is no `average_method` switch like the one in scikit-learn. With this normalizer, NMI is numerically identical to `v_measure_score`, the harmonic mean of homogeneity and completeness. This gives a rule of thumb. When the two clusterings have different numbers of clusters, use `adjusted_rand_index` or `adjusted_mutual_info`. Reserve NMI and V-measure for comparisons at a fixed K, where the uncorrected bias stays constant and cancels out.

The three functions share one signature. Only the meaning of the score differs.

```rust,ignore
pub fn adjusted_rand_index<S>(labels_true: &ArrayBase<S, Ix1>, labels_pred: &ArrayBase<S, Ix1>) -> f64
where S: Data<Elem = isize>;
// identical for adjusted_mutual_info and normalized_mutual_info
```

The degenerate cases differ between these functions. This matters when a clustering trivially puts every point in one cluster. `adjusted_rand_index` returns `1.0` when there are fewer than 2 samples, because there are no pairs to disagree on. It also returns `1.0` when its normalizer vanishes. `adjusted_mutual_info` returns `1.0` when its own normalizer is degenerate. `normalized_mutual_info` returns `0.0` whenever either partition has a single cluster, because zero entropy makes its denominator zero. Do not read these constants as quality judgments. They are only the defined values for inputs where the ratio is `0/0`.

```rust
use ndarray::array;
use rustyml::metrics::{adjusted_mutual_info, adjusted_rand_index, normalized_mutual_info};

fn main() {
    let labels_true = array![0isize, 0, 1, 1, 2, 2];
    let labels_pred = array![0isize, 0, 1, 2, 1, 2]; // one true cluster split in two

    // Verified against the closed forms: ARI = 1/6 ~= 0.167, AMI = 1/6, NMI ~= 0.579.
    println!("ARI = {:.4}", adjusted_rand_index(&labels_true, &labels_pred));
    println!("AMI = {:.4}", adjusted_mutual_info(&labels_true, &labels_pred));
    println!("NMI = {:.4}", normalized_mutual_info(&labels_true, &labels_pred));

    // Independent labelings: ARI and AMI hit their chance-corrected floor of about -0.5.
    // NMI also reaches 0 here, because MI is exactly 0 for independent partitions.
    let a = array![0isize, 0, 1, 1];
    let b = array![0isize, 1, 0, 1];
    assert!((adjusted_rand_index(&a, &b) - (-0.5)).abs() < 1e-9);
    assert!((adjusted_mutual_info(&a, &b) - (-0.5)).abs() < 1e-9);
    assert!(normalized_mutual_info(&a, &b).abs() < 1e-12);
}
```

The remaining external metrics add detail. They do not replace ARI or AMI. `homogeneity_score` asks whether each predicted cluster is pure, meaning it contains only one true class. `completeness_score` is its dual: it asks whether each true class stays in one cluster. `v_measure_score` is the harmonic mean of the two. `fowlkes_mallows_score` is the geometric mean of pairwise precision and recall. All 4 metrics range over `[0.0, 1.0]`. All 4 are based on entropy or pair ratios, so none of them are chance-corrected. The same fixed-K caveat from above applies to all 4.

## 5.3.4. The silhouette score

Clustering usually gives you no ground truth. The silhouette score is the main tool for this case. For each sample, it computes `a`, the mean distance to the other members of the sample's own cluster. It also computes `b`, the mean distance to the members of the nearest other cluster. It combines the two into this formula.

```rust,ignore
s = (b - a) / max(a, b)
```

The score `s` runs from `-1` to `+1`. A score of `-1` means the sample sits closer to a neighboring cluster than to its own, so it is likely misassigned. A score of `0` means the sample lies on the boundary between clusters. A score of `+1` means a tight own cluster and distant neighbors, so the sample is well clustered. `silhouette_score` returns the mean of `s` over all samples. Two edge cases matter in practice. A sample that is the only member of its cluster contributes `0`, because there is no `a` to compute. When every point coincides, so `a = b = 0` everywhere, the score is `0` rather than `NaN`.

```rust,ignore
pub fn silhouette_score<S1, S2>(
    x: &ArrayBase<S1, Ix2>,
    labels: &ArrayBase<S2, Ix1>,
    metric: DistanceCalculationMetric,
) -> f64
where S1: Data<Elem = f64> + Sync, S2: Data<Elem = isize>;
```

The `metric` argument goes through [`DistanceCalculationMetric`](../Chapter-06/6.1._Distance_Metrics.md), the same dispatch point the estimators use. `Euclidean`, `Manhattan`, and `Minkowski(p)` all work. The metric genuinely changes the result, not just its label. Pass `DistanceCalculationMetric::Euclidean` for the conventional silhouette. The enum has its own `Default` value, but this function does not use it. You must name the metric yourself.

`silhouette_score` panics in 4 cases. The row count of `x` differs from the length of `labels`. The input is empty. The number of distinct clusters falls outside `2..=n_samples - 1`. A single cluster has no `b` to compute, and a partition of all singletons has no `a` to compute. The fourth case is passing `Minkowski(p)` with `p < 1`. The internal `davies_bouldin_score` and `calinski_harabasz_score` enforce the same length, non-empty, and cluster-count bounds. Neither function takes a `metric` parameter, so the `Minkowski(p)` check does not apply to them.

```rust
use ndarray::array;
use rustyml::math::DistanceCalculationMetric;
use rustyml::metrics::silhouette_score;

fn main() {
    // 2 tight, well-separated 2-D clusters (not collinear, so the metric matters).
    let x = array![[0.0, 0.0], [0.0, 1.0], [10.0, 10.0], [10.0, 11.0]];
    let labels = array![0isize, 0, 1, 1];

    let euclidean = silhouette_score(&x, &labels, DistanceCalculationMetric::Euclidean);
    let manhattan = silhouette_score(&x, &labels, DistanceCalculationMetric::Manhattan);
    println!("euclidean silhouette = {euclidean:.4}");
    println!("manhattan silhouette = {manhattan:.4}");

    assert!(euclidean > 0.8 && euclidean <= 1.0);
    assert!(manhattan > 0.8 && manhattan <= 1.0);
    // Different metrics genuinely give different scores on these points.
    assert!((euclidean - manhattan).abs() > 1e-3);
}
```

## 5.3.5. Cost and the parallel fill

The silhouette score has a cost for its thoroughness. Computing every `a` and `b` needs the distance from each sample to every other sample. For `n` samples in `d` dimensions, this computation is inherently `O(n^2 * d)`, quadratic in the number of points. `davies_bouldin_score` and `calinski_harabasz_score` cost much less, because they only touch centroids. `davies_bouldin_score` runs an `O(n)` scan of each point against its own centroid, plus an `O(k^2)` loop over centroid pairs. `calinski_harabasz_score` runs a single `O(n)` scan, with no pairwise centroid loop at all. The silhouette's quadratic cost dominates on any dataset of real size. This is where RustyML spends most of its engineering effort.

Two techniques keep the silhouette's cost manageable. First, the implementation never builds the full `n x n` distance matrix. It accumulates a compact `dist_to_cluster[[i, c]]` table instead, holding the total distance from sample `i` to each cluster `c`. Memory use is `O(n * k)`, not `O(n^2)`. Second, the implementation uses the symmetry `d(i, j) = d(j, i)`. It scans only the upper triangle of the distance matrix, which halves the number of metric evaluations compared to a full scan. This halving matters more as the metric gets more expensive. `Manhattan` is cheapest. `Euclidean` adds a square root. `Minkowski(p)` adds a `powf` call and costs the most.

Above a work threshold, the upper-triangle fill runs in parallel. The gate is measured in scanned elements, `scan_work = n * n * d`. When `scan_work` reaches `SILHOUETTE_PARALLEL_MIN_ELEMS` (default `262_144`), the rows are dealt round-robin into `rayon`'s `current_num_threads()` buckets. Row `i` does `n - 1 - i` pair evaluations, so round-robin dealing balances the buckets better than a contiguous split would. Each bucket folds into its own accumulator, and the buckets are then summed in a fixed order. This fixed grouping makes the parallel result reproducible on the same machine across runs. The parallel result equals the serial fill numerically, though not always bit-for-bit. Below the gate, the serial path runs and matches a full scan bit-for-bit. A dedicated benchmark covers this fill.

```bash
cargo bench --bench silhouette
```

If the default crossover point is wrong for your hardware or data shape, change it at runtime. Use the tuning facade, `rustyml::tuning::metrics::set_silhouette(value)` and `get_silhouette()`. [7.3. Performance Tuning and Parallelism](../Chapter-07/7.3._Performance_Tuning_and_Parallelism.md) covers this mechanism and the reasoning behind parallel gates in general.

## 5.3.6. Picking K with a silhouette sweep

[KMeans](../Chapter-02/2.7._KMeans_Clustering.md) needs the number of clusters, K, chosen before it runs. The silhouette score is the best-known tool for choosing K. Fit the model for each candidate K. Score the resulting partition. Keep the K with the highest mean silhouette. KMeans returns its labels as `Array1<isize>`, which feeds `silhouette_score` directly, with no conversion.

```rust
use ndarray::array;
use rustyml::machine_learning::KMeans;
use rustyml::math::DistanceCalculationMetric;
use rustyml::metrics::silhouette_score;

fn main() {
    // 3 compact, well-separated blobs, 4 points each.
    let x = array![
        [0.0, 0.0], [0.2, 0.1], [0.1, 0.2], [0.0, 0.3],
        [5.0, 5.0], [5.2, 5.1], [5.1, 5.2], [5.0, 5.3],
        [0.0, 5.0], [0.2, 5.1], [0.1, 5.2], [0.0, 5.3],
    ];

    let mut best_k = 0usize;
    let mut best_score = f64::NEG_INFINITY;

    for k in 2..=5 {
        let labels = KMeans::new(k, 100, 1e-4)
            .unwrap()
            .with_random_state(42) // fixed seed for a reproducible sweep
            .fit_predict(&x)
            .unwrap();

        // silhouette_score needs 2..=n-1 distinct clusters. Skip a fit that collapsed a cluster.
        let mut distinct = labels.to_vec();
        distinct.sort_unstable();
        distinct.dedup();
        if distinct.len() < 2 {
            continue;
        }

        let s = silhouette_score(&x, &labels, DistanceCalculationMetric::Euclidean);
        println!("k = {k}: silhouette = {s:.4}");
        if s > best_score {
            best_score = s;
            best_k = k;
        }
    }

    println!("best k = {best_k} (silhouette = {best_score:.4})");
    assert_eq!(best_k, 3); // the 3 real blobs win
}
```

The fixed `with_random_state(42)` call makes the sweep reproducible. A different seed can move the KMeans initialization, and at the margin, it can change the winning K. See [7.1. Reproducibility and Random Seeds](../Chapter-07/7.1._Reproducibility_and_Random_Seeds.md). The `distinct.len() < 2` guard matters for a specific reason. KMeans can return an empty cluster, which would drop the distinct count below the silhouette's lower bound and cause a panic. Skipping such a fit is cheaper and clearer than catching the panic. For a faster sweep on large `n`, run `davies_bouldin_score` (lower is better) or `calinski_harabasz_score` (higher is better) instead. `davies_bouldin_score` costs an `O(n)` scan plus an `O(k^2)` loop over centroid pairs. `calinski_harabasz_score` costs a single `O(n)` scan. Both avoid the silhouette's quadratic cost, at the price of a coarser view of cluster shape that looks only at centroids.

## 5.3.7. Handling DBSCAN's noise sentinel

Labels from [DBSCAN](../Chapter-02/2.8._DBSCAN.md), or from [Mean Shift](../Chapter-02/2.9._Mean_Shift.md) with `cluster_all = false`, type-check against every metric on this page without a cast. Both the estimators and the metrics use `isize`. This removes the mechanical friction, but not the semantic problem. None of these functions has any notion of a noise sentinel. Each distinct label value counts as a full cluster, so a `-1` becomes an artificial noise cluster. For the silhouette score, this means the scattered noise points get scored as if they formed a real group. This is almost never what you want. (scikit-learn's silhouette score behaves the same way.)

The clean approach is to drop the noise rows before scoring. Subset both `x` and the labels down to the points DBSCAN actually clustered.

```rust
use ndarray::{array, Array1, Axis};
use rustyml::machine_learning::DBSCAN;
use rustyml::math::DistanceCalculationMetric;
use rustyml::metrics::silhouette_score;

fn main() {
    // 2 dense blobs plus one far-flung outlier.
    let x = array![
        [0.0, 0.0], [0.1, 0.0], [0.0, 0.1], [0.1, 0.1],
        [5.0, 5.0], [5.1, 5.0], [5.0, 5.1], [5.1, 5.1],
        [50.0, 50.0], // noise
    ];

    // eps = 1.0, min_samples = 3: each blob is a core cluster, the outlier is noise (-1).
    let labels = DBSCAN::new(1.0, 3).unwrap().fit_predict(&x).unwrap();

    // Keep only clustered rows (label >= 0). No cast: the labels are already isize.
    let keep: Vec<usize> = labels
        .iter()
        .enumerate()
        .filter(|&(_, &l)| l >= 0)
        .map(|(i, _)| i)
        .collect();

    let x_clustered = x.select(Axis(0), keep.as_slice());
    let labels_clustered = Array1::from_iter(keep.iter().map(|&i| labels[i]));

    let s = silhouette_score(
        &x_clustered,
        &labels_clustered,
        DistanceCalculationMetric::Euclidean,
    );
    println!("silhouette over non-noise points = {s:.4}"); // ~1.0, 2 clean blobs
    assert!(s > 0.9);
}
```

Filtering answers the usual question: how well separated are the points that formed real clusters. It also hides how much data was discarded, so report the noise fraction next to the score. An alternative is to keep the noise points under their own `-1` label. This is defensible only for the external metrics. ARI, AMI, and NMI will then compare a noise class against your ground truth like any other cluster. This is a coherent evaluation, though a strict one. The silhouette score does not support this alternative. A diffuse cloud of noise is not a cluster, and scoring it as one only muddies the result.