Sklearn-源码解析-书-v1-0-十-

Sklearn 源码解析（书）v1.0（十）

概述：predict 根据策略生成预测：most_frequent/prior 返回加权众数（argmax）；stratified 先调用 predict_proba 采样再取 argmax，保证两者一致；uniform 均匀随机采样；constant 广播用户常量。稀疏输出走 _random_choice_csc 生成 CSC 矩阵，稠密输出用 np.tile/np.vstack 广播。

策略行为对比：

下表总结了五种策略在预测行为、概率输出、随机性及样本权重影响上的差异，便于选择合适的基线策略：

| 策略 | predict 行为 | predict_proba 行为 | 随机性 | 样本权重影响 |

|------|----------------|----------------------|--------|--------------|

| most_frequent | 返回加权众数类别 | 对应类别概率 1.0，其余 0 | 无 | 是（影响众数判定） |

| prior | 返回加权众数类别 | 返回加权类别先验分布 | 无 | 是（影响先验分布） |

| stratified | 从先验分布采样类别 | 从先验分布多项式采样 one-hot | 有 | 是（影响采样分布） |

| uniform | 均匀随机采样类别 | 均匀分布 1/n_classes | 有 | 否 |

| constant | 返回用户指定常量 | 对应类别概率 1.0 | 无 | 否 |

源码路径：sklearn/dummy.py - DummyClassifier.predict_proba（第 242-310 行）

    def predict_proba(self, X):
        """Return probability estimates for the test vectors X."""
        check_is_fitted(self)

        n_samples = _num_samples(X)
        rs = check_random_state(self.random_state)

        n_classes_ = self.n_classes_
        classes_ = self.classes_
        class_prior_ = self.class_prior_
        constant = self.constant
        if self.n_outputs_ == 1:
            n_classes_ = [n_classes_]
            classes_ = [classes_]
            class_prior_ = [class_prior_]
            constant = [constant]

        P = []
        for k in range(self.n_outputs_):
            if self._strategy == "most_frequent":
                ind = class_prior_[k].argmax()
                out = np.zeros((n_samples, n_classes_[k]), dtype=np.float64)
                out[:, ind] = 1.0
            elif self._strategy == "prior":
                out = np.ones((n_samples, 1)) * class_prior_[k]  # 广播到所有样本

            elif self._strategy == "stratified":
                out = rs.multinomial(1, class_prior_[k], size=n_samples)  # 多项式采样
                out = out.astype(np.float64)

            elif self._strategy == "uniform":
                out = np.ones((n_samples, n_classes_[k]), dtype=np.float64)
                out /= n_classes_[k]

            elif self._strategy == "constant":
                ind = np.where(classes_[k] == constant[k])
                out = np.zeros((n_samples, n_classes_[k]), dtype=np.float64)
                out[:, ind] = 1.0

            P.append(out)

        if self.n_outputs_ == 1:
            P = P[0]

        return P

概述：生成各策略的概率输出：most_frequent 产生 one-hot；prior 广播先验分布到所有样本；stratified 用 multinomial(1, p) 采样 one-hot；uniform 生成均匀分布；constant 产生 one-hot。核心洞察是 most_frequent 和 prior 的 predict 相同（都取众数），但 predict_proba 截然不同：前者是 one-hot，后者是真实先验分布，这对校准评估至关重要。stratified 的 predict_proba 采样 one-hot，predict 取其 argmax，保证二者一致。

架构图：DummyClassifier.predict_proba 策略分发

graph TD A[predict_proba(X)] --> B[check_is_fitted] B --> C[_num_samples 获取 n_samples] C --> D[check_random_state] D --> E[单输出时包装为列表统一处理] E --> F[遍历 n_outputs_] F --> G{_strategy?} G -- most_frequent --> H[class_prior_.argmax 得到众数索引] H --> I[构造 one-hot: out[:, ind] = 1.0] G -- prior --> J[out = ones * class_prior_ 广播先验分布] G -- stratified --> K[rs.multinomial(1, class_prior_) 多项式采样] K --> L[astype float64] G -- uniform --> M[out = ones / n_classes_ 均匀分布] G -- constant --> N[np.where 定位常量索引] N --> O[构造 one-hot: out[:, ind] = 1.0] I --> P[加入 P 列表] J --> P L --> P M --> P O --> P P --> Q{单输出?} Q -- 是 --> R[返回 P[0]] Q -- 否 --> S[返回列表 P]

源码路径：sklearn/dummy.py - DummyClassifier.predict_log_proba（第 312-330 行）

    def predict_log_proba(self, X):
        """
        Return log probability estimates for the test vectors X.

        Parameters
        ----------
        X : {array-like, object with finite length or shape}
            Training data.

        Returns
        -------
        P : ndarray of shape (n_samples, n_classes) or list of such arrays
            Returns the log probability of the sample for each class in
            the model, where classes are ordered arithmetically for each
            output.
        """
        proba = self.predict_proba(X)
        if self.n_outputs_ == 1:
            return np.log(proba)
        else:
            return [np.log(p) for p in proba]

概述：直接对 predict_proba 结果取对数。单输出返回 np.log(proba)，多输出返回列表推导式 [np.log(p) for p in proba]。注意 most_frequent 和 constant 策略下 one-hot 概率中含 0，取对数后会产生 -inf，符合数学定义但需使用者注意处理。

源码路径：sklearn/dummy.py - DummyClassifier.score（第 360-390 行）

    def score(self, X, y, sample_weight=None):
        """Return the mean accuracy on the given test data and labels.

        In multi-label classification, this is the subset accuracy
        which is a harsh metric since you require for each sample that
        each label set be correctly predicted.

        Parameters
        ----------
        X : None or array-like of shape (n_samples, n_features)
            Test samples. Passing None as test samples gives the same result
            as passing real test samples, since DummyClassifier
            operates independently of the sampled observations.

        y : array-like of shape (n_samples,) or (n_samples, n_outputs)
            True labels for X.

        sample_weight : array-like of shape (n_samples,), default=None
            Sample weights.

        Returns
        -------
        score : float
            Mean accuracy of self.predict(X) w.r.t. y.
        """
        if X is None:
            X = np.zeros(shape=(len(y), 1))
        return super().score(X, y, sample_weight)

概述：实现分类器评分逻辑。关键设计是允许 X=None：因为预测完全不依赖输入特征，传入 None 等价于传入真实数据，内部构造全零矩阵仅用于获取样本数。最终调用父类 ClassifierMixin.score，默认使用 accuracy_score 计算准确率，并支持 sample_weight 加权。

21.5.3 DummyRegressor：回归基线的统计量策略

源码路径：sklearn/dummy.py - DummyRegressor.__init__（第 400-408 行）

    def __init__(self, *, strategy="mean", constant=None, quantile=None):
        self.strategy = strategy
        self.constant = constant
        self.quantile = quantile

概述：构造函数接收三个仅关键字参数：strategy 指定基线策略（默认 "mean"），constant 仅用于 "constant" 策略，quantile 仅用于 "quantile" 策略且需在 [0,1] 区间。

源码路径：sklearn/dummy.py - DummyRegressor.fit（第 410-480 行）

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y, sample_weight=None):
        """Fit the baseline regressor."""
        validate_data(self, X, skip_check_array=True)

        y = check_array(y, ensure_2d=False, input_name="y")
        if len(y) == 0:
            raise ValueError("y must not be empty.")

        if y.ndim == 1:
            y = np.reshape(y, (-1, 1))
        self.n_outputs_ = y.shape[1]

        check_consistent_length(X, y, sample_weight)

        if sample_weight is not None:
            sample_weight = _check_sample_weight(sample_weight, X)

        if self.strategy == "mean":
            self.constant_ = np.average(y, axis=0, weights=sample_weight)

        elif self.strategy == "median":
            if sample_weight is None:
                self.constant_ = np.median(y, axis=0)
            else:
                self.constant_ = _weighted_percentile(
                    y, sample_weight, percentile_rank=50.0
                )

        elif self.strategy == "quantile":
            if self.quantile is None:
                raise ValueError(
                    "When using `strategy='quantile', you have to specify the desired "
                    "quantile in the range [0, 1]."
                )
            percentile_rank = self.quantile * 100.0
            if sample_weight is None:
                self.constant_ = np.percentile(y, axis=0, q=percentile_rank)
            else:
                self.constant_ = _weighted_percentile(
                    y, sample_weight, percentile_rank=percentile_rank
                )

        elif self.strategy == "constant":
            if self.constant is None:
                raise TypeError(
                    "Constant target value has to be specified "
                    "when the constant strategy is used."
                )

            self.constant_ = check_array(
                self.constant,
                accept_sparse=["csr", "csc", "coo"],
                ensure_2d=False,
                ensure_min_samples=0,
            )

            if self.n_outputs_ != 1 and self.constant_.shape[0] != y.shape[1]:
                raise ValueError(
                    "Constant target value should have shape (%d, 1)." % y.shape[1]
                )

        self.constant_ = np.reshape(self.constant_, (1, -1))
        return self

概述：根据策略计算 constant_：mean 用加权均值；median 用加权中位数（_weighted_percentile 50%）；quantile 用加权分位数；constant 校验并保存用户常量。所有统计量均支持 sample_weight，体现了对加权样本的完整支持。最后将 constant_ 重塑为 (1, n_outputs) 便于广播。

策略行为对比：

下表对比了四种回归策略的学习目标、权重支持及典型用途，帮助理解不同基线的统计含义：

| 策略 | 学习的 constant_ | 样本权重支持 | 典型用途 |

|------|-------------------|--------------|----------|

| mean | 加权均值 np.average | 是 | 标准基线，对应 MSE 最优常数预测 |

| median | 加权中位数 _weighted_percentile(50%) | 是 | 鲁棒基线，对应 MAE 最优常数预测 |

| quantile | 加权分位数 _weighted_percentile(q*100) | 是 | 分位数回归基线，对应 Pinball Loss 最优 |

| constant | 用户指定常量 | 否 | 业务指定基准（如零预测） |

源码路径：sklearn/dummy.py - DummyRegressor.predict（第 482-515 行）

    def predict(self, X, return_std=False):
        """Perform classification on test vectors X."""
        check_is_fitted(self)
        n_samples = _num_samples(X)

        y = np.full(
            (n_samples, self.n_outputs_),
            self.constant_,
            dtype=np.array(self.constant_).dtype,
        )
        y_std = np.zeros((n_samples, self.n_outputs_))

        if self.n_outputs_ == 1:
            y = np.ravel(y)
            y_std = np.ravel(y_std)

        return (y, y_std) if return_std else y

概述：极简实现：np.full 广播 constant_ 到所有样本，return_std=True 时返回零方差（符合“确定性预测”语义）。constant_ 形状为 (1, n_outputs)，利用 NumPy 广播机制高效扩展。

源码路径：sklearn/dummy.py - DummyRegressor.score（第 520-550 行）

    def score(self, X, y, sample_weight=None):
        """Return the coefficient of determination R^2 of the prediction.

        The coefficient R^2 is defined as `(1 - u/v)`, where `u` is the
        residual sum of squares `((y_true - y_pred) ** 2).sum()` and `v` is the
        total sum of squares `((y_true - y_true.mean()) ** 2).sum()`. The best
        possible score is 1.0 and it can be negative (because the model can be
        arbitrarily worse). A constant model that always predicts the expected
        value of y, disregarding the input features, would get a R^2 score of
        0.0.

        Parameters
        ----------
        X : None or array-like of shape (n_samples, n_features)
            Test samples. Passing None as test samples gives the same result
            as passing real test samples, since `DummyRegressor`
            operates independently of the sampled observations.

        y : array-like of shape (n_samples,) or (n_samples, n_outputs)
            True values for X.

        sample_weight : array-like of shape (n_samples,), default=None
            Sample weights.

        Returns
        -------
        score : float
            R^2 of `self.predict(X)` w.r.t. y.
        """
        if X is None:
            X = np.zeros(shape=(len(y), 1))
        return super().score(X, y, sample_weight)

概述：与分类器类似，支持 X=None 并构造全零矩阵获取样本数。调用父类 RegressorMixin.score 计算 R²。常数预测若等于均值则 R²=0，若偏离则为负，直观反映“基线以上/以下”的性能。

架构图：DummyRegressor 核心流程

graph TD A[DummyRegressor] --> B[__init__: strategy/constant/quantile] B --> C[fit: validate_data 仅校验 X] C --> D[check_array 统一 y 为二维] D --> E{strategy?} E -- mean --> F[np.average 加权均值] E -- median --> G[_weighted_percentile 50%] E -- quantile --> H[_weighted_percentile q*100] E -- constant --> I[校验并保存用户常量] F --> J[reshape constant_ 为 (1, n_outputs)] G --> J H --> J I --> J J --> K[fit 完成] K --> L[predict: np.full 广播 constant_] L --> M{return_std?} M -- True --> N[返回 (y, zeros)] M -- False --> O[返回 y] K --> P[score: 支持 X=None] P --> Q[调用 RegressorMixin.score 计算 R²]

21.5.4 Dummy 估计器数据流图

