用 Rust 从头实现机器学习算法·聚类②：DBSCAN 与谱聚类——月牙形数据的两位奇兵
用 Rust 从头实现机器学习算法·聚类②：DBSCAN 与谱聚类——月牙形数据的两位奇兵
k 均值与 GMM 的簇都是”凸的蛋形”——遇到月牙、环形这类纠缠结构就束手无策。本篇的两位主角从完全不同的假设出发：DBSCAN 用密度（簇是高密度点的连通区域），谱聚类用图论（把聚类变成图切分）。数据是交错的月牙——上一篇的椭圆簇在这里只会让两簇被切得支离破碎。
核心思想
DBSCAN（Density-Based Spatial Clustering）的两个超参数：（邻域半径）与 MinPts（成为核心点所需的最小邻居数）。点的三种身份：核心点（邻居 ≥ MinPts）、边界点（在核心点邻域内但自己不达标）、噪声点（谁也管不着）。算法从每个未访问的核心点出发，把密度可达的点层层扩展成簇——无需预设簇数，噪声自动现形，任意形状的簇都能长出来。
谱聚类三步走（清单 4.6）：① 构图——每个点连 个最近邻，边权随距离指数衰减；② 谱嵌入——算图拉普拉斯 的前 个特征向量（ 是度矩阵），把样本嵌入到 维谱空间；③ 收尾——在嵌入空间里跑 k 均值。直觉：切图（Normalized Cut）想让”簇内边权大、簇间边权小”，其松弛解恰好是拉普拉斯的低维特征向量——月牙两弧之间的”桥”很窄，一切就开。
数学：密度可达与拉普拉斯
DBSCAN 没有目标函数，只有可达性定义： 从 密度直达若 在核心点 的 邻域内；密度可达是直达的传递闭包。簇 = 互相密度可达的最大点集。这解释了它的性格：对异常值免疫（噪声不入簇）、对密度不均敏感（一个 管不住胖瘦两簇）。
谱聚类的对象是归一化拉普拉斯 。清单说它的本质是 Normalized Cut 的松弛解——这里给出操作版： 的最小特征向量（对应 的最大特征向量）把图”劈”成联系最弱的两半；前两个特征向量组成的嵌入矩阵按行归一化后，普通 k 均值就能在谱空间里把纠缠的簇分开。
Rust 实现
// src/main.rs（段一：月牙数据、k 均值对照、DBSCAN、谱聚类与主程序）
// 聚类②：DBSCAN 与谱聚类——月牙形数据的两位奇兵
// 单文件、仅标准库。谱聚类用幂迭代求特征向量；绘图器与系列前篇一致。
// ===================== xorshift64 + Box-Muller（沿用 07 篇） =====================
struct XorShift(u64);
impl XorShift {
fn new(seed: u64) -> Self {
assert!(seed != 0, "种子不能为 0");
XorShift(seed)
}
fn next_u64(&mut self) -> u64 {
let mut x = self.0;
x ^= x << 13;
x ^= x >> 7;
x ^= x << 17;
self.0 = x;
x
}
fn next_f64(&mut self) -> f64 {
(self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
}
fn next_gauss(&mut self) -> f64 {
let u1 = 1.0 - self.next_f64();
let u2 = self.next_f64();
(-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
}
}
// ===================== 数据：两个月牙 =====================
/// 上弧圆心 (0, 0)，下弧圆心 (1, 0.3)，加小噪声——交错月牙。
fn make_moons(rng: &mut XorShift, n: usize, sigma: f64) -> (Vec<(f64, f64)>, Vec<usize>) {
let mut xs = Vec::with_capacity(n * 2);
let mut ys = Vec::with_capacity(n * 2);
for _ in 0..n {
let t = 0.2 + (std::f64::consts::PI - 0.4) * rng.next_f64();
xs.push((t.cos() + sigma * rng.next_gauss(), t.sin() + sigma * rng.next_gauss()));
ys.push(0);
}
for _ in 0..n {
let t = 0.2 + (std::f64::consts::PI - 0.4) * rng.next_f64();
xs.push((1.0 - t.cos() + sigma * rng.next_gauss(), 0.3 - t.sin() + sigma * rng.next_gauss()));
ys.push(1);
}
(xs, ys)
}
fn dist2(a: (f64, f64), b: (f64, f64)) -> f64 {
(a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)
}
// ===================== k 均值（对比用，沿用 15 篇） =====================
fn kmeans(xs: &[(f64, f64)], k: usize, iters: usize, rng: &mut XorShift) -> Vec<usize> {
let mut centroids: Vec<(f64, f64)> = vec![xs[(rng.next_f64() * xs.len() as f64) as usize]];
while centroids.len() < k {
let d2: Vec<f64> = xs
.iter()
.map(|x| centroids.iter().map(|c| dist2(*x, *c)).fold(f64::INFINITY, f64::min))
.collect();
let total: f64 = d2.iter().sum();
let mut u = rng.next_f64() * total;
let mut pick = 0;
for (i, d) in d2.iter().enumerate() {
u -= d;
if u <= 0.0 {
pick = i;
break;
}
}
centroids.push(xs[pick]);
}
let mut assign = vec![0usize; xs.len()];
for _ in 0..iters {
let mut moved = false;
for (i, x) in xs.iter().enumerate() {
let a = (0..k)
.min_by(|&a, &b| dist2(*x, centroids[a]).partial_cmp(&dist2(*x, centroids[b])).unwrap())
.unwrap();
if assign[i] != a {
assign[i] = a;
moved = true;
}
}
let mut sum = vec![(0.0, 0.0); k];
let mut cnt = vec![0usize; k];
for (x, &a) in xs.iter().zip(assign.iter()) {
sum[a].0 += x.0;
sum[a].1 += x.1;
cnt[a] += 1;
}
for c in 0..k {
if cnt[c] > 0 {
centroids[c] = (sum[c].0 / cnt[c] as f64, sum[c].1 / cnt[c] as f64);
}
}
if !moved {
break;
}
}
assign
}
// ===================== DBSCAN（密度连通） =====================
/// 返回每个点的簇标签：0..k-1 为簇，-1 为噪声。
fn dbscan(xs: &[(f64, f64)], eps: f64, min_pts: usize) -> Vec<i32> {
let n = xs.len();
let labels = vec![-1i32; n];
let mut labels = labels;
let mut cluster = 0i32;
for i in 0..n {
if labels[i] != -1 {
continue; // 已处理
}
// ε 邻域
let neighbors: Vec<usize> = (0..n).filter(|&j| dist2(xs[i], xs[j]) <= eps * eps).collect();
if neighbors.len() < min_pts {
labels[i] = -1; // 暂记噪声，可能被后来的核心点收留
continue;
}
// 扩展簇：BFS
labels[i] = cluster;
let mut queue: Vec<usize> = neighbors;
while let Some(p) = queue.pop() {
if labels[p] == -1 {
labels[p] = cluster;
}
if labels[p] != cluster {
continue; // 已属其它簇
}
let p_neighbors: Vec<usize> = (0..n).filter(|&j| dist2(xs[p], xs[j]) <= eps * eps).collect();
if p_neighbors.len() >= min_pts {
for &q in &p_neighbors {
if labels[q] == -1 || labels[q] == -2 {
queue.push(q);
}
}
}
}
cluster += 1;
}
labels
}
// ===================== 谱聚类（kNN 图 + 幂迭代） =====================
/// 幂迭代求算子 matvec 的主特征向量（|λ| 最大）。
/// 每次迭代做与已有特征向量的正交化（deflation），用于取第 2、3 个。
fn power_iteration(
n: usize,
matvec: &dyn Fn(&[f64]) -> Vec<f64>,
orthogonal_to: &[Vec<f64>],
rng: &mut XorShift,
iters: usize,
) -> (f64, Vec<f64>) {
let mut v: Vec<f64> = (0..n).map(|_| rng.next_gauss()).collect();
let norm = |v: &[f64]| v.iter().map(|x| x * x).sum::<f64>().sqrt();
let mut l = v.iter().map(|x| x * x).sum::<f64>().sqrt();
for _ in 0..iters {
let mut y = matvec(&v);
// deflation：与已知特征向量正交化
for u in orthogonal_to {
let dot: f64 = y.iter().zip(u.iter()).map(|(a, b)| a * b).sum();
for k in 0..n {
y[k] -= dot * u[k];
}
}
let ny = norm(&y);
if ny < 1e-12 {
break;
}
for k in 0..n {
v[k] = y[k] / ny;
}
l = ny;
}
(l, v)
}
/// 谱聚类三步骤：构图 → 拉普拉斯谱嵌入 → k 均值。
/// 返回簇分配。
fn spectral_cluster(xs: &[(f64, f64)], k: usize, knn: usize, sigma: f64, rng: &mut XorShift) -> Vec<usize> {
let n = xs.len();
// 1) k 近邻图（双向合并）：权重 w_ij = exp(−‖xi−xj‖² / 2σ²)
let mut adj: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
for i in 0..n {
let mut order: Vec<usize> = (0..n).filter(|&j| j != i).collect();
order.sort_by(|&a, &b| dist2(xs[i], xs[a]).partial_cmp(&dist2(xs[i], xs[b])).unwrap());
for &j in order.iter().take(knn) {
let w = (-dist2(xs[i], xs[j]) / (2.0 * sigma * sigma)).exp();
adj[i].push((j, w));
adj[j].push((i, w)); // 对称化
}
}
// 度 d_i = Σ_j w_ij；sqrt_d = D^{-1/2}
let mut sqrt_d = vec![0.0; n];
for i in 0..n {
let s: f64 = adj[i].iter().map(|&(_, w)| w).sum();
sqrt_d[i] = s.sqrt().max(1e-12);
}
// 2) M = D^{-1/2} W D^{-1/2} 的稀疏矩阵-向量乘
let matvec = |v: &[f64]| -> Vec<f64> {
let mut y = vec![0.0; n];
for i in 0..n {
let mut s = 0.0;
for &(j, w) in &adj[i] {
s += w * v[j] / sqrt_d[j];
}
y[i] = s / sqrt_d[i];
}
y
};
// M 的主特征向量 = L_sym = I − M 的最小特征向量
let (_, v1) = power_iteration(n, &matvec, &[], rng, 200);
let (_, v2) = power_iteration(n, &matvec, &[v1.clone()], rng, 200);
// 3) 嵌入 U = [v1, v2]，行归一化后 k 均值
let mut embedded: Vec<(f64, f64)> = (0..n).map(|i| (v1[i], v2[i])).collect();
for e in &mut embedded {
let nrm = (e.0 * e.0 + e.1 * e.1).sqrt();
if nrm > 1e-12 {
e.0 /= nrm;
e.1 /= nrm;
}
}
kmeans(&embedded, k, 100, rng)
}
// ===================== 主程序 =====================
fn main() {
let out_dir = "../../../frontend/public/images/series/rust-ml-16-density-spectral";
std::fs::create_dir_all(out_dir).expect("创建输出目录失败");
// ---- 1. 月牙数据 ----
let mut rng = XorShift::new(42);
let (xs, truth) = make_moons(&mut rng, 100, 0.06);
println!("数据：{} 个点的交错月牙（σ = 0.06）", xs.len());
// ---- 2. 对照组：k 均值在月牙上的失败 ----
let km = kmeans(&xs, 2, 100, &mut rng);
let km_acc = cluster_acc(&km, &truth);
println!("\n[对照] k 均值（球形假设）：调整兰德后 acc = {:.3}", km_acc);
// ---- 3. DBSCAN：ε 扫描 ----
println!("\n[DBSCAN] MinPts = 5，ε 扫描：");
println!(" ε 簇数 噪声点");
let mut best_eps = (0.0f64, i32::MIN, usize::MAX);
for eps in [0.10, 0.15, 0.20, 0.25, 0.30] {
let lb = dbscan(&xs, eps, 5);
let clusters = lb.iter().copied().filter(|&l| l >= 0).max().unwrap_or(-1) + 1;
let noise = lb.iter().filter(|&&l| l == -1).count();
println!(" {:.2} {:>4} {:>5}", eps, clusters, noise);
if clusters == 2 && noise < best_eps.2 {
best_eps = (eps, clusters, noise);
}
}
let db = dbscan(&xs, best_eps.0, 5);
let db_acc = cluster_acc(&db.iter().map(|&l| if l < 0 { 2usize } else { l as usize }).collect::<Vec<_>>(), &truth);
println!(" 选 ε = {:.2}：2 簇 {} 噪声，acc = {:.3}", best_eps.0, best_eps.2, db_acc);
// ---- 4. 谱聚类 ----
let sp = spectral_cluster(&xs, 2, 10, 0.15, &mut rng);
println!("\n[谱聚类] kNN = 10，σ = 0.15：acc = {:.3}", cluster_acc(&sp, &truth));
// ---- 5. 三张结果图 ----
// k 均值失败图
let mut c = Canvas::new(560.0, 420.0, -2.2, 2.8, -1.6, 1.8);
c.axes("x1", "x2");
for (k, color) in [(0usize, PALETTE[1]), (1, PALETTE[3])] {
let pts: Vec<(f64, f64)> = xs.iter().zip(km.iter()).filter(|(_, &a)| a == k).map(|(x, _)| *x).collect();
c.dots(&pts, color, 3.0);
}
c.legend(&[("k 均值簇 0", PALETTE[1]), ("k 均值簇 1", PALETTE[3])]);
let p1 = format!("{out_dir}/kmeans-moons.svg");
c.save(&p1);
// DBSCAN 结果图（噪声黑色）
let mut c = Canvas::new(560.0, 420.0, -2.2, 2.8, -1.6, 1.8);
c.axes("x1", "x2");
for (k, color) in [(0i32, PALETTE[1]), (1, PALETTE[3])] {
let pts: Vec<(f64, f64)> = xs.iter().zip(db.iter()).filter(|(_, &l)| l == k).map(|(x, _)| *x).collect();
c.dots(&pts, color, 3.0);
}
let noise: Vec<(f64, f64)> = xs.iter().zip(db.iter()).filter(|(_, &l)| l == -1).map(|(x, _)| *x).collect();
c.dots(&noise, "#16161d", 2.5);
c.legend(&[("簇 0", PALETTE[1]), ("簇 1", PALETTE[3]), ("噪声", "#16161d")]);
let p2 = format!("{out_dir}/dbscan.svg");
c.save(&p2);
// 谱聚类结果图
let mut c = Canvas::new(560.0, 420.0, -2.2, 2.8, -1.6, 1.8);
c.axes("x1", "x2");
for (k, color) in [(0usize, PALETTE[1]), (1, PALETTE[3])] {
let pts: Vec<(f64, f64)> = xs.iter().zip(sp.iter()).filter(|(_, &a)| a == k).map(|(x, _)| *x).collect();
c.dots(&pts, color, 3.0);
}
c.legend(&[("谱簇 0", PALETTE[1]), ("谱簇 1", PALETTE[3])]);
let p3 = format!("{out_dir}/spectral.svg");
c.save(&p3);
println!("图已写入：{p1}、{p2}、{p3}");
}
/// 聚类准确率：簇标签可任意置换，取与真标签最大重合（k=2 时直接两种对拍）。
fn cluster_acc(assign: &[usize], truth: &[usize]) -> f64 {
let direct = assign.iter().zip(truth.iter()).filter(|(a, t)| a == t).count();
let flipped = assign.len() - direct;
direct.max(flipped) as f64 / assign.len() as f64
}
// ===================== 迷你 SVG 绘图器（与系列前篇一致） =====================实现要点：
dbscan是 BFS 的变体：核心点的 邻域整体入队，弹出时若仍是核心点则继续扩展——labels[q] == -1 || -2的判断收留边界点、跳过已入队的。- 谱聚类的
matvec闭包把 作用到向量上，不建 矩阵——邻接表稀疏存储，一次乘法 。 power_iteration是系列第一个迭代式特征求解器：反复做 收敛到主特征向量；取第二个时每次迭代先做与 的正交化（deflation），把已找到的”劈法”从算子里剥掉。2×2 用闭式解， 用幂迭代——工具随问题尺寸升级。
// src/main.rs 段二：迷你 SVG 绘图器（与系列前篇同一份实现）
use std::fmt::Write as _;
/// 系列默认配色：朱橙 / 青绿 / 蓝 / 琥珀（取自站点设计令牌）。
const PALETTE: [&str; 4] = ["#d6491f", "#0f766e", "#2563eb", "#d97706"];
/// 一张图：坐标映射 + 已累积的 SVG 元素。
struct Canvas {
w: f64,
h: f64,
xmin: f64,
xmax: f64,
ymin: f64,
ymax: f64,
pad_l: f64,
pad_r: f64,
pad_t: f64,
pad_b: f64,
body: String,
}
impl Canvas {
/// 按数据范围自动留 5% 边距建图；xmin == xmax 时手动摊开避免除零。
fn new(w: f64, h: f64, xmin: f64, xmax: f64, ymin: f64, ymax: f64) -> Canvas {
let (xmin, xmax) = spread(xmin, xmax);
let (ymin, ymax) = spread(ymin, ymax);
Canvas {
w,
h,
xmin,
xmax,
ymin,
ymax,
pad_l: 52.0,
pad_r: 16.0,
pad_t: 16.0,
pad_b: 40.0,
body: String::new(),
}
}
fn px(&self, x: f64) -> f64 {
self.pad_l + (x - self.xmin) / (self.xmax - self.xmin) * (self.w - self.pad_l - self.pad_r)
}
fn py(&self, y: f64) -> f64 {
self.h - self.pad_b - (y - self.ymin) / (self.ymax - self.ymin) * (self.h - self.pad_t - self.pad_b)
}
/// 折线，点需按 x 升序传入。
#[allow(dead_code)] // 本篇未用，保持与系列其他篇的绘图器逐字一致
fn polyline(&mut self, pts: &[(f64, f64)], color: &str, width: f64) {
let mut d = String::new();
for (i, (x, y)) in pts.iter().enumerate() {
let _ = write!(d, "{}{:.1},{:.1}", if i == 0 { "M" } else { "L" }, self.px(*x), self.py(*y));
}
let _ = write!(
self.body,
r##"<path d="{d}" fill="none" stroke="{color}" stroke-width="{width}" stroke-linejoin="round"/>"##
);
}
/// 散点。
fn dots(&mut self, pts: &[(f64, f64)], color: &str, r: f64) {
for (x, y) in pts {
let _ = write!(
self.body,
r##"<circle cx="{:.1}" cy="{:.1}" r="{r}" fill="{color}" fill-opacity="0.65"/>"##,
self.px(*x),
self.py(*y)
);
}
}
/// 坐标轴 + 虚线网格 + 刻度标签。
fn axes(&mut self, xlabel: &str, ylabel: &str) {
let x0 = self.px(self.xmin);
let x1 = self.px(self.xmax);
let y0 = self.py(self.ymin);
let y1 = self.py(self.ymax);
for t in nice_ticks(self.xmin, self.xmax, 6) {
let x = self.px(t);
let _ = write!(
self.body,
r##"<line x1="{x:.1}" y1="{y0:.1}" x2="{x:.1}" y2="{y1:.1}" stroke="#16161d" stroke-opacity="0.08"/>"##
);
let _ = write!(
self.body,
r##"<text x="{x:.1}" y="{:.1}" text-anchor="middle" fill="#16161d" fill-opacity="0.45">{t:.6}</text>"##,
y0 + 16.0,
t = trim(t)
);
}
for t in nice_ticks(self.ymin, self.ymax, 6) {
let y = self.py(t);
let _ = write!(
self.body,
r##"<line x1="{x0:.1}" y1="{y:.1}" x2="{x1:.1}" y2="{y:.1}" stroke="#16161d" stroke-opacity="0.08"/>"##
);
let _ = write!(
self.body,
r##"<text x="{:.1}" y="{:.1}" text-anchor="end" fill="#16161d" fill-opacity="0.45">{t:.6}</text>"##,
x0 - 6.0,
y + 4.0,
t = trim(t)
);
}
let _ = write!(
self.body,
r##"<path d="M{x0:.1},{y0:.1}H{x1:.1}M{x0:.1},{y0:.1}V{y1:.1}" fill="none" stroke="#16161d" stroke-opacity="0.6"/>"##
);
let _ = write!(
self.body,
r##"<text x="{:.1}" y="{:.1}" text-anchor="middle" fill="#16161d" fill-opacity="0.7">{xlabel}</text>"##,
(x0 + x1) / 2.0,
self.h - 8.0
);
let _ = write!(
self.body,
r##"<text x="14" y="{:.1}" text-anchor="middle" transform="rotate(-90 14 {:.1})" fill="#16161d" fill-opacity="0.7">{ylabel}</text>"##,
(y0 + y1) / 2.0,
(y0 + y1) / 2.0
);
}
/// 置信带：upper/lower 两条曲线围成的区域，半透明填充。
#[allow(dead_code)] // 部分篇章不使用，保持各篇绘图器逐字一致
fn band(&mut self, upper: &[(f64, f64)], lower: &[(f64, f64)], color: &str) {
let mut d = String::new();
for (i, (x, y)) in upper.iter().enumerate() {
let _ = write!(d, "{}{:.1},{:.1}", if i == 0 { "M" } else { "L" }, self.px(*x), self.py(*y));
}
for (x, y) in lower.iter().rev() {
let _ = write!(d, "L{:.1},{:.1}", self.px(*x), self.py(*y));
}
let _ = write!(
self.body,
r##"<path d="{d}Z" fill="{color}" fill-opacity="0.15" stroke="none"/>"##
);
}
/// 图例：右上角依次画「色线 + 文字」。
fn legend(&mut self, entries: &[(&str, &str)]) { let sample_w = 22.0;
let line_h = 18.0;
let x_text = self.w - self.pad_r - 78.0 + sample_w + 6.0;
for (i, (label, color)) in entries.iter().enumerate() {
let y = self.pad_t + 14.0 + i as f64 * line_h;
let x0 = self.w - self.pad_r - 78.0;
let _ = write!(
self.body,
r##"<line x1="{x0:.1}" y1="{y:.1}" x2="{:.1}" y2="{y:.1}" stroke="{color}" stroke-width="2.5"/>"##,
x0 + sample_w
);
let _ = write!(
self.body,
r##"<text x="{x_text:.1}" y="{:.1}" fill="#16161d" fill-opacity="0.7">{label}</text>"##,
y + 4.0
);
}
}
fn save(&self, path: &str) {
let svg = format!(
r##"<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {w} {h}" font-family="'JetBrains Mono',monospace" font-size="11"><rect width="{w}" height="{h}" fill="#fbfaf7"/>{body}</svg>"##,
w = self.w,
h = self.h,
body = self.body
);
std::fs::write(path, svg).expect("写 SVG 失败");
}
}
/// 把 fmt 出来的 f64 末尾零去掉，让刻度标签短一点。
fn trim(v: f64) -> String {
let s = format!("{v:.6}");
s.trim_end_matches('0').trim_end_matches('.').to_string()
}
/// 数据范围退化成一点时向两侧摊开。
fn spread(min: f64, max: f64) -> (f64, f64) {
if min == max {
(min - 1.0, max + 1.0)
} else {
let pad = (max - min) * 0.05;
(min - pad, max + pad)
}
}
/// “好看”的刻度：步长取 1/2/2.5/5 × 10^k。
fn nice_ticks(min: f64, max: f64, n: usize) -> Vec<f64> {
let raw = (max - min) / n.max(1) as f64;
if raw <= 0.0 || !raw.is_finite() {
return vec![min];
}
let exp = raw.abs().log10().floor() as i32;
let base = 10f64.powi(exp);
let frac = raw / base;
let step = if frac <= 1.0 {
base
} else if frac <= 2.0 {
2.0 * base
} else if frac <= 2.5 {
2.5 * base
} else if frac <= 5.0 {
5.0 * base
} else {
10.0 * base
};
let mut ticks = Vec::new();
let mut t = (min / step).ceil() * step;
while t <= max + 1e-9 {
ticks.push(if t.abs() < step * 1e-9 { 0.0 } else { t });
t += step;
}
ticks
}// src/main.rs 末尾：单元测试（cargo test 运行）
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn dbscan_splits_two_blobs() {
let mut rng = XorShift::new(7);
let mut xs = Vec::new();
for _ in 0..30 {
xs.push((rng.next_gauss(), rng.next_gauss()));
}
for _ in 0..30 {
xs.push((10.0 + rng.next_gauss(), 10.0 + rng.next_gauss()));
}
let lb = dbscan(&xs, 1.5, 4);
let clusters = lb.iter().copied().filter(|&l| l >= 0).max().unwrap() + 1;
assert_eq!(clusters, 2);
assert!(lb.iter().filter(|&&l| l == -1).count() == 0);
}
#[test]
fn dbscan_flags_pure_noise() {
let mut rng = XorShift::new(9);
let xs: Vec<(f64, f64)> = (0..50).map(|_| (100.0 * rng.next_f64(), 100.0 * rng.next_f64())).collect();
let lb = dbscan(&xs, 0.5, 4);
assert!(lb.iter().all(|&l| l == -1));
}
#[test]
fn power_iteration_recovers_dominant_eigenvector() {
// 对角阵 diag(5, 1)：主特征值 5，特征向量 e0
let matvec = |v: &[f64]| vec![5.0 * v[0], 1.0 * v[1]];
let (l, v) = power_iteration(2, &matvec, &[], &mut XorShift::new(3), 100);
assert!((l - 5.0).abs() < 1e-6);
assert!(v[0].abs() > 0.999);
}
#[test]
fn spectral_splits_moons() {
let mut rng = XorShift::new(42);
let (xs, truth) = make_moons(&mut rng, 100, 0.05);
let sp = spectral_cluster(&xs, 2, 10, 0.12, &mut XorShift::new(5));
assert!(cluster_acc(&sp, &truth) > 0.95, "谱聚类应分开月牙: {}", cluster_acc(&sp, &truth));
}
}Rust 语法角：从函数指针到 trait 对象——dyn Fn 传算子
13 篇的核函数用裸函数指针（fn(...)），本篇的 power_iteration 升级成了 trait 对象：参数类型 matvec: &dyn Fn(&[f64]) -> Vec<f64>。区别在携带状态：裸函数指针只是地址，没有”捕获环境”的能力；dyn Fn 是”任何实现了 Fn trait 的东西”的胖指针——闭包、函数指针都行，闭包还能捕获外围变量（matvec 闭包捕获了 adj 和 sqrt_d）。代价是一次间接跳转。选择规则：无状态用 fn，有状态用 &dyn Fn 或泛型 impl Fn。见《Rust 程序设计语言》ch17-02（trait 对象）。
运行结果
cargo test（4 个用例：DBSCAN 分两 blob、纯噪声全标记、幂迭代主特征值精度、谱聚类分开月牙）全部通过后，cargo run：
数据：200 个点的交错月牙（σ = 0.06）
[对照] k 均值（球形假设）：调整兰德后 acc = 0.835
[DBSCAN] MinPts = 5，ε 扫描：
ε 簇数 噪声点
0.10 12 29
0.15 4 5
0.20 2 0
0.25 2 0
0.30 2 0
选 ε = 0.20：2 簇 0 噪声，acc = 1.000
[谱聚类] kNN = 10，σ = 0.15：acc = 1.000
图已写入：../../../frontend/public/images/series/rust-ml-16-density-spectral/kmeans-moons.svg、../../../frontend/public/images/series/rust-ml-16-density-spectral/dbscan.svg、../../../frontend/public/images/series/rust-ml-16-density-spectral/spectral.svg三张图：
怎么读这些数字和图
- k 均值对照组 0.835 分：球形簇假设逼着它横切两弧——图里能看到左上和右下被生硬地切成两半。这就是 15 篇结论的延伸：蛋形假设搞不定非凸结构。
- ε 扫描是 DBSCAN 的性格写照：ε=0.10 时碎成 12 簇（噪声 29）——邻域太小，弧被切成段；ε=0.20 起稳定 2 簇 0 噪声；acc 1.000。ε 的唯一性区间很宽（0.2~0.3），但数据密度不均时这个窗口会消失（清单缺点）。
- DBSCAN 与谱聚类殊途同归：两图都干净地分开月牙——前者靠”弧内部点密、弧之间桥窄”，后者靠”图上两弧之间的边权总和最小”。同一结构，两种数学目光。
- 选型速查（清单 4.7 落地版）：
| 数据特点 | 推荐算法 |
|---|---|
| 球形簇、大样本、已知簇数 | k 均值 |
| 需要概率 / 软归属 | GMM |
| 任意形状、有噪声、密度均匀 | DBSCAN / OPTICS |
| 任意形状、中小样本、要谱系 | 层次聚类 / 谱聚类 |
| 图结构 / 非流形结构 | 谱聚类 |
优缺点与适用场景
（抄清单 4.5–4.6 原文）
- DBSCAN：自动发现任意形状簇、自动识别噪声、无需预设簇数；但对 和 MinPts 敏感、密度不均匀时效果差、高维失效。扩展：OPTICS（用可达距离排序解决全局 敏感）、HDBSCAN（层次 + 密度自适应）。
- 谱聚类：能发现任意形状簇（月牙、环形等非凸结构）、有坚实图论基础；但构图和 k 近邻参数敏感、相似度图计算 、同样需要预设簇数。
- 适用场景：DBSCAN——异常检测、空间数据；谱聚类——图像分割、社交网络社区发现、形状不规则的数据。
小结
聚类模块收官。四种算法给出四种”相似”的定义：球形距离（k 均值）、概率密度（GMM）、密度连通（DBSCAN）、图切分（谱聚类）——选算法就是选相似性。无监督的两大主题（降维、聚类）至此走完，系列只剩终篇：集成学习——把前面所有弱学习器装进一个投票箱，看”三个臭皮匠”如何炼成工业界的默认王者。