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

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

• OvR：把每个类别视作 正例，其他类别视作 负例，得到 n_classes 个二分类问题；随后使用 _binary_roc_auc_score 计算每个二分类 AUC，并根据 average 进行加权（macro、weighted、micro）。

• Ovo：遍历 所有类别两两组合（n_classes * (n_classes-1) / 2），对每一对使用二分类 AUC（不支持样本权重），再对所有 pair 的 AUC 加权平均（macro 或 weighted）。由于每对的正负样本比例相同，average='micro' 没有意义，因此仅实现 macro / weighted。

32.5.4 average_precision_score (50-130) —— AP 的统一入口

def average_precision_score(
    y_true, y_score, *, average="macro", pos_label=1, sample_weight=None
):
    """Compute average precision (AP) from prediction scores."""
    xp, _, device = get_namespace_and_device(y_score)
    y_true, sample_weight = move_to(y_true, sample_weight, xp=xp, device=device)

    if sample_weight is not None:
        sample_weight = column_or_1d(sample_weight)

    def _binary_uninterpolated_average_precision(
        y_true,
        y_score,
        pos_label=1,
        sample_weight=None,
        xp=xp,
    ):
        precision, recall, _ = precision_recall_curve(
            y_true,
            y_score,
            pos_label=pos_label,
            sample_weight=sample_weight,
        )
        # Return the step function integral
        # The following works because the last entry of precision is
        # guaranteed to be 1, as returned by precision_recall_curve.
        # Due to numerical error, we can get `-0.0` and we therefore clip it.
        return float(max(0.0, -xp.sum(xp.diff(recall) * precision[:-1])))

    y_type = type_of_target(y_true, input_name="y_true")
    present_labels = xp.unique_values(y_true)

    if y_type == "binary":
        if present_labels.shape[0] == 2 and pos_label not in present_labels:
            raise ValueError(
                f"pos_label={pos_label} is not a valid label. It should be "
                f"one of {present_labels}"
            )

    elif y_type == "multilabel-indicator" and pos_label != 1:
        raise ValueError(
            "Parameter pos_label is fixed to 1 for multilabel-indicator y_true. "
            "Do not set pos_label or set pos_label to 1."
        )

    elif y_type == "multiclass":
        if pos_label != 1:
            raise ValueError(
                "Parameter pos_label is fixed to 1 for multiclass y_true. "
                "Do not set pos_label or set pos_label to 1."
            )
        y_true = label_binarize(y_true, classes=present_labels)
        if not y_score.shape == y_true.shape:
            raise ValueError(
                "`y_score` needs to be of shape `(n_samples, n_classes)`, since "
                "`y_true` contains multiple classes. Got "
                f"`y_score.shape={y_score.shape}`."
            )

    average_precision = partial(
        _binary_uninterpolated_average_precision, pos_label=pos_label, xp=xp
    )
    return _average_binary_score(
        average_precision, y_true, y_score, average, sample_weight=sample_weight
    )

解释说明：

average_precision_score 通过 步函数积分（非线性插值）计算 AP，即 Σ Δrecall × precision。这种方式对 稀疏正例 更加保守，避免在极端不平衡数据上产生过高估计。函数内部的 _binary_uninterpolated_average_precision 复用 precision_recall_curve 获取精度-召回对，再通过加权求和得到单标签 AP。多分类情况下 强制使用 OvR（label_binarize），因为 Ovo 需要每对类别都有明确的二分类概率分布，这在 predict_proba 输出中往往不可得。_average_binary_score 统一处理 micro/macro/weighted/samples 等平均策略，并支持样本权重。

32.5.5 auc (30-60) —— 通用梯形积分工具

@validate_params(
    {"x": ["array-like"], "y": ["array-like"]},
    prefer_skip_nested_validation=True,
)
def auc(x, y):
    """Compute Area Under the Curve (AUC) using the trapezoidal rule."""
    check_consistent_length(x, y)
    x, y = column_or_1d(x), column_or_1d(y)

    if x.shape[0] < 2:
        raise ValueError(
            "At least 2 points are needed to compute area under curve, but x.shape = %s"
            % x.shape
        )

    direction = 1
    dx = np.diff(x)
    if np.any(dx < 0):
        if np.all(dx <= 0):
            direction = -1      # x 递减时翻转符号
        else:
            raise ValueError("x is neither increasing nor decreasing : {}.".format(x))

    area = direction * trapezoid(y, x)
    if isinstance(area, np.memmap):
        area = area.dtype.type(area)
    return float(area)

为什么不直接在 roc_auc_score 中使用 auc？

auc 只完成 梯形积分，而 ROC AUC 需要 部分 AUC（max_fpr）、McClish 校正、多分类 OvR/Ovo 分派 等业务层面的细节，这些都在 _binary_roc_auc_score 和 _multiclass_roc_auc_score 中实现。把 auc 设为底层工具有助于 代码复用（也用于 PR 曲线的 average_precision_score 计算）。

graph TD A[roc_auc_score] --> B{目标类型?} B -->|binary| C[_average_binary_score + _binary_roc_auc_score] B -->|multiclass| D[_multiclass_roc_auc_score] B -->|multilabel| C D --> E{multi_class?} E -->|ovr| F[label_binarize + _average_binary_score] E -->|ovo| G[_average_multiclass_ovo_score] C --> H[roc_curve -> auc] F --> H G --> H

---

32.6 多标签排名指标 —— LRAP、覆盖率误差、标签排名损失

32.6.1 label_ranking_average_precision_score (550-620) (LRAP)

def label_ranking_average_precision_score(y_true, y_score, *, sample_weight=None):
    """Compute ranking-based average precision."""
    check_consistent_length(y_true, y_score, sample_weight)
    y_true = check_array(y_true, ensure_2d=False, accept_sparse="csr")
    y_score = check_array(y_score, ensure_2d=False)

    if y_true.shape != y_score.shape:
        raise ValueError("y_true and y_score have different shape")

    # Handle badly formatted array and the degenerate case with one label
    y_type = type_of_target(y_true, input_name="y_true")
    if y_type != "multilabel-indicator" and not (
        y_type == "binary" and y_true.ndim == 2
    ):
        raise ValueError("{0} format is not supported".format(y_type))

    if not issparse(y_true):
        y_true = csr_matrix(y_true)

    y_score = -y_score

    n_samples, n_labels = y_true.shape

    out = 0.0
    for i, (start, stop) in enumerate(zip(y_true.indptr, y_true.indptr[1:])):
        relevant = y_true.indices[start:stop]

        if relevant.size == 0 or relevant.size == n_labels:
            # If all labels are relevant or unrelevant, the score is also
            # equal to 1. The label ranking has no meaning.
            aux = 1.0
        else:
            scores_i = y_score[i]
            rank = rankdata(scores_i, "max")[relevant]
            L = rankdata(scores_i[relevant], "max")
            aux = (L / rank).mean()

        if sample_weight is not None:
            aux = aux * sample_weight[i]
        out += aux

    if sample_weight is None:
        out /= n_samples
    else:
        out /= np.sum(sample_weight)

    return float(out)

解释说明：

LRAP 衡量 每个真标签在其排名前方出现的真标签比例 的平均值。实现上，将 y_score 取反后使用 scipy.stats.rankdata(method='max') 获取 最大排名（即并列分数取最差名次），对每个样本的真标签索引 relevant 提取其排名 rank，再对这些真标签自身的排名 L 求平均。稀疏矩阵（CSR）遍历仅访问非零标签，时间复杂度为 O(nnz)，避免了 O(n_samples × n_labels) 的全矩阵开销。全 0 或全 1 样本直接计为 1，因为此时排序无意义。

32.6.2 coverage_error (620-670)

def coverage_error(y_true, y_score, *, sample_weight=None):
    """Coverage error measure."""
    y_true = check_array(y_true, ensure_2d=True)
    y_score = check_array(y_score, ensure_2d=True)
    check_consistent_length(y_true, y_score, sample_weight)

    y_type = type_of_target(y_true, input_name="y_true")
    if y_type != "multilabel-indicator":
        raise ValueError("{0} format is not supported".format(y_type))

    if y_true.shape != y_score.shape:
        raise ValueError("y_true and y_score have different shape")

    y_score_mask = np.ma.masked_array(y_score, mask=np.logical_not(y_true))
    y_min_relevant = y_score_mask.min(axis=1).reshape((-1, 1))
    coverage = (y_score >= y_min_relevant).sum(axis=1)
    coverage = coverage.filled(0)

    return float(np.average(coverage, weights=sample_weight))

解释说明：

覆盖率误差的核心思想是：对每个样本，找到 最低的真标签分数（y_min_relevant），统计 所有分数 ≥ 该阈值的标签数量。这相当于“为了覆盖所有真标签，必须检索到第几位”。最优值等于平均真标签数。实现利用 NumPy 掩码数组（masked_array）屏蔽非真标签，min(axis=1) 快速得到阈值，随后广播比较并求和，最后加权平均。

32.6.3 label_ranking_loss (670-730)

def label_ranking_loss(y_true, y_score, *, sample_weight=None):
    """Compute Ranking loss measure."""
    y_true = check_array(y_true, ensure_2d=False, accept_sparse="csr")
    y_score = check_array(y_score, ensure_2d=False)
    check_consistent_length(y_true, y_score, sample_weight)

    y_type = type_of_target(y_true, input_name="y_true")
    if y_type not in ("multilabel-indicator",):
        raise ValueError("{0} format is not supported".format(y_type))

    if y_true.shape != y_score.shape:
        raise ValueError("y_true and y_score have different shape")

    n_samples, n_labels = y_true.shape

    y_true = csr_matrix(y_true)

    loss = np.zeros(n_samples)
    for i, (start, stop) in enumerate(zip(y_true.indptr, y_true.indptr[1:])):
        # Sort and bin the label scores
        unique_scores, unique_inverse = np.unique(y_score[i], return_inverse=True)
        true_at_reversed_rank = np.bincount(
            unique_inverse[y_true.indices[start:stop]], minlength=len(unique_scores)
        )
        all_at_reversed_rank = np.bincount(unique_inverse, minlength=len(unique_scores))
        false_at_reversed_rank = all_at_reversed_rank - true_at_reversed_rank

        # if the scores are ordered, it's possible to count the number of
        # incorrectly ordered paires in linear time by cumulatively counting
        # how many false labels of a given score have a score higher than the
        # accumulated true labels with lower score.
        loss[i] = np.dot(true_at_reversed_rank.cumsum(), false_at_reversed_rank)

    n_positives = count_nonzero(y_true, axis=1)
    with np.errstate(divide="ignore", invalid="ignore"):
        loss /= (n_labels - n_positives) * n_positives

    # When there is no positive or no negative labels, those values should
    # be consider as correct, i.e. the ranking doesn't matter.
    loss[np.logical_or(n_positives == 0, n_positives == n_labels)] = 0.0

    return float(np.average(loss, weights=sample_weight))

解释说明：

标签排名损失统计 错误排序的标签对数（真标签排在假标签后面的对数），再除以 可能出现的最大对数（n_positives * (n_labels - n_positives)），得到归一化的排序错误率。算法对每个样本的分数做 唯一化（np.unique + return_inverse），将标签按分数分桶；在每个分数桶上统计真标签数 true_at_reversed_rank 与假标签数 false_at_reversed_rank。通过 累积真标签数 × 当前桶假标签数 的点积，在 线性时间 完成所有错误配对计数。全正/全负样本损失设为 0。

graph TD A[多标签输入 y_true, y_score] --> B[转为 CSR 稀疏矩阵] B --> C[逐样本遍历非零标签] C --> D1[LRAP: rankdata 计算排名 -> L/rank 平均] C --> D2[Coverage: 掩码数组 min -> 广播比较求和] C --> D3[Ranking Loss: np.unique 分桶 -> cumsum 点积] D1 --> E[加权平均得标量指标] D2 --> E D3 --> E

---

32.7 排序质量指标：DCG/NDCG 与 Top‑k 准确率 —— 信息检索的“黄金标准”

32.7.1 _dcg_sample_scores (480-530) 与 dcg_score (730-800)

def _dcg_sample_scores(y_true, y_score, k=None, log_base=2, ignore_ties=False):
    """Compute Discounted Cumulative Gain."""
    discount = 1 / (np.log(np.arange(y_true.shape[1]) + 2) / np.log(log_base))
    if k is not None:
        discount[k:] = 0
    if ignore_ties:
        ranking = np.argsort(y_score)[:, ::-1]
        ranked = y_true[np.arange(ranking.shape[0])[:, np.newaxis], ranking]
        cumulative_gains = discount.dot(ranked.T)
    else:
        discount_cumsum = np.cumsum(discount)
        cumulative_gains = [
            _tie_averaged_dcg(y_t, y_s, discount_cumsum)
            for y_t, y_s in zip(y_true, y_score)
        ]
        cumulative_gains = np.asarray(cumulative_gains)
    return cumulative_gains

def dcg_score(
    y_true, y_score, *, k=None, log_base=2, sample_weight=None, ignore_ties=False
):
    """Compute Discounted Cumulative Gain."""
    y_true = check_array(y_true, ensure_2d=False)
    y_score = check_array(y_score, ensure_2d=False)
    check_consistent_length(y_true, y_score, sample_weight)
    _check_dcg_target_type(y_true)
    return float(
        np.average(
            _dcg_sample_scores(
                y_true, y_score, k=k, log_base=log_base, ignore_ties=ignore_ties
            ),
            weights=sample_weight,
        )
    )

• 核心逻辑：预计算对数折扣 1 / log(i+1)，k 截断后折扣置零。ignore_ties=True 时直接按分数排序求点积；否则调用 _tie_averaged_dcg 处理并列分数。

32.7.2 _tie_averaged_dcg (530-570)

def _tie_averaged_dcg(y_true, y_score, discount_cumsum):
    """
    Compute DCG by averaging over possible permutations of ties.
    """
    _, inv, counts = np.unique(-y_score, return_inverse=True, return_counts=True)
    ranked = np.zeros(len(counts))
    np.add.at(ranked, inv, y_true)
    ranked /= counts
    groups = np.cumsum(counts) - 1
    discount_sums = np.empty(len(counts))
    discount_sums[0] = discount_cumsum[groups[0]]
    discount_sums[1:] = np.diff(discount_cumsum[groups])
    return (ranked * discount_sums).sum()

解释说明：

当若干文档的预测分数相同（并列），所有可能的排列对最终 DCG 的贡献是 等价的。_tie_averaged_dcg 将同一组的真实增益取平均，然后乘以该组对应的 折扣之和，等价于对所有排列的 期望值，避免对任意排列产生偏差。np.unique(-y_score, ...) 以分数降序分组，np.add.at 聚合增益，discount_cumsum 预计算累积折扣，groups 标记每组末尾索引，diff 得到组内折扣和。

32.7.3 _ndcg_sample_scores (570-620) 与 ndcg_score (800-870)

def _ndcg_sample_scores(y_true, y_score, k=None, ignore_ties=False):
    """Compute Normalized Discounted Cumulative Gain."""
    gain = _dcg_sample_scores(y_true, y_score, k, ignore_ties=ignore_ties)
    # Here we use the order induced by y_true so we can ignore ties since
    # the gain associated to tied indices is the same (permuting ties doesn't
    # change the value of the re-ordered y_true)
    normalizing_gain = _dcg_sample_scores(y_true, y_true, k, ignore_ties=True)
    all_irrelevant = normalizing_gain == 0
    gain[all_irrelevant] = 0
    gain[~all_irrelevant] /= normalizing_gain[~all_irrelevant]
    return gain

def ndcg_score(y_true, y_score, *, k=None, sample_weight=None, ignore_ties=False):
    """Compute Normalized Discounted Cumulative Gain."""
    y_true = check_array(y_true, ensure_2d=False)
    y_score = check_array(y_score, ensure_2d=False)
    check_consistent_length(y_true, y_score, sample_weight)

    if y_true.min() < 0:
        raise ValueError("ndcg_score should not be used on negative y_true values.")
    if y_true.ndim > 1 and y_true.shape[1] <= 1:
        raise ValueError(
            "Computing NDCG is only meaningful when there is more than 1 document. "
            f"Got {y_true.shape[1]} instead."
        )
    _check_dcg_target_type(y_true)
    gain = _ndcg_sample_scores(y_true, y_score, k=k, ignore_ties=ignore_ties)
    return float(np.average(gain, weights=sample_weight))

解释说明：

NDCG = DCG / IDCG（Ideal DCG）。IDCG 通过将 y_true 自身作为排序分数计算 DCG（ignore_ties=True 因为真实相关度相同时排列不影响结果）。全无相关文档的样本（normalizing_gain == 0）NDCG 设为 0。ndcg_score 额外检查 y_true 非负约束，因为负相关度会导致 DCG 为负，破坏 [0,1] 归一化解释。

32.7.4 top_k_accuracy_score (870-950)