graph TD A[fit(X, y, sample_weight)] --> B[validate_data 仅校验 X 形状] B --> C{strategy 类型} C -- 分类器 --> D[class_distribution 计算 classes_, n_classes_, class_prior_] D --> E{strategy 细分} E -- prior/most_frequent --> F[class_prior_.argmax() 得到众数类别] E -- stratified --> G[保存 class_prior_ 供采样用] E -- uniform --> H[记录 n_classes_ 供均匀采样用] E -- constant --> I[校验并保存用户常量] C -- 回归器 --> J{strategy 细分} J -- mean --> K[np.average 加权均值] J -- median --> L[_weighted_percentile 50%] J -- quantile --> M[_weighted_percentile q*100] J -- constant --> N[校验并保存用户常量] K --> O[存储 constant_ (1, n_outputs)] L --> O M --> O N --> O O --> P[fit 完成] P --> Q[predict(X)] Q --> R[_num_samples 获取 n_samples] R --> S{输出格式} S -- 稀疏 --> T[_random_choice_csc 生成稀疏预测] S -- 稠密 --> U[np.full 广播 constant_ / 采样生成预测] T --> V[返回预测] U --> V

21.6 设计中的取舍

为什么随机投影不学习数据相关的投影矩阵？

PCA 等方法通过 SVD 学习数据协方差的主成分方向，能以最少维度保留最大方差，但需要 O(n_features² × n_samples) 或 O(n_features³) 的计算，且必须看全量数据。随机投影完全无视数据分布，用随机矩阵“盲射”，理论上只需 O(n_components × n_features) 生成矩阵，transform 仅 O(n_samples × n_features × n_components)。JL 引理保证：只要目标维度足够高，任意数据的成对距离都能近似保持。权衡在于：优点是无需训练、可增量流式处理、极其适合超高维稀疏数据（如文本）、理论保证与数据无关；缺点是维度通常比 PCA 高得多（JL 界偏保守）、无法保留“最大方差方向”、解释性差。适用场景包括流式数据、极高维稀疏特征、作为下游任务的快速预处理、隐私保护（随机投影可视为一种扰动）。

为什么 Dummy 估计器不实现 partial_fit？

Dummy 估计器的 fit 仅计算全局统计量（众数、均值、分位数），这些统计量天然支持增量更新（如在线均值、在线分位数估计）。但 scikit-learn 故意不实现 partial_fit，原因是：定位清晰，Dummy 是“一次性基线”，用于实验开始前快速跑通流程、建立下限，不参与生产流水线；API 简洁，不引入增量学习的复杂性（如状态管理、批次合并逻辑）；替代方案是若真需流式基线，用户可自行维护统计量或用 SGDClassifier(loss='log_loss') 等真正的增量模型。

为什么 SparseRandomProjection 默认 dense_output=False？

稀疏矩阵乘稀疏矩阵结果往往更稀疏（非零元素数约为 nnz₁ × nnz₂ / n_features）。当 n_components 很大时，输出稀疏矩阵极其节省内存。但当 n_components 很小（如 < 100）时，输出几乎全非零，稠密格式反而更快（避免 CSR 索引开销）。因此暴露 dense_output 参数让用户根据实际维度权衡。

21.7 动手练习

1. 阅读随机投影实现

   • 阅读 sklearn/random_projection.py 中 GaussianRandomProjection 和 SparseRandomProjection 类的 fit 和 transform 方法，重点关注：

     1. 如何生成投影矩阵（高斯 vs 稀疏 Achlioptas 分布）

     2. fit 阶段完成什么工作（仅生成 components_）

     3. transform 如何进行矩阵乘法实现降维

   • 回答问题：

     • 两种投影在密集数据和稀疏数据上的计算效率有何不同？

     • 为什么 SparseRandomProjection 支持 dense_output 参数？

2. 探索 Dummy 估计器的策略行为

   • 阅读 sklearn/dummy.py 中 DummyClassifier 和 DummyRegressor 的 fit 和 predict 方法，重点关注：

     1. 不同 strategy 参数如何影响训练过程中学习到的常量（如 prior 是否考虑样本权重）

     2. predict 和 predict_proba 在分类器中的输出形式

     3. score 方法如何使用传入的 metric（默认 accuracy）评估预测

   • 回答问题：

     • DummyClassifier(strategy='uniform') 在多类别问题上会生成什么样的概率输出？

     • DummyRegressor(strategy='quantile', quantile=0.5) 的行为与什么特征等价？

3. 运行基线验证实验

   • 使用 sklearn 自带的 iris 或 boston 数据集，比较以下模型：

     1. LogisticRegression（或 LinearRegression）

     2. DummyClassifier(strategy='prior') / DummyRegressor(strategy='mean')

     3. GaussianRandomProjection + LogisticPipeline（可选）

   • 记录准确率或 R² 分数，回答：

     • 真实模型相较于 Dummy 基线的提升幅度是多少？

     • 添加随机投影是否显著改变了基线以上的模型性能？

21.8 本章小结

这一章中我们学习了随机投影与 Dummy 估计器这两个“化繁为简”的实用工具。首先，我们深入理解了 Johnson-Lindenstrauss 引理如何为随机投影提供理论保障，推导了最小投影维度的计算公式。其次，我们对比了 GaussianRandomProjection 与 SparseRandomProjection 的实现机制：前者用稠密高斯矩阵配合 BLAS 加速，适合中小规模稠密数据；后者用 CSR 稀疏矩阵配合 Achlioptas 分布，实现内存与计算的双重节省，适合高维稀疏或超大规模场景。接着，我们剖析了 BaseRandomProjection 基类如何统一 fit/transform/inverse_transform 流程，以及 n_components='auto' 与手动指定的两种模式。随后，我们详细拆解了 DummyClassifier 的五种策略在 fit/predict/predict_proba 中的行为差异，特别是 prior 与 most_frequent 在概率输出上的本质区别，以及 stratified 的采样一致性设计。最后，我们解读了 DummyRegressor 的四种策略对应的统计量（均值/中位数/分位数/常量）及其对样本权重的支持，理解了基线模型在模型评估中“参照系”的核心价值。

本章我们一起学习了以下概念：

下表汇总了本章核心概念及其解释，作为快速参考手册：

| 概念 | 解释 |

|------|------|

| Johnson-Lindenstrauss 引理 | 理论保障：低维嵌入可近似保持点间成对距离，是随机投影的核心依据 |

| GaussianRandomProjection | 使用高斯随机矩阵进行投影，适用于密集数据，投影矩阵由 N(0, 1) 生成 |

| SparseRandomProjection | 使用稀疏随机矩阵（Achlioptas 分布）进行投影，计算更快，内存占用更低 |

| johnson_lindenstrauss_min_dim | 计算满足 JL 引理所需的最小投影维度 n_components，是理论保障的工具函数 |

| DummyClassifier | 基线分类器，支持多种简单策略（如 prior、stratified、constant、uniform），用于建立模型性能下限 |

| DummyRegressor | 基线回归器，支持 mean、median、quantile、constant 等策略，提供回归任务的基准性能 |

下一章中，我们将学习 scikit-learn 概览 —— 数据科学家的“机器学习军火库”，从项目全貌、设计哲学与模块地图出发，带你鸟瞰 scikit-learn 的整体架构。

第 22 章 —— 聚类算法全景 —— 在“无标签世界”中寻找结构

22.1 学习目标

• 难度：★★★☆☆（3/5）

• 预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础

• 理解 K-Means 及其变体（MiniBatch、Bisecting）的核心流程、初始化策略与并行加速机制

• 掌握层次聚类的 Ward/Complete/Average/Single 链接准则、结构化聚类与堆优化合并

• 理解 DBSCAN/OPTICS 的邻域查询、核心样本扩展与可达距离计算

• 深入理解 HDBSCAN 的互达距离图、MST 构建、单链接树压缩、稳定性积分与 EOM 选择

• 掌握谱聚类的拉普拉斯特征嵌入、离散化分配与 QR 分解直接提标签

• 理解亲和传播的责任/可用性矩阵迭代与阻尼收敛机制

• 掌握 BIRCH 的 CF 树结构、增量更新、节点分裂与全局聚类两阶段流程

• 理解谱双聚类的三种归一化方法、奇异向量选择与行列协同划分

• 掌握均值漂移的带宽估计、网格分箱播种、并行爬山迭代与后处理去重

22.2 生活类比

想象聚类算法是一场大型社交派对的分组游戏：K-Means = 主持人预设桌数，客人按最近距离入座，主持人反复调整桌子中心位置直到稳定；MiniBatch 是分批次让客人入座，Bisecting 是先分两大桌再递归细分。Lloyd/Elkan = 两种找座位策略：Lloyd 遍历所有桌，Elkan 用三角不等式聪明地排除不可能的桌。层次聚类 = 从每人一桌开始合并，Ward 合并方差最小的桌，Complete 合并最远距离最近的桌，单链接是拉手链式合并。DBSCAN/OPTICS = 密度派对：核心客人（邻居多）拉人成团，边缘客人蹭团，噪音客人孤立；OPTICS 还记录每人被拉入团的难易度（可达距离）。HDBSCAN = 多层次密度派对：先算互达距离（考虑局部密度），建最小生成树，压缩树枝，算稳定性积分，最后选最稳分层。谱聚类 = 投影派对：把客人映射到低维特征空间，再用 K-Means 或离散化分组。亲和传播 = 传纸条选代表：互发责任/可用性纸条，收敛出 exemplar 代表。BIRCH = 流式派对：来一个客人合并到最近叶子簇，簇满则分裂，最后只对叶子簇中心做全局聚类。谱双聚类 = 行列同分派对：行客人和列客人同时分组，形成棋盘或块对角结构。均值漂移 = 爬山找峰派对：从种子点出发，往带宽内客人重心移动，直到到达密度峰顶。

22.3 源码地图

sklearn/cluster/_kmeans.py

├── kmeans_plusplus() # K-Means++ 初始化

├── _kmeans_plusplus() # 私有实现

├── _tolerance() # 数据依赖容差计算

├── k_means() # 函数式接口

├── _kmeans_single_elkan() # Elkan 单次运行

├── _kmeans_single_lloyd() # Lloyd 单次运行

├── _labels_inertia() # E步：标签与惯性

├── _BaseKMeans # 基类：参数验证、质心初始化、预测/变换/评分

│ ├── _check_params_vs_input()

│ ├── _init_centroids()

│ ├── fit_predict()

│ ├── predict()

│ ├── transform()

│ ├── fit_transform()

│ ├── score()

│ ├── _validate_center_shape()

│ ├── _check_test_data()

│ ├── _check_mkl_vcomp()

│ ├── _warn_mkl_vcomp()

│ ├── sklearn_tags()

│ └── _BaseKMeans

├── KMeans # 经典 K-Means

│ ├── _check_params_vs_input()

│ ├── _warn_mkl_vcomp()

│ └── fit()

├── _mini_batch_step() # MiniBatch 单步更新

└── MiniBatchKMeans # MiniBatch K-Means

├── _check_params_vs_input()

├── _warn_mkl_vcomp()

├── _mini_batch_convergence()

├── _random_reassign()

├── fit()

├── partial_fit()

└── init()

sklearn/cluster/_bisect_k_means.py

├── _BisectingTree # 层次树节点

│ ├── init()

│ ├── split()

│ ├── get_cluster_to_bisect()

│ ├── iter_leaves()

│ └── sklearn_tags()

└── BisectingKMeans # 二分 K-Means

├── _warn_mkl_vcomp()

├── _inertia_per_cluster()

├── _bisect()

├── fit()

├── predict()

├── _predict_recursive()

├── init()

└__sklearn_tags__()

sklearn/cluster/_k_means_lloyd.pyx

├── lloyd_iter_chunked_dense() # 稠密 Lloyd 并行迭代

├── _update_chunk_dense() # 稠密分块更新核心

├── lloyd_iter_chunked_sparse() # 稀疏 Lloyd 并行迭代

├── _update_chunk_sparse() # 稀疏分块更新核心

└── main

sklearn/cluster/_k_means_elkan.pyx

├── init_bounds_dense() # 稠密边界初始化

├── init_bounds_sparse() # 稀疏边界初始化

├── elkan_iter_chunked_dense() # 稠密 Elkan 并行迭代

├── _update_chunk_dense() # 稠密分块 Elkan 核心

├── elkan_iter_chunked_sparse() # 稀疏 Elkan 并行迭代

├── _update_chunk_sparse() # 稀疏分块 Elkan 核心

sklearn/cluster/_k_means_minibatch.pyx

├── _minibatch_update_dense() # 稠密 MiniBatch 中心更新

├── update_center_dense() # 单中心稠密更新

├── _minibatch_update_sparse() # 稀疏 MiniBatch 中心更新

├── update_center_sparse() # 单中心稀疏更新

sklearn/cluster/_k_means_common.pyx

├── CHUNK_SIZE # 分块常数

├── _euclidean_dense_dense() # 稠密-稠密距离

├── _euclidean_sparse_dense() # 稀疏-稠密距离

├── _inertia_dense() # 稠密惯性计算

├── _inertia_sparse() # 稀疏惯性计算

├── _relocate_empty_clusters_dense() # 稠密空簇重定位

├── _relocate_empty_clusters_sparse() # 稀疏空簇重定位

├── _average_centers() # 中心平均化

├── _center_shift() # 中心位移计算

├── _is_same_clustering() # 标签等价性检查

├── _euclidean_dense_dense_wrapper()

├── _euclidean_sparse_dense_wrapper()

└── main

sklearn/cluster/_k_means_common.pxd

├── _euclidean_dense_dense

├── _euclidean_sparse_dense

├── _relocate_empty_clusters_dense

├── _relocate_empty_clusters_sparse

├── _average_centers

├── _center_shift

sklearn/cluster/_agglomerative.py

