5.3. Clustering Metrics

Clustering has no residual to square and no confusion matrix, unlike regression and classification. 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. 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. 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, Mean Shift, and DBSCAN 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), 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.

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 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.

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.

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.

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.

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.

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, 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.

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.

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 covers this mechanism and the reasoning behind parallel gates in general.

5.3.6. Picking K with a silhouette sweep

KMeans 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.

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. 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, or from Mean Shift 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.

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.