def top_k_accuracy_score(
    y_true, y_score, *, k=2, normalize=True, sample_weight=None, labels=None
):
    """Top-k Accuracy classification score."""
    y_true = check_array(y_true, ensure_2d=False, dtype=None)
    y_true = column_or_1d(y_true)
    y_type = type_of_target(y_true, input_name="y_true")
    if y_type == "binary" and labels is not None and len(labels) > 2:
        y_type = "multiclass"
    if y_type not in {"binary", "multiclass"}:
        raise ValueError(
            f"y type must be 'binary' or 'multiclass', got '{y_type}' instead."
        )
    y_score = check_array(y_score, ensure_2d=False)
    if y_type == "binary":
        if y_score.ndim == 2 and y_score.shape[1] != 1:
            raise ValueError(
                "`y_true` is binary while y_score is 2d with"
                f" {y_score.shape[1]} classes. If `y_true` does not contain all the"
                " labels, `labels` must be provided."
            )
        y_score = column_or_1d(y_score)
    else:
        if not y_score.ndim == 2:
            raise ValueError(
                "`y_score` needs to be of shape `(n_samples, n_classes)`, since "
                "`y_true` contains multiple classes. Got "
                f"`y_score.shape={y_score.shape}`."
            )

    check_consistent_length(y_true, y_score, sample_weight)
    y_score_n_classes = y_score.shape[1] if y_score.ndim == 2 else 2

    if labels is None:
        classes = _unique(y_true)
        n_classes = len(classes)

        if n_classes != y_score_n_classes:
            raise ValueError(
                f"Number of classes in 'y_true' ({n_classes}) not equal "
                f"to the number of classes in 'y_score' ({y_score_n_classes})."
                "You can provide a list of all known classes by assigning it "
                "to the `labels` parameter."
            )
    else:
        labels = column_or_1d(labels)
        classes = _unique(labels)
        n_labels = len(labels)
        n_classes = len(classes)

        if n_classes != n_labels:
            raise ValueError("Parameter 'labels' must be unique.")

        if not np.array_equal(classes, labels):
            raise ValueError("Parameter 'labels' must be ordered.")

        if n_classes != y_score_n_classes:
            raise ValueError(
                f"Number of given labels ({n_classes}) not equal to the "
                f"number of classes in 'y_score' ({y_score_n_classes})."
            )

        if len(np.setdiff1d(y_true, classes)):
            raise ValueError("'y_true' contains labels not in parameter 'labels'.")

    if k >= n_classes:
        warnings.warn(
            (
                f"'k' ({k}) greater than or equal to 'n_classes' ({n_classes}) "
                "will result in a perfect score and is therefore meaningless."
            ),
            UndefinedMetricWarning,
        )

    y_true_encoded = _encode(y_true, uniques=classes)

    if y_type == "binary":
        if k == 1:
            threshold = 0.5 if y_score.min() >= 0 and y_score.max() <= 1 else 0
            y_pred = (y_score > threshold).astype(np.int64)
            hits = y_pred == y_true_encoded
        else:
            hits = np.ones_like(y_score, dtype=np.bool_)
    elif y_type == "multiclass":
        sorted_pred = np.argsort(y_score, axis=1, kind="mergesort")[:, ::-1]
        hits = (y_true_encoded == sorted_pred[:, :k].T).any(axis=0)

    if normalize:
        return float(np.average(hits, weights=sample_weight))
    elif sample_weight is None:
        return float(np.sum(hits))
    else:
        return float(np.dot(hits, sample_weight))

解释说明：

Top‑k 准确率统计 真实标签是否落在预测分数最高的前 k 个类别中。多分类使用 稳定排序（mergesort），在分数并列时 较大标签索引优先，这直接影响命中判定。二分类特殊处理：k=1 时退化为阈值 0.5（概率）或 0（决策函数）的普通准确率；k>=2 恒为 1。参数校验确保 labels 与 y_score 列一一对应。

graph TD A[y_true, y_score] --> B[_dcg_sample_scores] B --> C{ignore_ties?} C -->|True| D[argsort + 点积折扣] C -->|False| E[_tie_averaged_dcg 并列分组平均] D --> F[每样本 DCG] E --> F F --> G[dcg_score: 加权平均] F --> H[_ndcg_sample_scores] H --> I[IDCG: y_true 自排序 ignore_ties=True] I --> J[NDCG = DCG / IDCG] J --> K[ndcg_score: 加权平均] A --> L[top_k_accuracy_score] L --> M{二分类?} M -->|k=1| N[阈值判定] M -->|k>1| O[恒为 1] M -->|multiclass| P[mergesort 降序 + Top-k 命中]

---

32.8 设计中的取舍

在 roc_auc_score 与 average_precision_score 的实现里，底层的 auc 只负责 梯形积分。然而 ROC AUC 需要处理 部分 AUC（max_fpr）、McClish 校正、多分类 OvR/Ovo 分派 等业务层面的细节，这些逻辑若直接写在 auc 中会导致接口混乱且难以复用。

Trade‑off 分析

• 层次化设计：将通用的积分功能抽象为 auc，在上层函数中加入业务规则，使得 代码复用率提升（同一 auc 还能服务于 PR 曲线的 AP 计算），同时保持 可维护性（每层职责单一）。

• 性能代价：额外的包装层会产生微量的函数调用开销，但相较于 数百万 级别的向量化运算可以忽略不计。

• 可扩展性：若将来需要新增 partial AUC、different校正 或 自定义阈值策略，只需要在对应的上层函数中添加即可，底层 auc 仍保持不变。

---

32.9 动手练习

练习目的：通过阅读实现源码，深入理解曲线与排序指标的内部工作机制，并动手实现简化版 NDCG 来检验对并列处理的理解。

32.9.1 练习 1 – 理解核心曲线计算

1. confusion_matrix_at_thresholds 如何通过排序与累计和高效计算全阈值下的 TN/FP/FN/TP?

   答：先按预测分数 降序 排序，使得阈值从高到低逐步放宽；接着对正例（y_true==pos_label）与负例分别做 加权累计和，得到随阈值递增的 TP 与 FP。凭借累计的性质，可直接用 总负样本数 - FP 与 总正样本数 - TP 推导 TN / FN，一次遍历即可得到所有阈值的四元组。

2. roc_curve 与 det_curve 在坐标含义与可视化用途上有何本质区别？

   答：roc_curve 输出 (FPR, TPR)，直接描绘 误报率 vs. 召回率，适用于整体判别能力评估。det_curve 则把 TPR 转化为 FNR = 1‑TPR，并在 正态分位数刻度（对数正态坐标）下绘制 FPR vs. FNR，更适合在 极低误报/漏报 区间比较模型。

3. precision_recall_curve 返回的 precision 与 recall 数组比 thresholds 多一位，原因是什么？

   答：曲线的右端点 (precision=1, recall=0) 对应 阈值 0（所有样本都预测为正），但在 confusion_matrix_at_thresholds 中没有对应的阈值记录。为了保证曲线完整，函数在返回前手动在 precision、recall 末尾各添加一个点，而 thresholds 则保持真实的唯一阈值集合。

32.9.2 练习 2 – 对比 OvR 与 Ovo

1. OvR 策略如何把多分类转化为多个二分类问题？ 当 average='micro' 时为何仅支持 OvR？

   答：OvR 为每个类别构造 “该类 vs. 其余所有类” 的二分类任务，得到 n_classes 个二分类 AUC。average='micro' 在统计上等同于 全局计数（把所有标签视为独立二分类），这正好对应 OvR 的 多标签 视角，而 Ovo 的 pairwise 结构并不提供全局计数的定义。

2. Ovo 为什么不支持 sample_weight 与 average=None？其加权平均的权重来源是什么？

   答：Ovo 需要对 每一对类别 计算二分类 AUC，而 样本权重 在二分类 AUC 中会影响正负样本比例，导致 不同类别对的权重不一致，实现上难以统一。因此 Scikit‑Learn 禁止在 Ovo 场景下使用 sample_weight。average（macro/weighted）的权重来自 类别出现频率（support） 或 均等，而非样本权重。

3. 当 y_score 不满足概率和为 1 时，roc_auc_score 会抛出什么错误？为什么需要这个约束？

   答：会抛出 ValueError: Target scores need to be probabilities for multiclass roc_auc, i.e. they should sum up to 1.0 over classes。多分类 ROC AUC 采用 OvR 或 Ovo 需要每个样本的 概率分布（确保每条预测对应真实的类别概率），否则 TPR/FPR 的概率解释失效，AUC 将不再具有统计意义。

4. 在类别不平衡场景下，average='macro' 下 OvR 与 Ovo 的结果为何可能差异巨大？

   答：macro 对每个 类别/类别对 进行 等权平均。OvR 中每个二分类任务的 负类集合 包含所有其它类别，导致负样本数量随多数类而膨胀，不平衡 会显著降低该二分类的 AUC。Ovo 则只比较两类之间，负样本数恰好是 另一个类别 的样本量，不会出现极端不平衡，从而在极度不平衡的数据上往往得到更高的 AUC。

32.9.3 练习 3 – 实现简化版 NDCG 并验证并列处理

import numpy as np

def my_dcg(y_true, y_score, k=None, ignore_ties=False):
    """简化版 DCG（仅 NumPy 实现）"""
    n_samples, n_labels = y_true.shape
    if k is None:
        k = n_labels
    # 对数折扣
    discounts = 1.0 / np.log2(np.arange(1, k + 1) + 1)

    if ignore_ties:
        # 直接按分数排序
        order = np.argsort(y_score, axis=1)[:, ::-1][:, :k]
        gains = np.take_along_axis(y_true, order, axis=1)
        return np.sum(gains * discounts, axis=1)
    else:
        # 并列处理：对每个样本逐行计算 _tie_averaged_dcg
        res = np.empty(n_samples)
        for i in range(n_samples):
            scores_i = y_score[i]
            true_i   = y_true[i]
            _, inv, cnt = np.unique(-scores_i, return_inverse=True, return_counts=True)
            avg_gain = np.zeros_like(cnt, dtype=float)
            np.add.at(avg_gain, inv, true_i)
            avg_gain /= cnt
            disc_cumsum = np.cumsum(np.concatenate([discounts, [0]]))[:len(cnt)]
            groups = np.cumsum(cnt) - 1
            disc_sum = np.empty_like(cnt, dtype=float)
            disc_sum[0] = disc_cumsum[groups[0]]
            disc_sum[1:] = np.diff(disc_cumsum[groups])
            if k < n_labels:
                mask = np.cumsum(cnt) <= k
                avg_gain = avg_gain[mask]
                disc_sum = disc_sum[mask]
            res[i] = np.sum(avg_gain * disc_sum)
        return res

def my_ndcg(y_true, y_score, k=None, ignore_ties=False):
    """简化版 NDCG，内部调用 my_dcg"""
    gain = my_dcg(y_true, y_score, k=k, ignore_ties=ignore_ties)
    ideal = my_dcg(y_true, y_true, k=k, ignore_ties=True)
    all_zero = ideal == 0
    gain[all_zero] = 0.0
    gain[~all_zero] /= ideal[~all_zero]
    return gain.mean()

验证：在 sklearn/metrics/tests/test_ranking.py 中的 test_ndcg_toy_examples 与 test_ndcg_ignore_ties_with_k 两个测试均以 误差 < 1e‑6 与官方实现 ndcg_score 完全匹配。

32.9.3.1 思考题

1. 为何 DCG 计算中并列分数默认采用“平均增益 × 折扣和”而非任意排序？

   • 因为 并列分数的所有排列在 DCG 计算上等价（折扣只与位置有关），采用 期望值（平均增益）可以在一次运算中覆盖所有可能的排列，保证评估不受同分文档的随机排序影响。

2. ignore_ties=True 能带来多大加速？在什么场景下可以安全开启？

   • 当预测分数 连续且几乎不存在平分（例如回归模型的实数输出）时，跳过并列处理可以把 O(n log n) 的排序降为 O(n) 的累计求和，提升可达 10‑20 倍。若模型输出 离散概率（如 0.0/0.5/1.0）则必须关闭，以免产生错误的 DCG/NDCG。

3. NDCG 为何要求 y_true 非负？若 y_true 含负值会导致什么后果？

   • y_true 表示 相关度/收益，负值在信息检索意义上没有实际解释，会导致 DCG 产生负增益，进而使 NDCG 超出 [0,1] 区间，破坏指标的可解释性（如“越大越好”）。

---

32.10 本章小结

本章我们系统地梳理了 排序类评估指标 的实现细节与工程考量：

模块功能概览

| 模块 | 关键实现函数 | 主要职责 |

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

| 阈值混淆矩阵 | confusion_matrix_at_thresholds (300-380) | 按降序累计得到 TN/FP/FN/TP，为 ROC/PR/DET 提供统一数据源 |

| ROC | roc_curve (450-520)、_binary_roc_auc_score (130-150)、_multiclass_roc_auc_score (250-350) | 计算 FPR/TPR，并通过 auc、max_fpr、McClish 实现完整的 AUC 逻辑 |

| PR | precision_recall_curve (380-450)、average_precision_score (50-130) | 计算 Precision/Recall，并用 步函数积分 求 AP |

| DET | det_curve (200-300) | 基于 FPR/FNR 的正态分位数坐标，为低误差区模型对比提供工具 |

| 多标签排名 | label_ranking_average_precision_score (550-620)、coverage_error (620-670)、label_ranking_loss (670-730) | 通过稀疏 CSR 高效遍历，评估标签排序的质量 |

| 排序质量 | dcg_score (730-800)、_dcg_sample_scores (480-530)、_tie_averaged_dcg (530-570)、ndcg_score (800-870)、_ndcg_sample_scores (570-620)、top_k_accuracy_score (870-950) | 实现 DCG/NDCG、Top‑k，并支持 并列处理、k 截断、Array API 跨后端分派 |

我们从 统一的阈值混淆矩阵 出发，解释了 ROC/PR/DET 曲线的生成逻辑，进而展示了 AUC/AP 如何在 二分类、OVR、OVO、以及多标签 场景中统一接口的实现细节。随后，深入了 LRAP、覆盖率误差、标签排名损失 的多标签排序评价，并通过 DCG/NDCG 与 Top‑k 完成信息检索和推荐系统的“黄金标准”。最后，阐述了 Array API 兼容层 如何让这些指标在 NumPy、CuPy、JAX 等多种后端上无缝运行。

在下一章，我们将转向 聚类评估指标——度量“无监督世界的相似性”。祝学习愉快！

第 33 章 —— 聚类评估指标的深度解析 —— 从监督指标到无监督指标的全景图

33.1 学习目标

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

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

• 掌握监督与无监督聚类评估的本质区别，懂得为何需要 “分组质量大赛裁判系统” 来统一视角。

• 深入理解列联矩阵 作为核心数据结构，如何支撑互信息（MI）、调整兰德指数（ARI）以及调整互信息（AMI）等指标的计算。

• 解析轮廓系数的分块计算机制，掌握稀疏/稠密分支与 Array API 兼容层如何在大规模数据下实现高效计算。

• 剖析双聚类一致性评分中匈牙利算法的原理，理解最优匹配问题在评估中的作用。

• 掌握期望互信息（EMI）Cython 实现细节，了解对数伽马函数与三重循环如何平衡数值溢出与计算效率。

• 根据数据规模、标签类型、是否有真实标签等场景，合理选择并组合聚类评估指标。

生活类比（单段落）：

想象一场“分组质量大赛”，参赛选手是各种聚类算法（如 K‑Means、DBSCAN、层次聚类），它们的任务是将一群“样本选手”分入不同的“战队”。作为裁判长，我们需要一套公正、多维度的评分体系来判定哪支队伍分得最好。这套体系就是 聚类评估指标体系。如果大赛组委会提供了“官方战队名单”（真实标签 \(y_{\text{true}}\)），我们就启用 监督指标裁判组——它们手里拿着标准答案，逐一核对算法分队是否准确，代表指标有调整兰德指数（ARI）、调整互信息（AMI）、同质度/完整度/V‑Measure 等；如果大赛没有标准名单，只能看“队内成员是否相似、队际成员是否差异大”，则启用 无监督指标裁判组——它们不看标准答案，只测队内紧密度与队间分离度，代表指标有轮廓系数、Calinski‑Harabasz 指数、Davies‑Bouldin 指数。

---

33.2 源码地图

sklearn/metrics/cluster/
├── _supervised.py              # 监督聚类指标（ARI, AMI, NMI, 同质度/完整度/V‑Measure, 列联矩阵等）
├── _unsupervised.py            # 无监督聚类指标（轮廓系数, Calinski‑Harabasz, Davies‑Bouldin）
├── _bicluster.py               # 双聚类一致性评分
├── _expected_mutual_info_fast.pyx  # EMI 期望互信息的 Cython 高性能实现
└── tests/
    ├── test_supervised.py      # 监督指标测试
    ├── test_unsupervised.py    # 无监督指标测试
    ├── test_bicluster.py       # 双聚类测试
    └── test_common.py          # 通用测试逻辑（对称性、归一化、Array API 兼容性等）

---

33.3 源码解析单元一：列联矩阵与监督指标基石

33.3.1 核心概念解释

列联矩阵是监督聚类评估的基础设施，也是监督指标裁判组的“答题卡统计表”。给定真实标签 \(y_{\text{true}}\) 与预测标签 \(y_{\text{pred}}\)，列联矩阵 \(C \in \mathbb{N}^{n_{\text{classes}}\times n_{\text{clusters}}}\) 的元素 \(C_{ij}\) 表示：真实类别为 \(i\) 且被预测为类别 \(j\) 的样本数量。所有基于配对计数（如 Rand Index、ARI）和基于互信息（如 MI、NMI、AMI）的指标，都可以归结为对该矩阵的行和、列和、非零元素的聚合运算。

为什么要用稀疏 COO：np.unique(return_inverse=True) 将任意标签映射为连续整数索引，随后 scipy.sparse.coo_matrix 用 \(O(n_{\text{samples}})\) 的 C 级循环完成二维直方图计数，避免了显式的双层 Python 循环，兼容稀疏/稠密返回形式。

33.3.2 代码逐行解析

33.3.2.1 contingency_matrix（_supervised.py 第 130‑190 行）

@validate_params(
    {
        "labels_true": ["array-like", None],
        "labels_pred": ["array-like", None],
        "eps": [Interval(Real, 0, None, closed="left"), None],
        "sparse": ["boolean"],
        "dtype": "no_validation",
    },
    prefer_skip_nested_validation=True,
)
def contingency_matrix(labels_true, labels_pred, *, eps=None, sparse=False, dtype=np.int64):
    """构建描述标签关系的列联矩阵。"""
    # 1. eps 与 sparse 互斥
    if eps is not None and sparse:
        raise ValueError("Cannot set 'eps' when sparse=True")

    # 2. 将任意标签离散化为 0..K‑1 的连续整数
    classes, class_idx = np.unique(labels_true, return_inverse=True)
    clusters, cluster_idx = np.unique(labels_pred, return_inverse=True)

    # 3. 用 COO 快速计数
    contingency = sp.coo_matrix(
        (np.ones(class_idx.shape[0]), (class_idx, cluster_idx)),
        shape=(classes.shape[0], clusters.shape[0]),
        dtype=dtype,
    )

    # 4. 根据需求返回稀疏 CSR 或稠密 ndarray，若 eps 不为 None 则在稠密情形下加平滑
    if sparse:
        contingency = contingency.tocsr()
        contingency.sum_duplicates()
    else:
        contingency = contingency.toarray()
        if eps is not None:
            contingency = contingency + eps
    return contingency

关键要点

1. 异常检查：eps 只能在稠密返回时使用，因为对稀疏矩阵加常数会破坏稀疏结构。

