1. 聚类解决什么问题

给定没有标签的数据集 \(X=\{x_1,x_2,\ldots,x_n\},\quad x_i\in\mathbb R^d,\) 聚类希望找到一个分组函数 $z_i\in{1,\ldots,K}$，或者更一般的层次结构、密度区域与软归属概率，使得同一组内相似、不同组间有差异。这句话看似直观，却隐藏了三个必须先回答的问题：

  1. “相似”用什么度量？欧氏距离、余弦距离、马氏距离，还是图上的连通性？
  2. “组”是什么形状？球状、椭圆、任意密度连通，还是流形上的局部邻域？
  3. 结果如何验证？无监督场景没有唯一答案，指标只能回答某一种假设下的好坏。

因此聚类不是一个单一算法，而是一族带有不同归纳偏置（inductive bias）的模型。K-Means 假设簇近似球形且方差相近；GMM 假设数据由多个高斯分布混合生成；DBSCAN 假设簇是高密度区域并允许噪声；谱聚类则把数据看成相似度图上的划分问题。算法没有脱离假设的“最好”，只有和数据生成机制更匹配的选择1。

1.1 硬聚类、软聚类与重叠聚类

  • 硬聚类：每个样本只属于一个簇，例如 K-Means 的标签。
  • 软聚类：输出 $p(z_i=k\mid x_i)$，一个样本可同时对多个簇有较高概率，例如 GMM。
  • 重叠/模糊聚类：输出隶属度 $u_{ik}\in[0,1]$，且不一定满足概率归一，便于表达边界模糊的群组。
  • 层次聚类：不预先固定 $K$，输出一棵树（dendrogram），用户在不同高度切割得到不同粒度。
  • 密度聚类：簇由高密度区域定义，孤立点被标记为噪声。

1.2 统一视角：目标函数与归纳偏置

大多数算法都可写成“在假设空间中优化某个目标”：

\[\hat\theta=\arg\min_\theta \mathcal L(X;\theta) \quad\text{或}\quad \hat\theta=\arg\max_\theta \log p(X;\theta).\]

区别在于：

  • K-Means 最小化簇内平方误差（SSE）。
  • GMM 最大化观测数据似然。
  • 谱聚类最小化图割目标的松弛形式。
  • DBSCAN 不显式优化全局目标，而是按密度可达关系构造连通分量。
  • 深度聚类同时学习表示 $f_\phi(x)$ 与簇结构。

理解目标函数比记 API 更重要：当结果异常时，你可以回到目标函数判断是尺度、距离、初始化还是模型假设导致的。

2. 开始前的数学工具箱

2.1 距离与相似度

欧氏距离

\(d_2(x,y)=\sqrt{\sum_{j=1}^{d}(x_j-y_j)^2}.\) 平方欧氏距离去掉了开方，排序关系不变，适合 K-Means 与高斯模型。它对量纲极其敏感：年龄（0–100）和年收入（0–1,000,000）直接相加没有意义。

曼哈顿距离

\(d_1(x,y)=\sum_j|x_j-y_j|.\) 在高维稀疏或含异常值的数据上通常比欧氏距离更稳健，但对应的几何等值面是菱形，和 K-Means 的平方损失不匹配。

余弦相似度

\(\operatorname{cos}(x,y)=\frac{x^\top y}{\|x\|_2\|y\|_2}.\) 文本 TF-IDF、归一化 embedding 常用余弦距离 $1-\operatorname{cos}$。若向量方向比长度更重要，应先 L2 归一化再使用 K-Means；归一化后欧氏距离与余弦距离单调等价。

马氏距离

\(d_M(x,y)=\sqrt{(x-y)^\top\Sigma^{-1}(x-y)}.\) 它用协方差校正特征相关性。样本少、维度高时 $\Sigma$ 不可逆，需使用收缩估计 $\Sigma_\lambda=(1-\lambda)\Sigma+\lambda I$ 或 PCA 降维。

2.2 标准化、缺失值与异常值

常见预处理：

import numpy as np


def standardize(x: np.ndarray, eps: float = 1e-8) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """按列做 z-score 标准化。

    Args:
        x: 形状为 (n_samples, n_features) 的数值矩阵。
        eps: 防止常数列除零的稳定项。

    Returns:
        返回标准化矩阵、列均值和列标准差，便于线上复用训练期统计量。
    """
    mean = np.nanmean(x, axis=0)
    std = np.nanstd(x, axis=0)
    std = np.maximum(std, eps)
    return (x - mean) / std, mean, std

标准化统计量只能在训练集上计算，否则会产生数据泄漏。缺失值可以用中位数、KNN 或模型插补；对异常值，先区分“真正稀有样本”和“测量错误”，因为 DBSCAN 可能把前者正确标成噪声，而 K-Means 会被其质心拖动。

2.3 高维距离退化

