# 5.3. 聚类指标
聚类没有可以平方的残差，也没有混淆矩阵，这一点和[回归](./5.1._回归指标.md)、[分类](./5.2._分类指标.md)都不一样。聚类从不给自己的分组起名字，它只是把样本切分成若干部分。这一点决定了本页每一个指标的设计：外部指标拿你的划分和一份真值划分作对比，内部指标只看划分本身的几何形状打分，两者都必须无视簇叫什么名字。RustyML 的聚类指标位于 `rustyml::metrics::clustering`，被平铺重新导出到 `rustyml::metrics`，也可以通过 [prelude](../Chapter-01/1.5._Prelude与模块导入.md) 拿到。写一句 `use rustyml::metrics::*;`，就能把下面所有函数都引入作用域。
和 `metrics` 这个叶子模块的其余部分一样，这些函数在前置条件被违反时会 **panic**，例如长度不匹配、输入为空、簇数量越界，而不是返回 crate 的 `Error` 类型。这是有意为之：`metrics` 只依赖 `ndarray` 和 `ahash`，所以它照搬了 `ndarray` 自己"维度不匹配就 panic"的规则，而不去引入 [1.6. 错误处理](../Chapter-01/1.6._错误处理.md)里那套错误处理机制。如果一次 panic 会让程序崩掉，就在把标签交给指标之前先自行校验。
## 5.3.1. 两大类指标
外部指标接收两个标签数组，`labels_true` 和 `labels_pred`，都是 `isize` 类型，衡量预测划分复现参考划分的程度。本页每个指标都收 `isize` 标签，这正是 scikit-learn `labels_` 使用的类型。正因如此，crate 里任何一个聚类估计器都能喂给本页任何一个指标：[KMeans](../Chapter-02/2.7._KMeans聚类.md)、[MeanShift](../Chapter-02/2.9._MeanShift.md)、[DBSCAN](../Chapter-02/2.8._DBSCAN.md) 全都返回 `Array1<isize>`，永远不需要类型转换。DBSCAN 和 MeanShift 的 `-1` 噪声标签仍然需要你留意（见 [5.3.7](#537-处理-dbscan-的噪声哨兵值)），但已经不需要转换类型了。内部指标接收特征矩阵 `x` 加一个标签数组，拿划分和数据自身的几何形状比对打分，完全不需要真值。对无标签数据做聚类时（这是常态）用内部指标；手上有一份参考划分、想拿它给某个算法当基准时用外部指标。
| 函数 | 类别 | 输入 | 取值范围 | 完美得分 | 机会校正 |
|---|---|---|---|---|---|
| `adjusted_rand_index` | 外部 | 两个 `isize` 标签数组 | `[-0.5, 1.0]` | `1.0` | 是 |
| `adjusted_mutual_info` | 外部 | 两个 `isize` 标签数组 | 通常 `[-1.0, 1.0]` | `1.0` | 是 |
| `normalized_mutual_info` | 外部 | 两个 `isize` 标签数组 | `[0.0, 1.0]` | `1.0` | **否** |
| `v_measure_score` | 外部 | 两个 `isize` 标签数组 | `[0.0, 1.0]` | `1.0` | 否 |
| `homogeneity_score` | 外部 | 两个 `isize` 标签数组 | `[0.0, 1.0]` | `1.0` | 否 |
| `completeness_score` | 外部 | 两个 `isize` 标签数组 | `[0.0, 1.0]` | `1.0` | 否 |
| `fowlkes_mallows_score` | 外部 | 两个 `isize` 标签数组 | `[0.0, 1.0]` | `1.0` | 否 |
| `silhouette_score` | 内部 | `x`、标签、metric | `[-1.0, 1.0]` | `1.0`（越高越好） | 不适用 |
| `davies_bouldin_score` | 内部 | `x`、标签 | `>= 0.0` | `0.0`（越低越好） | 不适用 |
| `calinski_harabasz_score` | 内部 | `x`、标签 | `>= 0.0` | 越高越好 | 不适用 |
外部指标里除 2 个之外都是对称的：交换 `labels_true` 和 `labels_pred`，其余指标的得分都不变。这 2 个例外是 `homogeneity_score` 和 `completeness_score`，它们互为对偶，交换参数会互换两个得分。它们的调和平均 `v_measure_score` 又是对称的。
## 5.3.2. 为什么标签匹配是错误的工具
从分类指标那边带来的直觉，是把两个标签向量逐元素对齐、数一数有多少个对上了。这套直觉在聚类上不成立，因为簇的编号本来就是任意的。设想一次运行把某个分组标成 `0`，另一次运行把同一个分组标成 `7`：两者描述的是同一个划分，可逐元素的"准确率"依然会报告完全不一致。这不是什么罕见的边角情况，每次换个种子重跑 [KMeans](../Chapter-02/2.7._KMeans聚类.md) 都会发生，因为 KMeans 是按初始化顺序给簇编号的。
本页的指标避开了这个问题。它们基于两组标签的**列联表**（contingency table）：有多少样本落进每个（真值簇，预测簇）的格子里。它们也基于由列联表派生出的成对一致性计数。这两样东西在簇被任意重命名后都保持不变。下面的例子把一个划分里两个簇的名字对调，可以看到朴素准确率直接掉到 0，而 ARI 和 NMI 稳稳钉在 `1.0`。
```rust
use ndarray::array;
use rustyml::metrics::{adjusted_rand_index, normalized_mutual_info};
fn main() {
let truth = array![0isize, 0, 1, 1];
// 相同的划分，只是簇名对调：{0,1} -> 标签 1，{2,3} -> 标签 0。
let relabeled = array![1isize, 1, 0, 0];
// 分组完全相同，逐元素的"准确率"却掉到了 0。
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);
}
```
正因为有这种置换不变性，打分之前永远不需要先解一个指派问题，比如用匈牙利算法把预测簇和真值类别配对。指标本身已经处理了这一步。
## 5.3.3. 外部指标与机会校正
还有一个更隐蔽的陷阱：一个指标可以是置换不变的，却依然会误导你，因为随机标签并不会得 0 分。兰德指数（Rand index）统计的是两个划分在全部 `C(n, 2)` 个样本对里达成一致的比例，一致的意思是两个划分都把这对样本分到一起，或者都把它们分开。兰德指数在随机标签下的期望值不是 0，而且随着簇数增多还会升高，所以孤零零一个 0.7 的兰德指数本身说明不了任何问题。**调整**兰德指数减掉这个期望值再重新缩放，于是相互独立的标签得分约为 `0.0`，完美匹配得 `1.0`，而系统性地比随机还差的一致程度可以为负，低到约 `-0.5`。
互信息（mutual information）有同样的问题，而且更严重：MI 会随着数据被切成更多簇而不断上升，等到每个样本自成一簇时便达到完整的熵。因此拿不同 K 候选值下的原始 MI 互相比较毫无意义。`adjusted_mutual_info` 减去了**期望互信息**（EMI）。RustyML 在一个簇大小相同的随机划分超几何模型下精确计算 EMI，用一张共享的对数阶乘表在对数空间里求出每一个二项式系数。因此 AMI 就是 ARI 在互信息上的对应物：对相互独立的标签约为 `0.0`，对完全相同的标签为 `1.0`，偶尔会略微为负。
`normalized_mutual_info` **不**做机会校正，它只是把 MI 缩放到 `[0.0, 1.0]`。做法是把 MI 除以两个聚类结果熵的算术平均，`(H_true + H_pred) / 2`。这个归一化方式在 RustyML 里是固定的约定，不像 scikit-learn 那样提供 `average_method` 开关。采用这个归一化因子后，NMI 在数值上与 `v_measure_score`（同质性与完整性的调和平均）完全相等。由此得到一条经验法则：当两个聚类结果的簇数量不同时，用 `adjusted_rand_index` 或 `adjusted_mutual_info`；把 NMI 和 V-measure 留给固定 K 的比较，那种情况下未经校正的偏差是个常数，会相互抵消。
这三个函数共用同一套签名，区别只在得分的含义。
```rust,ignore
pub fn adjusted_rand_index<S>(labels_true: &ArrayBase<S, Ix1>, labels_pred: &ArrayBase<S, Ix1>) -> f64
where S: Data<Elem = isize>;
// adjusted_mutual_info 和 normalized_mutual_info 完全相同
```
这几个函数的退化情形各不相同，而当某个聚类平凡地把所有点都塞进一个簇时，这种差异就很关键。`adjusted_rand_index` 在样本少于 2 个时返回 `1.0`，因为没有样本对可供分歧；它在归一化因子消失时也返回 `1.0`。`adjusted_mutual_info` 在自己的归一化因子退化时返回 `1.0`。`normalized_mutual_info` 只要任一划分只有一个簇就返回 `0.0`，因为零熵会让分母为零。不要把这些常数当成质量判断，它们只是比值为 `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]; // 一个真值簇被拆成了两个
// 对照闭式解验证过：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));
// 相互独立的标签：ARI 和 AMI 触到了各自经机会校正后的下限，约为 -0.5。
// 这里 NMI 也归零，因为对相互独立的划分，MI 精确为 0。
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);
}
```
其余的外部指标是在补充细节，不是替代 ARI 或 AMI。`homogeneity_score` 问的是每个预测簇是否纯净，也就是只包含一个真值类别。`completeness_score` 是它的对偶：问的是每个真值类别是否都待在同一个簇里。`v_measure_score` 是这两者的调和平均。`fowlkes_mallows_score` 是成对精确率与召回率的几何平均。这 4 个指标都落在 `[0.0, 1.0]` 区间，而且都基于熵或成对比值，因此都不做机会校正，前面固定 K 的告诫对这 4 个同样适用。
## 5.3.4. 轮廓系数
聚类通常没有真值可用，轮廓系数（silhouette score）就是应对这种情况的主力工具。对每个样本，它计算 `a`，即该样本到自己簇内其他成员的平均距离；再计算 `b`，即到最近的那个其他簇成员的平均距离；然后把两者组合成下面的公式。
```rust,ignore
s = (b - a) / max(a, b)
```
`s` 的取值从 `-1` 到 `+1`。得分 `-1` 表示样本离相邻簇比离自己簇还近，多半是分错了。得分 `0` 表示样本正好落在簇与簇的边界上。得分 `+1` 表示自己簇紧凑、邻簇遥远，样本聚得很好。`silhouette_score` 返回的是 `s` 在所有样本上的均值。有两种边界情况在实践中很关键：某个样本若是自己簇里唯一的成员，它贡献 `0`，因为根本没有 `a` 可算；当所有点重合、处处 `a = b = 0` 时，得分是 `0` 而不是 `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>;
```
`metric` 参数走的是 [`DistanceCalculationMetric`](../Chapter-06/6.1._距离度量.md)，和各个估计器用的是同一个调度点。`Euclidean`、`Manhattan`、`Minkowski(p)` 都能用，而且度量方式实实在在地改变结果，不只是换个标签而已。想要常规的轮廓系数就传 `DistanceCalculationMetric::Euclidean`。这个枚举自带 `Default` 值，但本函数不使用它，你必须自己指明度量方式。
`silhouette_score` 会在 4 种情况下 panic：`x` 的行数和 `labels` 的长度不一致；输入为空；不同簇的数量落在 `2..=n_samples - 1` 之外（只有一个簇就没有 `b` 可算，全是单点簇的划分就没有 `a` 可算）；第四种情况是传了 `p < 1` 的 `Minkowski(p)`。内部的 `davies_bouldin_score` 和 `calinski_harabasz_score` 强制执行同样的长度、非空、簇数量边界，但这两个函数都不接收 `metric` 参数，所以 `Minkowski(p)` 那条检查对它们不适用。
```rust
use ndarray::array;
use rustyml::math::DistanceCalculationMetric;
use rustyml::metrics::silhouette_score;
fn main() {
// 2 个紧凑、彼此分得很开的二维簇（不共线，所以度量方式会起作用）。
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);
// 不同的度量方式在这些点上确实会给出不同的得分。
assert!((euclidean - manhattan).abs() > 1e-3);
}
```
## 5.3.5. 计算成本与并行填充
轮廓系数的全面是有代价的。要算出每一个 `a` 和 `b`，就需要每个样本到所有其他样本的距离。对 `d` 维空间中的 `n` 个样本来说，这项计算本质上是 `O(n^2 * d)`，关于点数是二次的。`davies_bouldin_score` 和 `calinski_harabasz_score` 的代价小得多，因为它们只碰质心：`davies_bouldin_score` 对每个点相对自己质心做一次 `O(n)` 扫描，再加一个遍历质心两两组合的 `O(k^2)` 循环；`calinski_harabasz_score` 只做一次 `O(n)` 扫描，完全没有质心两两配对的循环。在任何有点规模的数据集上，轮廓系数的二次代价都占主导，RustyML 的大部分工程功夫也正花在这里。
有两项技术让轮廓系数的代价保持在可控范围。其一，实现从不构建完整的 `n x n` 距离矩阵，而是累积一张紧凑的 `dist_to_cluster[[i, c]]` 表，记录样本 `i` 到每个簇 `c` 的距离总和，内存占用是 `O(n * k)`，不是 `O(n^2)`。其二，实现利用了 `d(i, j) = d(j, i)` 这一对称性，只扫描距离矩阵的上三角，相比全量扫描把度量计算次数砍掉一半。度量越贵，这一半就越划算：`Manhattan` 最便宜，`Euclidean` 多一次开方，`Minkowski(p)` 多一次 `powf` 调用，代价最高。
超过某个工作量阈值后，上三角的填充就会并行执行。这道闸门以扫描元素数来衡量，`scan_work = n * n * d`。一旦 `scan_work` 达到 `SILHOUETTE_PARALLEL_MIN_ELEMS`（默认 `262_144`），各行就会以轮转（round-robin）方式分发到 `rayon` 的 `current_num_threads()` 个桶里。第 `i` 行要做 `n - 1 - i` 次样本对计算，所以轮转分配比连续切分更能让各个桶负载均衡。每个桶折叠进自己的累加器，随后各桶按固定顺序求和。这种固定的分组让并行结果在同一台机器上多次运行都可复现。并行结果在数值上等于串行填充的结果，但未必逐比特相同。低于这道闸门时走串行路径，结果与全量扫描逐比特一致。这项填充有专门的基准测试。
```bash
cargo bench --bench silhouette
```
如果这个默认的切换点不适合你的硬件或数据形状，可以在运行时调整它。使用调优门面，`rustyml::tuning::metrics::set_silhouette(value)` 和 `get_silhouette()`。[7.3. 性能调优与并行](../Chapter-07/7.3._性能调优与并行.md)讲了这套机制，以及关于并行闸门的通用原理。
## 5.3.6. 用轮廓系数扫描来选 K
[KMeans](../Chapter-02/2.7._KMeans聚类.md) 需要在运行之前先定好簇的数量 K。轮廓系数是挑选 K 最广为人知的工具。对每个候选 K 拟合一遍模型，给得到的划分打分，留下平均轮廓系数最高的那个 K。KMeans 把标签以 `Array1<isize>` 返回，无需任何转换就能直接喂给 `silhouette_score`。
```rust
use ndarray::array;
use rustyml::machine_learning::KMeans;
use rustyml::math::DistanceCalculationMetric;
use rustyml::metrics::silhouette_score;
fn main() {
// 3 个紧凑、彼此分得很开的团簇，每个 4 个点。
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) // 固定种子，扫描可复现
.fit_predict(&x)
.unwrap();
// silhouette_score 需要 2..=n-1 个不同的簇，跳过那种把某个簇塌缩掉的拟合。
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); // 3 个真实的团簇胜出
}
```
固定的 `with_random_state(42)` 调用让这次扫描可复现，换个种子可能改变 KMeans 的初始化，在临界处甚至会改变胜出的 K。参见 [7.1. 可复现性与随机种子](../Chapter-07/7.1._可复现性与随机种子.md)。`distinct.len() < 2` 这道守卫有其具体理由：KMeans 可能返回一个空簇，那会让不同簇的数量掉到轮廓系数的下界以下，从而引发 panic。跳过这样的拟合，比去捕获 panic 更省事也更清晰。要在大 `n` 上跑得更快，就改用 `davies_bouldin_score`（越低越好）或 `calinski_harabasz_score`（越高越好）。`davies_bouldin_score` 的代价是一次 `O(n)` 扫描加一个遍历质心两两组合的 `O(k^2)` 循环，`calinski_harabasz_score` 的代价只是一次 `O(n)` 扫描。两者都避开了轮廓系数的二次代价，代价是对簇形状只有一个更粗糙、仅看质心的视角。
## 5.3.7. 处理 DBSCAN 的噪声哨兵值
来自 [DBSCAN](../Chapter-02/2.8._DBSCAN.md)，或者 `cluster_all = false` 的 [MeanShift](../Chapter-02/2.9._MeanShift.md) 的标签，不需要任何转换就能通过本页每个指标的类型检查，因为估计器和指标两边都用 `isize`。这消掉了机械上的摩擦，却没有消掉语义上的问题：**这些函数没有任何一个懂得噪声哨兵值的概念**。每一个不同的标签值都被当成一个完整的簇，于是 `-1` 就成了一个人造的噪声簇。对轮廓系数而言，这意味着零散的噪声点会被当作一个真实分组来打分，而这几乎从不是你想要的结果。（scikit-learn 的轮廓系数行为完全一样。）
干净的做法是在打分之前把噪声行丢掉。把 `x` 和标签都取子集，只留下 DBSCAN 真正聚成簇的那些点。
```rust
use ndarray::{array, Array1, Axis};
use rustyml::machine_learning::DBSCAN;
use rustyml::math::DistanceCalculationMetric;
use rustyml::metrics::silhouette_score;
fn main() {
// 2 个稠密的团簇，外加一个远远甩出去的离群点。
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], // 噪声
];
// eps = 1.0，min_samples = 3：每个团簇是一个核心簇，离群点是噪声（-1）。
let labels = DBSCAN::new(1.0, 3).unwrap().fit_predict(&x).unwrap();
// 只保留成簇的行（标签 >= 0）。不需要转换：标签本来就是 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 个干净的团簇
assert!(s > 0.9);
}
```
过滤回答的是通常真正要问的问题：那些确实聚成了簇的点，彼此分得有多开。但它也掩盖了到底丢弃了多少数据，所以要把噪声占比和得分一起报告出来。另一种做法，是把噪声点保留在它们自己的 `-1` 标签下，这只对外部指标站得住脚：ARI、AMI、NMI 会把这个噪声类别当成任意其他簇一样和你的真值比对，这是一种自洽的评估，虽然严格。轮廓系数不支持这种做法，因为一团弥散的噪声根本不是簇，把它当簇来打分只会把结果搅浑。