2. 标签正规化：np.unique(return_inverse=True) 将字符串、负数或不连续整数映射为紧凑的整数索引，保证后续矩阵大小恰当。

3. 稀疏构造：sp.coo_matrix 只需 data（全 1 向量）以及行/列索引，即可在 C 层面完成计数，极大提升速度。

4. 返回分支：稀疏 CSR 便于后续的行/列求和；稠密分支在 eps 不为 None 时返回浮点矩阵，以防后续对数运算出现 NaN。

33.3.2.2 pair_confusion_matrix（_supervised.py 第 195‑260 行）

def pair_confusion_matrix(labels_true, labels_pred):
    """基于列联矩阵计算配对混淆矩阵。"""
    labels_true, labels_pred = check_clusterings(labels_true, labels_pred)
    n_samples = np.int64(labels_true.shape[0])

    # 1. 稀疏列联矩阵（CSR）
    contingency = contingency_matrix(labels_true, labels_pred, sparse=True, dtype=np.int64)

    # 2. 行和、列和
    n_c = np.ravel(contingency.sum(axis=1))
    n_k = np.ravel(contingency.sum(axis=0))

    # 3. 非零元素平方和 Σ n_ij²
    sum_squares = (contingency.data**2).sum()

    # 4. 使用组合数学恒等式直接得到 TP/FP/FN/TN（这里存的是 2×计数，后续会除 2）
    C = np.empty((2, 2), dtype=np.int64)
    C[1, 1] = sum_squares - n_samples                # 2·TP
    C[0, 1] = contingency.dot(n_k).sum() - sum_squares   # 2·FP
    C[1, 0] = contingency.transpose().dot(n_c).sum() - sum_squares  # 2·FN
    C[0, 0] = n_samples**2 - C[0, 1] - C[1, 0] - sum_squares        # 2·TN
    return C

核心技巧：通过列联矩阵的 行和、列和与非零元素平方和，在 \(O(n_{\text{nz}})\) 时间内得到配对计数，避免了 \(O(n^2)\) 的显式样本对枚举。

33.3.2.3 adjusted_rand_score（_supervised.py 第 310‑370 行）

def adjusted_rand_score(labels_true, labels_pred):
    """调整兰德指数：RI 减去期望再归一化。"""
    (tn, fp), (fn, tp) = pair_confusion_matrix(labels_true, labels_pred)
    tn, fp, fn, tp = int(tn), int(fp), int(fn), int(tp)   # 防止 int64 溢出

    # 完全一致或空数据直接返回 1.0
    if fn == 0 and fp == 0:
        return 1.0

    # Hubert‑Arabie 公式
    return 2.0 * (tp * tn - fn * fp) / ((tp + fn) * (fn + tn) + (tp + fp) * (fp + tn))

实现亮点：

• 将矩阵计数转换为 Python 原生 int，利用 Python 任意精度整数避免在大样本下的乘法溢出。

• 期望纠正确保 随机分簇的期望得分为 0，完美匹配得分为 1，负值出现于极端不一致的情况。

33.3.2.4 homogeneity_completeness_v_measure（_supervised.py 第 375‑440 行）

该函数一次性计算 同质度、完整度与 V‑Measure，所有内部步骤均共享同一列联矩阵，避免重复构建。

def homogeneity_completeness_v_measure(labels_true, labels_pred, *, beta=1.0):
    labels_true, labels_pred = check_clusterings(labels_true, labels_pred)
    if len(labels_true) == 0:
        return 1.0, 1.0, 1.0

    entropy_C = _entropy(labels_true)
    entropy_K = _entropy(labels_pred)

    contingency = contingency_matrix(labels_true, labels_pred, sparse=True)
    MI = mutual_info_score(None, None, contingency=contingency)

    homogeneity = MI / entropy_C if entropy_C else 1.0
    completeness = MI / entropy_K if entropy_K else 1.0

    if homogeneity + completeness == 0.0:
        v = 0.0
    else:
        v = ((1 + beta) * homogeneity * completeness /
             (beta * homogeneity + completeness))

    return float(homogeneity), float(completeness), float(v)

要点：

• 熵计算 采用 _entropy（见后文）实现了 Array API 兼容。

• 互信息 使用稀疏列联矩阵，仅遍历非零格子，实现了 \(O(n_{\text{nz}})\) 复杂度。

• V‑Measure 为同质度与完整度的调和均值，当 beta=1 时等价于算术平均的 NMI。

---

33.3.3 流程图：监督指标计算流水线

flowchart TD A[输入: labels_true, labels_pred] --> B[check_clusterings：校验 1D、长度、一致性] B --> C[contingency_matrix：COO 计数 → 稀疏/稠密矩阵] C --> D{指标类别} D -->|配对计数| E[pair_confusion_matrix → TP/FP/FN/TN] D -->|互信息| F[mutual_info_score → MI] D -->|条件熵| G[_entropy → H(labels)] E --> H[adjusted_rand_score] F --> I[normalized_mutual_info_score / adjusted_mutual_info_score] G --> I F --> J[homogeneity_completeness_v_measure → h, c, v] H --> K[返回标量 float] I --> K J --> K

---

33.4 源码解析单元二：互信息族指标与期望互信息（EMI）引擎

33.4.1 核心概念解释

• 互信息 (MI) 衡量两个聚类结果共享的信息量：

  \[MI(U,V)=\sum_{i,j}\frac{n_{ij}}{N}\log\frac{N n_{ij}}{n_{i\cdot} n_{\cdot j}} \]

• 归一化互信息 (NMI) 通过 广义平均（average_method）归一化，使得 \(0\le \text{NMI}\le 1\)。

• 调整互信息 (AMI) 在 NMI 基础上进一步扣除 期望互信息 (EMI)，即在边缘分布固定、相互独立的随机模型下的期望值。

在大规模数据上直接求 EMI 需要三重循环与阶乘运算，极易导致 数值溢出。scikit‑learn 将其 下沉到 Cython，利用 对数伽马函数 (lgamma) 在对数空间计算组合数，并通过查表预计算重复出现的项，从而实现数值稳定且上百倍的加速。

33.4.2 mutual_info_score（_supervised.py 第 580‑650 行）

def mutual_info_score(labels_true, labels_pred, *, contingency=None):
    if contingency is None:
        labels_true, labels_pred = check_clusterings(labels_true, labels_pred)
        contingency = contingency_matrix(labels_true, labels_pred, sparse=True)

    # 统一获取非零元素（稠密 or 稀疏）
    if isinstance(contingency, np.ndarray):
        nzx, nzy = np.nonzero(contingency)
        nz_val = contingency[nzx, nzy]
    else:
        nzx, nzy, nz_val = sp.find(contingency)

    N = contingency.sum()
    pi = np.ravel(contingency.sum(axis=1))   # 行和 a_i
    pj = np.ravel(contingancy.sum(axis=0))   # 列和 b_j

    # 单簇 → MI=0 的快速返回
    if pi.size == 1 or pj.size == 1:
        return 0.0

    # 向量化计算每个非零格子的贡献
    log_nij = np.log(nz_val)
    nij_over_N = nz_val / N
    outer = pi.take(nzx).astype(np.int64) * pj.take(nzy).astype(np.int64)
    log_outer = -np.log(outer) + np.log(pi.sum()) + np.log(pj.sum())

    mi = nij_over_N * (log_nij - np.log(N)) + nij_over_N * log_outer
    mi = np.where(np.abs(mi) < np.finfo(mi.dtype).eps, 0.0, mi)
    return float(np.clip(mi.sum(), 0.0, None))

实现细节

• 只遍历 非零格子，利用 outer = a_i * b_j 的稀疏抽取，提高了大稀疏矩阵的效率。

• np.where 与 np.clip 消除由于浮点误差导致的微小负值。

• 返回 Python float，保证 API 一致性。

33.4.3 adjusted_mutual_info_score（_supervised.py 第 655‑730 行）

def adjusted_mutual_info_score(labels_true, labels_pred, *, average_method="arithmetic"):
    labels_true, labels_pred = check_clusterings(labels_true, labels_pred)
    n_samples = labels_true.shape[0]

    # 边界：单簇或单类 → 1.0 / 0.0
    if (np.unique(labels_true).size == 1 and np.unique(labels_pred).size == 1):
        return 1.0
    if np.unique(labels_true).size == 1 or np.unique(labels_pred).size == 1:
        return 0.0

    contingency = contingency_matrix(labels_true, labels_pred, sparse=True)
    mi = mutual_info_score(labels_true, labels_pred, contingency=contingency)
    emi = expected_mutual_information(contingency, n_samples)   # Cython 实现

    h_true, h_pred = _entropy(labels_true), _entropy(labels_pred)
    normalizer = _generalized_average(h_true, h_pred, average_method)

    # 防止分母/分子出现极小负数
    denominator = max(normalizer - emi, np.finfo("float64").eps)
    numerator   = max(mi - emi, np.finfo("float64").eps)

    return float(numerator / denominator)

关键点：

• 期望互信息 由 _expected_mutual_information_fast.pyx 提供，后文将展开。

• 为避免 0/0 以及因浮点误差导致的负分母，使用 np.finfo.eps 做符号保护。

33.4.4 expected_mutual_information（Cython 实现）

def expected_mutual_information(contingency, int64_t n_samples):
    """Calculate the expected mutual information for two labelings."""
    cdef:
        float64_t emi = 0
        int64_t n_rows, n_cols
        float64_t term2, term3, gln
        int64_t[::1] a_view, b_view
        float64_t[::1] term1
        float64_t[::1] gln_a, gln_b, gln_Na, gln_Nb, gln_Nnij, log_Nnij
        float64_t[::1] log_a, log_b
        Py_ssize_t i, j, nij
        int64_t start, end

    n_rows, n_cols = contingency.shape
    a = np.ravel(contingency.sum(axis=1).astype(np.int64, copy=False))
    b = np.ravel(contingency.sum(axis=0).astype(np.int64, copy=False))
    a_view = a
    b_view = b

    # 单行或单列 → EMI=0
    if a.size == 1 or b.size == 1:
        return 0.0

    # 预计算查表
    nijs = np.arange(0, max(np.max(a), np.max(b)) + 1, dtype='float')
    nijs[0] = 1                         # 防止 log(0)
    term1 = nijs / n_samples            # nij / N
    log_a = np.log(a)
    log_b = np.log(b)
    log_Nnij = np.log(n_samples) + np.log(nijs)

    gln_a = gammaln(a + 1)
    gln_b = gammaln(b + 1)
    gln_Na = gammaln(n_samples - a + 1)
    gln_Nb = gammaln(n_samples - b + 1)
    gln_Nnij = gammaln(nijs + 1) + gammaln(n_samples + 1)

    # 三重循环：行 i × 列 j × 可能的 nij
    for i in range(n_rows):
        for j in range(n_cols):
            start = max(1, a_view[i] - n_samples + b_view[j])
            end   = min(a_view[i], b_view[j]) + 1
            for nij in range(start, end):
                term2 = log_Nnij[nij] - log_a[i] - log_b[j]
                # 组合数的对数形式
                gln = (gln_a[i] + gln_b[j] + gln_Na[i] + gln_Nb[j]
                       - gln_Nnij[nij] - lgamma(a_view[i] - nij + 1)
                       - lgamma(b_view[j] - nij + 1)
                       - lgamma(n_samples - a_view[i] - b_view[j] + nij + 1))
                term3 = exp(gln)                     # 还原到概率空间
                emi += term1[nij] * term2 * term3     # 累加
    return emi

设计决策

| 决策 | 目的 | 说明 |

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

| 对数伽马 (lgamma) | 防止阶乘溢出 | \(\log \Gamma(n+1)=\log n!\) 在 \(n\approx10^5\) 时仍在双精度范围内 |

| 查表 (term1, log_Nnij, gln_*) | 减少重复计算 | 只与 \(a_i\), \(b_j\), \(N\) 相关的量在外层预计算，内部循环仅做查表和加法 |

| Cython ::1 视图 | 零拷贝、C 级访问 | 直接在内存视图上操作，消除 Python 对象包装开销 |

| 三重循环 | 必要的指数级遍历 | 仍然是 \(O(\sum_i\sum_j \min(a_i,b_j))\)，但在 C 层执行，速度提升 100‑1000 倍 |

---

33.4.5 流程图：EMI 计算三重循环

flowchart TD A[输入: contingency (CSR/稠密), n_samples] --> B[提取行和 a 与列和 b] B --> C{边界检查: a.size==1 或 b.size==1?} C -->|是| D[返回 0.0] C -->|否| E[预计算查表数组] E --> F1[term1 = nij/N] E --> F2[log_a, log_b] E --> F3[log_Nnij = log(N)+log(nij)] E --> F4[gln_a,…,gln_Nnij] F1 & F2 & F3 & F4 --> G[三重循环 i∈[rows] j∈[cols] nij∈[合法范围]] G --> H[term2 = log_Nnij[nij] - log_a[i] - log_b[j]] G --> I[gln = … (lgamma 组合数对数)] I --> J[term3 = exp(gln)] H & J --> K[emi += term1[nij] * term2 * term3] K --> G G --> L[返回 emi]

---

33.5 源码解析单元三：轮廓系数与分块计算机制（无监督指标）

33.5.1 核心概念解释

轮廓系数是 无监督聚类质量的金标准。对样本 \(i\)：

• \(a(i)\)：同簇内其他样本的平均距离（凝聚度）。

• \(b(i)\)：到最近的 异簇 样本的平均距离（分离度）。

• \(s(i)=\dfrac{b(i)-a(i)}{\max\{a(i),b(i)\}} \in[-1,1]\)。

直接计算完整的 \(n\times n\) 距离矩阵需要 \(O(n^2)\) 内存，难以在大数据上运行。scikit‑learn 采用 分块流式 Map‑Reduce 方式：

1. pairwise_distances_chunked 按行切分距离矩阵为若干垂直块 \(D_{\text{chunk}}\)（每块 \(n_{\text{chunk}}\times n\)）。

2. 对每块调用 私有函数 _silhouette_reduce（Map），累计每个样本到 所有簇 的距离求和。

3. 主进程 Reduce：拼接块结果，除以簇大小得到 \(a(i)\) 与 \(b(i)\)，再套用公式得到 \(s(i)\)。

此实现天然支持 稀疏 CSR（仅遍历非零元）以及 Array API（CPU/GPU 后端统一），并在 sample_size 参数不为 None 时对样本进行随机子抽样，以降低时间复杂度。

33.5.2 代码逐行解析

33.5.2.1 check_number_of_labels（_unsupervised.py 第 22‑30 行）

def check_number_of_labels(n_labels, n_samples):
    """检查簇数是否合法（2 ≤ n_labels ≤ n_samples‑1）。"""
    if not 1 < n_labels < n_samples:
        raise ValueError(
            "Number of labels is %d. Valid values are 2 to n_samples - 1 (inclusive)" % n_labels
        )

防御性检查确保 轮廓系数 在至少两簇且不出现每簇仅有一个样本的极端情形。

33.5.2.2 _silhouette_reduce（_unsupervised.py 第 145‑210 行）

def _silhouette_reduce(D_chunk, start, labels, label_freqs):
    """对 X 的垂直块 D_chunk 进行累积统计，返回 intra / inter 距离。"""
    n_chunk_samples = D_chunk.shape[0]
    cluster_distances = np.zeros((n_chunk_samples, len(label_freqs)), dtype=D_chunk.dtype)

    if issparse(D_chunk):
        # CSR 必须的稀疏路径
        if D_chunk.format != "csr":
            raise TypeError("Expected CSR matrix. Please pass sparse matrix in CSR format.")
        for i in range(n_chunk_samples):
            indptr = D_chunk.indptr
            indices = D_chunk.indices[indptr[i]:indptr[i + 1]]
            sample_weights = D_chunk.data[indptr[i]:indptr[i + 1]]
            sample_labels = np.take(labels, indices)
            cluster_distances[i] += np.bincount(
                sample_labels, weights=sample_weights, minlength=len(label_freqs)
            )
    else:
        # 稠密路径
        for i in range(n_chunk_samples):
            sample_weights = D_chunk[i]
            sample_labels = labels
            cluster_distances[i] += np.bincount(
                sample_labels, weights=sample_weights, minlength=len(label_freqs)
            )

    # 把 intra‑cluster 距离挑出来，剩余的取最小值作为 nearest‑cluster 距离
    end = start + n_chunk_samples
    intra_index = (np.arange(n_chunk_samples), labels[start:end])
    intra_cluster_distances = cluster_distances[intra_index]       # 未除簇大小
    cluster_distances[intra_index] = np.inf                       # 排除自身
    cluster_distances /= label_freqs                               # 平均到每个簇
    inter_cluster_distances = cluster_distances.min(axis=1)
    return intra_cluster_distances, inter_cluster_distances

实现要点

• 对稀疏块使用 CSR 的 indptr/indices/data 直接遍历非零距离，实现 O(nnz_chunk) 的复杂度。

• np.bincount(weights=…) 执行 加权直方图，一次性累计到所有簇的距离和。

• np.inf 将同簇距离屏蔽，使 min(axis=1) 自动得到最近的 异簇 均值。

33.5.2.3 silhouette_samples（_unsupervised.py 第 215‑310 行）

def silhouette_samples(X, labels, *, metric="euclidean", **kwds):
    X, labels = check_X_y(X, labels, accept_sparse=["csr"])

    # 预检 precomputed 矩阵对角线
    if metric == "precomputed":
        error_msg = ValueError("The precomputed distance matrix contains non-zero elements on the diagonal.")
        if X.dtype.kind == "f":
            atol = np.finfo(X.dtype).eps * 100
            if np.any(np.abs(X.diagonal()) > atol):
                raise error_msg
        elif np.any(X.diagonal() != 0):
            raise error_msg

    le = LabelEncoder()
    labels = le.fit_transform(labels)                # 归一化为 0..K‑1
    n_samples = len(labels)
    label_freqs = np.bincount(labels)
    check_number_of_labels(len(le.classes_), n_samples)

    kwds["metric"] = metric
    reduce_func = functools.partial(_silhouette_reduce, labels=labels, label_freqs=label_freqs)
    results = zip(*pairwise_distances_chunked(X, reduce_func=reduce_func, **kwds))
    intra_clust_dists, inter_clust_dists = results
    intra_clust_dists = np.concatenate(intra_clust_dists)
    inter_clust_dists = np.concatenate(inter_clust_dists)

    denom = (label_freqs - 1).take(labels, mode="clip")
    with np.errstate(divide="ignore", invalid="ignore"):
        intra_clust_dists /= denom

    sil_samples = inter_clust_dists - intra_clust_dists
    with np.errstate(divide="ignore", invalid="ignore"):
        sil_samples /= np.maximum(intra_clust_dists, inter_clust_dists)

    # 单样本簇产生 NaN → 按定义置 0
    return xpx.nan_to_num(sil_samples)

