# 2.7. KMeans聚类
`KMeans` 把样本划分到固定数量的簇里。它交替执行两步：先把每个点分配给最近的质心，再把每个质心重新计算为其成员的均值。RustyML 的实现是一套并行的 Lloyd 算法，配上 k-means++ 初始化和「重启若干次、取最优」的策略。如果你熟悉 scikit-learn 的 `sklearn.cluster.KMeans`，那套心智模型的大部分都能直接搬过来用。唯一一处刻意不同的默认值是 `n_init`，原因见 [2.7.1](#271-估计器的工作原理)。
## 2.7.1. 估计器的工作原理
一次 `fit` 调用会把整套流程跑 `n_init` 遍，保留 inertia 最低的那一遍。每一遍里，质心都从 **k-means++** 初始化开始。第一个中心是从样本中均匀随机抽取的一个点。
之后每个中心都从剩余的点里抽取，抽中的概率正比于该点到最近的已选中心的平方距离（经典的 D^2 轮盘赌）。这样能把初始中心尽量拉开，给 Lloyd 迭代一个比均匀随机初始化更好的起点。如果所有候选点的平方距离都为零，实现会转而对这个中心做均匀随机挑选。这种情况只有在所有点都与已选中心重复时才会发生。
初始化之后，每轮迭代做两件事：分配步把每个点归到最近的质心，更新步把每个质心替换成分配给它的那些点的均值。
当质心不再移动时，迭代停止。收敛判据比较的是所有质心的**平方位移之和**与一个按方差缩放的容差。这个阈值等于数据每维总体方差的均值，再乘以 `tol`。所以 `tol` 是*相对*容差，而不是绝对距离。这与 scikit-learn 的约定一致。收敛解是一个不动点：每个质心都恰好等于分配给它的点的均值。
一次因为用光预算、而不是因为收敛而结束的拟合，在返回之前会多跑一趟分配。Lloyd 循环先按*当前*质心给点打标签，然后才装入更新后的质心。所以停在 `max_iterations` 上，会让 `labels` 和 `inertia` 描述的是旧质心，而 `get_centroids` 报告的却是新质心。这样一来，`predict(x)` 就会与 `get_labels()` 对不上。最后这一趟按模型真正存下来的质心重新分配。scikit-learn 出于同样的原因重跑它最后那一步 E-step。
与 scikit-learn 唯一一处刻意的差别是重启次数。**这里 `n_init` 默认是 10。scikit-learn 在 k-means++ 下的 `n_init='auto'` 是 1。** 这不是疏漏。scikit-learn 的默认值建立在它的*贪心* k-means++ 之上，每个中心会抽 `2 + ln(k)` 个候选并留下最好的那个。所以那里单次初始化本身方差就已经很低。
RustyML 用的是朴素 k-means++：每个中心只做一次 D^2 抽取。重启机制存在的意义，正是为了弥补这种更高方差的初始化。想要严格对齐 scikit-learn，就传 `with_n_init(1)`。这样做之后，预期拟合结果会随之改变。各次重启的种子由 `random_state` 确定性地派生，所以带种子的拟合仍然可复现。
## 2.7.2. 构造模型
构造函数接收 3 个位置参数，并会即时校验。配置不合法会在构造时就失败，而不是拖到 `fit`。
```rust,ignore
pub fn new(n_clusters: usize, max_iterations: usize, tolerance: f64) -> Result<Self, Error>
```
| 参数 | 位置 | 类型 | 含义 |
| --- | --- | --- | --- |
| `n_clusters` | 第 1 个 | `usize` | 要形成的簇数 *k*。必须大于 0。 |
| `max_iterations` | 第 2 个 | `usize` | 每次重启内部 Lloyd 迭代的上限。必须大于 0。 |
| `tolerance` | 第 3 个 | `f64` | 相对收敛容差（按特征方差缩放）。必须为正且有限。 |
任何一条不满足都会返回 [`Error::InvalidParameter`](../Chapter-01/1.6._错误处理.md)：`n_clusters == 0`、`max_iterations == 0`，或者 `tolerance` 为零、为负、`NaN` 或无穷。注意它和数据校验错误之间的不对称。越界的*超参数*会返回 `InvalidParameter`。*数据*里的非有限值则返回 `NonFinite`。这个区分是全 crate 通用的约定。
`KMeans::default()` 给你 `n_clusters = 8`、`max_iterations = 300`、`tolerance = 1e-4`、`n_init = 10`，且不带种子。这份配置永远合法。
余下两项设置是 builder 步骤：它们各自拿走实例的所有权、再把它交还，所以可以接在 `new` 后面链式调用。`with_random_state(seed)` 不会失败。`with_n_init(n)` 返回 `Result`，因为 `0` 次重启会触发 `Error::InvalidParameter`：
```rust,ignore
let mut km = KMeans::new(3, 300, 1e-4)
.unwrap()
.with_n_init(1) // 对齐 scikit-learn；返回 Result
.unwrap()
.with_random_state(42); // 返回 Self
```
不带种子时，k-means++ 从系统熵取随机数，每次 `fit` 都会得到不同的划分。带上种子，`fit` 就可复现，详见 [2.7.6](#276-可复现的聚类局部种子与全局种子)。这种可复现性也覆盖每一次重启：每次重启都会从你设的种子确定性地派生出自己的子种子。
## 2.7.3. 训练、预测与读取结果
驱动模型的方法有 3 个，它们在返回值和改动的状态上各不相同。
`fit(&mut self, data)` 原地训练模型，返回 `Result<&mut Self, Error>`。它会算出并存下质心、训练集标签、inertia 和迭代次数。这 4 个值都描述胜出的那次重启。
`predict(&self, data)` 把新矩阵的每一行归到最近的已拟合质心，返回一个持有所有权的 `Array1<isize>`。它以不可变方式借用 `self`，不改动任何已存状态，所以拟合一次之后，想调用多少次都行。
标签用的是有符号类型，而非无符号类型。这样才能让 crate 里每个聚类估计器共享同一种标签类型。[DBSCAN](./2.8._DBSCAN.md) 和 [MeanShift](./2.9._MeanShift.md) 都用 `-1` 表示噪声，共享的标签类型让它们不加转换就能喂给 [5.3. 聚类指标](../Chapter-05/5.3._聚类指标.md) 里的各项指标。k-means 自己永远不会返回负标签。
`fit_predict(&mut self, data)` 会调用 `fit`，然后返回训练标签的一份克隆。只想要刚训练那批数据的标签时，用它。想给新样本打分时，改用 `fit` 加 `predict`。
已拟合的状态通过一组 getter 暴露出来，`fit` 之前它们都返回 `None`：
| Getter | 返回 | 说明 |
| --- | --- | --- |
| `get_centroids()` | `Option<&Array2<f64>>` | 形状 `(n_clusters, n_features)`。对应 scikit-learn 的 `cluster_centers_`。 |
| `get_labels()` | `Option<&Array1<isize>>` | 训练集的归属。对应 scikit-learn 的 `labels_`。 |
| `get_inertia()` | `Option<f64>` | 各点到最近质心的平方距离之和。对应 scikit-learn 的 `inertia_`。无论是否收敛，都与 `get_centroids()` 保持一致。 |
| `get_actual_iterations()` | `Option<usize>` | 胜出那次重启实际跑的迭代次数，`1..=max_iterations`。对应 scikit-learn 的 `n_iter_`。 |
| `get_n_clusters()` / `get_max_iterations()` / `get_tolerance()` / `get_n_init()` / `get_random_state()` | 普通值 / `Option<u64>` | 回读配置本身。 |
`fit` 会报 3 类数据错误。矩阵行数为零时，返回 [`Error::EmptyInput`](../Chapter-01/1.6._错误处理.md)。任何值为 `NaN` 或无穷时，返回 [`Error::NonFinite`](../Chapter-01/1.6._错误处理.md)。样本数少于簇数时，返回 [`Error::InvalidInput`](../Chapter-01/1.6._错误处理.md)，因为点数不到 *k* 个，就摆不下 *k* 个质心。
`predict` 还会另外报 2 种错误。在 `fit` 之前调用，返回 [`Error::NotFitted`](../Chapter-01/1.6._错误处理.md)。特征数与训练数据对不上，返回 [`Error::DimensionMismatch`](../Chapter-01/1.6._错误处理.md)。它还会做和 `fit` 一样的空输入与非有限值检查。
## 2.7.4. 端到端聚类三个数据簇
3 个紧凑、彼此分得很开的数据簇是最经典的正确性自检。任何正确的 k=3 拟合都必须把每个数据簇单独分成一类。这个例子用了一份小而确定的数据集，数据本身不含随机数，再加一个固定种子。它瞬间就能跑完，结果的结构也在意料之中。
```rust
use rustyml::machine_learning::KMeans;
use ndarray::{array, Array2};
fn main() {
// 三个数据簇，各 5 个点，中心在 (0,0)、(10,0)、(5,10)。
let data: Array2<f64> = array![
[-0.05, 0.03], [ 0.04, -0.02], [ 0.01, 0.05], [-0.03, -0.04], [ 0.02, 0.01],
[ 9.95, 0.03], [10.04, -0.02], [10.01, 0.05], [ 9.97, -0.04], [10.02, 0.01],
[ 4.95, 10.03], [ 5.04, 9.98], [ 5.01, 10.05], [ 4.97, 9.96], [ 5.02, 10.01],
];
let mut km = KMeans::new(3, 300, 1e-4).unwrap().with_random_state(42);
km.fit(&data).unwrap();
let labels = km.get_labels().unwrap();
let centroids = km.get_centroids().unwrap();
println!("labels: {:?}", labels);
println!("centroids: {:?}", centroids);
println!("inertia: {:.6}", km.get_inertia().unwrap());
println!("iters: {}", km.get_actual_iterations().unwrap());
// 给落在各数据簇真实中心上的新样本打分。
let new_points = array![[0.0, 0.0], [10.0, 0.0], [5.0, 10.0]];
let predicted = km.predict(&new_points).unwrap();
println!("new-point labels: {:?}", predicted);
}
```
具体的簇编号是任意的：k-means 按发现顺序给簇编号，所以哪个数据簇成为 0 号取决于种子。但*结构*是固定的：
```text
labels: 每个数据簇的 5 个点共享同一个编号。3 个数据簇会拿到 3 个互不相同的编号
（0、1、2 的某种排列）。
centroids: 形状 (3, 2)。三行按某种顺序分别落在 (0,0)、(10,0)、(5,10) 附近约 0.1 的范围内。
inertia: 一个很小的正 f64（各点到质心距离的平方和）。
iters: 若干次，远低于 max_iterations。
new-point labels: 每个测试点都映射到它所在数据簇对应的那个编号。
```
簇编号取决于排列方式，所以永远不要用相等去比较 2 个分别拟合的模型的标签。应该改用排列不变的指标，比如 [5.3. 聚类指标](../Chapter-05/5.3._聚类指标.md) 里的调整兰德指数。
## 2.7.5. 如何选择 k
k-means 没法告诉你 *k* 是多少，得你自己给出。有 2 种方法能帮你缩小范围，RustyML 为两者都备好了原料。
**肘部法**把 inertia 对 *k* 画成曲线。inertia 会随着 *k* 增大而单调下降，因为质心越多，平方距离之和只会更小。到 *k = n* 时它会降到零。要找的是「肘部」，也就是下降变平的那个点。inertia 可以直接从 `get_inertia` 读出来。
**轮廓系数**更有决断力，因为它有一个内部最优值，而不是单调趋势。对每个点，它衡量的是这个点离自己所在的簇比离最近的其他簇近多少。它把这些值平均成 `[-1, 1]` 里的一个分数，越高越好。metrics 模块把它作为 [`silhouette_score`](../Chapter-05/5.3._聚类指标.md) 提供。它接收特征矩阵、标签和一个 [`DistanceCalculationMetric`](../Chapter-06/6.1._距离度量.md)。它只对 `2..=n-1` 个不同的簇有定义，所以没法给 `k = 1` 打分。
```rust
use rustyml::machine_learning::KMeans;
use rustyml::metrics::silhouette_score;
use rustyml::math::DistanceCalculationMetric;
use ndarray::{array, Array2};
fn main() {
let data: Array2<f64> = array![
[-0.05, 0.03], [ 0.04, -0.02], [ 0.01, 0.05], [-0.03, -0.04], [ 0.02, 0.01],
[ 9.95, 0.03], [10.04, -0.02], [10.01, 0.05], [ 9.97, -0.04], [10.02, 0.01],
[ 4.95, 10.03], [ 5.04, 9.98], [ 5.01, 10.05], [ 4.97, 9.96], [ 5.02, 10.01],
];
println!(" k inertia silhouette");
for k in 1..=5usize {
let mut km = KMeans::new(k, 300, 1e-4).unwrap().with_random_state(42);
let labels = km.fit_predict(&data).unwrap();
let inertia = km.get_inertia().unwrap();
if k >= 2 {
let s = silhouette_score(&data, &labels, DistanceCalculationMetric::Euclidean);
println!("{k:2} {inertia:8.4} {s:7.4}");
} else {
println!("{k:2} {inertia:8.4} (n/a)");
}
}
}
```
在这 3 个干净的数据簇上，inertia 从 `k = 1` 到 `k = 3` 骤降，然后趋平。轮廓系数在 `k = 3` 达到峰值。两者都指向真实结构。真实数据要浑浊得多，所以轮廓系数的内部极大值往往能给出更清晰的信号。完整的指标清单见 [5.3. 聚类指标](../Chapter-05/5.3._聚类指标.md)，其中还有 Davies-Bouldin 和 Calinski-Harabasz，能给出独立的第二意见。
## 2.7.6. 可复现的聚类：局部种子与全局种子
k-means++ 带随机性，所以不设种子的模型每次 `fit` 都会产出不同的划分。有 2 种办法能把它钉死，而且它们组合起来的行为是可预料的。
**局部种子**用 `with_random_state(seed)` 设定。它是自足的：直接为该估计器的初始化 RNG 播种。它无视任何全局状态，也绝不影响其他组件。
`n_init` 次重启各自的子种子，是通过一个固定的混合步骤从它派生出来的，而不是靠推进同一个共享 RNG。所以某次重启的初始化不会取决于之前几次重启消耗了多少随机数。用同一个局部种子构建、并在同一份数据上拟合的 2 个模型，会产出完全一致的质心、标签和 inertia。迭代次数也一样一致，精确到最后一位。
**全局种子**通过 `rustyml::random::set_global_seed(seed)` 设定，能只用一次调用，就固定此后在同一线程上构建并拟合的每一个*未设种子*的随机化组件。这包括 k-means、神经网络的初始化器、`train_test_split` 等等。这照搬了 Keras 的全局种子行为。未设种子的 k-means 模型会在 `fit` 时从全局流里取一个独立的子种子。要复现某次运行，就在重新拟合之前再设一遍全局种子。
这里有 2 点需要留意。全局种子是**线程局部**的，所以要在拟合模型的那个线程上设它。未设种子的组件还会按*拟合顺序*消费这个流，所以它们的可复现性对顺序敏感。`with_random_state` 种子能绕开这两个问题，这正是它适合用来钉死单个估计器的原因。
```rust
use rustyml::machine_learning::KMeans;
use rustyml::random::{set_global_seed, clear_global_seed};
use ndarray::{array, Array2};
fn main() {
let data: Array2<f64> = array![
[0.0, 0.0], [0.1, 0.0], [0.0, 0.1],
[10.0, 0.0], [10.1, 0.0], [10.0, 0.1],
[5.0, 10.0], [5.1, 10.0], [5.0, 10.1],
];
// 局部种子：结果完全一致，与任何全局状态无关。
let mut a = KMeans::new(3, 300, 1e-4).unwrap().with_random_state(42);
let mut b = KMeans::new(3, 300, 1e-4).unwrap().with_random_state(42);
a.fit(&data).unwrap();
b.fit(&data).unwrap();
assert_eq!(a.get_labels().unwrap(), b.get_labels().unwrap());
// 全局种子：每次拟合前重设一遍，就能复现未设种子的运行。
set_global_seed(7);
let mut c = KMeans::new(3, 300, 1e-4).unwrap();
c.fit(&data).unwrap();
let labels_c = c.get_labels().unwrap().clone();
set_global_seed(7);
let mut d = KMeans::new(3, 300, 1e-4).unwrap();
d.fit(&data).unwrap();
assert_eq!(&labels_c, d.get_labels().unwrap());
clear_global_seed();
}
```
这种确定性之所以成立，是因为并行算术从设计上就追求可复现，而不只是追求速度。具体做法见 [2.7.8](#278-并行与性能)。整个 crate 的完整播种模型，见 [7.1. 可复现性与随机种子](../Chapter-07/7.1._可复现性与随机种子.md)。
## 2.7.7. 失败模式与常见陷阱
**空簇。**当 Lloyd 的分配让某个质心一个点都没分到时，RustyML 不会丢掉这个簇。它也不会留下一个陈旧的质心。每一轮迭代都会给每个空簇重新播种，用的是离它自己所属质心最远的那个点，也就是对 inertia 贡献最大的那个点。这样模型始终保持恰好 `n_clusters` 个质心，并在下一轮把它从退化状态里推出来。所以 `get_centroids` 永远会返回 `n_clusters` 行。在病态数据上，比如大量重复、或者 `k` 接近 `n`，*标签*仍可能解析出少于 `k` 个不同的值。
**对特征尺度敏感。**k-means 最小化的是欧几里得距离。所以一个以千计的特征会压过一个以零点几计的特征，聚类实际上就把小尺度那个特征忽略了。它没有内置的标准化。只要你的特征处在不同量级，就先用 [4.2. 标准化与归一化](../Chapter-04/4.2._标准化与归一化.md) 里的工具做标准化。k-means 结果看起来不对，最常见的原因就是这一个。
**对离群点敏感。**质心只是普通均值，所以几个极端点就能把它从簇的稠密区域拽走，并抬高 inertia。数据里有离群点时，要么把它们裁掉，要么改用把它们当作噪声处理的基于密度的方法。[2.8. DBSCAN](./2.8._DBSCAN.md) 会显式地标出离群点。[2.9. MeanShift](./2.9._MeanShift.md) 无需预设 *k* 就能找到密度峰值。
**非球状簇。**k-means 把空间划分成质心周围的 Voronoi 元胞。所以它只能还原大致凸形、大小相近的簇。同心圆环或狭长的流形，无论 *k* 取多少都能难倒它。这种几何结构正是 DBSCAN 和 [2.12. t-SNE](./2.12._t-SNE.md) 里的流形方法要对付的。
**样本太少。**要求的簇数比点数还多，会返回 `Error::InvalidInput`。它不会悄悄截断 `k`。如果 `k` 取决于数据，就在拟合前先检查 `n_clusters <= data.nrows()`。
## 2.7.8. 并行与性能
分配步是开销的大头。每轮迭代它都要拿每个点去比对每个质心，并行也就落在这里。每轮迭代把所有点到质心的投影当成一次并行矩阵乘积算出来（见 [6.2. 矩阵乘法](../Chapter-06/6.2._矩阵乘法.md)）。再用一次简短的 arg-min 扫描找出每个点最近的质心。
质心更新会把每个点的那一行累加进它所属簇的累积和里，做法是一次**确定性分块归约**（见 [6.3. 并行归约](../Chapter-06/6.3._并行归约.md)）。k-means++ 初始化的距离计算也用同样的方式并行。
“确定性”这个词在这里很关键。上面每个阶段都只在超过一个标定过的工作量阈值时才并行。归约始终按固定的分块顺序求和，而不是按线程到达的顺序。正因如此，并行路径才会给出和串行路径一样的答案。也正因如此，带种子的拟合在同一台机器上才可复现，不管 rayon 用了多少线程。
浮点加法不满足结合律。按线程恰好完成的顺序求和，会把线程数泄漏进结果里。跨不同机器或不同构建目标时，算术后端最后一位上的微小差异仍有可能出现。但在同一台机器内，固定种子就给出固定答案。
输入较小时，一切都保持单线程，这样能避免在用不着的数据上白付 rayon 的协调开销。这里有 3 个阈值控制着切换点。arg-min 扫描那一步用的是 f64 扫描门限。累加那一步用的是 f64 求和门限。把每个质心除以其簇大小那一步，用的是 f64 cheap-map 门限。
扫描门限和求和门限都在 `rustyml::tuning::reduction` 里。cheap-map 门限在 `rustyml::tuning::elementwise` 里。如果你在聚类一些不寻常的形状、想挪动某个阈值，就调它们。多数用户从不碰它们。[7.3. 性能调优与并行](../Chapter-07/7.3._性能调优与并行.md) 讲解了这些旋钮。
开启 `show_progress` feature 构建时，`fit` 过程中会画一条实时的 inertia 与迭代进度条。这在大数据集上很趁手，feature 关掉时它不产生任何开销。
## 2.7.9. 保存与加载已训练的模型
`KMeans` 派生了 serde 的 `Serialize` 和 `Deserialize`。配套生成的 `save_to_path` 和 `load_from_path` 方法，会把整个模型持久化成一段紧凑的 postcard 二进制块，其中包括质心、标签、超参数和元数据。文件扩展名并不重要，因为格式始终是二进制 postcard。加载回来的模型与原模型预测完全一致，连质心的最后一位都不差。
引入重启机制时，序列化布局多了 `n_init` 这个字段。所以旧版本写出的二进制块加载不了。请重新拟合模型、再重新保存。
```rust
use rustyml::machine_learning::KMeans;
use ndarray::{array, Array2};
use std::fs::remove_file;
fn main() {
let data: Array2<f64> = array![
[0.0, 0.0], [0.1, 0.0], [0.0, 0.1],
[10.0, 0.0], [10.1, 0.0], [10.0, 0.1],
[5.0, 10.0], [5.1, 10.0], [5.0, 10.1],
];
let mut km = KMeans::new(3, 300, 1e-4).unwrap().with_random_state(42);
km.fit(&data).unwrap();
km.save_to_path("kmeans_model.bin").unwrap();
let loaded = KMeans::load_from_path("kmeans_model.bin").unwrap();
let original = km.predict(&data).unwrap();
let restored = loaded.predict(&data).unwrap();
assert_eq!(original, restored);
remove_file("kmeans_model.bin").unwrap();
}
```
I/O 和反序列化的问题会以 [`Error::Io`](../Chapter-01/1.6._错误处理.md) 的形式出现。加载时文件不存在是最常见的情形。持久化格式、版本兼容方面的考量，以及它如何与 crate 其余部分配合，见 [7.2. 深入模型持久化](../Chapter-07/7.2._深入模型持久化.md)。