├── _fix_connectivity() # 连接矩阵修复

├── _single_linkage_tree() # 单链接树构建

├── ward_tree() # Ward 结构化聚类

├── linkage_tree() # Average/Complete/Single 链接

├── _hc_cut() # 树切割生成标筓

├── _complete_linkage()

├── _average_linkage()

├── _single_linkage()

├── AgglomerativeClustering # 层次聚类估计器

│ ├── init()

│ ├── fit()

│ ├── _fit()

│ └── fit_predict()

└── FeatureAgglomeration # 特征凝聚

├── init()

├── fit()

└── fit_predict()

sklearn/cluster/_feature_agglomeration.py

└── AgglomerationTransform # 变换混入

├── transform()

└── inverse_transform()

sklearn/cluster/_hierarchical_fast.pyx

├── compute_ward_dist() # Ward 距离增量计算

├── _hc_get_descendent() # 获取后代叶子

├── hc_get_heads() # 获取树根标筥

├── _get_parents() # 并查集找祖先

├── max_merge() # Complete 合并策略

├── average_merge() # Average 合并策略

├── WeightedEdge # 堆边对象

│ ├── init()

│ ├── richcmp()

│ └── repr()

├── UnionFind # 并查集结构

│ ├── init()

│ ├── union()

│ └── fast_find()

├── _single_linkage_label() # MST 转单链接树

├── single_linkage_label() # 公开接口

└── mst_linkage_core() # MST-LINKAGE-CORE 算法

sklearn/cluster/_hierarchical_fast.pxd

└── UnionFind # 并查集结构声明

├── union

└── fast_find

sklearn/cluster/_dbscan.py

├── dbscan() # 函数式接口

└── DBSCAN # DBSCAN 估计器

├── init()

├── fit()

├── fit_predict()

└── sklearn_tags()

sklearn/cluster/_dbscan_inner.pyx

└── dbscan_inner() # Cython 核心 DFS 扩展

sklearn/cluster/_optics.py

├── _validate_size() # 参数大小校验

├── compute_core_distances() # 核心距离分块计算

├── compute_optics_graph() # 可达图主循环

├── _set_reach_dist() # 更新可达距离

├── cluster_optics_dbscan() # DBSCAN 风格提取

├── cluster_optics_xi() # Xi 方法提取

├── _xi_cluster() # Xi 核心逻辑

├── _extend_region() # 陡峭区域延伸

├── _update_filter_sdas() # 陡峭下降区过滤

├── _correct_predecessor() # 前驱校正

├── _extract_xi_labels() # 标签提取

└── OPTICS # OPTICS 估计器

├── init()

└── fit()

sklearn/cluster/_hdbscan/hdbscan.py

├── _OUTLIER_ENCODING # 异常标签编码

├── _brute_mst() # 稠密/稀疏 MST 构建

├── _process_mst() # MST 转单链接树

├── _hdbscan_brute() # 暴力模式完整流程

├── _hdbscan_prims() # Prim 算法模式

├── remap_single_linkage_tree() # 非有限值重映射

├── _get_finite_row_indices() # 有限行索引

└── HDBSCAN # HDBSCAN 估计器

├── init()

├── fit()

├── fit_predict()

├── _weighted_cluster_center()

├── dbscan_clustering()

└── sklearn_tags()

sklearn/cluster/_hdbscan/_linkage.pyx

├── MST_edge_dtype # MST 边结构体

├── mst_from_mutual_reachability() # 互达图 Prim MST

├── mst_from_data_matrix() # 数据矩阵 Prim MST

├── make_single_linkage() # MST 转层级树

sklearn/cluster/_hdbscan/_reachability.pyx

├── mutual_reachability_graph() # 互达图入口

├── _dense_mutual_reachability_graph() # 稠密实现

└── _sparse_mutual_reachability_graph() # 稀疏实现

sklearn/cluster/_hdbscan/_tree.pxd

├── HIERARCHY_t # 层级树节点结构

└── CONDENSED_t # 压缩树节点结构

sklearn/cluster/_hdbscan/_tree.pyx

├── HIERARCHY_dtype / CONDENSED_dtype

├── tree_to_labels() # 树转标签主流程

├── bfs_from_hierarchy() # 层级树广度优先搜索

├── _condense_tree() # 树压缩剪枝

├── _compute_stability() # 簇稳定性积分计算

├── bfs_from_cluster_tree() # 压缩树广度优先搜索

├── max_lambdas() # 最大 lambda 计算

├── TreeUnionFind # 标签分配并查集

│ ├── init()

│ ├── union()

│ └── find()

├── labelling_at_cut() # DBSCAN* 平切标签

├── _do_labelling() # 标签分配核心

├── get_probabilities() # 成员概率计算

├── recurse_leaf_dfs() # 叶子簇 DFS

├── get_cluster_tree_leaves() # 获取叶子簇

├── traverse_upwards() # 向上遍历合并

├── epsilon_search() # Epsilon 搜索合并

└── _get_clusters() # EOM/Leaf 选择主逻辑

sklearn/cluster/_spectral.py

├── cluster_qr() # QR 分解离散化

├── discretize() # 迭代离散化搜索

├── spectral_clustering() # 函数式接口

└── SpectralClustering # 谱聚类估计器

├── init()

├── fit()

├── fit_predict()

└── sklearn_tags()

sklearn/cluster/_affinity_propagation.py

├── _equal_similarities_and_preferences()

├── _affinity_propagation() # 核心消息传递循环

├── affinity_propagation() # 函数式接口

└── AffinityPropagation # 亲和传播估计器

├── init()

├── fit()

├── predict()

├── fit_predict()

└── sklearn_tags()

sklearn/cluster/_birch.py

├── _iterate_sparse_X()

├── _split_node() # CF节点分裂

├── _CFNode # CF 树节点

│ ├── init()

│ ├── append_subcluster()

│ ├── update_split_subclusters()

│ └── insert_cf_subcluster()

├── _CFSubcluster # CF 子簇

│ ├── init()

│ ├── update()

│ ├── merge_subcluster()

│ └── radius

└── Birch # BIRCH 估计器

├── init()

├── _fit() / partial_fit()

├── _get_leaves()

├── _global_clustering()

├── _predict()

├── predict()

├── transform()

└── sklearn_tags()

sklearn/cluster/_bicluster.py

├── _scale_normalize() # 行列独立归一化

├── _bistochastic_normalize() # Sinkhorn 双随机归一化

├── _log_normalize() # 对数交互归一化

├── BaseSpectral # 光谱双聚类基类

│ ├── init()

│ ├── _check_parameters()

│ ├── _svd()

│ ├── _k_means()

│ ├── fit()

│ └── sklearn_tags()

├── SpectralCoclustering # 共聚类

│ ├── init()

│ ├── _check_parameters()

│ └── _fit()

└── SpectralBiclustering # 双聚类

├── init()

├── _check_parameters()

├── _fit()

├── _fit_best_piecewise()

├── _project_and_cluster()

└── sklearn_tags()

sklearn/cluster/_mean_shift.py

├── estimate_bandwidth() # 带宽估计

├── _mean_shift_single_seed() # 单种子爬山

├── mean_shift() # 函数式接口

├── get_bin_seeds() # 网格分箱播种

└── MeanShift # 均值漂移估计器

├── init()

└── fit() / predict()

22.4 K-Means 及其变体 —— 聚类的“瑞士军刀”

K-Means 聚类算法是机器学习中最经典的基于中心的划分方法，通过迭代优化簇中心来最小化簇内方差（惯性）。其核心流程包括参数验证、数据预处理、质心初始化、多次运行选择最优惯性等步骤。

KMeans 类的核心流程与参数验证

源码路径：sklearn/cluster/_kmeans.py - KMeans.fit（600-700行）

def fit(self, X, y=None, sample_weight=None):
    """Compute k-means clustering.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Training instances to cluster. It must be noted that the data
        will be converted to C ordering, which will cause a memory
        copy if the given data is not C-contiguous.
        If a sparse matrix is passed, a copy will be made if it's not in
        CSR format.

    y : Ignored
        Not used, present here for API consistency by convention.

    sample_weight : array-like of shape (n_samples,), default=None
        The weights for each observation in X. If None, all observations
        are assigned equal weight. `sample_weight` is not used during
        initialization if `init` is a callable or a user provided array.

        .. versionadded:: 0.20

    Returns
    -------
    self : object
        Fitted estimator.
    """
    X = validate_data(
        self,
        X,
        accept_sparse="csr",
        dtype=[np.float64, np.float32],
        order="C",
        copy=self.copy_x,
        accept_large_sparse=False,
    )

    self._check_params_vs_input(X)

    random_state = check_random_state(self.random_state)
    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
    self._n_threads = _openmp_effective_n_threads()

    # Validate init array
    init = self.init
    init_is_array_like = _is_arraylike_not_scalar(init)
    if init_is_array_like:
        init = check_array(init, dtype=X.dtype, copy=True, order="C")
        self._validate_center_shape(X, init)

    # subtract of mean of x for more accurate distance computations
    if not sp.issparse(X):
        X_mean = X.mean(axis=0)
        # The copy was already done above
        X -= X_mean

        if init_is_array_like:
            init -= X_mean

    # precompute squared norms of data points
    x_squared_norms = row_norms(X, squared=True)

    if self._algorithm == "elkan":
        kmeans_single = _kmeans_single_elkan
    else:
        kmeans_single = _kmeans_single_lloyd
        self._check_mkl_vcomp(X, X.shape[0])

    best_inertia, best_labels = None, None

    for i in range(self._n_init):
        # Initialize centers
        centers_init = self._init_centroids(
            X,
            x_squared_norms=x_squared_norms,
            init=init,
            random_state=random_state,
            sample_weight=sample_weight,
        )
        if self.verbose:
            print("Initialization complete")

        # run a k-means once
        labels, inertia, centers, n_iter_ = kmeans_single(
            X,
            sample_weight,
            centers_init,
            max_iter=self.max_iter,
            verbose=self.verbose,
            tol=self._tol,
            n_threads=self._n_threads,
        )

        # determine if these results are the best so far
        # we chose a new run if it has a better inertia and the clustering is
        # different from the best so far (it's possible that the inertia is
        # slightly better even if the clustering is the same with potentially
        # permuted labels, due to rounding errors)
        if best_inertia is None or (
            inertia < best_inertia
            and not _is_same_clustering(labels, best_labels, self.n_clusters)
        ):
            best_labels = labels
            best_centers = centers
            best_inertia = inertia
            best_n_iter = n_iter_

    if not sp.issparse(X):
        if not self.copy_x:
            X += X_mean
        best_centers += X_mean

    distinct_clusters = len(set(best_labels))
    if distinct_clusters < self.n_clusters:
        warnings.warn(
            "Number of distinct clusters ({}) found smaller than "
            "n_clusters ({}). Possibly due to duplicate points "
            "in X.".format(distinct_clusters, self.n_clusters),
            ConvergenceWarning,
            stacklevel=2,
        )

    self.cluster_centers_ = best_centers
    self._n_features_out = self.cluster_centers_.shape[0]
    self.labels_ = best_labels
    self.inertia_ = best_inertia
    self.n_iter_ = best_n_iter
    return self

这段代码定义了 KMeans 类的主训练流程。它首先验证输入数据，然后进行参数检查（包括动态设置 _n_init 和容差 _tol），接着根据初始化策略（k-means++ 或 random）生成初始质心，并通过多次运行选择惯性最小的结果。值得注意的是，它支持 lloyd 和 elkan 两种算法，并在处理稀疏数据时仅接受 CSR 格式，因为稀疏矩阵的中心化操作会破坏其稀疏性。

MiniBatchKMeans 增量学习机制

源码路径：sklearn/cluster/_kmeans.py - MiniBatchKMeans.fit（800-950行）