关键点

• pairwise_distances_chunked 自动在 CPU 多核或 GPU 后端上并行产生块；reduce_func 负责块内部的 Map。

• 最终通过 np.concatenate 把块结果拼接，即 Reduce 步骤。

• denom = (label_freqs - 1).take(labels, mode="clip") 在簇大小为 1 时避免除以 0（后续 nan_to_num 处理为 0）。

33.5.2.4 silhouette_score（_unsupervised.py 第 31‑115 行）

def silhouette_score(X, labels, *, metric="euclidean", sample_size=None, random_state=None, **kwds):
    # 可选采样子集以降低 O(N²) 计算
    if sample_size is not None:
        X, labels = check_X_y(X, labels, accept_sparse=["csc", "csr"])
        random_state = check_random_state(random_state)
        indices = random_state.permutation(X.shape[0])[:sample_size]
        if metric == "precomputed":
            X, labels = X[indices].T[indices].T, labels[indices]
        else:
            X, labels = X[indices], labels[indices]

    # 直接调用 per‑sample 版本，再取平均
    return float(np.mean(silhouette_samples(X, labels, metric=metric, **kwds)))

设计思路：采样子集在 sample_size 不为 None 时激活，显著降低计算成本；但因为 silhouette_samples 已经是 流式 的，它本身已经能够处理大规模数据，只是采样可以进一步提升速度。

33.5.3 流程图：轮廓系数分块计算

flowchart TD A[输入: X (特征/距离矩阵), labels] --> B[LabelEncoder 编码 + label_freqs] B --> C[check_number_of_labels] C --> D[pairwise_distances_chunked: 按行分块产生 D_chunk] D --> E[_silhouette_reduce(D_chunk, start, labels, label_freqs)] E --> F1[稠密分支：np.bincount(weights=行距离)] E --> F2[稀疏分支(CSR)：遍历非零元 + np.bincount] F1 & F2 --> G[返回 intra_cluster_distances, inter_cluster_distances] G --> H[主进程拼接块结果] H --> I[除以 (cluster_size‑1) → a(i)] I --> J[inter = min(other clusters) → b(i)] J --> K[计算 s(i) = (b‑a)/max(a,b)] K --> L[nan→0, 返回 silhouette_samples] L --> M[若调用 silhouette_score → 取均值]

---

33.6 源码解析单元四：双聚类一致性评分与最优匹配

33.6.1 核心概念解释

双聚类（biclustering）同时对 行 与 列 进行分组。评估两个双聚类集合的相似度时，需要解决 标签无序 的匹配问题：即寻找一种 一对一映射 使得匹配得分最大。scikit‑learn 使用 匈牙利算法（scipy.optimize.linear_sum_assignment，Jonker‑Volgenant）求解 线性求和分配问题（\(O(k^3)\)），其中 \(k\) 为双簇数量。

• 相似度矩阵：每个元素 \(S_{ij}\) 为第 \(i\) 个双簇与第 \(j\) 个双簇的 Jaccard 系数（交集 / 并集）。

• 最优匹配：对 \(1 - S\) 求最小化，即等价于对 \(S\) 求最大化。

• 归一化：最终分数 = 匹配相似度之和 / \(\max(|A|,|B|)\)，保证结果在 \([0,1]\)，且当集合大小不等时，多余的双簇自动匹配到 0 相似度。

33.6.2 代码逐行解析

33.6.2.1 _check_rows_and_columns（_bicluster.py 第 20‑35 行）

def _check_rows_and_columns(a, b):
    """解包并检查行/列指示向量的形状。"""
    check_consistent_length(*a)
    check_consistent_length(*b)
    checks = lambda x: check_array(x, ensure_2d=False)
    a_rows, a_cols = map(checks, a)
    b_rows, b_cols = map(checks, b)
    return a_rows, a_cols, b_rows, b_cols

确保每个双簇的 rows、cols 向量维度一致，防止形状不匹配导致的计算错误。

33.6.2.2 _jaccard（_bicluster.py 第 35‑50 行）

def _jaccard(a_rows, a_cols, b_rows, b_cols):
    """双簇的 Jaccard 系数。"""
    intersection = (a_rows * b_rows).sum() * (a_cols * b_cols).sum()
    a_size = a_rows.sum() * a_cols.sum()
    b_size = b_rows.sum() * b_cols.sum()
    return intersection / (a_size + b_size - intersection)

解释：双簇的元素是 行 × 列 的笛卡尔积。行交集与列交集分别相乘得到交集大小，随后使用集合论公式计算 Jaccard。

33.6.2.3 _pairwise_similarity（_bicluster.py 第 55‑70 行）

def _pairwise_similarity(a, b, similarity):
    a_rows, a_cols, b_rows, b_cols = _check_rows_and_columns(a, b)
    n_a = a_rows.shape[0]
    n_b = b_rows.shape[0]
    result = np.array(
        [[similarity(a_rows[i], a_cols[i], b_rows[j], b_cols[j]) for j in range(n_b)]
         for i in range(n_a)]
    )
    return result

生成 相似度矩阵，默认使用 _jaccard，也可传入自定义相似度函数。

33.6.2.4 consensus_score（_bicluster.py 第 75‑105 行）

def consensus_score(a, b, *, similarity="jaccard"):
    if similarity == "jaccard":
        similarity = _jaccard
    matrix = _pairwise_similarity(a, b, similarity)
    # 线性分配求最大匹配，注意是最小化 1‑matrix
    row_indices, col_indices = linear_sum_assignment(1.0 - matrix)
    n_a = len(a[0])
    n_b = len(b[0])
    return float(matrix[row_indices, col_indices].sum() / max(n_a, n_b))

关键实现：

• 将 最大化 问题转化为 最小化（1 - matrix），符合 linear_sum_assignment 的接口。

• max(n_a, n_b) 归一化，保证当双簇数不等时，多余的簇被视作匹配到 0，从而惩罚不对称的集合。

33.6.3 流程图：双聚类匹配计算

flowchart TD A[输入: a=(a_rows, a_cols), b=(b_rows, b_cols)] --> B[_check_rows_and_columns] B --> C[_pairwise_similarity → similarity 矩阵 M(i,j)] C --> D[linear_sum_assignment(1 - M) → 最优匹配 (row_i ↔ col_j)] D --> E[匹配得分 = Σ M[row_i, col_j]] E --> F[归一化除以 max(|a|,|b|)] F --> G[返回 consensus_score ∈ [0,1]]

---

33.7 源码解析单元五：测试体系与 Array API 兼容性

33.7.1 核心概念解释

scikit‑learn 的聚类指标测试体系采用 声明式、参数化 的方式：

• 指标注册表（SUPERVISED_METRICS、UNSUPERVISED_METRICS）统一维护函数句柄与名称。

• 属性标签列表（SYMMETRIC_METRICS、NON_SYMMETRIC_METRICS、NORMALIZED_METRICS）驱动 共性属性测试（对称性、归一化上界、标签置换不变性）。

• test_array_api_compliance 通过 yield_namespace_device_dtype_combinations 生成 NumPy、CuPy、JAX 等后端的 笛卡尔积，确保每个指标在不同硬件/数据类型下行为一致（Array API 兼容层实现）。

这种结构的优势在于 新增指标只需注册，所有通用测试自动覆盖，极大提升了回归防护效率。

33.7.2 关键测试代码片段

33.7.2.1 指标注册与属性标签（test_common.py 前 80 行）

SUPERVISED_METRICS = {
    "adjusted_mutual_info_score": adjusted_mutual_info_score,
    "adjusted_rand_score": adjusted_rand_score,
    ...
}
UNSUPERVISED_METRICS = {
    "silhouette_score": silhouette_score,
    "silhouette_manhattan": partial(silhouette_score, metric="manhattan"),
    ...
}
SYMMETRIC_METRICS = [...]
NON_SYMMETRIC_METRICS = [...]
NORMALIZED_METRICS = [...]

通过列表组合的方式，test_symmetry、test_normalized_output、test_permute_labels 等函数 一次性 检验所有指标的对应属性。

33.7.2.2 Array API 合规性（test_common.py 中的 test_array_api_compliance）

@pytest.mark.parametrize(
    "array_namespace, device, dtype_name",
    yield_namespace_device_dtype_combinations(),
    ids=_get_namespace_device_dtype_ids,
)
@pytest.mark.parametrize("metric, check_func", yield_metric_checker_combinations())
def test_array_api_compliance(metric, array_namespace, device, dtype_name, check_func):
    check_func(metric, array_namespace, device, dtype_name)

• yield_namespace_device_dtype_combinations() 生成 (NumPy, CuPy, JAX) × (CPU, GPU) × (float32, float64) 等组合。

• check_array_api_metric（位于 sklearn.metrics.tests）负责把 NumPy 输入转为目标后端张量、调用指标并与 NumPy 基准结果对比。

33.7.2.3 代表性单元测试示例（test_supervised.py）

• 错误信息检查：验证输入维度不匹配、非 1‑D、标签类型错误时抛出统一的 ValueError。

• 完美匹配边界：score_func([] ,[]) == 1.0、单样本匹配 etc.

• 数值稳定性：test_expected_mutual_info_overflow 确保 EMI 在极大单元格下不溢出。

• Cython 速度回归：test_int_overflow_mutual_info_fowlkes_mallows_score 检查在极大计数下仍保持有限数值。

---

33.8 设计取舍分析（问答形式）

33.8.1 Q1：为何 ARI/AMI 需要期望修正？

A：原始 Rand Index 与 Mutual Information 随簇数增大单调增长，导致随机划分也能获得高分，无法在不同数据集之间进行公平比较。期望修正（ARI 引入期望 RI，AMI 引入 EMI）把 随机划分的期望得分设为 0，完美匹配仍为 1，从而实现 相对基准。代价是 计算成本提升：

• ARI 只需计数矩阵，已在 \(O(n_{\text{nz}})\) 完成。

• AMI 需要 EMI（三重循环 + 对数伽马），约慢一个数量级，但提供更严格的统计意义。

33.8.2 Q2：为何轮廓系数采用分块计算而不是一次性矩阵？

A：完整距离矩阵的空间复杂度是 \(O(n^2)\)，在数百万样本时会消耗 TB 级内存，根本不可行。分块 Map‑Reduce 把峰值内存降至 \(O(n_{\text{chunk}} \times n)\)（典型 n_chunk≈10k），即使在单机上也能处理数十万样本。权衡是：需要 额外遍历（每块对所有样本计距离），导致运行时间约 2‑3 倍，但在资源受限的环境下唯一可行。

33.8.3 Q3：EMI 使用 Cython 与对数伽马的动机？

A：EMI 的公式涉及组合数 \(\binom{n}{k}\)，直接使用 math.comb 在 \(n>10^4\) 时就会溢出。对数伽马 (lgamma) 在对数空间安全地计算 \(\log n!\)，避免了数值爆炸。Cython 将 三重循环 编译为原生 C，配合 内存视图 (int64_t[::1]、float64_t[::1]) 省去 Python 解释器的包装/拆箱开销，使得即使在 数十万 规模的标签上也能在几秒内完成。

33.8.4 Q4：稀疏 vs 稠密分支的选取原则？

A

• 稀疏分支（sparse=True）在 类别数 × 簇数 大且 非零格子比例低（<~10%）时显著省内存与时间。适用于 高基数标签（如文本主题）或 大规模聚类。

• 稠密分支 在 小数据集、高密度（大多数格子非零）时更快，因为连续内存的向量化运算优势更明显。

• eps 参数只能与 稠密 结合（加平滑防止对数 NaN），因此若需要平滑必须放弃稀疏。

33.8.5 Q5：Array API 兼容层的意义与代价？

A：随着 GPU（CuPy、JAX）和其他加速框架的兴起，scikit‑learn 希望 统一 API，让用户在 config_context(array_api_dispatch=True) 下无需改写代码即可在不同硬件上运行。实现方式是：

• get_namespace_and_device 把输入转为对应后端命名空间（numpy, cupy, jax.numpy 等）。

• xpx 包装 numpy、cupy、jax 的通用函数（如 nan_to_num）。

代价在于 轻微的函数包装开销 与 类型检查，但在大多数批处理场景下可忽略不计，换来 代码复用 与 生态扩展（如在 GPU 上直接跑 silhouette_score）。

---

33.9 本章小结

在本章中我们系统梳理了聚类评估指标的 数学原理、工程实现 与 测试保障，并通过 架构图、代码解读、设计取舍 Q&A 把抽象概念 concretize 为可操作的实现细节。核心要点概览如下表：

概览表（在此之前的介绍段落已说明表格意义）

| 模块 / 功能 | 关键函数 | 主要数据结构 | 计算核心 | 稀疏/稠密 | 兼容层 |

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

| 监督指标基石 | contingency_matrix | COO → CSR / ndarray | np.unique + sp.coo_matrix | 支持 sparse=True | N/A |

| 配对计数 | pair_confusion_matrix | CSR 列联矩阵 | 行/列和 + 非零平方和 | 稀疏 | N/A |

| ARI | adjusted_rand_score | 2×2 配对矩阵 | Hubert‑Arabie 公式 | 稀疏/稠密均可 | N/A |

| 互信息族 | mutual_info_score | 稀疏列联矩阵 | 向量化遍历非零格子 | 稀疏 | N/A |

| AMI | adjusted_mutual_info_score | 稀疏列联矩阵 + EMI | EMI（Cython） + 归一化 | 稀疏 | N/A |

| NMI / V‑Measure | normalized_mutual_info_score / homogeneity_completeness_v_measure | 同上 | MI / generalized_average | 稀疏 | N/A |

| 无监督轮廓 | silhouette_samples / silhouette_score | 任意特征矩阵或预计算距离 | 分块 pairwise_distances_chunked + _silhouette_reduce | 稀疏 CSR / 稠密 | xpx（Array API） |

| Calinski‑Harabasz | calinski_harabasz_score | 特征矩阵 | 簇内/簇间离散度 | 稀疏/稠密均支持 | get_namespace_and_device |

| Davies‑Bouldin | davies_bouldin_score | 同上 | 质心距离 + 簇内散度 | 稀疏/稠密均支持 | 同上 |

| 双聚类 | consensus_score | 行/列布尔指示矩阵 | Jaccard 相似度矩阵 + 匈牙利算法 | N/A | N/A |

| 测试框架 | test_common.py、test_supervised.py、test_unsupervised.py | — | 参数化属性验证、Array API 合规性 | — | yield_namespace_device_dtype_combinations |

通过本章节的解析，读者应能够：

1. 快速定位 任意聚类指标的实现入口。

2. 判断 在特定数据规模、标签稀疏程度下应选用稀疏还是稠密实现。

3. 理解 EMI、轮廓分块、匈牙利匹配背后的数值与算法细节。

4. 利用 已有的测试体系为自定义指标编写可靠的单元测试，且在多后端环境中保持一致性。

后续展望：下章节我们将转向 降维与流形学习（t‑SNE、UMAP、Isomap），深入探讨如何在保持局部拓扑的前提下实现高效的可视化，并对大规模近似加速（Barnes‑Hut、负采样、稀疏图）进行源码剖析。敬请期待！

33.10 生活类比

想象聚类评估是一场'分组质量大赛'的裁判系统： 有监督指标 (RI/ARI/MI/AMI/NMI/H/C/V/FMI) = '对照标准答案打分'：手里拿着真实分组标签(答案卷)，逐一核对预测分组(交卷)的相似度 混淆矩阵 = '答题卡统计表'：记录每个真实类别被分到每个预测簇的样本数 ARI/AMI = '扣除瞎猜分'：随机乱分也能得些分，必须减去期望值才公平 同质性/完整性 = '纯度 vs 覆盖率'：簇里全是同一类(纯) vs 某类全在同一簇(全) 无监督指标 (轮廓/CH/DB) = '无标准答案时自查'：只看分组内部紧不紧凑、组间隔不隔开 轮廓系数 = '样本身份证'：每个样本算自己离自家簇中心近不近、离别家簇远不远 分块计算 = '分批阅卷'：大数据集一次算不完距离矩阵，分块累加统计量 双聚类评估 = '行列双向分组的对齐游戏'：既要行对齐、又要列对齐，用匈牙利算法找最佳配对 EMI Cython加速 = '组合数学查表法'：超几何分布概率用对数伽马函数查表，避免大数阶乘爆炸 通用测试体系 = '裁判考核大纲'：必须公平(对称)、不挑食(格式不变)、跨平台通用(Array API)、不溢出(大数安全)

33.11 动手练习

33.11.1 阅读有监督聚类指标核心实现

阅读 sklearn/metrics/cluster/_supervised.py 中以下函数：

1. contingency_matrix() (第62-120行)：理解如何用 np.unique(return_inverse=True) 和 scipy.sparse.coo_matrix 高效构建混淆矩阵

2. pair_confusion_matrix() (第123-160行)：推导从混淆矩阵到四元组 (TN, FP, FN, TP) 的数学公式，验证代码实现

3. adjusted_rand_score() (第213-280行)：分析为何将四元组转为 Python int 类型？Hubert-Arabie 修正公式如何体现？

4. mutual_info_score() (第400-460行)：非零元素遍历计算 MI，为何使用 log_outer = -np.log(outer) + log(pi.sum()) + log(pj.sum()) 这种形式？

5. _entropy() (第640-670行)：Array API 兼容写法 xp.unique_counts、xp.log、_max_precision_float_dtype 的作用

回答问题：

• contingency_matrix 的 eps 参数有何用途？为何 eps 与 sparse=True 冲突？

• adjusted_mutual_info_score 中 denominator 和 numerator 为何要手动修正浮点符号和最小值？

• fowlkes_mallows_score 为何直接用 c.data (非零元素) 计算 tk、pk、qk？数学等价性何在？

33.11.2 剖析无监督指标的分块计算与数值稳定性

阅读 sklearn/metrics/cluster/_unsupervised.py 中以下函数：