维度增大时，最近邻和最远邻的距离相对差异会变小，称为距离集中（concentration）。直觉上，所有点都“差不多远”，基于距离的聚类会失去分辨率。应对方法包括：

  • 删除业务上无意义的特征；
  • PCA/随机投影/自编码器降维；
  • 对文本和图像使用预训练 embedding；
  • 采用余弦、相关系数或学习到的度量；
  • 在原空间与低维空间分别聚类并比较稳定性。

3. K-Means：最重要的基线

3.1 目标函数与 Lloyd 算法

帮助理解：K-Means 可以理解为“先猜 K 个代表点，再把每个样本交给最近的代表点，然后把代表点移动到组内平均位置”。反复执行后，得到的是硬标签：每个样本只能属于一个簇。

K-Means 最小化簇内平方和：

\[\min_{\{\mu_k\},z}\;J=\sum_{i=1}^n\|x_i-\mu_{z_i}\|_2^2.\]

Lloyd 算法交替执行两步：

  1. 分配（assignment）：$z_i\leftarrow\arg\min_k|x_i-\mu_k|^2$。
  2. 更新（update）：$\mu_k\leftarrow\frac{1}{|C_k|}\sum_{i:z_i=k}x_i$。

每一步都不会增加 $J$，因此有限步后收敛到局部最优；但问题本身是 NP-hard，局部最优不保证全局最优2。

3.2 NumPy 手写实现

import numpy as np