def fit(self, X, y=None, sample_weight=None):
    """Compute the centroids on X by chunking it into mini-batches.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Training instances to cluster. It must be noted that the data
        will be converted to C ordering, which will cause a memory copy
        if the given data is not C-contiguous.
        If a sparse matrix is passed, a copy will be made if it's not in
        CSR format.

    y : Ignored
        Not used, present here for API consistency by convention.

    sample_weight : array-like of shape (n_samples,), default=None
        The weights for each observation in X. If None, all observations
        are assigned equal weight. `sample_weight` is not used during
        initialization if `init` is a callable or a user provided array.

        .. versionadded:: 0.20

    Returns
    -------
    self : object
        Fitted estimator.
    """
    X = validate_data(
        self,
        X,
        accept_sparse="csr",
        dtype=[np.float64, np.float32],
        order="C",
        accept_large_sparse=False,
    )

    self._check_params_vs_input(X)
    random_state = check_random_state(self.random_state)
    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
    self._n_threads = _openmp_effective_n_threads()
    n_samples, n_features = X.shape

    # Validate init array
    init = self.init
    if _is_arraylike_not_scalar(init):
        init = check_array(init, dtype=X.dtype, copy=True, order="C")
        self._validate_center_shape(X, init)

    self._check_mkl_vcomp(X, self._batch_size)

    # precompute squared norms of data points
    x_squared_norms = row_norms(X, squared=True)

    # Validation set for the init
    validation_indices = random_state.randint(0, n_samples, self._init_size)
    X_valid = X[validation_indices]
    sample_weight_valid = sample_weight[validation_indices]

    # perform several inits with random subsets
    best_inertia = None
    for init_idx in range(self._n_init):
        if self.verbose:
            print(f"Init {init_idx + 1}/{self._n_init} with method {init}")

        # Initialize the centers using only a fraction of the data as we
        # expect n_samples to be very large when using MiniBatchKMeans.
        cluster_centers = self._init_centroids(
            X,
            x_squared_norms=x_squared_norms,
            init=init,
            random_state=random_state,
            init_size=self._init_size,
            sample_weight=sample_weight,
            )

        # Compute inertia on a validation set.
        _, inertia = _labels_inertia_threadpool_limit(
            X_valid,
            sample_weight_valid,
            cluster_centers,
            n_threads=self._n_threads,
        )

        if self.verbose:
            print(f"Inertia for init {init_idx + 1}/{self._n_init}: {inertia}")
        if best_inertia is None or inertia < best_inertia:
            init_centers = cluster_centers
            best_inertia = inertia

    centers = init_centers
    centers_new = np.empty_like(centers)

    # Initialize counts
    self._counts = np.zeros(self.n_clusters, dtype=X.dtype)

    # Attributes to monitor the convergence
    self._ewa_inertia = None
    self._ewa_inertia_min = None
    self._no_improvement = 0

    # Initialize number of samples seen since last reassignment
    self._n_since_last_reassign = 0

    sum_of_weights = np.sum(sample_weight)

    n_steps = (self.max_iter * n_samples) // self._batch_size
    normalized_sample_weight = sample_weight / sum_of_weights
    unit_sample_weight = np.ones_like(sample_weight, shape=(self._batch_size,))

    with _get_threadpool_controller().limit(limits=1, user_api="blas"):
        # Perform the iterative optimization until convergence
        for i in range(n_steps):
            # Sample a minibatch from the full dataset
            minibatch_indices = random_state.choice(
                n_samples,
                self._batch_size,
                p=normalized_sample_weight,
                replace=True,
            )
            # Perform the actual update step on the minibatch data
            # Note: since the sampling of the minibatch is sample_weight aware,
            # we pass fixed unit weights to the `_mini_batch_step` call to avoid
            # accounting for the weights twice. Also note that `_mini_batch_step`
            # can be called with non-unit weights when the caller constructs
            # the batches explicitly by calling the public `partial_fit` method
            # instead.
            batch_inertia = _mini_batch_step(
                X=X[minibatch_indices],
                sample_weight=unit_sample_weight,
                centers=centers,
                centers_new=centers_new,
                weight_sums=self._counts,
                random_state=random_state,
                random_reassign=self._random_reassign(),
                reassignment_ratio=self.reassignment_ratio,
                verbose=self.verbose,
                n_threads=self._n_threads,
            )

            if self._tol > 0.0:
                centers_squared_diff = np.sum((centers_new - centers) ** 2)
            else:
                centers_squared_diff = 0

            centers, centers_new = centers_new, centers

            # Monitor convergence and do early stopping if necessary
            if self._mini_batch_convergence(
                i, n_steps, n_samples, centers_squared_diff, batch_inertia
            ):
                break

    self.cluster_centers_ = centers
    self._n_features_out = self.cluster_centers_.shape[0]

    self.n_steps_ = i + 1
    self.n_iter_ = int(np.ceil(((i + 1) * self._batch_size) / n_samples))

    if self.compute_labels:
        self.labels_, self.inertia_ = _labels_inertia_threadpool_limit(
            X,
            sample_weight,
            self.cluster_centers_,
            n_threads=self._n_threads,
        )
    else:
        self.inertia_ = self._ewa_inertia * sum_of_weights

    return self

这段代码实现了 MiniBatchKMeans 的增量学习机制。它通过将数据分成小批次（mini-batch）来处理大规模数据集，避免一次性加载全部数据到内存。核心是 _mini_batch_step 函数，它在每个小批次上执行标签分配和中心更新，并支持通过 reassignment_ratio 控制低计数簇的随机重分配，防止过早收敛。早停策略基于指数加权平均惯性（EWA）和中心位移的双重判断，使算法在收敛前就能停止，提高效率。

BisectingKMeans 分层二分策略

源码路径：sklearn/cluster/_bisect_k_means.py - BisectingKMeans.fit（200-300行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None, sample_weight=None):
    """Compute bisecting k-means clustering.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)

        Training instances to cluster.

        .. note:: The data will be converted to C ordering,
            which will cause a memory copy
            if the given data is not C-contiguous.

    y : Ignored
        Not used, present here for API consistency by convention.

    sample_weight : array-like of shape (n_samples,), default=None
        The weights for each observation in X. If None, all observations
        are assigned equal weight. `sample_weight` is not used during
        initialization if `init` is a callable.

    Returns
        self
            Fitted estimator.
    """
    X = validate_data(
        self,
        X,
        accept_sparse="csr",
        dtype=[np.float64, np.float32],
        order="C",
        copy=self.copy_x,
        accept_large_sparse=False,
    )

    self._check_params_vs_input(X)

    self._random_state = check_random_state(self.random_state)
    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
    self._n_threads = _openmp_effective_n_threads()

    if self.algorithm == "lloyd" or self.n_clusters == 1:
        self._kmeans_single = _kmeans_single_lloyd
        self._check_mkl_vcomp(X, X.shape[0])
    else:
        self._kmeans_single = _kmeans_single_elkan

    # Subtract of mean of X for more accurate distance computations
    if not sp.issparse(X):
        self._X_mean = X.mean(axis=0)
        X -= self._X_mean

    # Initialize the hierarchical clusters tree
    self._bisecting_tree = _BisectingTree(
        indices=np.arange(X.shape[0]),
        center=X.mean(axis=0),
        score=0,
    )

    x_squared_norms = row_norms(X, squared=True)

    for _ in range(self.n_clusters - 1):
        # Chose cluster to bisect
        cluster_to_bisect = self._bisecting_tree.get_cluster_to_bisect()

        # Split this cluster into 2 subclusters
        self._bisect(X, x_squared_norms, sample_weight, cluster_to_bisect)

    # Aggregate final labels and centers from the bisecting tree
    self.labels_ = np.full(X.shape[0], -1, dtype=np.int32)
    self.cluster_centers_ = np.empty((self.n_clusters, X.shape[1]), dtype=X.dtype)

    for i, cluster_node in enumerate(self._bisecting_tree.iter_leaves()):
        self.labels_[cluster_node.indices] = i
        self.cluster_centers_[i] = cluster_node.center
        cluster_node.label = i  # label final clusters for future prediction
        cluster_node.indices = None  # release memory

    # Restore original data
    if not sp.issparse(X):
        X += self._X_mean
        self.cluster_centers_ += self._X_mean

    _inertia = _inertia_sparse if sp.issparse(X) else _inertia_dense
    self.inertia_ = _inertia(
        X, sample_weight, self.cluster_centers_, self.labels_, self._n_threads
    )

    self._n_features_out = self.cluster_centers_.shape[0]

    return self

这段代码实现了 BisectingKMeans 的分层二分策略。它不是直接将所有点分成 k 个簇，而是采用自底向上的层次方法：开始时所有点在一个簇中，然后反复选择“最值得分裂”的簇（基于 inertia 或样本数）并将其二分，直到达到所需的簇数。这种方法在处理大数据集时更高效，因为每次只需在一个子集上运行 K-Means（二分），而非全量数据。预测时，它通过递归下沉树结构（_predict_recursive）来分配标签，而不是直接使用全局质心，这确保了预测与训练时的层次结构一致。

22.5 Cython 加速核心 —— K-Means 的“引擎室”

K-Means 的计算瓶颈主要在于距离计算和中心更新。为了提升性能，scikit-learn 使用 Cython 实现了 Lloyd 和 Elkan 算法的并行版本，并通过分块处理和 BLAS 加速来优化内存访问和计算效率。

Lloyd 算法并行实现

源码路径：sklearn/cluster/_k_means_lloyd.pyx - lloyd_iter_chunked_dense（20-120行）

def lloyd_iter_chunked_dense(
        const floating[:, ::1] X,            # IN
        const floating[::1] sample_weight,   # IN
        const floating[:, ::1] centers_old,  # IN
        floating[:, ::1] centers_new,        # OUT
        floating[::1] weight_in_clusters,    # OUT
        int[::1] labels,                     # OUT
        floating[::1] center_shift,          # OUT
        int n_threads,
        bint update_centers=True):
    """Single iteration of K-means lloyd algorithm with dense input.

    Update labels and centers (inplace), for one iteration, distributed
    over data chunks.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features), dtype=floating
        The observations to cluster.

    sample_weight : ndarray of shape (n_samples,), dtype=floating
        The weights for each observation in X.

    centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
        Centers before previous iteration, placeholder for the centers after
        previous iteration.

    centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
        Centers after previous iteration, placeholder for the new centers
        computed during this iteration. `centers_new` can be `None` if
        `update_centers` is False.

    weight_in_clusters : ndarray of shape (n_clusters,), dtype=floating
        Placeholder for the sums of the weights of every observation assigned
        to each center. `weight_in_clusters` can be `None` if `update_centers`
        is False.

    labels : ndarray of shape (n_samples,), dtype=int
        labels assignment.

    center_shift : ndarray of shape (n_clusters,), dtype=floating
        Distance between old and new centers.

    n_threads : int
        The number of threads to be used by openmp.

    update_centers : bool
        - If True, the labels and the new centers will be computed, i.e. runs
          the E-step and the M-step of the algorithm.
        - If False, only the labels will be computed, i.e runs the E-step of
          the algorithm. This is useful especially when calling predict on a
          fitted model.
    """
    cdef:
        int n_samples = X.shape[0]
        int n_features = X.shape[1]
        int n_clusters = centers_old.shape[0]

    if n_samples == 0:
        # An empty array was passed, do nothing and return early (before
        # attempting to compute n_chunks). This can typically happen when
        # calling the prediction function of a bisecting k-means model with a
        # large fraction of outliers.
        return

    cdef:
        # hard-coded number of samples per chunk. Appeared to be close to
        # optimal in all situations.
        int n_samples_chunk = CHUNK_SIZE if n_samples > CHUNK_SIZE else n_samples
        int n_chunks = n_samples // n_samples_chunk
        int n_samples_rem = n_samples % n_samples_chunk
        int chunk_idx
        int start, end

        int j, k

        floating[::1] centers_squared_norms = row_norms(centers_old, squared=True)

        floating *centers_new_chunk
        floating *weight_in_clusters_chunk
        floating *pairwise_distances_chunk

        omp_lock_t lock

    # count remainder chunk in total number of chunks
    n_chunks += n_samples != n_chunks * n_samples_chunk

    # number of threads should not be bigger than number of chunks
    n_threads = min(n_threads, n_chunks)

    if update_centers:
        memset(&centers_new[0, 0], 0, n_clusters * n_features * sizeof(floating))
        memset(&weight_in_clusters[0], 0, n_clusters * sizeof(floating))
        omp_init_lock(&lock)

    with nogil, parallel(num_threads=n_threads):
        # thread local buffers
        centers_new_chunk = <floating*> calloc(n_clusters * n_features, sizeof(floating))
        weight_in_clusters_chunk = <floating*> calloc(n_clusters, sizeof(floating))
        pairwise_distances_chunk = <floating*> malloc(n_samples_chunk * n_clusters * sizeof(floating))

        for chunk_idx in prange(n_chunks, schedule='static'):
            start = chunk_idx * n_samples_chunk
            if chunk_idx == n_chunks - 1 and n_samples_rem > 0:
                end = start + n_samples_rem
            else:
                end = start + n_samples_chunk

            _update_chunk_dense(
                X[start: end],
                sample_weight[start: end],
                centers_old,
                centers_squared_norms,
                labels[start: end],
                centers_new_chunk,
                weight_in_clusters_chunk,
                pairwise_distances_chunk,
                update_centers)

        # reduction from local buffers.
        if update_centers:
            # The lock is necessary to avoid race conditions when aggregating
            # info from different thread-local buffers.
            omp_set_lock(&lock)
            for j in range(n_clusters):
                weight_in_clusters[j] += weight_in_clusters_chunk[j]
                for k in range(n_features):
                    centers_new[j, k] += centers_new_chunk[j * n_features + k]

            omp_unset_lock(&lock)

        free(centers_new_chunk)
        free(weight_in_clusters_chunk)
        free(pairwise_distances_chunk)

    if update_centers:
        omp_destroy_lock(&lock)
        _relocate_empty_clusters_dense(
            X, sample_weight, centers_old, centers_new, weight_in_clusters, labels
        )

        _average_centers(centers_new, weight_in_clusters)
        _center_shift(centers_old, centers_new, center_shift)

这段代码实现了 Lloyd 算法的稠密数据并行版本。它将数据按 CHUNK_SIZE=256 分块，使用 OpenMP 的 prange 实现样本级并行。每个线程处理自己的数据块，计算局部的标签分配和中心更新贡献，然后通过锁机制将结果归约（reduction）到全局数组。关键优化在于 _update_chunk_dense 函数中使用 BLAS 的 _gemm 来高效计算 -2 * X @ C.T 项，避免了显式的三重循环，显著提升了矩阵乘法的性能。稀疏版本则直接遍历 CSR 格式的非零元素，避免了不必要的零值计算。

Elkan 算法三角不等式加速

源码路径：sklearn/cluster/_k_means_elkan.pyx - elkan_iter_chunked_dense（80-200行）