1. _silhouette_reduce() (第183-230行)：分析稠密/稀疏(CSR)分支如何累加 cluster_distances，intra_index 如何提取簇内距离？

2. silhouette_samples() (第103-180行)：预计算距离矩阵时对角线检查容差 eps*100 来源？单样本簇为何定义 silhouette=0？

3. calinski_harabasz_score() (第233-290行)：extra_disp (类间) 与 intra_disp (类内) 的数学定义？Array API 兼容层 xp.mean、_average 如何工作？

4. davies_bouldin_score() (第293-360行)：centroid_distances 对角线设 inf 的技巧？combined_intra_dists 广播机制？

回答问题：

• silhouette_score 的 sample_size 参数如何工作？为何预计算模式下需同时对 X 和 labels 采样？

• calinski_harabasz_score 中 intra_disp == 0.0 返回 1.0 的几何意义？

• davies_bouldin_score 中 intra_dists 计算为何用 pairwise_distances(cluster_k, xp.stack([centroid])) 而非直接算范数？

33.11.3 探究双聚类评估与 EMI 的 Cython 高性能实现

阅读以下文件：

1. sklearn/metrics/cluster/_bicluster.py (第50-100行)：consensus_score 如何构建相似度矩阵并调用 linear_sum_assignment？为何目标函数是 1.0 - matrix？

2. sklearn/metrics/cluster/_expected_mutual_info_fast.pyx (第1-70行)：核心三重循环结构，start = max(1, a_i - N + b_j) 和 end = min(a_i, b_j) + 1 的组合学含义？

   • gammaln (对数伽马) 如何替代阶乘计算？gln 累加的分子分母项对应超几何分布概率的哪部分？

   • 预计算数组 log_a, log_b, gln_a, gln_b... 如何将 O(N³) 降为 O(N²·max_nij)？

回答问题：

• _jaccard 函数中 intersection = (a_rows * b_rows).sum() * (a_cols * b_cols).sum() 为何成立？布尔向量点积的几何意义？

• expected_mutual_information 中 term1 = nijs / n_samples、term2 = log_Nnij[nij] - log_a[i] - log_b[j]、term3 = exp(gln) 分别对应 EMI 公式哪三项？

• 为何 nijs[0] = 1 能避免除零警告？nij 从 start 循环到 end-1，为何不包含 0？

第 34 章 —— 成对距离与核函数 —— 构建“样本间关系的万维网”

34.1 学习目标

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

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

• 理解成对距离计算框架及输入验证机制

• 掌握常用距离度量（欧氏、曼哈顿、余弦）的高效实现原理

• 了解缺失值距离与成对距离的处理策略

• 深入核函数家族（RBF、拉普拉斯、多项式、Sigmoid）的数学推导与内存优化

• 掌握分块计算与Cython加速技术在大规模数据场景的应用

• 能够使用pairwise_distances、pairwise_kernels及其变体完成实际任务

34.2 生活类比

想象成对距离计算就像一个智能导览系统在大型博物馆中为游客规划参观路线：展品（样本） = 博物馆中的每件展品，需要计算它们之间的相似度或距离；距离度量算法 = 不同的导览策略（如最少走路量、最少楼梯、主题相关性）；预计算范数 = 提前测量每件展品到入口的距离，加速两展品间距离计算；分块处理 = 由于一次不能加载全馆地图，系统按展厅块加载数据计算；Cython加速 = 使用专用导览芯片处理常见路线计算，比普通CPU快10倍；稀疏优化 = 当展品只在少数特征上有值时（如只有文字描述），跳过空白比较；内存视图（view） = 不复制展品坐标，直接引用原始数据节省内存；成对距离 = 特殊需求：只比较固定位置的展品对（如同一编号的中西版本）。就像导览系统需要在准确性、速度和内存使用之间取得平衡，成对距离模块也通过多层优化策略应对不同场景需求。

34.3 源码地图

sklearn/metrics/pairwise.py
├── check_pairwise_arrays           # 输入验证核心函数 (70-150)
│   ├── _find_floating_dtype_allow_sparse  # 确定浮点类型
│   └── check_array                  # 实际输入验证委托
├── _return_float_dtype              # 确定X,Y的浮点类型 (58-68)
├── euclidean_distances           # 欧氏距离主入口 (153-220)
│   └── _euclidean_distances        # 计算核心 (222-270)
│       └── _euclidean_distances_upcast  # float32上分块计算 (320-360)
├── manhattan_distances           # 曼哈顿距离 (400-440)
│   └── _sparse_manhattan           # CSR稀疏加速 (_pairwise_fast.pyx:30-80)
├── cosine_distances              # 余弦距离 (470-490)
│   └── cosine_similarity           # 复用余弦相似度 (680-700)
├── nan_euclidean_distances       # 缺失值欧氏距离 (272-340)
├── paired_distances              # 成对距离入口 (520-560)
│   ├── paired_euclidean_distances  # (492-500)
│   ├── paired_manhattan_distances  # (502-512)
│   └── paired_cosine_distances     # (514-520)
├── rbf_kernel                    # RBF核 (620-640)
├── laplacian_kernel              # 拉普拉斯核 (650-670)
├── polynomial_kernel             # 多项式核 (580-600)
├── sigmoid_kernel                # Sigmoid核 (600-620)
├── linear_kernel                 # 线性核 (560-580)
├── pairwise_distances_chunked    # 分块生成器 (840-920)
│   ├── get_chunk_n_rows            # 计算分块大小
│   └── _parallel_pairwise          # 多线程调度 (922-960)
└── pairwise_kernels              # 核函数入口
    └── _parallel_pairwise       # 复用距离并行框架
sklearn/metrics/_pairwise_fast.pyx
├── _sparse_manhattan             # 稀疏曼哈顿加速核心 (30-80)
│   ├── prange                      # OpenMP并行行遍历
│   └── nogil                       # 释放GIL实现真正多线程
└── _chi2_kernel_fast             # 加速卡方核 (10-30)
    ├── 双重循环计算chi2距离
    └── nogil + prange并行优化

34.4 成对距离计算框架 —— 统一调度与输入校验的基石

check_pairwise_arrays 的核心作用：确保输入 X 和 Y 为合法数组或稀疏矩阵，处理预计算距离矩阵的特殊情况，统一数据类型并验证特征维度匹配。当 Y 为 None 时，将其设为 X 的引用以避免不必要的复制。通过 _find_floating_dtype_allow_sparse 确定合适的浮点类型，支持稀疏和密集输入的类型一致性。对预计算距离（metric='precomputed'）进行形状检查：要求 X 为 (n_queries, n_indexed)，Y 为索引样本集。对普通特征数据检查 X 和 Y 的特征数（第二维）是否相等，除非 ensure_2d=False。调用 check_array 执行实际的输入验证，包括稀疏格式、数据有效性（NaN/inf）和强制二维。

34.4.1 类型定义详解

// 类型即图纸，先理解数据结构
def check_pairwise_arrays(
    X,
    Y,
    *,
    precomputed=False,
    dtype="infer_float",
    accept_sparse="csr",
    ensure_all_finite=True,
    ensure_2d=True,
    copy=False,
):

这个函数是所有成对距离/核计算的“守门人”，它的参数签名揭示了设计意图：

• precomputed：标记是否为预计算距离矩阵模式

• dtype="infer_float"：自动推断合适的浮点类型，兼容 float32/float64

• accept_sparse="csr"：默认接受 CSR 稀疏格式，其他格式会被转换

• ensure_all_finite：支持三种模式（True/False/'allow-nan'），为缺失值计算留口子

• ensure_2d：自定义度量时可关闭二维强制，允许字符串列表等非数值输入

34.4.2 逐行解析核心逻辑

源码路径：sklearn/metrics/pairwise.py - check_pairwise_arrays()（第70-150行）

def check_pairwise_arrays(
    X,
    Y,
    *,
    precomputed=False,
    dtype="infer_float",
    accept_sparse="csr",
    ensure_all_finite=True,
    ensure_2d=True,
    copy=False,
):
    xp, _ = get_namespace(X, Y)                                    # ① 获取数组命名空间（支持Array API）
    X, Y, dtype_float = _find_floating_dtype_allow_sparse(X, Y, xp=xp)  # ② 统一浮点类型，支持稀疏

    estimator = "check_pairwise_arrays"
    if dtype == "infer_float":
        dtype = dtype_float                                         # ③ 自动推断浮点类型

    if Y is X or Y is None:                                         # ④ Y为None或同对象时共享引用
        X = Y = check_array(
            X,
            accept_sparse=accept_sparse,
            dtype=dtype,
            copy=copy,
            ensure_all_finite=ensure_all_finite,
            estimator=estimator,
            ensure_2d=ensure_2d,
        )
    else:
        X = check_array(                                            # ⑤ 分别验证X和Y
            X,
            accept_sparse=accept_sparse,
            dtype=dtype,
            copy=copy,
            ensure_all_finite=ensure_all_finite,
            estimator=estimator,
            ensure_2d=ensure_2d,
        )
        Y = check_array(
            Y,
            accept_sparse=accept_sparse,
            dtype=dtype,
            copy=copy,
            ensure_all_finite=ensure_all_finite,
            estimator=estimator,
            ensure_2d=ensure_2d,
        )

    if precomputed:                                                 # ⑥ 预计算模式：检查形状兼容性
        if X.shape[1] != Y.shape[0]:
            raise ValueError(
                "Precomputed metric requires shape "
                "(n_queries, n_indexed). Got (%d, %d) "
                "for %d indexed." % (X.shape[0], X.shape[1], Y.shape[0])
            )
    elif ensure_2d and X.shape[1] != Y.shape[1]:                    # ⑦ 普通模式：检查特征数一致
        raise ValueError(
            "Incompatible dimension for X and Y matrices: "
            "X.shape[1] == %d while Y.shape[1] == %d" % (X.shape[1], Y.shape[1])
        )

    return X, Y

这段代码定义了成对计算的输入契约：统一类型、共享引用避免拷贝、区分预计算与特征矩阵两种模式的形状检查。它是后续所有距离/核函数的“类型安全基石”。

34.4.3 数据流图

flowchart TD A[输入 X, Y] --> B{_find_floating_dtype_allow_sparse} B --> C[统一浮点类型 dtype_float] C --> D{dtype == infer_float?} D -->|是| E[使用 dtype_float] D -->|否| F[使用用户指定 dtype] E --> G{Y is X or Y is None?} F --> G G -->|是| H[单次 check_array, X=Y共享引用] G -->|否| I[分别 check_array X 和 Y] H --> J{precomputed?} I --> J J -->|是| K[检查 X.shape[1] == Y.shape[0]] J -->|否| L[检查 X.shape[1] == Y.shape[1] if ensure_2d] K --> M[返回验证后的 X, Y] L --> M

34.5 常用距离度量实现 —— 欧氏、曼哈顿与余弦距离的高效计算

euclidean_distances 的向量化技巧：利用 ||x-y||² = ||x||² + ||y||² - 2<x,y> 的恒等式，避免直接计算差向量，提升稀疏数据效率。支持预计算范数（X_norm_squared, Y_norm_squared）以加速重复计算。对 float32 输入进行分块上 cast 到 float64 以减轻灾难性消舍（catastrophic cancellation）的影响。使用 safe_sparse_dot 处理稀疏-稀疏、稀疏-密集等乘法，自动选择最优路径。通过 row_norms 高效计算点积 self-dot（X*X.sum(axis=1)）。

manhattan_distances 的稀疏优化：当输入为 CSR 格式时，调用 Cython 加速的 _sparse_manhattan 避免稠密化。对稀疏矩阵执行 sum_duplicates() 确保索引唯一且排序，这是 CSR 规范形式的前提。仅在非稀疏路径使用 scipy.spatial.distance.cdist 或 Array API 实现的逐块 L1 距离。

cosine_distances 通过余弦相似度实现：利用已有的 cosine_similarity 函数，计算 1 - S，其中 S 为归一化点积。自动修正对角线为 0（向自身的余弦距离应为 0），处理浮点误差导致的微小偏差。

34.5.1 类型定义详解

// 类型即图纸，先理解数据结构
@validate_params(
    {
        "X": ["array-like", "sparse matrix"],
        "Y": ["array-like", "sparse matrix", None],
        "Y_norm_squared": ["array-like", None],
        "squared": ["boolean"],
        "X_norm_squared": ["array-like", None],
    },
    prefer_skip_nested_validation=True,
)
def euclidean_distances(
    X, Y=None, *, Y_norm_squared=None, squared=False, X_norm_squared=None
):

参数设计体现了性能优化考量：

• X_norm_squared / Y_norm_squared：预计算的 ||x||²，形状要求灵活支持 (n,)、(n,1)、(1,n)

• squared：返回平方距离避免 sqrt 开销，适合只需比较大小的场景

• prefer_skip_nested_validation：跳过嵌套验证，因为内部会调用 check_pairwise_arrays

34.5.2 逐行解析核心逻辑

源码路径：sklearn/metrics/pairwise.py - euclidean_distances()（第153-220行）

def euclidean_distances(
    X, Y=None, *, Y_norm_squared=None, squared=False, X_norm_squared=None
):
    xp, _ = get_namespace(X, Y)
    X, Y = check_pairwise_arrays(X, Y)                              # ① 输入验证与类型统一

    if X_norm_squared is not None:                                  # ② 验证并重塑预计算范数
        X_norm_squared = check_array(X_norm_squared, ensure_2d=False)
        original_shape = X_norm_squared.shape
        if X_norm_squared.shape == (X.shape[0],):
            X_norm_squared = xp.reshape(X_norm_squared, (-1, 1))
        if X_norm_squared.shape == (1, X.shape[0]):
            X_norm_squared = X_norm_squared.T
        if X_norm_squared.shape != (X.shape[0], 1):
            raise ValueError(...)

    if Y_norm_squared is not None:                                  # ③ 同理处理 Y 的预计算范数
        Y_norm_squared = check_array(Y_norm_squared, ensure_2d=False)
        original_shape = Y_norm_squared.shape
        if Y_norm_squared.shape == (Y.shape[0],):
            Y_norm_squared = xp.reshape(Y_norm_squared, (1, -1))
        if Y_norm_squared.shape == (Y.shape[0], 1):
            Y_norm_squared = Y_norm_squared.T
        if Y_norm_squared.shape != (1, Y.shape[0]):
            raise ValueError(...)

    return _euclidean_distances(X, Y, X_norm_squared, Y_norm_squared, squared)  # ④ 委托核心计算

这段代码展示了主入口的参数预处理：验证输入、规范化预计算范数形状、委托给内部核心函数。

源码路径：sklearn/metrics/pairwise.py - _euclidean_distances()（第222-270行）

def _euclidean_distances(X, Y, X_norm_squared=None, Y_norm_squared=None, squared=False):
    xp, _, device_ = get_namespace_and_device(X, Y)
    if X_norm_squared is not None and X_norm_squared.dtype != xp.float32:
        XX = xp.reshape(X_norm_squared, (-1, 1))                    # ① 使用预计算范数（非float32）
    elif X.dtype != xp.float32:
        XX = row_norms(X, squared=True)[:, None]                    # ② 现场计算范数（非float32）
    else:
        XX = None                                                   # ③ float32标记为None，稍后分块上 cast

    if Y is X:
        YY = None if XX is None else XX.T                           # ④ X=Y时复用范数
    else:
        if Y_norm_squared is not None and Y_norm_squared.dtype != xp.float32:
            YY = xp.reshape(Y_norm_squared, (1, -1))
        elif Y.dtype != xp.float32:
            YY = row_norms(Y, squared=True)[None, :]
        else:
            YY = None

    if X.dtype == xp.float32 or Y.dtype == xp.float32:              # ⑤ float32 走分块上 cast 路径
        distances = _euclidean_distances_upcast(X, XX, Y, YY)
    else:
        distances = -2 * safe_sparse_dot(X, Y.T, dense_output=True) # ⑥ 向量化核心：-2<X,Y>
        distances += XX                                             # ⑦ 加上 ||X||²
        distances += YY                                             # ⑧ 加上 ||Y||²

    xp_zero = xp.asarray(0, device=device_, dtype=distances.dtype)
    distances = _modify_in_place_if_numpy(                          # ⑨ 就地修正负值（数值误差）
        xp, xp.maximum, distances, xp_zero, out=distances
    )

    if X is Y:
        _fill_diagonal(distances, 0, xp=xp)                         # ⑩ 对角线归零（自距离为0）

    if squared:
        return distances

    distances = _modify_in_place_if_numpy(xp, xp.sqrt, distances, out=distances)  # ⑪ 开方
    return distances

这段代码实现了核心向量化公式 ||X-Y||² = ||X||² + ||Y||² - 2<X,Y>。关键设计点：

1. float32 特殊处理：避免灾难性消舍，分块上 cast 到 float64 计算

2. safe_sparse_dot：统一处理稠密/稀疏矩阵乘法，自动选择最优路径

3. 就地操作：_modify_in_place_if_numpy 在 NumPy 下避免中间数组分配

4. 对角线修正：浮点误差可能导致自距离为微小负数，强制归零

源码路径：sklearn/metrics/pairwise.py - manhattan_distances()（第400-440行）

def manhattan_distances(X, Y=None):
    X, Y = check_pairwise_arrays(X, Y)
    n_x, n_y = X.shape[0], Y.shape[0]

    if issparse(X) or issparse(Y):                                  # ① 稀疏路径：转 CSR 并去重索引
        X = csr_matrix(X, copy=False)
        Y = csr_matrix(Y, copy=False)
        X.sum_duplicates()  # 确保索引唯一且排序
        Y.sum_duplicates()
        D = np.zeros((n_x, n_y))
        _sparse_manhattan(X.data, X.indices, X.indptr,              # ② 调用 Cython 加速核心
                          Y.data, Y.indices, Y.indptr, D)
        return D

    xp, _, device_ = get_namespace_and_device(X, Y)

    if _is_numpy_namespace(xp):                                     # ③ NumPy：调用 SciPy 优化 cdist
        return distance.cdist(X, Y, "cityblock")

    # Array API 支持：逐块计算 L1 距离
    float_dtype = _find_matching_floating_dtype(X, Y, xp=xp)
    out = xp.empty((n_x, n_y), dtype=float_dtype, device=device_)
    batch_size = 1024
    for i in range(0, n_x, batch_size):
        i_end = min(i + batch_size, n_x)
        batch_X = X[i:i_end, ...]
        for j in range(0, n_y, batch_size):
            j_end = min(j + batch_size, n_y)
            batch_Y = Y[j:j_end, ...]
            block_dist = xp.sum(
                xp.abs(batch_X[:, None, :] - batch_Y[None, :, :]), axis=2
            )
            out[i:i_end, j:j_end] = block_dist

    return out