def kmeans_numpy(
    x: np.ndarray,
    n_clusters: int,
    n_init: int = 10,
    max_iter: int = 300,
    tol: float = 1e-4,
    seed: int = 0,
) -> tuple[np.ndarray, np.ndarray, float]:
    """使用 NumPy 实现多次随机初始化的 K-Means。

    Args:
        x: 输入数据，形状为 (n_samples, n_features)。
        n_clusters: 簇数量 K，必须小于等于样本数。
        n_init: 独立初始化次数，取目标函数最小的一次。
        max_iter: 单次运行的最大迭代次数。
        tol: 质心移动的平方范数阈值。
        seed: 随机数种子。

    Returns:
        labels: 每个样本的簇标签。
        centers: 形状为 (K, d) 的质心。
        inertia: 最终 SSE。

    Raises:
        ValueError: 输入维度、K 或 n_init 不合法时抛出。
    """
    if x.ndim != 2 or n_clusters < 1 or n_clusters > len(x) or n_init < 1:
        raise ValueError("x 必须为二维矩阵，且 1 <= n_clusters <= n_samples")
    rng = np.random.default_rng(seed)
    # 多次随机初始化，最后选择 SSE 最小的一次，降低局部最优风险。
    best = None
    for _ in range(n_init):
        centers = x[rng.choice(len(x), size=n_clusters, replace=False)].copy()
        for _ in range(max_iter):
            # 广播计算“每个样本到每个质心”的平方距离矩阵。
            distances = ((x[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
            labels = distances.argmin(axis=1)
            # 先复制旧质心，逐簇写入均值，避免更新顺序影响本轮分配。
            new_centers = centers.copy()
            for k in range(n_clusters):
                members = x[labels == k]
                if len(members) > 0:
                    new_centers[k] = members.mean(axis=0)
                else:
                    # 空簇重置为当前误差最大的样本，避免 NaN 传播。
                    worst = distances.min(axis=1).argmax()
                    new_centers[k] = x[worst]
            shift = ((new_centers - centers) ** 2).sum()
            centers = new_centers
            if shift <= tol:
                break
        # 质心收敛后重新分配一次，确保返回的标签与最终质心一致。
        distances = ((x[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
        labels = distances.argmin(axis=1)
        final_dist = distances[np.arange(len(x)), labels]
        inertia = float(final_dist.sum())
        if best is None or inertia < best[2]:
            best = (labels.copy(), centers.copy(), inertia)
    return best

复杂度约为 $O(nKdT)$，其中 $T$ 是迭代次数；内存中的距离矩阵是 $O(nK)$。当 $n$ 或 $K$ 很大时，应使用分块计算或 Mini-Batch K-Means。

3.3 K-Means++ 初始化

随机初始化可能把多个质心放在同一簇。K-Means++ 先随机选一个质心，之后按样本到最近质心的平方距离 $D(x)^2$ 作为抽样权重，重复选点。该策略在期望意义上给出了 $O(\log K)$ 近似保证，并显著减少坏局部最优3。


def kmeans_plus_plus(x: np.ndarray, k: int, rng: np.random.Generator) -> np.ndarray:
    """按 K-Means++ 分布选择初始质心。"""
    centers = [x[rng.integers(len(x))]]
    closest_sq = ((x - centers[0]) ** 2).sum(axis=1)
    for _ in range(1, k):
        total = closest_sq.sum()
        if total <= 0:
            centers.append(x[rng.integers(len(x))])
            continue
        index = rng.choice(len(x), p=closest_sq / total)
        centers.append(x[index])
        closest_sq = np.minimum(closest_sq, ((x - x[index]) ** 2).sum(axis=1))
    return np.asarray(centers)

工程默认应优先 k-means++ 加少量 n_init，并固定随机种子记录可复现结果。n_init 不是越大越好：它线性增加成本，应该和数据规模及结果不稳定程度匹配。

3.4 Mini-Batch K-Means 与 PyTorch Tensor

帮助理解：Mini-Batch K-Means 是 K-Means 的“省内存版”：每次只看一小批样本，马上更新质心。它速度快、适合大数据，但结果会有更多随机波动。

Mini-Batch K-Means 每次只取一个小批次，用指数或计数加权更新质心。它把每轮成本从 $O(nKd)$ 降到 $O(bKd)$，适合流式数据，但会引入随机噪声，需要学习率衰减和多轮遍历。

import torch


def minibatch_kmeans_tensor(
    x: torch.Tensor,
    k: int,
    batch_size: int = 1024,
    steps: int = 2000,
    seed: int = 0,
) -> tuple[torch.Tensor, torch.Tensor]:
    """在 PyTorch Tensor 上运行 Mini-Batch K-Means。

    Args:
        x: 浮点 Tensor，形状为 (n_samples, n_features)，可在 CPU 或 GPU 上。
        k: 簇数量。
        batch_size: 每步随机抽取的样本数。
        steps: 更新步数。
        seed: PyTorch 随机种子。

    Returns:
        返回整型标签和最终质心；标签在函数末尾对全量样本计算。
    """
    if x.ndim != 2 or k < 1 or k > x.shape[0]:
        raise ValueError("x 必须为二维 Tensor，且 k 合法")
    generator = torch.Generator(device=x.device).manual_seed(seed)
    initial = torch.randperm(x.shape[0], generator=generator, device=x.device)[:k]
    centers = x[initial].clone()
    counts = torch.zeros(k, device=x.device, dtype=x.dtype)
    for _ in range(steps):
        # 随机抽取一个小批次，只用它近似全量数据的更新方向。
        indices = torch.randint(x.shape[0], (min(batch_size, x.shape[0]),), generator=generator, device=x.device)
        batch = x[indices]
        # 计算批次样本到当前质心的平方距离。
        dist = torch.cdist(batch, centers).square()
        labels = dist.argmin(dim=1)
        for cluster in range(k):
            points = batch[labels == cluster]
            if points.numel() == 0:
                continue
            count = points.shape[0]
            old = counts[cluster]
            counts[cluster] += count
            rate = count / counts[cluster]
            centers[cluster] = (1 - rate) * centers[cluster] + rate * points.mean(dim=0)
    labels = torch.cdist(x, centers).square().argmin(dim=1)
    return labels, centers

3.5 K-Means 的结构性盲点

  • 非球形簇：两个同心环会被切成扇区。
  • 尺度不同：未标准化时大尺度特征支配距离。
  • 簇大小/密度不同：大簇可能吞并小簇。
  • 离群点敏感：平方损失放大极端值。
  • K 必须预先给定：需要指标、领域知识或模型选择。

不要只凭二维散点图判断“看起来分得开”。高维中可分离性、业务可解释性和线上稳定性必须单独验证。

4. 如何选择 K 与评估聚类

4.1 内部指标

SSE 与肘部法

\(\operatorname{SSE}(K)=\sum_i\|x_i-\mu_{z_i}\|^2.\) 随着 K 增大 SSE 必然下降，肘部点只能提供启发式建议，不能当成统计检验。

轮廓系数

对样本 $i$，令 $a(i)$ 为同簇平均距离，$b(i)$ 为最近其他簇平均距离：

\(s(i)=\frac{b(i)-a(i)}{\max(a(i),b(i))}.\) 取值范围为 [-1,1]。它偏好凸、均匀的簇，对密度变化和非球形结构不一定公平；实现时可直接使用 silhouette_score，但要确保距离度量与训练聚类一致45。

Calinski-Harabasz 与 Davies-Bouldin

CH 比较簇间离散度和簇内离散度；DB 计算每个簇与最相似簇的平均相似度，越低越好，具体公式和实现细节可参考评估工具文档5。指标之间冲突时，不要简单平均，应该回到业务目标与数据假设。

4.2 有标签时的外部指标

如果事后拥有真实标签，可以使用 Adjusted Rand Index（ARI）、Normalized Mutual Information（NMI）、V-measure、Fowlkes-Mallows 等。ARI 对随机一致性做校正，但标签编号置换不影响结果。外部指标只评价和给定标签的一致性，不代表聚类揭示了更有用的业务结构。

4.3 稳定性与可复现性

对样本重采样、随机种子、时间窗口和特征子集重复聚类，比较标签的一致性（ARI/NMI）或共识矩阵。一个轮廓系数高但重采样后完全改变的结果，通常不适合生产。应记录：数据版本、预处理统计量、距离、随机种子、算法版本、超参数、指标和可视化。

5. Gaussian Mixture Model：从硬标签到概率模型

5.1 生成式假设

帮助理解：GMM 不强迫样本做“非黑即白”的选择，而是认为数据由多个高斯分布混合生成，并给出“属于每个组件的概率”。椭圆形、重叠的簇通常更适合它。

GMM 假设每个样本先按混合权重 $\pi_k$ 选择组件，再从高斯分布采样：

\[p(x)=\sum_{k=1}^{K}\pi_k\mathcal N(x\mid\mu_k,\Sigma_k),\quad \pi_k\ge0,\ \sum_k\pi_k=1.\]

后验责任度（responsibility）为：

\[\gamma_{ik}=p(z_i=k\mid x_i)=\frac{\pi_k\mathcal N(x_i\mid\mu_k,\Sigma_k)}{\sum_j\pi_j\mathcal N(x_i\mid\mu_j,\Sigma_j)}.\]

5.2 EM 推导

引入隐变量 $z$ 后，直接最大化似然困难。EM 交替执行：

  • E 步：固定参数，计算 $\gamma_{ik}$。
  • M 步：固定责任度，更新 \(N_k=\sum_i\gamma_{ik},\quad \mu_k=\frac{1}{N_k}\sum_i\gamma_{ik}x_i,\) \(\Sigma_k=\frac{1}{N_k}\sum_i\gamma_{ik}(x_i-\mu_k)(x_i-\mu_k)^\top,\quad \pi_k=N_k/n.\)

EM 保证观测对数似然不下降，但仍可能陷入局部最优。协方差矩阵接近奇异时需加 reg_covar；高维全协方差参数量为 $O(Kd^2)$，可选 diagonal、tied 或 spherical 结构。

5.3 PyTorch Tensor 的对数域实现

import torch


def gmm_log_likelihood(x: torch.Tensor, means: torch.Tensor, log_vars: torch.Tensor, logits: torch.Tensor) -> torch.Tensor:
    """计算对角协方差 GMM 的逐样本对数似然。

    Args:
        x: (n, d) 输入样本。
        means: (k, d) 组件均值。
        log_vars: (k, d) 对角方差的对数，避免直接优化负方差。
        logits: (k,) 未归一化混合权重。

    Returns:
        (n,) 每个样本的 log p(x)，使用 logsumexp 保持数值稳定。
    """
    # 在对数方差参数化下保证方差为正，并避免极小值导致数值爆炸。
    var = log_vars.exp().clamp_min(1e-6)
    diff = x[:, None, :] - means[None, :, :]
    log_prob = -0.5 * ((diff.square() / var) + log_vars + torch.log(torch.tensor(2 * torch.pi, device=x.device))).sum(dim=-1)
    # logsumexp 聚合各组件，避免先转回概率后发生下溢。
    return torch.logsumexp(log_prob + torch.log_softmax(logits, dim=0), dim=1)

与 K-Means 相比，GMM 能表示椭圆簇并给出不确定性；但它对高斯假设、协方差估计和初始化更敏感。组件数可用 BIC/AIC 选择：

\(\operatorname{BIC}=-2\log L+p\log n,\) 其中 $p$ 是自由参数数量。BIC 惩罚复杂模型更强，适合在候选模型中做相对比较。

6. 层次聚类：从相似度树观察多尺度结构

6.1 凝聚式算法

凝聚式（agglomerative）从每个样本一个簇开始，每次合并最相近的两个簇，直到只剩一个簇。簇间距离由 linkage 定义：

  • single：两簇最近点距离，容易出现“链式效应”。
  • complete：两簇最远点距离，倾向紧凑簇。
  • average：所有跨簇点对平均距离。
  • Ward：合并后 SSE 增量最小，和欧氏空间 K-Means 目标一致。

帮助理解：层次聚类像是在“搭积木”：开始时每个样本都是独立小组，算法不断合并最相近的两组，最后形成一棵树。你可以在树的不同高度切一刀，得到不同数量的簇。

结果是一棵 dendrogram，切割高度即可得到不同 K。朴素实现需要维护 $O(n^2)$ 距离矩阵，样本数很大时内存先成为瓶颈。

import numpy as np


def ward_merge_cost(a: np.ndarray, b: np.ndarray) -> float:
    """返回 Ward 合并两个簇带来的 SSE 增量。"""
    return len(a) * len(b) / (len(a) + len(b)) * float(((a.mean(0) - b.mean(0)) ** 2).sum())

层次聚类的优点是无需预先固定 K、可解释多粒度关系；缺点是早期错误合并不可撤销，噪声会污染树，且大样本复杂度高。若需要大规模层次结构，可先用 BIRCH 或抽样，再对子簇中心做层次聚类。

7. 密度聚类：DBSCAN、OPTICS 与 HDBSCAN

7.1 DBSCAN 的定义

帮助理解：DBSCAN 先找“足够拥挤”的区域，再把相互连通的拥挤区域合并成簇；落在稀疏区域的点会被标记为噪声，因此不需要预先指定簇数量。

给定半径 $\varepsilon$ 和最小点数 min_samples：

  • 核心点：$\varepsilon$-邻域内至少有 min_samples 个点。
  • 边界点：自身不是核心点，但在某个核心点邻域内。
  • 噪声点：既非核心点也不从属于任何核心点。

从任一核心点出发，按密度可达关系扩展，即得到一个簇。DBSCAN 不需要 K，能发现任意形状并显式识别噪声6。

$\varepsilon$ 的选择可用 k-distance 图：计算每个点到第 min_samples 个近邻的距离并排序，在曲线拐点附近取值。特征尺度和距离度量必须先处理好，否则参数没有可迁移性。

7.2 OPTICS 与 HDBSCAN

DBSCAN 使用单一全局 $\varepsilon$，在不同密度的簇上表现受限。OPTICS 输出有序样本及可达距离，从一次运行中得到多个密度层次7。HDBSCAN 对层次密度树做凝聚和稳定性选择，能自动选择簇并为样本给出成员概率，通常比 DBSCAN 更适合密度差异明显的数据8。

实战建议：

  • 低维、噪声明显且密度近似均匀：DBSCAN。
  • 想探索多个密度阈值：OPTICS。
  • 簇密度不同、需要软成员概率：HDBSCAN。
  • 高维原始特征：先学习或压缩表示，再做密度聚类；直接在几百维上用欧氏距离通常不可靠。

8. 谱聚类：把聚类变成图切分

8.1 相似度图与拉普拉斯矩阵

帮助理解：谱聚类把样本看成一张图：相似的样本之间连一条权重较大的边。它先在图上找容易切开的区域，再对新的坐标做 K-Means，所以能处理月牙、同心圆等非凸形状。

构造相似度矩阵 $W$，例如 k-nearest-neighbor 图或高斯核：

\[W_{ij}=\exp(-\|x_i-x_j\|^2/(2\sigma^2)).\]

度矩阵 $D_{ii}=\sum_jW_{ij}$，未归一化拉普拉斯为 $L=D-W$，归一化版本为 $L_{sym}=I-D^{-1/2}WD^{-1/2}$。取最小的前 K 个特征向量组成矩阵 $U$，对 U 的行再运行 K-Means。这个“先嵌入后聚类”能把同心圆、月牙等非凸簇变成线性可分结构9。

8.2 PyTorch 计算特征向量

import torch


def spectral_embedding(x: torch.Tensor, k_neighbors: int = 10, n_components: int = 2) -> torch.Tensor:
    """构造 kNN 图并返回归一化拉普拉斯的低频特征向量。"""
    # 先得到全体样本的两两距离；大数据时应改用稀疏近邻搜索。
    distances = torch.cdist(x, x).square()
    knn = distances.topk(k_neighbors + 1, largest=False).indices[:, 1:]
    # 用邻居权重填充相似度图，其余边保持为 0。
    w = torch.zeros_like(distances)
    sigma = distances.gather(1, knn).median().clamp_min(1e-6)
    weights = torch.exp(-distances.gather(1, knn) / (2 * sigma))
    rows = torch.arange(x.shape[0], device=x.device)[:, None].expand_as(knn)
    w[rows, knn] = weights
    w = torch.maximum(w, w.T)
    degree = w.sum(dim=1).clamp_min(1e-8)
    inv_sqrt = degree.rsqrt()
    laplacian = torch.eye(x.shape[0], device=x.device) - inv_sqrt[:, None] * w * inv_sqrt[None, :]
    # 取拉普拉斯矩阵最小特征值对应的低频方向作为新表示。
    eigenvalues, eigenvectors = torch.linalg.eigh(laplacian)
    return eigenvectors[:, :n_components]

复杂度主要来自 $n\times n$ 相似度矩阵和特征分解，内存为 $O(n^2)$，不适合超大样本。可用 Nyström 近似、稀疏图、近似近邻和随机特征分解扩展规模。谱聚类对图构造超参数敏感，kNN 太小会断图，太大则抹平局部结构。

9. Mean Shift、Affinity Propagation 与 BIRCH

9.1 Mean Shift

帮助理解：Mean Shift 不先猜簇数，而是把每个点反复推向“周围最拥挤的地方”。最后落到同一个密度峰的点被分到一起；带宽越大，看到的结构越粗。

Mean Shift 在核密度估计的梯度方向上迭代：

\(m(x)=\frac{\sum_i K((x-x_i)/h)x_i}{\sum_i K((x-x_i)/h)}-x.\) 所有点向密度峰移动，收敛到同一峰的点归为一簇。它无需指定 K，能发现任意形状，但每次迭代都要访问大量样本，带宽 $h$ 决定簇粒度，样本量大时成本很高。

9.2 Affinity Propagation

算法在样本间传递 responsibility 和 availability 两类消息，自动选出 exemplars（代表样本）10。相似度矩阵为 $O(n^2)$，不适合大数据；preference 参数控制簇数，通常需要通过分位数或验证集调节。

9.3 BIRCH

BIRCH 用聚类特征（CF）树压缩大规模数据：每个节点保存样本数 N、线性和 LS、平方和 SS，从而可快速计算质心和半径11。先构建 CF 树，再对叶节点子簇运行 K-Means 或层次聚类。它适合内存受限的流式场景，但依赖欧氏空间，阈值半径对结果影响很大。

10. 模糊 C 均值与不确定性

帮助理解：模糊 C 均值允许一个样本同时属于多个簇，例如“更像 A，也有一点像 B”。隶属度越平均，说明样本越接近边界，决策时应保留不确定性。

模糊 C 均值把硬标签放宽为隶属度矩阵 U，目标为：

\(J_m=\sum_{i=1}^n\sum_{k=1}^K u_{ik}^m\|x_i-c_k\|^2, \quad \sum_k u_{ik}=1.\) 更新公式：

\[u_{ik}=\left[\sum_j\left(\frac{\|x_i-c_k\|}{\|x_i-c_j\|}\right)^{2/(m-1)}\right]^{-1}.\]

m 越大，边界越模糊。它适合一个样本可能同时属于多个状态的场景，如医学表型，但仍继承球状簇和欧氏距离假设，并且对异常点敏感。生产中应同时报告最大隶属度与熵：

\(H_i=-\sum_k u_{ik}\log u_{ik},\) 高熵样本应进入人工复核或后续决策的“不确定”分支。

11. 深度聚类：同时学习表示与簇

帮助理解：深度聚类先学习“什么样的表示更适合分组”，再在表示空间里聚类。它能处理图像、文本等原始距离不代表语义的高维数据，但训练不稳定和塌缩风险也更高。

经典算法通常在固定特征空间中聚类。图像、语音、文本的原始空间距离未必表达语义，因此深度聚类把编码器 $z=f_\phi(x)$ 与聚类目标联合优化。

11.1 自编码器 + K-Means 基线

先训练自编码器重建输入，再取瓶颈向量做 K-Means：

import torch
from torch import nn


class AutoEncoder(nn.Module):
    """用于聚类前表示学习的最小自编码器。"""

    def __init__(self, input_dim: int, latent_dim: int) -> None:
        super().__init__()
        self.encoder = nn.Sequential(nn.Linear(input_dim, 128), nn.ReLU(), nn.Linear(128, latent_dim))
        self.decoder = nn.Sequential(nn.Linear(latent_dim, 128), nn.ReLU(), nn.Linear(128, input_dim))

    def forward(self, x: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor]:
        """返回潜变量和重建结果。"""
        # 编码器压缩输入，潜变量将作为后续聚类的特征。
        latent = self.encoder(x)
        return latent, self.decoder(latent)

只优化重建损失会保留与聚类无关的细节；只优化聚类损失又可能学到退化表示。因此通常采用预训练 + 聚类微调，并监控重建、聚类稳定性和下游任务三类信号。

11.2 DEC：自训练的软分配

Deep Embedded Clustering（DEC）在潜空间引入 Student-t 核得到软分配：

\[q_{ij}=\frac{(1+\|z_i-\mu_j\|^2/\alpha)^{-(\alpha+1)/2}}{\sum_{j'}(1+\|z_i-\mu_{j'}\|^2/\alpha)^{-(\alpha+1)/2}}.\]

再构造强化高置信度样本的目标分布：

\(p_{ij}=\frac{q_{ij}^2/f_j}{\sum_{j'}q_{ij'}^2/f_{j'}},\quad f_j=\sum_iq_{ij},\) 最小化 KL 散度 $\sum_i\sum_j p_{ij}\log(p_{ij}/q_{ij})$12。平方会放大高置信度分配，但也可能造成确认偏差；应使用目标分布更新间隔、早停和稳定性检查。


def dec_kl(q: torch.Tensor, eps: float = 1e-8) -> torch.Tensor:
    """根据 DEC 公式计算 KL(p || q) 聚类损失。"""
    # 统计每个簇当前吸收了多少样本，用于校正大簇偏置。
    frequency = q.sum(dim=0, keepdim=True)
    # 平方会强化高置信度分配，再按簇频率做归一化。
    target = (q.square() / frequency.clamp_min(eps))
    target = target / target.sum(dim=1, keepdim=True).clamp_min(eps)
    return (target * (target.add(eps).log() - q.clamp_min(eps).log())).sum(dim=1).mean()

11.3 对比学习与深度聚类的现代组合

近年的方法通常把以下组件组合起来：

  1. 增强一致性：同一样本的两种增强应得到相近表示；不同样本通过 InfoNCE 或 supervised contrastive loss 分离。
  2. 原型/中心学习：维护可学习原型，使用 Sinkhorn-Knopp 或均衡约束避免所有样本塌缩到一个簇。
  3. 伪标签自训练：高置信度样本生成伪标签，低置信度样本延迟决策。
  4. 图结构正则：在 kNN 图上约束邻居表示或标签一致。
  5. 最优传输（OT）：给定簇容量先验，求样本到原型的软匹配，缓解簇不平衡与塌缩。

一个可训练的联合损失可以写成：

\[\mathcal L=\lambda_{rec}\mathcal L_{rec}+\lambda_{con}\mathcal L_{InfoNCE}+\lambda_{clu}\mathcal L_{KL}+\lambda_{bal}\mathcal L_{balance}.\]

权重不是越多越好；先固定表示学习和聚类的最小闭环，再逐项加入正则并做消融。对于文本/图像，优先使用领域预训练模型的 embedding，再比较简单 K-Means、GMM、HDBSCAN 与深度微调，深度方法只有在表示确实不足时才值得承担复杂度。

11.4 现代深度聚类方法谱系

“深度聚类”不是单一模型，而是几条逐渐融合的路线：

  • 生成式路线：VaDE 将 VAE 的先验设为高斯混合，在重建与 KL 正则之间学习可采样的簇结构；适合需要生成和不确定性估计的场景，但对先验和解码器容量敏感。
  • 重建增强路线：IDEC 在 DEC 的 KL 目标之外保留重建损失，减少微调阶段表示漂移；当原始特征细节对簇有用时通常更稳。
  • 自监督原型路线：DeepCluster/DeepCluster-v2 交替执行表示网络训练与 K-Means 伪标签更新；SwAV 用在线聚类和最优传输在不同增强视图间交换分配，避免显式构造所有负样本。
  • 邻域一致路线：SCAN 先用自监督表示构建近邻，再最大化邻居预测一致性，同时用熵正则避免塌缩；适合类别边界由局部邻域决定的图像数据。
  • 图融合路线：SDCN、DCRN 等把自编码器与 GCN/图正则联合，利用样本关系传播聚类信号；必须防止图构造中的标签泄漏和过平滑。

这些方法的共同难点是簇塌缩、伪标签确认偏差和目标冲突。实验时至少记录每个 epoch 的簇占比、潜变量协方差的有效秩、原型间最小距离和不同随机种子的 ARI；只看训练损失无法发现退化解。现代模型带来的提升通常首先来自更好的表示和增强，而不是更复杂的聚类头，故应把“冻结预训练 embedding + 简单聚类”作为不可省略的对照组。

11.5 图神经网络与多视图聚类

当样本天然构成图（社交网络、分子、知识图谱）时，可用 GCN/GAT 编码邻居信息，再在节点 embedding 上聚类。多视图聚类则为同一对象的文本、图像、行为分别编码，通过协同正则或共享潜空间融合。关键风险是图泄漏：如果边包含未来信息，离线聚类会高估效果；时间切分和边遮蔽必须与下游任务一致。

12. 统一的工程选型流程

12.1 先做数据诊断

  1. 统计样本数、维度、稀疏度、缺失率、异常值比例。
  2. 明确特征的尺度、单位、时间切分和业务语义。
  3. 抽样计算距离分布，比较欧氏、余弦、曼哈顿的集中程度。
  4. 用 PCA/UMAP 仅作可视化，不把二维图当成真实结构证明。

12.2 建立可比较的基线

推荐顺序：

  • 标准化 + K-Means++（多个 K）。
  • 标准化 + GMM（多个协方差类型，用 BIC 比较）。
  • kNN 图 + 谱聚类（样本量允许时）。
  • PCA/embedding + HDBSCAN（需要噪声与任意形状时）。

所有候选使用同一数据切分、同一预处理和同一评估协议。不要拿原空间 K-Means 的轮廓系数和 embedding 空间 HDBSCAN 的业务指标直接混为一谈。

12.3 选择算法的快速决策表

数据特征/目标 首选 备选与注意事项
大规模、球状、需快速分群 Mini-Batch K-Means BIRCH；先标准化与抽样估计 K
椭圆簇、需要概率与不确定性 GMM 协方差正则、BIC、组件退化检查
任意形状且有噪声 DBSCAN/HDBSCAN 先处理尺度；检查密度差异
多尺度层次关系 凝聚式层次聚类 Ward 适合欧氏；大样本先压缩
非凸簇、相似度图可靠 谱聚类 稀疏图、近似特征分解
无法给 K、密度峰明显 Mean Shift 带宽决定粒度，成本较高
高维语义数据 预训练 embedding + K-Means/HDBSCAN 深度微调需防塌缩与伪标签偏差
图或多视图对象 GNN/多视图深度聚类 防止时间泄漏，做视图消融

13. 常见失败模式与排查顺序

  1. 所有点挤在一个簇：检查特征是否未标准化、距离是否溢出、深度模型是否塌缩；查看每簇样本数和中心间距离。
  2. K 改一点结果就完全不同：增加初始化次数，做重采样稳定性；确认数据中是否没有清晰簇。
  3. DBSCAN 全是噪声：检查 $\varepsilon$ 与尺度，绘制 k-distance；高维先降维或换度量。
  4. GMM 出现极大似然但簇无意义：协方差奇异导致单点组件；增大正则、限制最小方差或改用 diagonal/tied。
  5. 深度聚类训练损失下降但标签退化：监控簇占比、表示方差、原型间距离；加入均衡约束和早停。
  6. 指标很好但业务不可用：轮廓系数只反映几何分离，必须做画像、稳定性、下游 uplift 和人工抽检。
  7. 线上漂移：记录新样本到中心/密度模型的距离分布、簇占比、未知/噪声比例；设定重训和回滚门槛。

14. 从实验到生产的最小闭环

一个可审计的聚类流水线至少包含：

数据快照 -> 训练期预处理 -> 表示学习/降维 -> 聚类模型 -> 版本化标签
        -> 内部指标与稳定性 -> 业务画像与抽检 -> 线上漂移监控 -> 重训/回滚

保存对象不仅是模型参数，还包括特征顺序、均值方差、距离定义、K 或密度参数、随机种子、训练数据指纹和评估报告。标签编号本身没有语义，模型更新后簇 0 不一定仍然代表原来的群体；需要用中心相似度、最优匹配或业务规则对齐新旧簇。

在推荐、风控、营销等高风险场景，聚类不应直接成为自动决策唯一依据。应把簇作为特征或候选分层，并检查不同群体的误差、覆盖率和资源分配是否产生不公平结果。无监督不等于无责任：没有标签只意味着验证更困难。

15. 一页纸总结

  • 聚类结果由距离/相似度 + 模型假设 + 超参数 + 随机性共同决定。
  • K-Means 是最好的第一基线，但不是通用答案；K-Means++、多次初始化和标准化是最低配置。
  • GMM 提供概率与椭圆簇，DBSCAN/HDBSCAN 处理噪声与任意形状，层次聚类揭示多尺度结构，谱聚类利用图关系。
  • 深度聚类的关键不在“网络更深”，而在于学习一个适合聚类的表示，并防止塌缩、确认偏差和数据泄漏。
  • 轮廓系数、BIC、ARI 只能回答局部问题；稳定性、业务解释、下游收益和线上漂移才决定能否落地。
  • 最可靠的工程路径是：先诊断数据，再建立简单基线，逐项增加模型能力，每一步都保留可复现的验证证据。

参考文献

  1. Rui Xu, Donald Wunsch. Survey of Clustering Algorithms. IEEE Transactions on Neural Networks, 2005. ↩

  2. Stuart Lloyd. Least Squares Quantization in PCM. IEEE Transactions on Information Theory, 1982. ↩

  3. David Arthur, Sergei Vassilvitskii. k-means++: The Advantages of Careful Seeding. SODA, 2007. ↩

  4. Jean-Bastien Rousseeuw. Silhouettes: A Graphical Aid to the Interpretation and Validation of Cluster Analysis. Journal of Computational and Applied Mathematics, 1987. ↩

  5. MathWorks/Scikit-learn documentation. Clustering performance evaluation and estimator guides. 2024. ↩ ↩2

  6. Martin Ester, Hans-Peter Kriegel, Jörg Sander, Xiaowei Xu. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. KDD, 1996. ↩

  7. Mihael Ankerst, Markus Breunig, Jörg Sander, Hans-Peter Kriegel. OPTICS: Ordering Points To Identify the Clustering Structure. SIGMOD, 1999. ↩

  8. Ricardo J. G. Campello, Davoud Moulavi, Jörg Sander. Density-Based Clustering Based on Hierarchical Density Estimates. PAKDD, 2013. ↩

  9. Andrew Y. Ng, Michael I. Jordan, Yair Weiss. On Spectral Clustering: Analysis and an Algorithm. NIPS, 2002. ↩

  10. Brendan J. Frey, Delbert Dueck. Clustering by Passing Messages Between Data Points. Science, 2007. ↩

  11. Tian Zhang, Raghu Ramakrishnan, Miron Livny. BIRCH: An Efficient Data Clustering Method for Very Large Databases. SIGMOD, 1996. ↩

  12. Junjie Xie, Ross Girshick, Ali Farhadi. Unsupervised Deep Embedding for Clustering Analysis. ICML, 2016. ↩