def elkan_iter_chunked_dense(
        const floating[:, ::1] X,                      # IN
        const floating[::1] sample_weight,             # IN
        const floating[:, ::1] centers_old,            # IN
        floating[:, ::1] centers_new,                  # OUT
        floating[::1] weight_in_clusters,              # OUT
        const floating[:, ::1] center_half_distances,  # IN
        const floating[::1] distance_next_center,      # IN
        floating[::1] upper_bounds,                    # INOUT
        floating[:, ::1] lower_bounds,                 # INOUT
        int[::1] labels,                               # INOUT
        floating[::1] center_shift,                    # OUT
        int n_threads,
        bint update_centers=True):
    """Single iteration of K-means Elkan algorithm with dense input.

    Update labels and centers (inplace), for one iteration, distributed
    over data chunks.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features), dtype=floating
        The observations to cluster.

    sample_weight : ndarray of shape (n_samples,), dtype=floating
        The weights for each observation in X.

    centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
        Centers before previous iteration, placeholder for the centers after
        previous iteration.

    centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
        Centers after previous iteration, placeholder for the new centers
        computed during this iteration.

    weight_in_clusters : ndarray of shape (n_clusters,), dtype=floating
        Placeholder for the sums of the weights of every observation assigned
        to each center.

    center_half_distances : ndarray of shape (n_clusters, n_clusters), \
            dtype=floating
        Half pairwise distances between centers.

    distance_next_center : ndarray of shape (n_clusters,), dtype=floating
        Distance between each center its closest center.

    upper_bounds : ndarray of shape (n_samples,), dtype=floating
        Upper bound for the distance between each sample and its center,
        updated inplace.

    lower_bounds : ndarray of shape (n_samples, n_clusters), dtype=floating
        Lower bound for the distance between each sample and each center,
        updated inplace.

    labels : ndarray of shape (n_samples,), dtype=int
        labels assignment.

    center_shift : ndarray of shape (n_clusters,), dtype=floating
        Distance between old and new centers.

    n_threads : int
        The number of threads to be used by openmp.

    update_centers : bool
        - If True, the labels and the new centers will be computed, i.e. runs
          the E-step and the M-step of the algorithm.
        - If False, only the labels will be computed, i.e runs the E-step of
          the algorithm. This is useful especially when calling predict on a
          fitted model.
    """
    cdef:
        int n_samples = X.shape[0]
        int n_features = X.shape[1]
        int n_clusters = centers_new.shape[0]

    if n_samples == 0:
        # An empty array was passed, do nothing and return early (before
        # attempting to compute n_chunks). This can typically happen when
        # calling the prediction function of a bisecting k-means model with a
        # large fraction of outliers.
        return

    cdef:
        # hard-coded number of samples per chunk. Splitting in chunks is
        # necessary to get parallelism. Chunk size chosen to be same as lloyd's
        int n_samples_chunk = CHUNK_SIZE if n_samples > CHUNK_SIZE else n_samples
        int n_chunks = n_samples // n_samples_chunk
        int n_samples_rem = n_samples % n_samples_chunk
        int chunk_idx
        int start, end

        int i, j, k

        floating *centers_new_chunk
        floating *weight_in_clusters_chunk

        omp_lock_t lock

    # count remainder chunk in total number of chunks
    n_chunks += n_samples != n_chunks * n_samples_chunk

    # number of threads should not be bigger than number of chunks
    n_threads = min(n_threads, n_chunks)

    if update_centers:
        memset(&centers_new[0, 0], 0, n_clusters * n_features * sizeof(floating))
        memset(&weight_in_clusters[0], 0, n_clusters * sizeof(floating))
        omp_init_lock(&lock)

    with nogil, parallel(num_threads=n_threads):
        # thread local buffers
        centers_new_chunk = <floating*> calloc(n_clusters * n_features, sizeof(floating))
        weight_in_clusters_chunk = <floating*> calloc(n_clusters, sizeof(floating))

        for chunk_idx in prange(n_chunks, schedule='static'):
            start = chunk_idx * n_samples_chunk
            if chunk_idx == n_chunks - 1 and n_samples_rem > 0:
                end = start + n_samples_rem
            else:
                end = start + n_samples_chunk

            _update_chunk_dense(
                X[start: end],
                sample_weight[start: end],
                centers_old,
                center_half_distances,
                distance_next_center,
                labels[start: end],
                upper_bounds[start: end],
                lower_bounds[start: end],
                centers_new_chunk,
                weight_in_clusters_chunk,
                update_centers)

        # reduction from local buffers.
        if update_centers:
            # The lock is necessary to avoid race conditions when aggregating
            # info from different thread-local buffers.
            omp_set_lock(&lock)
            for j in range(n_clusters):
                weight_in_clusters[j] += weight_in_clusters_chunk[j]
                for k in range(n_features):
                    centers_new[j, k] += centers_new_chunk[j * n_features + k]
            omp_unset_lock(&lock)

        free(centers_new_chunk)
        free(weight_in_clusters_chunk)

    if update_centers:
        omp_destroy_lock(&lock)
        _relocate_empty_clusters_dense(X, sample_weight, centers_old,
                                       centers_new, weight_in_clusters, labels)

        _average_centers(centers_new, weight_in_clusters)
        _center_shift(centers_old, centers_new, center_shift)

        # update lower and upper bounds
        for i in range(n_samples):
            upper_bounds[i] += center_shift[labels[i]]

            for j in range(n_clusters):
                lower_bounds[i, j] -= center_shift[j]
                if lower_bounds[i, j] < 0:
                    lower_bounds[i, j] = 0

这段代码实现了 Elkan 算法的核心加速机制。Elkan 算法通过维护两个边界数组来避免不必要的距离计算：upper_bounds[i] 记录样本 i 到其最近质心的距离上界，lower_bounds[i, j] 记录样本 i 到质心 j 的距离下界。在每次迭代中，只有当上界大于下界且上界大于质心间半距离时，才需要重新计算真实距离。这种基于三角不等式的剪枝策略可以显著减少距离计算次数，特别是在簇间距离较大、簇内紧凑的数据集上。然而，它需要额外的内存来存储 lower_bounds 数组（大小为 n_samples × n_clusters），当簇数很多时，这会成为内存瓶颈。_update_chunk_dense 中使用 _gemm 计算 -2 * X @ C.T 是为了高效构建距离矩阵的一部分，利用 BLAS 库的优化矩阵乘法，比显式循环快得多。

22.6 层次聚类与特征凝聚 —— 构建数据的“系谱树”

层次聚类通过自底向上（或自顶向下）合并或分裂簇来构建簇的层次结构（dendrogram），无需预先指定簇数。其核心在于链接准则的确定哪两个簇应该先合并，以及如何高效地更新簇间距离。

AgglomerativeClustering 树构建策略

源码路径：sklearn/cluster/_agglomerative.py - ward_tree（100-250行）