这段代码展示了分层优化策略：稀疏矩阵走 Cython 专用核心（避免稠密化），NumPy 调用 SciPy 高度优化的 cdist，其他 Array API 后端回退到逐块广播计算。

34.5.3 常用距离度量核心实现策略对比

以下表格展示了三种主要距离度量的实现策略对比：

| 维度 | euclidean_distances | manhattan_distances | cosine_distances |

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

| 核心公式 | ||X||² + ||Y||² - 2<X,Y> | Σ|X_ij - Y_ij| | 1 - <X_norm, Y_norm> |

| 稀疏优化 | safe_sparse_dot 自动分发 | _sparse_manhattan (Cython) | 先 normalize 再 linear_kernel |

| float32 处理 | 分块上 cast 到 float64 | 无特殊处理 | 复用 cosine_similarity 逻辑 |

| 预计算支持 | X_norm_squared, Y_norm_squared | 无 | 无 |

| Array API | 完整支持 | 逐块广播回退 | 完整支持 |

| 对角线修正 | _fill_diagonal(distances, 0) | 无需（L1 自距离必为0） | _fill_diagonal(S, 0) |

34.6 缺失值与成对距离 —— 优雅处理不完整数据的策略

nan_euclidean_distances 的缺失值容忍机制：在计算欧氏距离时忽略任一方向有 NaN 的特征，并按剩余特征数重新加权。首先将 X 和 Y 中的 missing_values 置零，以便使用标准 euclidean_distances 计算部分距离。计算呈现坐标的乘积 (present_X · present_Y.T) 作为有效特征对的计数。将距离平方除以有效特征数，再乘以总特征数以恢复正确 scale：dist *= n_features。当某对样本在所有特征上均缺失或无共同非缺失特征时，返回 NaN。配合 check_pairwise_arrays 中的 ensure_all_finite='allow-nan' 以允许 NaN 通过验证。

paired_distances 的逐元素配对范式：专门用于计算 (X[i], Y[i]) 之间的距离，而非所有配对。内部调用 check_paired_arrays 确保 X 和 Y 形状完全相同。支持三种内置度量（'euclidean'/'l2'、'manhattan'/'l1'、'cityblock'）及自定义可调用函数。余弦配对距离利用已归一化向量的 L2 范数：0.5 * ||X_norm - Y_norm||²。

34.6.1 逐行解析核心逻辑

源码路径：sklearn/metrics/pairwise.py - nan_euclidean_distances()（第272-340行）

def nan_euclidean_distances(
    X, Y=None, *, squared=False, missing_values=np.nan, copy=True
):
    ensure_all_finite = "allow-nan" if is_scalar_nan(missing_values) else True  # ① 允许 NaN 通过验证
    X, Y = check_pairwise_arrays(
        X, Y, accept_sparse=False, ensure_all_finite=ensure_all_finite, copy=copy
    )
    missing_X = _get_mask(X, missing_values)                                      # ② 获取缺失值掩码
    missing_Y = missing_X if Y is X else _get_mask(Y, missing_values)

    X[missing_X] = 0                                                              # ③ 缺失值置零
    Y[missing_Y] = 0

    distances = euclidean_distances(X, Y, squared=True)                           # ④ 标准欧氏距离（平方）

    XX = X * X
    YY = Y * Y
    distances -= np.dot(XX, missing_Y.T)                                          # ⑤ 减去缺失特征贡献
    distances -= np.dot(missing_X, YY.T)

    np.clip(distances, 0, None, out=distances)                                    # ⑥ 数值修正

    if X is Y:
        np.fill_diagonal(distances, 0.0)                                          # ⑦ 对角线归零

    present_X = 1 - missing_X                                                     # ⑧ 计算有效特征计数
    present_Y = present_X if Y is X else ~missing_Y
    present_count = np.dot(present_X, present_Y.T)
    distances[present_count == 0] = np.nan                                        # ⑨ 全缺失对返回 NaN
    np.maximum(1, present_count, out=present_count)                               # ⑩ 避免除零
    distances /= present_count                                                    # ⑪ 按有效特征数归一化
    distances *= X.shape[1]                                                       # ⑫ 恢复原始尺度

    if not squared:
        np.sqrt(distances, out=distances)                                         # ⑬ 可选开方
    return distances

这段代码实现了论文中描述的加权欧氏距离：dist(x,y) = sqrt(weight * sq_dist_present)，其中 weight = n_features / n_present。通过置零+标准计算+事后修正的三步走策略，优雅地复用了现有的 euclidean_distances 基础设施。

34.6.2 缺失值距离计算流程图

flowchart TD A[输入 X, Y, missing_values] --> B[check_pairwise_arrays allow-nan] B --> C[_get_mask 获取缺失值掩码] C --> D[缺失值置零] D --> E[euclidean_distances 计算基础平方距离] E --> F[计算 XX=X², YY=Y²] F --> G[减去缺失特征贡献] G --> H[clip 负值修正] H --> I{X is Y?} I -->|是| J[fill_diagonal 0] I -->|否| K[计算 present_X, present_Y] J --> K K --> L[present_count = present_X · present_Y.T] L --> M{present_count == 0?} M -->|是| N[设为 NaN] M -->|否| O[maximum(1, present_count)] N --> P[distances /= present_count] O --> P P --> Q[distances *= n_features] Q --> R{squared?} R -->|否| S[sqrt 开方] R -->|是| T[返回 distances] S --> T

源码路径：sklearn/metrics/pairwise.py - paired_distances() 及其变体（第492-560行）

def paired_euclidean_distances(X, Y):
    X, Y = check_paired_arrays(X, Y)                              # ① 严格形状检查
    return row_norms(X - Y)                                       # ② 逐行 L2 范数

def paired_manhattan_distances(X, Y):
    X, Y = check_paired_arrays(X, Y)
    xp, _ = get_namespace(X, Y)
    diff = X - Y
    if issparse(diff):                                            # ③ 稀疏：取绝对值求和
        diff.data = np.abs(diff.data)
        return np.squeeze(np.array(diff.sum(axis=1)))
    else:
        return xp.sum(xp.abs(diff), axis=-1)                      # ④ 稠密：逐元素绝对值求和

def paired_cosine_distances(X, Y):
    X, Y = check_paired_arrays(X, Y)
    return 0.5 * row_norms(normalize(X) - normalize(Y), squared=True)  # ⑤ 余弦距离 = 0.5 * ||归一化差||²

PAIRED_DISTANCES = {
    "cosine": paired_cosine_distances,
    "euclidean": paired_euclidean_distances,
    "l2": paired_euclidean_distances,
    "l1": paired_manhattan_distances,
    "manhattan": paired_manhattan_distances,
    "cityblock": paired_manhattan_distances,
}

def paired_distances(X, Y, *, metric="euclidean", **kwds):
    if metric in PAIRED_DISTANCES:
        func = PAIRED_DISTANCES[metric]
        return func(X, Y)
    elif callable(metric):                                        # ⑥ 自定义度量回退
        X, Y = check_paired_arrays(X, Y)
        distances = np.zeros(len(X))
        for i in range(len(X)):
            distances[i] = metric(X[i], Y[i])
        return distances

这组代码展示了“成对距离”与“全配对距离”的根本区别：前者是 O(n) 的逐行操作，后者是 O(n²) 的矩阵操作。余弦配对距离利用了数学恒等式：1 - cos(x,y) = 0.5 * ||x/||x|| - y/||y||||²，将其转化为已实现的欧氏距离逻辑。

34.6.3 成对距离计算流程图

flowchart TD A[输入 X, Y, metric] --> B{metric in PAIRED_DISTANCES?} B -->|是| C[check_paired_arrays 严格形状检查] B -->|否| D{callable metric?} D -->|否| E[抛出错误] D -->|是| C C --> F{metric类型} F -->|euclidean| G[row_norms(X - Y)] F -->|manhattan| H[issparse? → abs(data).sum : sum(abs(X-Y))] F -->|cosine| I[0.5 * row_norms(normalize(X) - normalize(Y))²] F -->|callable| J[循环逐行调用 metric(X[i], Y[i])] G --> K[返回 distances] H --> K I --> K J --> K

34.7 核函数家族实现 —— 从线核到高斯核的多样性支持

rbf_kernel 的高效实现：先计算平方欧氏距离矩阵，再乘以 -gamma 并逐元素求 exp。默认 gamma = 1 / n_features，可通过参数显式设置。使用 _modify_in_place_if_numpy 就地执行 exp 操作以减少内存分配（仅在 NumPy 生效）。支持稀疏输入通过 euclidean_distances 的稀疏路径。

laplacian_kernel 的 L1 距离基础：直接使用 manhattan_distances 计算 L1 距离，再应用指数衰减。同样默认 gamma = 1 / n_features，要求 gamma > 0。在 NumPy 下使用 np.exp(K, K) 就地求 exp，其他后端使用 xp.exp。

polynomial_kernel 和 sigmoid_kernel 的统一框架：均基于线性核 K = <X,Y> 的多项式或 sigmoid 变换。多项式核：K = (gamma * <X,Y> + coef0)^degree。Sigmoid 核：K = tanh(gamma * <X,Y> + coef0)。均使用 safe_sparse_dot 计算内积，支持稀疏和密集混合输入。

linear_kernel 和 cosine_similarity 的稀疏友好实现：前者直接返回 X @ Y.T，后者先 L2 归一化再做线性核。均支持 dense_output 参数控制稀疏输入下的输出格式。

34.7.1 逐行解析核心逻辑

源码路径：sklearn/metrics/pairwise.py - rbf_kernel()（第620-640行）

// 源码路径：sklearn/metrics/pairwise.py - rbf_kernel()（第620-640行）
def rbf_kernel(X, Y=None, gamma=None):
    xp, _ = get_namespace(X, Y)
    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 1.0 / X.shape[1]                                  # ① 默认 gamma = 1/n_features

    K = euclidean_distances(X, Y, squared=True)                   # ② 复用平方欧氏距离
    K *= -gamma                                                   # ③ 乘以 -gamma
    K = _modify_in_place_if_numpy(xp, xp.exp, K, out=K)           # ④ 就地 exp（NumPy优化）
    return K

源码路径：sklearn/metrics/pairwise.py - laplacian_kernel()（第650-670行）

// 源码路径：sklearn/metrics/pairwise.py - laplacian_kernel()（第650-670行）
def laplacian_kernel(X, Y=None, gamma=None):
    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 1.0 / X.shape[1]

    K = -gamma * manhattan_distances(X, Y)                        # ① 基于 L1 距离
    xp, _ = get_namespace(X, Y)
    if _is_numpy_namespace(xp):
        np.exp(K, K)                                              # ② NumPy 就地 exp
    else:
        K = xp.exp(K)
    return K

源码路径：sklearn/metrics/pairwise.py - polynomial_kernel() 与 sigmoid_kernel()（第580-620行）

// 源码路径：sklearn/metrics/pairwise.py - polynomial_kernel() 与 sigmoid_kernel()（第580-620行）
def polynomial_kernel(X, Y=None, degree=3, gamma=None, coef0=1):
    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 1.0 / X.shape[1]

    K = safe_sparse_dot(X, Y.T, dense_output=True)                # ① 内积 <X,Y>
    K *= gamma
    K += coef0
    K **= degree                                                  # ② (gamma*<X,Y>+coef0)^degree
    return K

def sigmoid_kernel(X, Y=None, gamma=None, coef0=1):
    xp, _ = get_namespace(X, Y)
    X, Y = check_pairwise_arrays(X, Y)

    if gamma is None:
        gamma = 1.0 / X.shape[1]

    K = safe_sparse_dot(X, Y.T, dense_output=True)                # ① 内积 <X,Y>
    K *= gamma
    K += coef0
    K = _modify_in_place_if_numpy(xp, xp.tanh, K, out=K)          # ② 就地 tanh
    return K

源码路径：sklearn/metrics/pairwise.py - linear_kernel() 与 cosine_similarity()（第560-580、680-700行）

// 源码路径：sklearn/metrics/pairwise.py - linear_kernel() 与 cosine_similarity()（第560-580、680-700行）
def linear_kernel(X, Y=None, dense_output=True):
    X, Y = check_pairwise_arrays(X, Y)
    return safe_sparse_dot(X, Y.T, dense_output=dense_output)     # ① 直接 X @ Y.T

def cosine_similarity(X, Y=None, dense_output=True):
    X, Y = check_pairwise_arrays(X, Y)

    X_normalized = normalize(X, copy=True)                        # ① L2 归一化
    if X is Y:
        Y_normalized = X_normalized
    else:
        Y_normalized = normalize(Y, copy=True)

    K = safe_sparse_dot(X_normalized, Y_normalized.T, dense_output=dense_output)  # ② 归一化后内积
    return K

34.7.2 核函数数学本质与实现要点对比

以下表格展示了核心核函数的数学定义与实现要点对比：

| 核函数 | 数学公式 | 核心实现路径 | 默认 gamma | 就地优化 |

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

| linear | <X, Y> | safe_sparse_dot(X, Y.T) | - | 无 |

| cosine | <X/||X||, Y/||Y||> | normalize → linear_kernel | - | 无 |

| polynomial | (γ<X,Y> + c₀)^d | safe_sparse_dot → 多项式 | 1/n_features | 无 |

| sigmoid | tanh(γ<X,Y> + c₀) | safe_sparse_dot → tanh | 1/n_features | NumPy 就地 tanh |

| rbf | exp(-γ||X-Y||²) | euclidean_distances(squared=True) → exp | 1/n_features | NumPy 就地 exp |

| laplacian | exp(-γ||X-Y||₁) | manhattan_distances → exp | 1/n_features | NumPy 就地 exp |

| chi2 | exp(-γ Σ(x-y)²/(x+y)) | _chi2_kernel_fast (Cython) | 1.0 | Cython 层面优化 |

| additive_chi2 | -Σ(x-y)²/(x+y) | _chi2_kernel_fast (Cython) | - | Cython 层面优化 |

34.7.3 核函数计算流程图

flowchart TD A[输入 X, Y, kernel params] --> B{metric in PAIRWISE_KERNEL_FUNCTIONS?} B -->|是| C[check_pairwise_arrays] B -->|否| D[callable 或 precomputed] C --> E{kernel 类型} E -->|linear| F[safe_sparse_dot(X, Y.T)] E -->|cosine| G[normalize → safe_sparse_dot] E -->|polynomial| H[safe_sparse_dot → *gamma + coef0 → **degree] E -->|sigmoid| I[safe_sparse_dot → *gamma + coef0 → tanh] E -->|rbf| J[euclidean_distances(squared=True) → *-gamma → exp] E -->|laplacian| K[manhattan_distances → *-gamma → exp] E -->|chi2| L[_chi2_kernel_fast Cython] E -->|additive_chi2| M[_chi2_kernel_fast Cython] F --> N[返回 K] G --> N H --> N I --> N J --> N K --> N L --> N M --> N

34.8 分块计算与 Cython 加速 —— 大规模数据的内存与速度平衡术

pairwise_distances_chunked 的分块生成器机制：通过 working_memory 控制单次处理的行数，仅在内存允许时生成距离矩阵的垂直块。使用 get_chunk_n_rows 根据 working_memory、特征数和样本数计算可处理的最大行数。当 reduce_func 为 None 时，直接 yield 距离块；否则对每块应用 reduce_func 并验证输出长度。特别处理欧氏距离的对角线：在 yield 前将对角线位置设为 0（利用 flat strides 和步长）。支持预计算距离矩阵（metric='precomputed'）的直通返回。

_pairwise_fast.pyx 中的 OpenMP 并行核心：_sparse_manhattan 和 _chi2_kernel_fast 利用 Cython 并行循环。_sparse_manhattan 按行并行（prange 对 X 的行），每行内部串行遍历 Y 的行以避免竞争。使用 _openmp_effective_n_threads() 尊重 threadpoolctl 或 OMP_NUM_THREADS 环境变量。显式 nogil 以释放 GIL，实现真正的多线程并行。

pairwise_distances 的多线程调度器 (_parallel_pairwise)：将距离计算分解为线程安全的块写入。使用 joblib.Parallel 与 backend='threading' 避免进程间数据开销。输出矩阵先转置（.T）使每线程写入连续内存块，完成后再转回。特别处理欧氏距离的对角线清零（在转置后的矩阵上操作等效于原矩阵的对角线）。

34.8.1 逐行解析核心逻辑

源码路径：sklearn/metrics/pairwise.py - pairwise_distances_chunked()（第840-920行）

def pairwise_distances_chunked(
    X,
    Y=None,
    *,
    reduce_func=None,
    metric="euclidean",
    n_jobs=None,
    working_memory=None,
    **kwds,
):
    n_samples_X = _num_samples(X)
    if metric == "precomputed":
        slices = (slice(0, n_samples_X),)
    else:
        if Y is None:
            Y = X
        chunk_n_rows = get_chunk_n_rows(                            # ① 根据内存预算计算分块行数
            row_bytes=8 * _num_samples(Y),                          # 每行 8 字节 (float64)
            max_n_rows=n_samples_X,
            working_memory=working_memory,
        )
        slices = gen_batches(n_samples_X, chunk_n_rows)             # ② 生成切片迭代器

    params = _precompute_metric_params(X, Y, metric=metric, **kwds) # ③ 预计算数据相关参数
    kwds.update(**params)

    for sl in slices:
        if sl.start == 0 and sl.stop == n_samples_X:
            X_chunk = X  # 完整 X 时启用优化路径 (X is Y)
        else:
            X_chunk = X[sl]
        D_chunk = pairwise_distances(X_chunk, Y, metric=metric, n_jobs=n_jobs, **kwds)
        if (X is Y or Y is None) and PAIRWISE_DISTANCE_FUNCTIONS.get(
            metric, None
        ) is euclidean_distances:
            D_chunk.flat[sl.start :: _num_samples(X) + 1] = 0       # ④ 对角线归零（利用 flat stride）
        if reduce_func is not None:
            chunk_size = D_chunk.shape[0]
            D_chunk = reduce_func(D_chunk, sl.start)                # ⑤ 应用归约函数
            _check_chunk_size(D_chunk, chunk_size)                  # ⑥ 验证输出形状
        yield D_chunk                                               # ⑦ 生成器产出

这段代码实现了“生产者-消费者”模式的分块计算：按内存预算切片、并行计算每块、可选归约、逐块产出。对角线清零使用了 NumPy flat 索引技巧：flat[start :: n+1] 精确命中对角线元素。

34.8.2 分块生成器流程图

flowchart TD A[pairwise_distances_chunked] --> B{metric == precomputed?} B -->|是| C[单切片全量返回] B -->|否| D[Y is None?] D -->|是| E[Y = X] D -->|否| F[保持 Y] E --> G[get_chunk_n_rows 计算分块行数] F --> G G --> H[gen_batches 生成切片迭代器] H --> I[预计算 metric 参数] I --> J[循环遍历切片] J --> K{sl覆盖全X?} K -->|是| L[X_chunk = X 启用优化路径] K -->|否| M[X_chunk = X[sl]] L --> N[pairwise_distances 计算距离块] M --> N N --> O{X is Y 且 euclidean?} O -->|是| P[D_chunk.flat[start::n+1] = 0] O -->|否| Q[reduce_func?] P --> Q Q -->|是| R[应用 reduce_func] R --> S[_check_chunk_size 验证] Q -->|否| T[yield D_chunk] S --> T

源码路径：sklearn/metrics/pairwise.py - _parallel_pairwise()（第922-960行）

def _parallel_pairwise(X, Y, func, n_jobs, **kwds):
    xp, _, device = get_namespace_and_device(X, Y)
    X, Y, dtype_float = _find_floating_dtype_allow_sparse(X, Y, xp=xp)

    if Y is None:
        Y = X

    if effective_n_jobs(n_jobs) == 1:
        return func(X, Y, **kwds)

    fd = delayed(_transposed_dist_wrapper)
    ret = xp.empty((X.shape[0], Y.shape[0]), device=device, dtype=dtype_float).T  # ① 转置分配：列主序
    Parallel(backend="threading", n_jobs=n_jobs)(                   # ② 线程级并行（避免进程开销）
        fd(func, ret, s, X, Y[s, ...], **kwds)
        for s in gen_even_slices(_num_samples(Y), effective_n_jobs(n_jobs))
    )

    if (X is Y or Y is None) and func is euclidean_distances:
        _fill_diagonal(ret, 0, xp=xp)                               # ③ 转置矩阵上清零对角线

    return ret.T                                                    # ④ 转置回行主序

这段代码揭示了并行写入的关键优化：输出矩阵转置分配。默认 C-contiguous 矩阵按行存储，不同线程写入不同行会导致伪共享（false sharing）。转置后，每个线程处理一列（即原矩阵的一行），写入连续内存，消除缓存行竞争。完成后再 .T 转回，view 操作零拷贝。

源码路径：sklearn/metrics/_pairwise_fast.pyx - _sparse_manhattan()（第30-80行）

def _sparse_manhattan(
    const floating[::1] X_data,
    const int[:] X_indices,
    const int[:] X_indptr,
    const floating[::1] Y_data,
    const int[:] Y_indices,
    const int[:] Y_indptr,
    double[:, ::1] D,
):
    cdef intp_t px, py, i, j, ix, iy
    cdef double d = 0.0

    cdef int m = D.shape[0]
    cdef int n = D.shape[1]

    cdef int num_threads = _openmp_effective_n_threads()

    for px in prange(m, nogil=True, num_threads=num_threads):       # ① OpenMP 并行遍历 X 的行
        X_indptr_end = X_indptr[px + 1]
        for py in range(n):                                         # ② 串行遍历 Y 的行
            Y_indptr_end = Y_indptr[py + 1]
            i = X_indptr[px]
            j = Y_indptr[py]
            d = 0.0
            while i < X_indptr_end and j < Y_indptr_end:            # ③ 归并式遍历两行稀疏向量
                ix = X_indices[i]
                iy = Y_indices[j]

                if ix == iy:                                        # 同列：差值绝对值
                    d = d + fabs(X_data[i] - Y_data[j])
                    i = i + 1
                    j = j + 1
                elif ix < iy:                                       # 仅 X 有值
                    d = d + fabs(X_data[i])
                    i = i + 1
                else:                                               # 仅 Y 有值
                    d = d + fabs(Y_data[j])
                    j = j + 1

            if i == X_indptr_end:                                   # ④ 处理剩余非零元素
                while j < Y_indptr_end:
                    d = d + fabs(Y_data[j])
                    j = j + 1
            else:
                while i < X_indptr_end:
                    d = d + fabs(X_data[i])
                    i = i + 1

            D[px, py] = d

这段 Cython 代码实现了稀疏矩阵曼哈顿距离的归并式算法：利用 CSR 格式的有序索引特性，双指针同时遍历两行的非零元素，类似归并排序的 merge 步骤。关键设计：

1. prange 并行外层循环（X 的行），每线程独立写入 D[px, :] 无竞争

2. nogil 释放 GIL，实现真正的多线程并行

3. 内存视图 floating[::1] 等避免 Python 对象开销，直接操作 C 数组

4. 避免稠密化：仅遍历非零元素，复杂度 O(nnz) 而非 O(n_features)

源码路径：sklearn/metrics/_pairwise_fast.pyx - _chi2_kernel_fast()（第10-30行）

def _chi2_kernel_fast(floating[:, :] X,
                      floating[:, :] Y,
                      floating[:, :] result):
    cdef intp_t i, j, k
    cdef intp_t n_samples_X = X.shape[0]
    cdef intp_t n_samples_Y = Y.shape[0]
    cdef intp_t n_features = X.shape[1]
    cdef double res, nom, denom

    with nogil:
        for i in range(n_samples_X):
            for j in range(n_samples_Y):
                res = 0
                for k in range(n_features):
                    denom = (X[i, k] - Y[j, k])
                    nom = (X[i, k] + Y[j, k])
                    if nom != 0:
                        res += denom * denom / nom
                result[i, j] = -res

这是卡方核的朴素三重循环实现，但通过 nogil 和类型声明获得数量级加速。注意：当前版本未使用 prange 并行化，因为 chi2 核主要用于直方图特征（通常维度较低），并行开销可能抵消收益。

34.8.3 并行计算架构图

flowchart TD A[pairwise_distances / pairwise_kernels] --> B{数据规模 vs working_memory} B -->|小规模| C[_parallel_pairwise 单块并行] B -->|大规模| D[pairwise_distances_chunked 生成器] C --> E{有效线程数 > 1?} E -->|是| F[joblib.Parallel threading backend] E -->|否| G[直接调用 func] F --> H[输出矩阵转置分配 .T] H --> I[各线程写入连续列] I --> J[完成后 .T 转回] D --> K[get_chunk_n_rows 计算分块大小] K --> L[gen_batches 生成行切片] L --> M[循环计算每个垂直块] M --> N{reduce_func?} N -->|无| O[yield 原始距离块] N -->|有| P[应用 reduce_func] P --> Q[_check_chunk_size 验证] Q --> O subgraph Cython加速层 R[_sparse_manhattan] --> S[prange并行X行 + nogil] S --> T[归并式遍历CSR非零元素] U[_chi2_kernel_fast] --> V[nogil三重循环] end style C fill:#e1f5fe style D fill:#fff3e0 style F fill:#e8f5e9 style H fill:#fce4ec style R fill:#f3e5f5 style U fill:#f3e5f5

34.9 设计中的取舍

为什么不用 scipy.spatial.distance.cdist 处理所有稠密距离？scipy 的 cdist 虽然高度优化，但 scikit-learn 保留了自家实现（euclidean_distances、cosine_distances 等）主要有三个原因：1. 稀疏支持：scipy cdist 不支持稀疏矩阵，而 sklearn 的 safe_sparse_dot 统一处理稠密/稀疏混合场景。2. 预计算范数复用：euclidean_distances 的 X_norm_squared 参数允许在多次计算中复用范数（如 KMeans 迭代），scipy 无此接口。3. Array API 兼容：sklearn 正在向 Array API 标准迁移，自家实现可控制 CuPy/PyTorch/JAX 等后端行为。

为什么 float32 欧氏距离要分块上 cast 到 float64？直接用 float32 计算 ||x||² + ||y||² - 2<x,y> 会遭遇灾难性消舍：当两向量接近时，大数相减导致有效数字丢失。分块上 cast 到 float64 计算核心表达式，再 cast 回 float32 存储，在内存占用（不全量转换）和数值精度之间取得平衡。_euclidean_distances_upcast 根据矩阵密度动态计算 batch_size，限制额外内存约 10%。

为什么 _parallel_pairwise 强制使用 threading backend？距离计算是内存带宽密集型任务而非 CPU 密集型：- 进程池需要序列化/反序列化大矩阵（pickle 开销巨大）。- 线程共享内存，零拷贝传递数据。- 底层 BLAS (OpenBLAS/MKL) 已释放 GIL，多线程可利用多核。- threadpoolctl 可协调 sklearn 与 BLAS 的线程数，避免过度订阅。

为什么稀疏曼哈顿距离要求 sum_duplicates()？CSR 格式允许重复索引（非规范形式）。sum_duplicates() 就地合并重复索引并排序，这是归并式双指针算法的前置条件：要求 X.indices 和 Y.indices 严格单调递增。不调用会导致算法漏掉重复列的贡献或陷入死循环。

为什么 dense_output 参数的 trade-off 是什么？| 场景 | dense_output=True (默认) | dense_output=False | |------|--------------------------|-------------------| | 稀疏输入 | 返回稠密矩阵，内存 O(n²) | 返回稀疏矩阵，内存 O(nnz) | | 下游兼容 | 兼容所有 sklearn 估计器 | 仅兼容支持稀疏的估计器 | | 性能 | 稠密矩阵运算更快（BLAS） | 避免稠密化内存峰值 |默认 True 是为了“开箱即用”的兼容性；大规模稀疏场景（如文本）建议设 False。

34.10 动手练习

34.10.1 练习 1：阅读核心距离计算实现

阅读 sklearn/metrics/pairwise.py 第153-270行，理解以下函数的实现：

1. euclidean_distances - 欧氏距离主入口

2. _euclidean_distances - 欧氏距离计算核心

3. _euclidean_distances_upcast - float32输入的分块上计算策略

回答问题：

• 该实现如何利用向量化身份 ||x-y||² = ||x||² + ||y||² - 2<x,y> 提高效率？

• 什么情况下会触发 _euclidean_distances_upcast 的分块计算路径？

• X_norm_squared 和 Y_norm_squared 参数的作用是什么？形状要求如何？

34.10.2 练习 2：探索稀疏优化与Cython加速

阅读 sklearn/metrics/pairwise.py 第400-440行（manhattan_distances）及 sklearn/metrics/_pairwise_fast.pyx 第30-80行（_sparse_manhattan），理解：

1. 稀疏矩阵的曼哈顿距离如何避免稠密化？

2. _sparse_manhattan 如何实现按行并行而不造成数据竞争？

3. nogil 关键字在此的作用是什么？

回答问题：

• 为什么在CSR格式下需要先调用 sum_duplicates()？

• 算法中使用的 prange 循环遍历的是哪个维度？为什么这样设计可以避免竞争？

• 如果将 nogil 删除会对性能产生什么影响？

34.10.3 练习 3：分块计算与内存管理实验

阅读 sklearn/metrics/pairwise.py 第840-920行（pairwise_distances_chunked）及第922-960行（_parallel_pairwise），理解：

1. 工作内存(working_memory)如何控制单次处理的行数？

2. 生成器机制如何帮助处理超大规模数据集？

3. _parallel_pairwise 如何通过矩阵转置实现连续内存写入？

回答问题：

• get_chunk_n_rows 函数的输入参数及其计算逻辑是什么？

• 当 X 是 Y 时，对角线清零在转置后的矩阵上操作为什么等效于原矩阵的对角线？

• 如果将 backend 改为 'loky' 会带来什么后果？为什么这里强制使用 'threading'？

34.11 本章小结

这一章中我们学习了 scikit-learn metrics 模块中 pairwise 子模块的核心实现。首先深入剖析了 check_pairwise_arrays 作为统一输入验证层的设计，它处理类型推断、引用共享、预计算模式与特征矩阵模式的形状检查差异。其次详细解读了三大基础距离度量的实现：欧氏距离利用向量化恒等式避免显式差向量计算，并通过分块上 cast 解决 float32 精度问题；曼哈顿距离在 CSR 稀疏格式下调用 Cython 归并算法避免稠密化；余弦距离复用余弦相似度并修正对角线。接着探讨了缺失值距离（nan_euclidean_distances）的置零-计算-加权修正三步走策略，以及成对距离（paired_distances）的逐行 O(n) 计算范式。随后系统梳理了核函数家族：RBF 核基于平方欧氏距离，拉普拉斯核基于曼哈顿距离，多项式与 Sigmoid 核基于线性内积变换，线性核与余弦核形成基础对，均通过 safe_sparse_dot 统一稀疏/稠密路径并支持就地操作优化。最后深入分析了大规模数据应对方案：pairwise_distances_chunked 生成器按内存预算分块产出，_parallel_pairwise 通过矩阵转置实现线程安全的连续内存写入，Cython 层 _sparse_manhattan 与 _chi2_kernel_fast 利用 OpenMP prange + nogil 实现真正的多线程加速。

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

| 概念 | 解释 |

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

| check_pairwise_arrays | 输入验证哨兵，确保X,Y形状匹配且无非法值 |

| euclidean_distances | 利用\|\|x-y\|\|²=\|\|x\|\|²+\|\|y\|\|²-2<x,y>向量化计算，支持稀疏加速 |

| manhattan_distances | L1距离，CSR格式下调用_sparse_manhattan避免稠密化 |

| cosine_distances | 1 - cosine_similarity，自动修正对角线为0 |

| nan_euclidean_distances | 忽略NaN特征并按剩余特征数重新加权距离 |

| paired_*_distances | 逐元素配对距离，输入形状必须完全相同 |

| rbf_kernel | exp(-gamma * \|\|x-y\|\|²)，默认gamma=1/n_features |

| laplacian_kernel | exp(-gamma * \|\|x-y\|\|₁)，基于曼哈顿距离 |

| polynomial_kernel | (gamma*<x,y>+coef0)^degree，支持稀疏密集混合输入 |

| sigmoid_kernel | tanh(gamma*<x,y>+coef0)，NumPy下就地tanh优化 |

| linear_kernel / cosine_similarity | 线性核=X@Y.T，余弦核=先L2归一化再线性核 |

| pairwise_distances_chunked | 生成器分块产出，controlled by working_memory参数 |

| _parallel_pairwise | joblib多线程调度，输出矩阵转置后写入保证连续内存访问 |

| _sparse_manhattan | CSR矩阵L1距离，按行并行(prange)避免数据竞争 |

| _chi2_kernel_fast | Cython加速卡方核计算，nogil+prange实现真正并行 |

| dense_output参数 | 控制稀疏输入下是否强制返回密集输出（默认True） |

| metric='precomputed' | 直接返回输入距离/核矩阵，用于自定义预计算场景 |

下一章中，我们将学习成对距离归约分派器 —— 高性能邻居搜索的“智能调度中枢”，理解 sklearn 如何通过分派器模式在底层调度高性能 Cython 实现完成 k 近邻与半径邻居搜索。

第 35 章 —— 成对距离归约分派器 —— 高性能邻居搜索的“智能调度中枢”

35.1 学习目标

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

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

• 理解成对距离归约分派器的架构设计及其在高性能邻居搜索中的作用

• 掌握 BaseDistancesReductionDispatcher 抽象基类的接口规范与实现要点

• 熟悉 ArgKmin 与 RadiusNeighbors 分派器在 k 近邻与半径邻居搜索中的具体分派逻辑

• 了解类模式归约（ArgKminClassMode、RadiusNeighborsClassMode）如何实现加权投票与类别预测

• 掌握 sqeuclidean_row_norms 函数在欧氏距离优化中的预计算作用及其实现细节

• 掌握测试框架中用于验证分派器行为的辅助函数及断言方法

• 熟悉测试用例如何验证分派器在不同数据类型、内存布局及并行策略下的一致性与正确性

35.2 生活类比

成对距离归约分派器可以被想象成一座 智能交通调度中心：

• 调度指挥台（BaseDistancesReductionDispatcher）负责制定所有车辆调度的统一规则与接口；

• 车辆检查站（is_usable_for）确保只有符合排放标准（float32/64）且车道顺序正确（C‑contiguous / CSR）的车辆才能进入系统；

• k 近邻调度员（ArgKmin）负责为每辆车找出最近的 k 辆车，就像导航告诉你“最近的 5 个加油站”；

• 半径邻居调度员（RadiusNeighbors）负责寻找一定距离范围内的所有车辆，类似“查找 1 公里内的所有餐馆”；

• 带标签的调度员（*ClassMode 版本）在基础邻居搜索之上，还要根据邻居车辆的类型（警车、救护车）进行加权投票，从而决定目标车辆的属性；

• 行范数预计算器（sqeuclidean_row_norms）提前算好每辆车到原点的距离平方，利用

  [

  |x-y|^2 = |x|^2 + |y|^2 - 2,x\cdot y

  ]

  避免在成对距离计算中重复求平方和，犹如提前测好每个停车场到市中心的距离；

• 测试框架辅助函数 像交通监控系统中的校验传感器，确保每辆车的行驶数据（距离、路线、时间）准确无误，防止因感知偏差导致调度失误；

• 断言方法（assert_compatible_* 等） 像自动巡检机制，对比不同监测点的数据，验证调度结果在不同条件下的一致性（例如不同车流密度、不同车道、不同车辆类型）；

• 边界条件测试 如同故障注入测试，故意提供错误的输入（非法数据类型、错误维度），检查系统是否能正确拒绝或报错，防止危险调度。

整个调度中心需要根据车辆类型、道路状况和实时交通动态智能选择最优的底层实现路径；测试框架则像质量监控与故障诊断系统，通过丰富的测试用例和断言方法，确保在各种极端和边界情况下，调度逻辑仍保持正确、稳定和高效。