def ward_tree(X, *, connectivity=None, n_clusters=None, return_distance=False):
    """Ward clustering based on a Feature matrix.

    Recursively merges the pair of clusters that minimally increases
    within-cluster variance.

    The inertia matrix uses a Heapq-based representation.

    This is the structured version, that takes into account some topological
    structure between samples.

    Read more in the :ref:`User Guide <hierarchical_clustering>`.

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Feature matrix representing `n_samples` samples to be clustered.

    connectivity : {array-like, sparse matrix}, default=None
        Connectivity matrix. Defines for each sample the neighboring samples
        following a given structure of the data. The matrix is assumed to
        be symmetric and only the upper triangular half is used.
        Default is None, i.e, the Ward algorithm is unstructured.

    n_clusters : int, default=None
        `n_clusters` should be less than `n_samples`.  Stop early the
        construction of the tree at `n_clusters.` This is useful to decrease
        computation time if the number of clusters is not small compared to the
        number of samples. In this case, the complete tree is not computed, thus
        the 'children' output is of limited use, and the 'parents' output should
        rather be used. This option is valid only when specifying a connectivity
        matrix.

    return_distance : bool, default=False
        If `True`, return the distance between the clusters.

    Returns
    -------
    children : ndarray of shape (n_nodes-1, 2)
        The children of each non-leaf node. Values less than `n_samples`
        correspond to leaves of the tree which are the original samples.
        A node `i` greater than or equal to `n_samples` is a non-leaf
        node and has children `children_[i - n_samples]`. Alternatively
        at the i-th iteration, children[i][0] and children[i][1]
        are merged to form node `n_samples + i`.

    n_connected_components : int
        The number of connected components in the graph.

    n_leaves : int
        The number of leaves in the tree.

    parents : ndarray of shape (n_nodes,) or None
        The parent of each node. Only returned when a connectivity matrix
        is specified, elsewhere 'None' is returned.

    distances : ndarray of shape (n_nodes-1,)
        Only returned if `return_distance` is set to `True` (for compatibility).
        The distances between the centers of the nodes. `distances[i]`
        corresponds to a weighted Euclidean distance between
        the nodes `children[i, 1]` and `children[i, 2]`. If the nodes refer to
        leaves of the tree, then `distances[i]` is their unweighted Euclidean
        distance. Distances are updated in the following way
        (from scipy.hierarchy.linkage):

        The new entry :math:`d(u,v)` is computed as follows,

        .. math::

           d(u,v) = \\sqrt{\\frac{|v|+|s|}
                               {T}d(v,s)^2
                        + \\frac{|v|+|t|}
                               {T}d(v,t)^2
                        - \\frac{|v|}
                               {T}d(s,t)^2}

        where :math:`u` is the newly joined cluster consisting of
        clusters :math:`s` and :math:`t`, :math:`v` is an unused
        cluster in the forest, :math:`T=|v|+|s|+|t|`, and
        :math:`|*|` is the cardinality of its argument. This is also
        known as the incremental algorithm.

    Examples
    --------
    >>> import numpy as np
    >>> from sklearn.cluster import ward_tree
    >>> X = np.array([[1, 2], [1, 4], [1, 0],
    ...               [4, 2], [4, 4], [4, 0]])
    >>> children, n_connected_components, n_leaves, parents = ward_tree(X)
    >>> children
    array([[0, 1],
           [3, 5],
           [2, 6],
           [4, 7],
           [8, 9]])
    >>> n_connected_components
    1
    >>> n_leaves
    6
    """
    X = np.asarray(X)
    if X.ndim == 1:
        X = np.reshape(X, (-1, 1))
    n_samples, n_features = X.shape

    if connectivity is None:
        from scipy.cluster import hierarchy  # imports PIL

        if n_clusters is not None:
            warnings.warn(
                (
                    "Partial build of the tree is implemented "
                    "only for structured clustering (i.e. with "
                    "explicit connectivity). The algorithm "
                    "will build the full tree and only "
                    "retain the lower branches required "
                    "for the specified number of clusters"
                ),
                stacklevel=2,
            )
        X = np.require(X, requirements="W")
        out = hierarchy.ward(X)
        children_ = out[:, :2].astype(np.intp)

        if return_distance:
            distances = out[:, 2]
            return children_, 1, n_samples, None, distances
        else:
            return children_, 1, n_samples, None

    connectivity, n_connected_components = _fix_connectivity(
        X, connectivity, affinity="euclidean"
    )
    if n_clusters is None:
        n_nodes = 2 * n_samples - 1
    else:
        if n_clusters > n_samples:
            raise ValueError(
                "Cannot provide more clusters than samples. "
                "%i n_clusters was asked, and there are %i "
                "samples." % (n_clusters, n_samples)
            )
        n_nodes = 2 * n_samples - n_clusters

    # create inertia matrix
    coord_row = []
    coord_col = []
    A = []
    for ind, row in enumerate(connectivity.rows):
        A.append(row)
        # We keep only the upper triangular for the moments
        # Generator expressions are faster than arrays on the following
        row = [i for i in row if i < ind]
        coord_row.extend(
            len(row)
            * [
                ind,
            ]
        )
        coord_col.extend(row)

    coord_row = np.array(coord_row, dtype=np.intp, order="C")
    coord_col = np.array(coord_col, dtype=np.intp, order="C")

    # build moments as a list
    moments_1 = np.zeros(n_nodes, order="C")
    moments_1[:n_samples] = 1
    moments_2 = np.zeros((n_nodes, n_features), order="C")
    moments_2[:n_samples] = X
    inertia = np.empty(len(coord_row), dtype=np.float64, order="C")
    _hierarchical.compute_ward_dist(moments_1, moments_2, coord_row, coord_col, inertia)
    inertia = list(zip(inertia, coord_row, coord_col))
    heapify(inertia)

    # prepare the main fields
    parent = np.arange(n_nodes, dtype=np.intp)
    used_node = np.ones(n_nodes, dtype=bool)
    children = []
    if return_distance:
        distances = np.empty(n_nodes - n_samples)

    not_visited = np.empty(n_nodes, dtype=bool, order="C"

    # recursive merge loop
    for k in range(n_samples, n_nodes):
        # identify the merge
        while True:
            inert, i, j = heappop(inertia)
            if used_node[i] and used_node[j]:
                break
        parent[i], parent[j] = k, k
        children.append((i, j))
        used_node[i] = used_node[j] = False
        if return_distance:  # store inertia value
            distances[k - n_samples] = inert

        # update the moments
        moments_1[k] = moments_1[i] + moments_1[j]
        moments_2[k] = moments_2[i] + moments_2[j]

        # update the structure matrix A and the inertia matrix
        coord_col = []
        not_visited.fill(1)
        not_visited[k] = 0
        _hierarchical._get_parents(A[i], coord_col, parent, not_visited)
        _hierarchical._get_parents(A[j], coord_col, parent, not_visited)
        # List comprehension is faster than a for loop
        [A[col].append(k) for col in coord_col]
        A.append(coord_col)
        coord_col = np.array(coord_col, dtype=np.intp, order="C")
        coord_row = np.empty(coord_col.shape, dtype=np.intp, order="C")
        coord_row.fill(k)
        n_additions = len(coord_row)
        ini = np.empty(n_additions, dtype=np.float64, order="C")

        _hierarchical.compute_ward_dist(moments_1, moments_2, coord_row, coord_col, ini)

        # List comprehension is faster than a for loop
        [heappush(inertia, (ini[idx], k, coord_col[idx])) for idx in range(n_additions)]

    # Separate leaves in children (empty lists up to now)
    n_leaves = n_samples
    # sort children to get consistent output with unstructured version
    children = [c[::-1] for c in children]
    children = np.array(children)  # return numpy array for efficient caching

    if return_distance:
        # 2 is scaling factor to compare w/ unstructured version
        distances = np.sqrt(2.0 * distances)
        return children, n_connected_components, n_leaves, parent, distances
    else:
        return children, n_connected_components, n_leaves, parent

这段代码实现了 Ward 链接的层次聚类。Ward 方法通过最小化簇内方差增加来选择合并对，这是一种基于方差的准则，� Aid 于保持簇的紧凑性和球形。关键创新在于它维护了一阶矩（样本和）和二阶矩（平方和），从而可以增量更新簇间距离，而无需重新计算所有 pairewise 距离。当存在连接矩阵（结构化约束）时，它使用堆优化的策略来只考虑相连的簇对；否则，它回退到 SciPy 的实现。树切割函数 _hc_cut 使用最大堆（通过负索引实现）迭代移除距离最大的内部节点，直到达到目标簇数，这种方法比重新构建整棵树更高效。

单链接优化：最小生成树 + 并查集

源码路径：sklearn/cluster/_agglomerative.py - _single_linkage_tree（50-100行）

def _single_linkage_tree(
    connectivity,
    n_samples,
    n_nodes,
    n_clusters,
    n_connected_components,
    return_distance,
):
    """
    Perform single linkage clustering on sparse data via the minimum
    spanning tree from scipy.sparse.csgraph, then using union-find to label.
    The parent array is then generated by walking through the tree.
    """
    from scipy.sparse.csgraph import minimum_spanning_tree

    # explicitly cast connectivity to ensure safety
    connectivity = connectivity.astype(np.float64, copy=False)

    # Ensure zero distances aren't ignored by setting them to "epsilon"
    epsilon_value = np.finfo(dtype=connectivity.data.dtype).eps
    connectivity.data[connectivity.data == 0] = epsilon_value

    # Use scipy.sparse.csgraph to generate a minimum spanning tree
    mst = minimum_spanning_tree(connectivity.tocsr())

    # Convert the graph to scipy.cluster.hierarchy array format
    mst = mst.tocoo()

    # Undo the epsilon values
    mst.data[mst.data == epsilon_value] = 0

    mst_array = np.vstack([mst.row, mst.col, mst.data]).T

    # Sort edges of the min_spanning_tree by weight
    mst_array = mst_array[np.argsort(mst_array.T[2], kind="mergesort"), :]

    # Convert edge list into standard hierarchical clustering format
    single_linkage_tree = _hierarchical._single_linkage_label(mst_array)
    children_ = single_linkage_tree[:, :2].astype(int)

    # Compute parents
    parent = np.arange(n_nodes, dtype=np.intp)
    for i, (left, right) in enumerate(children_, n_samples):
        if n_clusters is not None and i >= n_nodes:
            break
        if left < n_nodes:
            parent[left] = i
        if right < n_nodes:
            parent[right] = i

    if return_distance:
        distances = single_linkage_tree[:, 2]
        return children_, n_connected_components, n_samples, parent, distances
    return children_, n_connected_components, n_samples, parent

这段代码展示了单链接聚类的高效实现。单链接定义为两个簇之间的最小距离，因此整个聚类过程等价于构建数据点的最小生成树（MST），然后通过断开 MST 中最长的边来形成簇。通过使用 scipy.sparse.csgraph.minimum_spanning_tree 来构建 MST（仅在非零边上操作，适用于稀疏数据），并结合并查集（Union-Find）数据结构来高效标记簇，该算法 achieves O(n log n) 时间复杂度（对于稀疏图），远优于 naive 的 O(n²) 实现。并查集通过路径压缩和按秩合并，使得查找和 union 操作近乎常数时间。

树切割与标签生成

源码路径：sklearn/cluster/_agglomerative.py - _hc_cut（450-480行）

def _hc_cut(n_clusters, children, n_leaves):
    """Function cutting the ward tree for a given number of clusters.

    Parameters
    ----------
    n_clusters : int or ndarray
        The number of clusters to form.

    children : ndarray of shape (n_nodes-1, 2)
        The children of each non-leaf node. Values less than `n_samples`
        correspond to leaves of the tree which are the original samples.
        A node `i` greater than or equal to `n_samples` is a non-leaf
        node and has children `children_[i - n_samples]`. Alternatively
        at the i-th iteration, children[i][0] and children[i][1]
        are merged to form node `n_samples + i`.

    n_leaves : int
        Number of leaves of the tree.

    Returns
    -------
    labels : array [n_samples]
        Cluster labels for each point.
    """
    if n_clusters > n_leaves:
        raise ValueError(
            "Cannot extract more clusters than samples: "
            f"{n_clusters} clusters were given for a tree with {n_leaves} leaves."
        )
    # In this function, we store nodes as a heap to avoid recomputing
    # the max of the nodes: the first element is always the smallest
    # We use negated indices as heaps work on smallest elements, and we
    # are interested in largest elements
    # children[-1] is the root of the tree
    nodes = [-(max(children[-1]) + 1)]
    for _ in range(n_clusters - 1):
        # As we have a heap, nodes[0] is the smallest element
        these_children = children[-nodes[0] - n_leaves]
        # Insert the 2 children and remove the largest node
        heappush(nodes, -these_children[0])
        heappushpop(nodes, -these_children[1])
    label = np.zeros(n_leaves, dtype=np.intp)
    for i, node in enumerate(nodes):
        label[_hierarchical._hc_get_descendent(-node, children, n_leaves)] = i
    return label

这段代码实现了层次聚类树的切割逻辑。它不从根开始构建，而是从根开始，使用最大堆（通过存储负索引来模拟）来追踪当前的内部节点。在每次迭代中，它弹出堆顶（当前最大的节点），将其两个子节点压入堆中，直到恰好剩下 n_clusters 个节点（即达到目标簇数）。这种方法避免了对整棵树进行遍历，只需处理 O(n_clusters log n) 个节点。辅助函数 _hc_get_descendent 则用于将内部节点映射回原始样本的标签，通过递归收集所有后代叶子节点。

22.7 DBSCAN 与 OPTICS —— 密度驱动聚类的“双子星”

DBSCAN 和 OPTICS 都是基于密度的聚类算法，它们能够发现任意形状的簇并将低密度区域的点标记为噪声。DBSCAN 使用全局密度阈值（eps 和 min_samples），而 OPTICS 则通过构建可达性图来提供多尺度的聚类结果。

DBSCAN 邻域查询与核心样本扩展

源码路径：sklearn/cluster/_dbscan.py - DBSCAN.fit（150-220行）

@_fit_context(
    # DBSCAN.metric is not validated yet
    prefer_skip_nested_validation=False
)
def fit(self, X, y=None, sample_weight=None):
    """Perform DBSCAN clustering from features, or distance matrix.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features), or \
            (n_samples, n_samples)
        Training instances to cluster, or distances between instances if
        ``metric='precomputed'``. If a sparse matrix is provided, it will
        be converted into a sparse ``csr_matrix``.

    y : Ignored
        Not used, present here for API consistency by convention.

    sample_weight : array-like of shape (n_samples,), default=None
        Weight of each sample, such that a sample with a weight of at least
        ``min_samples`` is by itself a core sample; a sample with a
        negative weight may inhibit its eps-neighbor from being core.
        Note that weights are absolute, and default to 1.

    Returns
    -------
    self : object
            Returns a fitted instance of self.
    """
    X = validate_data(self, X, accept_sparse="csr")

    if sample_weight is not None:
        sample_weight = _check_sample_weight(sample_weight, X)

    # Calculate neighborhood for all samples. This leaves the original
    # point in, which needs to be considered later (i.e. point i is in the
    # neighborhood of point i. While True, its useless information)
    if self.metric == "precomputed" and sparse.issparse(X):
        # set the diagonal to explicit values, as a point is its own
        # neighbor
        X = X.copy()  # copy to avoid in-place modification
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", sparse.SparseEfficiencyWarning)
            X.setdiag(X.diagonal())

    neighbors_model = NearestNeighbors(
        radius=self.eps,
        algorithm=self.algorithm,
        leaf_size=self.leaf_size,
        metric=self.metric,
        metric_params=self.metric_params,
        p=self.p,
        n_jobs=self.n_jobs,
    )
    neighbors_model.fit(X)
    # This has worst case O(n^2) memory complexity
    neighborhoods = neighbors_model.radius_neighbors(X, return_distance=False)

    if sample_weight is None:
        n_neighbors = np.array([len(neighbors) for neighbors in neighborhoods])
    else:
        n_neighbors = np.array(
            [np.sum(sample_weight[neighbors]) for neighbors in neighborhoods]
            )

    # Initially, all samples are noise.
    labels = np.full(X.shape[0], -1, dtype=np.intp)

    # A list of all core samples found.
    core_samples = np.asarray(n_neighbors >= self.min_samples, dtype=np.uint8)
    dbscan_inner(core_samples, neighborhoods, labels)

    self.core_sample_indices_ = np.where(core_samples)[0]
    self.labels_ = labels

    if len(self.core_sample_indices_):
        # fix for scipy sparse indexing issue
        self.components_ = X[self.core_sample_indices_].copy()
    else:
        # no core samples
        self.components_ = np.empty((0, X.shape[1]))
    return self

这段代码实现了 DBSCAN 的核心逻辑。它首先为每个样本查询其 eps 邻域（使用 NearestNeighbors.radius_neighbors），这一步在最坏情况下会导致 O(n²) 的内存复杂度，因为它可能需要存储所有样本之间的距离。然而，通过使用 sample_weight 支持加权核心判定（其中邻域的“大小”由权重和决定），它能够处理不同样本的重要性。然后，它将邻域点数（或权重和）大于等于 min_samples 的样本标记为核心样本。最后，它调用 Cython 实现的 dbscan_inner 函数来执行深度优先搜索（DFS）：从每个未访问的核心样本开始，递归地将所有可达的核心样本（以及直接连接的非核心边界点）标记为同一簇，而噪声点（非核心且不可达任何核心样本）保持为 -1。这种方法能够发现任意形状的簇，并且对噪声具有鲁棒性。

OPTICS 多尺度密度排序

源码路径：sklearn/cluster/_optics.py - compute_optics_graph（300-400行）

@validate_params(
    {
        "X": [np.ndarray, "sparse matrix"],
        "min_samples": [
            Interval(Integral, 2, None, closed="left"),
            Interval(RealNotInt, 0, 1, closed="both"),
        ],
        "max_eps": [Interval(Real, 0, None, closed="both")],
        "metric": [StrOptions(set(_VALID_METRICS) | {"precomputed"}), callable],
        "p": [Interval(Real, 0, None, closed="right"), None],
        "metric_params": [dict, None],
        "algorithm": [StrOptions({"auto", "brute", "ball_tree", "kd_tree"})],
        "leaf_size": [Interval(Integral, 1, None, closed="left")],
        "n_jobs": [Integral, None],
    },
    prefer_skip_nested_validation=False,  # metric is not validated yet
)
def compute_optics_graph(
    X, *, min_samples, max_eps, metric, p, metric_params, algorithm, leaf_size, n_jobs
):
    """Compute the OPTICS reachability graph.

    Read more in the :ref:`User Guide <optics>`.

    Parameters
    ----------
    X : {ndarray, sparse matrix} of shape (n_samples, n_features), or \
            (n_samples, n_samples) if metric='precomputed'
        A feature array, or array of distances between samples if
        metric='precomputed'.

    min_samples : int > 1 or float between 0 and 1
        The number of samples in a neighborhood for a point to be considered
        as a core point. Expressed as an absolute number or a fraction of the
        number of samples (rounded to be at least 2).

    max_eps : float, default=np.inf
        The maximum distance between two samples for one to be considered as
        in the neighborhood of the other. Default value of ``np.inf`` will
        identify clusters across all scales; reducing ``max_eps`` will result
        in shorter run times.

    metric : str or callable, default='minkowski'
        Metric to use for distance computation. Any metric from scikit-learn
        or scipy.spatial.distance can be used.

        If metric is a callable function, it is called on each
        pair of instances (rows) and the resulting value recorded. The callable
        should take two arrays as input and return one value indicating the
        distance between them. This works for Scipy's metrics, but is less
        efficient than passing the metric name as a string. If metric is
        "precomputed", X is assumed to be a distance matrix and must be square.

        Valid values for metric are:

        - from scikit-learn: ['cityblock', 'cosine', 'euclidean', 'l1', 'l2',
          'manhattan']

        - from scipy.spatial.distance: ['braycurtis', 'canberra', 'chebyshev',
          'correlation', 'dice', 'hamming', 'jaccard', 'kulsinski',
          'mahalanobis', 'minkowski', 'rogerstanimoto', 'russellrao',
          'seuclidean', 'sokalmichener', 'sokalsneath', 'sqeuclidean',
          'yule']

        See the documentation for scipy.spatial.distance for details on these
        metrics.

        .. note::
           `'kulsinski'` is deprecated from SciPy 1.9 and will be removed in SciPy 1.11.

    p : float, default=2
        Parameter for the Minkowski metric from
        :class:`~sklearn.metrics.pairwise_distances`. When p = 1, this is
        equivalent to using manhattan_distance (l1), and euclidean_distance
        (l2) for p = 2. For arbitrary p, minkowski_distance (l_p) is used.

    metric_params : dict, default=None
        Additional keyword arguments for the metric function.

    algorithm : {'auto', 'ball_tree', 'kd_tree', 'brute'}, default='auto'
        Algorithm used to compute the nearest neighbors:

        - 'ball_tree' will use :class:`~sklearn.neighbors.BallTree`.
        - 'kd_tree' will use :class:`~sklearn.neighbors.KDTree`.
        - 'brute' will use a brute-force search.
        - 'auto' (default) will attempt to decide the most appropriate
          algorithm based on the values passed to `fit` method. (default)

        Note: fitting on sparse input will override the setting of
        this parameter, using brute force.

    leaf_size : int, default=30
        Leaf size passed to :class:`~sklearn.neighbors.BallTree` or
        :class:`~sklearn.neighbors.KDTree`. This can affect the speed of the
        construction and query, as well as the memory required to store the
        tree. The optimal value depends on the nature of the problem.

    n_jobs : int, default=None
        The number of parallel jobs to run for neighbors search.
        ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
        ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
        for more details.

    Returns
    -------
    ordering_ : array of shape (n_samples,)
        The cluster ordered list of sample indices.

    core_distances_ : array of shape (n_samples,)
        Distance at which each sample becomes a core point, indexed by object
        order. Points which will never be core have a distance of inf. Use
        ``clust.core_distances_[clust.ordering_]`` to access in cluster order.

    reachability_ : array of shape (n_samples,)
        Reachability distances per sample, indexed by object order. Use
        ``clust.reachability_[clust.ordering_]`` to access in cluster order.

    predecessor_ : array of shape (n_samples,)
        Point that a sample was reached from, indexed by object order.
        Seed points have a predecessor of -1.

    References
    ----------
    .. [1] Ankerst, Mihael, Markus M. Breunig, Hans-Peter Kriegel,
       and Jörg Sander. "OPTICS: ordering points to identify the clustering
       structure." ACM SIGMOD Record 28, no. 2 (1999): 49-60.

    Examples
    --------
    >>> import numpy as np
    >>> from sklearn.cluster import compute_optics_graph
    >>> X = np.array([[1, 2], [2, 5], [3, 6],
    ...               [8, 7], [8, 8], [7, 3]])
    >>> ordering, core_distances, reachability, predecessor = compute_optics_graph(
    ...     X,
    ...     min_samples=2,
    ...     max_eps=np.inf,
    ...     metric="minkowski",
    ...     p=2,
    ...     metric_params=None,
    ...     algorithm="auto",
    ...     leaf_size=30,
    ...     n_jobs=None,
    ... )
    >>> ordering
    array([0, 1, 2, 5, 3, 4])
    >>> core_distances
    array([3.16, 1.41, 1.41, 1.        , 1.        ,
           4.12])
    >>> reachability
    array([       inf, 3.16, 1.41, 4.12, 1.        ,
           5.        ])
    >>> predecessor
    array([-1,  0,  1,  5,  3,  2])
    """
    n_samples = X.shape[0]
    _validate_size(min_samples, n_samples, "min_samples")
    if min_samples <= 1:
        min_samples = max(2, int(min_samples * n_samples))

    # Start all points as 'unprocessed' ##
    reachability_ = np.empty(n_samples)
    reachability_.fill(np.inf)
    predecessor_ = np.empty(n_samples, dtype=int)
    predecessor_.fill(-1)

    nbrs = NearestNeighbors(
        n_neighbors=min_samples,
        algorithm=algorithm,
        leaf_size=leaf_size,
        metric=metric,
        metric_params=metric_params,
        p=p,
        n_jobs=n_jobs,
    )

    nbrs.fit(X)
    # Here we first do a kNN query for each point, this differs from
    # the original OPTICS that only used epsilon range queries.
    # TODO: handle working_memory somehow?
    core_distances_ = _compute_core_distances_(
        X=X, neighbors=nbrs, min_samples=min_samples, working_memory=None
    )
    # OPTICS puts an upper limit on these, use inf for undefined.
    core_distances_[core_distances_ > max_eps] = np.inf
    np.around(
        core_distances_,
        decimals=np.finfo(core_distances_.dtype).precision,
        out=core_distances_,
    )

    # Main OPTICS loop. Not parallelizable. The order that entries are
    # written to the 'ordering_' list is important!
    # Note that this implementation is O(n^2) theoretically, but
    # supposedly with very low constant factors.
    processed = np.zeros(X.shape[0], dtype=bool)
    ordering = np.zeros(X.shape[0], dtype=int)
    for ordering_idx in range(X.shape[0]):
        # Choose next based on smallest reachability distance
        # (And prefer smaller ids on ties, possibly np.inf!)
        index = np.where(processed == 0)[0]
        point = index[np.argmin(reachability_[index])]

        processed[point] = True
        ordering[ordering_idx] = point
        if core_distances_[point] != np.inf:
            _set_reach_dist(
                core_distances_=core_distances_,
                reachability_=reachability_,
                predecessor_=predecessor_,
                point_index=point,
                processed=processed,
                X=X,
                nbrs=nbrs,
                metric=metric,
                metric_params=metric_params,
                p=p,
                max_eps=max_eps,
            )
    if np.all(np.isinf(reachability_)):
        warnings.warn(
            (
                "All reachability values are inf. Set a larger"
                " max_eps or all data will be considered outliers."
            ),
            UserWarning,
        )
    return ordering, core_distances_, reachability_, predecessor_

这段代码实现了 OPTICS 算法的核心——计算可达性图。与 DBSCAN 不同，OPTICS 不依赖于全局的 eps 参数，而是首先为每个样本计算其核心距离（第 min_samples 近邻的距离），这通过 kNN 查询高效完成。然后，它按照可达距离（reachability distance）从小到大的顺序处理点：对于每个点，它更新所有未处理邻域的可达距离为 max(核心距离, 实际距离)。这一步是按需计算的（仅当邻域未被处理时才计算距离），这使得算法在稀疏数据或需要精细控制时更具灵活性，但也导致了难以并行化的瓶颈，因为每个点的处理依赖于之前点的状态。可达距离的最终序列捕捉了数据的密度层次结构，后续可以通过不同的提取方法（如 DBSCAN 风格或 Xi 方法）获得不同样数目的聚类结果。

22.8 HDBSCAN —— 层次密度聚类的“巅峰之作”

HDBSCAN 将层次聚类的思想与基于密度的聚类相结合，通过构建互达距离图和最小生成树来处理不同密度的簇，并利用簇的稳定性（持续性）来选择有意义的簇，使其能够在变密度数据上表现出色。

互达距离与最小生成树构建

源码路径：sklearn/cluster/_hdbscan/_reachability.pyx - mutual_reachability_graph（30-100行）

def mutual_reachability_graph(
    distance_matrix, min_samples=5, max_distance=0.0
):
    """Compute the weighted adjacency matrix of the mutual reachability graph.

    The mutual reachability distance used to build the graph is defined as::

        max(d_core(x_p), d_core(x_q), d(x_p, x_q))

    and the core distance `d_core` is defined as the distance between a point
    `x_p` and its k-th nearest neighbor.

    Note that all computations are done in-place.

    Parameters
    ----------
    distance_matrix : {ndarray, sparse matrix} of shape (n_samples, n_samples)
        Array of distances between samples. If sparse, the array must be in
        `CSR` format.

    min_samples : int, default=5
        The parameter `k` used to calculate the distance between a point
        `x_p` and its k-th nearest neighbor.

    max_distance : float, default=0.0
        The distance which `np.inf` is replaced with. When the true mutual-
        reachability distance is measured to be infinite, it is instead
        truncated to `max_dist`. Only used when `distance_matrix` is a sparse
        matrix.

    Returns
    -------
    mututal_reachability_graph: {ndarray, sparse matrix} of shape \
            (n_samples, n_samples)
        Weighted adjacency matrix of the mutual reachability graph.

    References
    ----------
    .. [1] Campello, R. J., Moulavi, D., & Sander, J. (2013, April).
       Density-based clustering based on hierarchical density estimates.
       In Pacific-Asia Conference on Knowledge Discovery and Data Mining
       (pp. 160-172). Springer Berlin Heidelberg.
    """
    further_neighbor_idx = min_samples - 1
    if issparse(distance_matrix):
        if distance_matrix.format != "csr":
            raise ValueError(
                "Only sparse CSR matrices are supported for `distance_matrix`."
            )
        _sparse_mutual_reachability_graph(
            distance_matrix.data,
            distance_matrix.indices,
            distance_matrix.indptr,
            distance_matrix.shape[0],
            further_neighbor_idx=further_neighbor_idx,
            max_distance=max_distance,
        )
    else:
        _dense_mutual_reachability_graph(
            distance_matrix, further_neighbor_idx=further_neighbor_idx
        )
    return distance_matrix

这段代码实现了互达距离图的构建，这是 HDBSCAN 的核心创新。互达距离定义为两点之间距离与它们各自核心距离（到第 k 近邻的距离）的最大值：d_mutual-reachability(p, q) = max(核心距离(p), 核心距离(q), 欧氏距离(p, q))。这个度量在某种程度上“平滑”了局部密度的影响：在高密度区域，核心距离较小，因此互达距离更接近实际欧氏距离；在低密度区域，核心距离较大，因此即使两点实际很近，如果它们属于不同的低密度区域，它们的互达距离也会被放大。这种设计使得算法对密度变化不敏感，能够在同一运行中发现不同密度的簇。计算过程中，它首先为每个点计算核心距离（到第 k 近邻的距离），然后对距离矩阵的每个元素应用上述最大值规则。对于稀疏数据，它仅遍历 CSR 格式的非零元素以保持效率。

层级树压缩与稳定性评估

源码路径：sklearn/cluster/_hdbscan/_tree.pyx - _condense_tree（50-150行）

cpdef cnp.ndarray[CONDENSED_t, ndim=1, mode='c'] _condense_tree(
    const HIERARCHY_t[::1] hierarchy,
    cnp.intp_t min_cluster_size=10
):
    """Condense a tree according to a minimum cluster size. This is akin
    to the runt pruning procedure of Stuetzle. The result is a much simpler
    tree that is easier to visualize. We include extra information on the
    lambda value at which individual points depart clusters for later
    analysis and computation.

    Parameters
    ----------
    hierarchy : ndarray of shape (n_samples,), dtype=HIERARCHY_dtype
        A single linkage hierarchy in scipy.cluster.hierarchy format.

    min_cluster_size : int, optional (default 10)
        The minimum size of clusters to consider. Clusters smaller than this
        are pruned from the tree.

    Returns
    -------
    condensed_tree : ndarray of shape (n_samples,), dtype=CONDENSED_dtype
        Effectively an edgelist encoding a parent/child pair, along with a
        value and the corresponding cluster_size in each row providing a tree
        structure.
    """

    cdef:
        cnp.intp_t root = 2 * hierarchy.shape[0]
        cnp.intp_t n_samples = hierarchy.shape[0] + 1
        cnp.intp_t next_label = n_samples + 1
        list result_list, node_list = bfs_from_hierarchy(hierarchy, root)

        cnp.intp_t[::1] relabel
        cnp.uint8_t[::1] ignore

        cnp.intp_t node, sub_node, left, right
        cnp.float64_t lambda_value, distance
        cnp.intp_t left_count, right_count
        HIERARCHY_t children

    relabel = np.empty(root + 1, dtype=np.intp)
    relabel[root] = n_samples
    result_list = []
    ignore = np.zeros(len(node_list), dtype=bool)

    for node in node_list:
        if ignore[node] or node < n_samples:
            continue

        children = hierarchy[node - n_samples]
        left = children.left_node
        right = children.right_node
        distance = children.value
        if distance > 0.0:
            lambda_value = 1.0 / distance
        else:
            lambda_value = INFTY

        if left >= n_samples:
            left_count = hierarchy[left - n_samples].cluster_size
        else:
            left_count = 1

        if right >= n_samples:
            right_count = hierarchy[right - n_samples].cluster_size
        else:
            right_count = 1

        if left_count >= min_cluster_size and right_count >= min_cluster_size:
            relabel[left] = next_label
            next_label += 1
            result_list.append(
                (relabel[node], relabel[left], lambda_value, left_count)
            )

            relabel[right] = next_label
            next_label += 1
            result_list.append(
                (relabel[node], relabel[right], lambda_value, right_count)
            )

        elif left_count < min_cluster_size and right_count < min_cluster_size:
            for sub_node in bfs_from_hierarchy(hierarchy, left):
                if sub_node < n_samples:
                    result_list.append(
                        (relabel[node], sub_node, lambda_value, 1)
                    )
                ignore[sub_node] = True

            for sub_node in bfs_from_hierarchy(hierarchy, right):
                if sub_node < n_samples:
                    result_list.append(
                        (relabel[node], sub_node, lambda_value, 1)
                    )
                ignore[sub_node] = True

        elif left_count < min_cluster_size:
            relabel[right] = relabel[node]
            for sub_node in bfs_from_hierarchy(hierarchy, left):
                if sub_node < n_samples:
                    result_list.append(
                        (relabel[node], sub_node, lambda_value, 1)
                    )
                ignore[sub_node] = True

        else:
            relabel[left] = relabel[node]
            for sub_node in bfs_from_hierarchy(hierarchy, right):
                if sub_node < n_samples:
                    result_list.append(
                        (relabel[node], sub_node, lambda_value, 1)
                    )
                ignore[sub_node] = True

    return np.array(result_list, dtype=CONDENSED_dtype)

这段代码实现了层级树的压缩（condensing）过程，这是 HDBSCAN 稳定性分析的关键步骤。它从单链接树（由 MST 构建）开始，广度优先遍历树结构。对于每个内部节点（代表一个潜在的簇合并事件），它检查其两个子簇的大小：只有当两个子簇都大于等于 min_cluster_size 时，才保留这个合并事件；否则，它将较小的子簇的所有叶子节点直接连接到父节点上，实际上是跳过了这个无意义的合并（因为产生的簇太小而不可靠）。在保留的合并事件中，它记录了一个 lambda 值（定义为 1/合并距离），这个值反过来代表了簇在层级结构中的“诞生时间” — — 越小的合并距离意味着在更高的密度层次（更小的 lambda 值）就形成了簇。这个 lambda 值后来被用于计算簇的稳定性。

源码路径：sklearn/cluster/_hdbscan/_tree.pyx - _compute_stability（150-200行）

cdef dict _compute_stability(
    cnp.ndarray[CONDENSED_t, ndim=1, mode='c'] condensed_tree
):

    cdef:
        cnp.float64_t[::1] result, births
        cnp.intp_t[:] parents = condensed_tree['parent']

        cnp.intp_t parent, cluster_size, result_index, idx
        cnp.float64_t lambda_val
        CONDENSED_t condensed_node
        cnp.intp_t largest_child = condensed_tree['child'].max()
        cnp.intp_t smallest_cluster = np.min(parents)
        cnp.intp_t num_clusters = np.max(parents) - smallest_cluster + 1
        dict stability_dict = {}

    largest_child = max(largest_child, smallest_cluster)
    births = np.full(largest_child + 1, np.nan, dtype=np.float64)

    for idx in range(PyArray_SHAPE(<cnp.PyArrayObject*> condensed_tree)[0]):
        condensed_node = condensed_tree[idx]
        births[condensed_node.child] = condensed_node.value

    births[smallest_cluster] = 0.0

    result = np.zeros(num_clusters, dtype=np.float64)
    for idx in range(PyArray_SHAPE(<cnp.PyArrayObject*> condensed_tree)[0]):
        condensed_node = condensed_tree[idx]
        parent = condensed_node.parent
        lambda_val = condensed_node.value
        cluster_size = condensed_node.cluster_size

        result_index = parent - smallest_cluster
        result[result_index] += (lambda_val - births[parent]) * cluster_size

    for idx in range(num_clusters):
        stability_dict[idx + smallest_cluster] = result[idx]

    return stability_dict

这段代码计算每个簇的稳定性，这是 HDBSCAN 选择最终簇的依据。稳定性定义为簇在层级树中“存活”的时间长度，通过在簇的整个生命周期上积分 lambda 值（即 1/距离）来衡量。具体来说，对于一个簇，其稳定性等于所有在它之下（在树中更深层）子簇的 lambda 值之和，减去父簇的 lambda 值乘以子簇大小，再乘以簇自身的大小：稳定性 = Σ (lambda_子簇 - lambda_父簇) × 大小_子簇。在代码中，这是通过遍历压缩树中的每个节点（代表一个簇），并�加 (当前节点的 lambda - 父节点的 lambda) × 当前簇大小 来实现的。直观地说，这衡量了簇在不同密度阈值下的持久性：一个真正稳定的簇应该在较大的密度变化范围内（对应较大的 lambda 值范围）保持存在，而临时的或噪声引起的簇则只在很窄的密度范围内出现。这种基于持续性的簇选择比简单的大小或密度阈值更具理论依据和鲁棒性。

22.9 谱聚类与亲和传播 —— 图视角下的“社区发现”

谱聚类通过将数据嵌入到低维特征空间（基于相似度图的拉普拉斯矩阵），然后在这些空间中应用传统的聚类算法（如 K-Means）来发现非凸簇。亲和传播则通过在样本之间传递“责任”和“可用性”信息来自动识别 exemplar（簇的代表点）。

谱聚类嵌入与离散化

源码路径：sklearn/cluster/_spectral.py - SpectralClustering.fit（400-500行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
    """Perform spectral clustering from features, or affinity matrix.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features) or \
                (n_samples, n_samples)
            Training instances to cluster, similarities / affinities between
            instances if ``affinity='precomputed'``, or distances between
            instances if ``affinity='precomputed_nearest_neighbors``. If a
            sparse matrix is provided in a format other than ``csr_matrix``,
            ``csc_matrix``, or ``coo_matrix``, it will be converted into a
            sparse ``csr_matrix``.

    y : Ignored
        Not used, present here for API consistency by convention.

    Returns
    -------
    self : object
        A fitted instance of the estimator.
    """
    X = validate_data(
        self,
        X,
        accept_sparse=["csr", "csc", "coo"],
        dtype=np.float64,
        ensure_min_samples=2,
    )
    allow_squared = self.affinity in [
        "precomputed",
        "precomputed_nearest_neighbors",
    ]
    if X.shape[0] == X.shape[1] and not allow_squared:
        warnings.warn(
            "The spectral clustering API has changed. ``fit``"
            "now constructs an affinity matrix from data. To use"
            " a custom affinity matrix, "
            "set ``affinity=precomputed``."
        )

    if self.affinity == "nearest_neighbors":
        connectivity = kneighbors_graph(
            X, n_neighbors=self.n_neighbors, include_self=True, n_jobs=self.n_jobs
            )
        self.affinity_matrix_ = 0.5 * (connectivity + connectivity.T)
    elif self.affinity == "precomputed_nearest_neighbors":
        estimator = NearestNeighbors(
            n_neighbors=self.n_neighbors, n_jobs=self.n_jobs, metric="precomputed"
            ).fit(X)
        connectivity = estimator.kneighbors_graph(X=X, mode="connectivity")
        self.affinity_matrix_ = 0.5 * (connectivity + connectivity.T)
    elif self.affinity == "precomputed":
        self.affinity_matrix_ = X
    else:
        params = self.kernel_params
        if params is None:
            params = {}
        if not callable(self.affinity):
            params["gamma"] = self.gamma
            params["degree"] = self.degree
            params["coef0"] = self.coef0
        self.affinity_matrix_ = pairwise_kernels(
            X, metric=self.affinity, filter_params=True, **params
        )

    random_state = check_random_state(self.random_state)
    n_components = (
        self.n_clusters if self.n_components is None else self.n_components
    )
    # We now obtain the real valued solution matrix to the
    # relaxed Ncut problem, solving the eigenvalue problem
    # L_sym x = lambda x  and recovering u = D^-1/2 x.
    # The first eigenvector is constant only for fully connected graphs
    # and should be kept for spectral clustering (drop_first = False)
    # See spectral_embedding documentation.
    maps = _spectral_embedding(
        self.affinity_matrix_,
        n_components=n_components,
        eigen_solver=self.eigen_solver,
        random_state=random_state,
        eigen_tol=self.eigen_tol,
        drop_first=False,
    )
    if self.verbose:
        print(f"Computing label assignment using {self.assign_labels}")

    if self.assign_labels == "kmeans":
        _, self.labels_, _ = k_means(
            maps,
            self.n_clusters,
            random_state=random_state,
            n_init=self.n_init,
            verbose=self.verbose,
        )
    elif self.assign_labels == "cluster_qr":
        self.labels_ = cluster_qr(maps)
    else:
        self.labels_ = discretize(maps, random_state=random_state)

    return self

这段代码实现了谱聚类的完整流程。首先，它根据指定的亲和度构建方法（如 RBF 核、最近邻或预计算矩阵）构建亲和度矩阵。然后，它从这个亲和度矩阵计算归一化拉普拉斯矩阵（通过 _spectral_embedding 函数），并提取其前 k 个特征向量（其中 k 是期望的簇数），这 k 维嵌入空间试图将非凸簇线性可分。最后，它在这些特征向量上应用聚类算法：可以是 K-Means（默认）、离散化迭代搜索（通过旋转特征向量寻找接近离散分区的旋转矩阵），或直接使用 QR 分解（cluster_qr 函数）从特征向量中提取离散标签。QR 方法特别有趣，因为它利用了特征向量矩阵的 QR 分解，其中列主元（pivoting）的位置直接给出了簇的分配，无需迭代，这使得它在某些情况下比 K-Means 更快且更鲁棒于初始化。

亲和传播消息传递机制

源码路径：sklearn/cluster/_affinity_propagation.py - _affinity_propagation（30-120行）

def _affinity_propagation(
    S,
    *,
    preference,
    convergence_iter,
    max_iter,
    damping,
    verbose,
    return_n_iter,
    random_state,
):
    """Main affinity propagation algorithm."""
    n_samples = S.shape[0]
    if n_samples == 1 or _equal_similarities_and_preferences(S, preference):
        # It makes no sense to run the algorithm in this case, so return 1 or
        # n_samples clusters, depending on preferences
        warnings.warn(
            "All samples have mutually equal similarities. "
            "Returning arbitrary cluster center(s)."
        )
        if preference.flat[0] > S.flat[n_samples - 1]:
            return (
                (np.arange(n_samples), np.arange(n_samples), 0)
                if return_n_iter
                else (np.arange(n_samples), np.arange(n_samples))
            )
        else:
            return (
                (np.array([0]), np.array([0] * n_samples), 0)
                if return_n_iter
                else (np.array([0]), np.array([0] * n_samples))
            )

    # Place preference on the diagonal of S
    S.flat[:: (n_samples + 1)] = preference

    A = np.zeros((n_samples, n_samples))
    R = np.zeros((n_samples, n_samples))  # Initialize messages
    # Intermediate results
    tmp = np.zeros((n_samples, n_samples))

    # Remove degeneracies
    S += (
        np.finfo(S.dtype).eps * S + np.finfo(S.dtype).tiny * 100
    ) * random_state.standard_normal(size=(n_samples, n_samples))

    # Execute parallel affinity propagation updates
    e = np.zeros((n_samples, convergence_iter))

    ind = np.arange(n_samples)

    for it in range(max_iter):
        # tmp = A + S; compute responsibilities
        np.add(A, S, tmp)
        I = np.argmax(tmp, axis=1)
        Y = tmp[ind, I]  # np.max(A + S, axis=1)
        tmp[ind, I] = -np.inf
        Y2 = np.max(tmp, axis=1)

        # tmp = Rnew
        np.subtract(S, Y[:, None], tmp)
        tmp[ind, I] = S[ind, I] - Y2

        # Damping
        tmp *= 1 - damping
        R *= damping
        R += tmp

        # tmp = Rp; compute availabilities
        np.maximum(R, 0, out=tmp)
        tmp.flat[:: n_samples + 1] = R.flat[:: n_samples + 1]

        # tmp = -Anew
        tmp -= np.sum(tmp, axis=0)
        dA = np.diag(tmp).copy()
        tmp.clip(0, np.inf, tmp)
        tmp.flat[:: n_samples + 1] = dA

        # Damping
        tmp *= 1 - damping
        A *= damping
        A -= tmp

        # Check for convergence
        E = (np.diag(A) + np.diag(R)) > 0
        e[:, it % convergence_iter] = E
        K = np.sum(E, axis=0)

        if it >= convergence_iter:
            se = np.sum(e, axis,axis=1)
            unconverged = np.sum((se == convergence_iter) + (se == 0)) != n_samples
            if (not unconverged and (K > 0)) or (it == max_iter):
                never_converged = False
                if verbose:
                    print("Converged after %d iterations." % it)
                break
    else:
        never_converged = True
        if verbose:
            print("Did not converge")

    I = np.flatnonzero(E)
    K = I.size  # Identify exemplars

    if K > 0:
        if never_converged:
            warnings.warn(
                (
                    "Affinity propagation did not converge, this model "
                    "may return degenerate cluster centers and labels."
                ),
                ConvergenceWarning,
            )
        c = np.argmax(S[:, I], axis=1)
        c[I] = np.arange(K)  # Identify clusters
        # Refine the final set of exemplars and clusters and return results
        for k in range(K):
            ii = np.asarray(c == k).nonzero()[0]
            j = np.argmax(np.sum(S[ii[:, np.newaxis], ii], axis=0))
            I[k] = ii[j]

        c = np.argmax(S[:, I], axis=1)
        c[I] = np.arange(K)
        labels = I[c]
        # Reduce labels to a sorted, gapless, list
        cluster_centers_indices = np.unique(labels)
        labels = np.searchsorted(cluster_centers_indices, labels)
    else:
        warnings.warn(
            (
                "Affinity propagation did not converge and this model "
                "will not have any cluster centers."
            ),
            ConvergenceWarning,
        )
        labels = np.array([-1] * n_samples)
        cluster_centers_indices = []

    if return_n_iter:
        return cluster_centers_indices, labels, it + 1
    else:
        return cluster_centers_indices, labels