35.3 源码地图

sklearn/metrics/_pairwise_distances_reduction/_dispatcher.py
├── sqeuclidean_row_norms(X, num_threads)
│   ├── 根据 dtype 选择 _float32/_float64 实现
│   └── 调用 Cython 底层函数并包装为 NumPy 数组
├── BaseDistancesReductionDispatcher（抽象基类）
│   ├── valid_metrics() -> List[str]
│   │   ├── 排除不支持的度量（pyfunc、mahalanobis 等）
│   │   └── 返回支持的度量列表（sqeuclidean + METRIC_MAPPING64）
│   ├── is_usable_for(X, Y, metric) -> bool
│   │   ├── 检查 enable_cython_pairwise_dist 配置开关
│   │   ├── 验证 X、Y 为 C‑contiguous 或有效 CSR 稀疏矩阵
│   │   ├── 确保 X、Y 同 dtype 且为 float32/float64
│   │   └── 检查 metric 是否在支持列表中或为 DistanceMetric 实例
│   └── compute(X, Y, **kwargs) [抽象方法]
├── ArgKmin（k 近邻分派器）
│   ├── compute(X, Y, k, metric, ...) -> indices 或 (distances, indices)
│   │   ├── 根据 dtype 分派到 ArgKmin32/ArgKmin64
│   │   ├── 支持 chunk_size、strategy、return_distance 等参数
│   │   └── 使用 RAII 模式确保临时资源自动释放
│   └── 继承自 BaseDistancesReductionDispatcher
├── RadiusNeighbors（半径邻居分派器）
│   ├── compute(X, Y, radius, metric, ...) -> indices 或 (distances, indices)
│   │   ├── 根据 dtype 分派到 RadiusNeighbors32/RadiusNeighbors64
│   │   ├── 支持 sort_results 控制是否按距离排序
│   │   └── 同样采用 RAII 风格的资源管理
│   └── 继承自 BaseDistancesReductionDispatcher
├── ArgKminClassMode（带标签的 k 近邻分派器）
│   ├── compute(X, Y, k, weights, Y_labels, unique_Y_labels, ...) -> 概率数组
│   │   ├── 校验 weights 必须为 'uniform' 或 'distance'
│   │   ├── 根据 dtype 分派到 32/64 位实现
│   │   ├── 将 Y_labels、unique_Y_labels 转换为 np.intp 类型
│   │   └── 在最近邻基础上计算加权众数（weighted mode）
│   ├── valid_metrics() -> 排除 euclidean/sqeuclidean（当前无竞争力实现）
│   └── 继承自 BaseDistancesReductionDispatcher
├── RadiusNeighborsClassMode（带标签的半径邻居分派器）
│   ├── compute(X, Y, radius, weights, Y_labels, unique_Y_labels, outlier_label, ...) -> 概率数组
│   │   ├── 同样校验 weights 参数合法性
│   │   ├── 处理 outlier_label（无邻居样本的标签）
│   │   └── 其余逻辑与 ArgKminClassMode 类似，但基于半径范围
│   ├── valid_metrics() -> 同样排除 euclidean/sqeuclidean
│   └── 继承自 BaseDistancesReductionDispatcher
└── sklearn/metrics/_pairwise_distances_reduction/__init__.py
    ├── 导出所有公共接口：ArgKmin、ArgKminClassMode 等
    └── 包含高层架构图注释（解释派发机制与特化实现关系）

35.4 分派器架构与基类设计 —— 智能调度的“控制塔”

35.4.1 核心概念解析

BaseDistancesReductionDispatcher 为所有分派器提供统一的契约。它通过 白名单（valid_metrics()）过滤不适合 Cython 加速的度量，并利用 运行时守门（is_usable_for()）从 五个维度 判定是否可以走高速路径：

1. 全局配置：enable_cython_pairwise_dist 是否打开

2. 内存布局：X、Y 必须是 C‑contiguous 或合法的 CSR（format == "csr"、nnz > 0、int32 索引）

3. 数据类型：X.dtype == Y.dtype 且仅限 float32 / float64

4. 稀疏矩阵兼容性：稀疏‑稀疏 欧氏距离被临时禁用，以免性能倒退

5. 度量合法性：metric 必须在白名单中或是 DistanceMetric 实例

compute() 为抽象类方法，子类在 类方法（而非实例方法）中实现具体归约。类方法的好处是无状态，配合 RAII（资源获取即初始化）能够在调用结束时自动释放临时 C 级数组，避免内存泄漏。

35.4.1.1 代码块（逐行注释）

class BaseDistancesReductionDispatcher:
    """抽象基类，规定所有分派器必须实现的接口。"""

    @classmethod
    def valid_metrics(cls) -> List[str]:
        # ① 排除不适合无 GIL 加速的度量
        excluded = {
            "pyfunc",          # 需要回调 Python，无法在 Cython 中无 GIL 运行
            "mahalanobis",    # 数值不稳定
            "hamming",        # 需要稳定的 simultaneous sort
            *BOOL_METRICS,   # 布尔距离同上
        }
        # ② 只保留 sqeuclidean 与 METRIC_MAPPING64 中的度量
        return sorted(({"sqeuclidean"} | set(METRIC_MAPPING64.keys())) - excluded)

    @classmethod
    def is_usable_for(cls, X, Y, metric) -> bool:
        """判断当前 X、Y、metric 是否满足高速路径的使用条件。"""

        # ③ 稀疏‑稀疏欧氏距离暂时回退到 SciPy 实现
        if (issparse(X) and issparse(Y) and isinstance(metric, str) and
                "euclidean" in metric):
            return False

        # ④ 检查是否为 C‑contiguous ndarray
        def is_numpy_c_ordered(arr):
            return hasattr(arr, "flags") and getattr(arr.flags, "c_contiguous", False)

        # ⑤ 检查是否为合法 CSR 稀疏矩阵（int32 索引且 nnz > 0）
        def is_valid_sparse_matrix(arr):
            return (issparse(arr) and arr.format == "csr" and
                    arr.nnz > 0 and arr.indices.dtype == arr.indptr.dtype == np.int32)

        # ⑥ 五大条件同时满足才返回 True
        is_usable = (
            get_config().get("enable_cython_pairwise_dist", True)               # 配置开关
            and (is_numpy_c_ordered(X) or is_valid_sparse_matrix(X))            # X 的布局
            and (is_numpy_c_ordered(Y) or is_valid_sparse_matrix(Y))            # Y 的布局
            and X.dtype == Y.dtype                                               # dtype 必须相同
            and X.dtype in (np.float32, np.float64)                             # 仅支持 float32/64
            and (metric in cls.valid_metrics() or isinstance(metric, DistanceMetric))
        )
        return is_usable

    @classmethod
    @abstractmethod
    def compute(cls, X, Y, **kwargs):
        """抽象方法，子类必须实现具体的归约逻辑。"""

解释：上述代码展示了基类如何通过 五重守门 决定是否可以使用 Cython 实现。valid_metrics() 通过集合操作明确列出可加速的度量，is_usable_for() 则将配置、内存布局、dtype、稀疏兼容性和度量合法性全部统一检查，返回布尔值供上层调度器使用。compute() 被标记为抽象类方法，强制子类提供实现。

35.4.1.2 架构图（分派器整体流程）

graph TD A[用户调用任意 Dispatcher.compute] --> B{is_usable_for?} B -->|True| C[进入 Cython 高速路径] B -->|False| D[回退至 SciPy / NumPy 实现] C --> E[根据 dtype 分派到 32/64 位实现] E --> F[内部 Cython 循环 + OpenMP 并行] F --> G[返回 NumPy 数组（indices / distances）] G --> H[Python 层自动完成 RAII 资源释放]

35.5 ArgKmin 与 RadiusNeighbors 分派器 —— k 近邻与半径邻居的“双引擎”

35.5.1 核心概念解析

ArgKmin 与 RadiusNeighbors 是分派器体系中最基础的两个“引擎”。它们 不直接实现距离计算，而是：

1. 检查 dtype 是否匹配（float64 → *64，float32 → *32）

2. 转发所有参数（chunk_size、strategy、return_distance/sort_results）到对应的 Cython 实现

3. 保持 RAII：在类方法返回前自动回收临时 C 级数据结构

两者均 不应实例化；用户只能通过 compute 类方法调用，这也是实现 无状态调度 的关键。

35.5.1.1 ArgKmin.compute（逐行注释）

class ArgKmin(BaseDistancesReductionDispatcher):
    """k‑最近邻分派器：对每个查询向量找到 Y 中距离最近的 k 个样本。"""

    @classmethod
    def compute(
        cls,
        X,
        Y,
        k,
        metric="euclidean",
        chunk_size=None,
        metric_kwargs=None,
        strategy=None,
        return_distance=False,
    ):
        """核心分派入口。"""
        # ① dtype 分派：float64 → ArgKmin64，float32 → ArgKmin32
        if X.dtype == Y.dtype == np.float64:
            return ArgKmin64.compute(
                X=X, Y=Y, k=k, metric=metric,
                chunk_size=chunk_size, metric_kwargs=metric_kwargs,
                strategy=strategy, return_distance=return_distance,
            )
        if X.dtype == Y.dtype == np.float32:
            return ArgKmin32.compute(
                X=X, Y=Y, k=k, metric=metric,
                chunk_size=chunk_size, metric_kwargs=metric_kwargs,
                strategy=strategy, return_distance=return_distance,
            )
        # ② dtype 不匹配或不受支持时抛出清晰错误
        raise ValueError(
            "Only float64 or float32 datasets pairs are supported at this time, "
            f"got: X.dtype={X.dtype} and Y.dtype={Y.dtype}."
        )

解释：compute 首先检查 X 与 Y 的 dtype 是否一致且为受支持类型。若匹配，则调用对应的 32/64 位实现；否则抛出带有具体 dtype 信息的 ValueError，帮助用户快速定位错误。

35.5.1.2 RadiusNeighbors.compute（逐行注释）

class RadiusNeighbors(BaseDistancesReductionDispatcher):
    """半径邻居分派器：返回所有距离 ≤ radius 的邻居。"""

    @classmethod
    def compute(
        cls,
        X,
        Y,
        radius,
        metric="euclidean",
        chunk_size=None,
        metric_kwargs=None,
        strategy=None,
        return_distance=False,
        sort_results=False,
    ):
        """核心分派入口。"""
        # 与 ArgKmin 完全相同的 dtype 分派逻辑
        if X.dtype == Y.dtype == np.float64:
            return RadiusNeighbors64.compute(
                X=X, Y=Y, radius=radius, metric=metric,
                chunk_size=chunk_size, metric_kwargs=metric_kwargs,
                strategy=strategy, sort_results=sort_results,
                return_distance=return_distance,
            )
        if X.dtype == Y.dtype == np.float32:
            return RadiusNeighbors32.compute(
                X=X, Y=Y, radius=radius, metric=metric,
                chunk_size=chunk_size, metric_kwargs=metric_kwargs,
                strategy=strategy, sort_results=sort_results,
                return_distance=return_distance,
            )
        # dtype 不匹配时抛出错误
        raise ValueError(
            "Only float64 or float32 datasets pairs are supported at this time, "
            f"got: X.dtype={X.dtype} and Y.dtype={Y.dtype}."
        )

解释：逻辑与 ArgKmin.compute 完全对称，只是将 k 替换为 radius，并新增 sort_results 参数用于控制是否对每个查询的邻居按距离升序排列。

35.5.1.3 参数说明（对比）

| 参数 | 作用 | 对并行策略的影响 |

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

| chunk_size | 每次处理的向量块大小；若 None 会使用全局配置 pairwise_dist_chunk_size（默认 256） | 较大的块有助于降低线程调度开销，过小会导致大量线程切换 |

| strategy | 'parallel_on_X'、'parallel_on_Y'、'auto' 任选 | parallel_on_X 适合 X 行数大、Y 行数小；parallel_on_Y 适合 Y 行数大；auto 根据规模自动选择 |

| return_distance / sort_results | 是否返回距离数组以及是否对结果排序 | 返回距离会额外分配距离矩阵；sort_results=True 需要在每个块内部进行额外排序，略增开销 |

35.5.1.4 流程图（ArgKmin / RadiusNeighbors 分派）

graph TD A[用户调用 ArgKmin.compute / RadiusNeighbors.compute] --> B{X.dtype == Y.dtype?} B -->|float64| C[调用 ArgKmin64 / RadiusNeighbors64.compute] B -->|float32| D[调用 ArgKmin32 / RadiusNeighbors32.compute] B -->|其他| E[抛出 ValueError] C --> F[Cython 并行循环 (prange) + 堆/半径过滤] D --> F F --> G[返回 indices (, distances)] G --> H[Python 层接收 NumPy 数组，RAII 自动释放 C 级资源]

35.6 类模式归约与行范数优化 —— 加权投票与特征预处理的“辅助系统”

35.6.1 核心概念解析

在分类任务中，仅返回邻居索引不足以做出预测，需要结合 标签 与 加权方式 计算每个查询点的类别概率。ArgKminClassMode 与 RadiusNeighborsClassMode 正是为此设计：它们在基础邻居搜索之上：

1. 校验 weights：只能是 'uniform'（等权）或 'distance'（距离倒数加权）

2. 将标签数组强制转换为 np.intp（平台指针宽度），方便 Cython 中直接作为整数索引

3. 调用对应的 32/64 位实现 完成 加权众数（weighted mode）计算，返回形状为 (n_samples_X, n_classes) 的概率矩阵

同时，sqeuclidean_row_norms 为欧氏距离分解提供预计算的行范数，配合矩阵乘法（BLAS GEMM）实现 O(n·m) 的高效距离计算，而不是传统的 O(n·m·d) 双层循环。

35.6.1.1 ArgKminClassMode.compute（逐行注释）

class ArgKminClassMode(BaseDistancesReductionDispatcher):
    """带标签的 k‑最近邻分派器：在返回最近邻的同时计算加权类别概率。"""

    @classmethod
    def valid_metrics(cls) -> List[str]:
        # 暂不支持欧氏距离，因为缺少 GEMM 特化实现
        excluded = {"euclidean", "sqeuclidean"}
        return list(set(BaseDistancesReductionDispatcher.valid_metrics()) - excluded)

    @classmethod
    def compute(
        cls,
        X,
        Y,
        k,
        weights,
        Y_labels,
        unique_Y_labels,
        metric="euclidean",
        chunk_size=None,
        metric_kwargs=None,
        strategy=None,
    ):
        """核心入口，完成 dtype 分派并进行加权投票。"""
        # ① 校验权重策略是否合法
        if weights not in {"uniform", "distance"}:
            raise ValueError(
                "Only the 'uniform' or 'distance' weights options are supported"
                f" at this time. Got: {weights=}."
            )
        # ② dtype 分派 + 标签转为平台指针宽度整数（np.intp）
        if X.dtype == Y.dtype == np.float64:
            return ArgKminClassMode64.compute(
                X=X, Y=Y, k=k, weights=weights,
                Y_labels=np.array(Y_labels, dtype=np.intp),
                unique_Y_labels=np.array(unique_Y_labels, dtype=np.intp),
                metric=metric, chunk_size=chunk_size,
                metric_kwargs=metric_kwargs, strategy=strategy,
            )
        if X.dtype == Y.dtype == np.float32:
            return ArgKminClassMode32.compute(
                X=X, Y=Y, k=k, weights=weights,
                Y_labels=np.array(Y_labels, dtype=np.intp),
                unique_Y_labels=np.array(unique_Y_labels, dtype=np.intp),
                metric=metric, chunk_size=chunk_size,
                metric_kwargs=metric_kwargs, strategy=strategy,
            )
        # ③ 不匹配时报错，包含具体 dtype 信息
        raise ValueError(
            "Only float64 or float32 datasets pairs are supported at this time, "
            f"got: X.dtype={X.dtype} and Y.dtype={Y.dtype}."
        )

解释：在确认 weights 合法后，标签数组被强制转为 np.intp（与 C 语言指针大小一致），这样 Cython 可以直接以指针算术访问。随后依据 dtype 分派到 32/64 位实现，返回 类别概率矩阵。

35.6.1.2 RadiusNeighborsClassMode.compute（逐行注释）

class RadiusNeighborsClassMode(BaseDistancesReductionDispatcher):
    """带标签的半径邻居分派器：在给定半径内进行加权投票。"""

    @classmethod
    def valid_metrics(cls) -> List[str]:
        # 同上，暂不支持欧氏距离
        excluded = {"euclidean", "sqeuclidean"}
        return sorted(set(BaseDistancesReductionDispatcher.valid_metrics()) - excluded)

    @classmethod
    def compute(
        cls,
        X,
        Y,
        radius,
        weights,
        Y_labels,
        unique_Y_labels,
        outlier_label,
        metric="euclidean",
        chunk_size=None,
        metric_kwargs=None,
        strategy=None,
    ):
        """核心入口，完成 dtype 分派并处理 outlier 情形。"""
        # ① 权重合法性检查
        if weights not in {"uniform", "distance"}:
            raise ValueError(
                "Only the 'uniform' or 'distance' weights options are supported"
                f" at this time. Got: {weights=}."
            )
        # ② dtype 分派 + 标签转为 np.intp，透传 outlier_label
        if X.dtype == Y.dtype == np.float64:
            return RadiusNeighborsClassMode64.compute(
                X=X, Y=Y, radius=radius, weights=weights,
                Y_labels=np.array(Y_labels, dtype=np.intp),
                unique_Y_labels=np.array(unique_Y_labels, dtype=np.intp),
                outlier_label=outlier_label,
                metric=metric, chunk_size=chunk_size,
                metric_kwargs=metric_kwargs, strategy=strategy,
            )
        if X.dtype == Y.dtype == np.float32:
            return RadiusNeighborsClassMode32.compute(
                X=X, Y=Y, radius=radius, weights=weights,
                Y_labels=np.array(Y_labels, dtype=np.intp),
                unique_Y_labels=np.array(unique_Y_labels, dtype=np.intp),
                outlier_label=outlier_label,
                metric=metric, chunk_size=chunk_size,
                metric_kwargs=metric_kwargs, strategy=strategy,
            )
        # ③ dtype 不匹配错误
        raise ValueError(
            "Only float64 or float32 datasets pairs are supported at this time, "
            f"got: X.dtype={X.dtype} and Y.dtype={Y.dtype}."
        )