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

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

X=[[0.5]] 无法高置信预测，首轮无新增样本，触发 no_change 终止。

---

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_zero_iterations()（第129-148行）

@pytest.mark.parametrize(
    "estimator",
    [
        KNeighborsClassifier(),
        CalibratedClassifierCV(SVC(gamma="scale", random_state=0), ensemble=False),
    ],
)
@pytest.mark.parametrize("y", [y_train_missing_labels, y_train_missing_strings])
def test_zero_iterations(estimator, y):
    estimator = clone(estimator)
    clf1 = SelfTrainingClassifier(estimator, max_iter=0)
    clf1.fit(X_train, y)
    clf2 = estimator.fit(X_train[:n_labeled_samples], y[:n_labeled_samples])
    assert_array_equal(clf1.predict(X_test), clf2.predict(X_test))
    assert clf1.termination_condition_ == "max_iter"

max_iter=0 直接跳过循环，仅在初始标注集上拟合，等价纯监督。

---

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_no_unlabeled()（第178-192行）

def test_no_unlabeled():
    knn = KNeighborsClassifier()
    knn.fit(X_train, y_train)
    st = SelfTrainingClassifier(knn)
    with pytest.warns(UserWarning, match="y contains no unlabeled samples"):
        st.fit(X_train, y_train)
    assert_array_equal(knn.predict(X_test), st.predict(X_test))
    assert np.all(st.labeled_iter_ == 0)
    assert st.termination_condition_ == "all_labeled"

全标注数据触发警告，但结果等价纯监督，所有样本 labeled_iter_=0。

---

17.14.6 异常处理与元估计器兼容性

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_warns_k_best()（第36-42行）

def test_warns_k_best():
    st = SelfTrainingClassifier(KNeighborsClassifier(), criterion="k_best", k_best=1000)
    with pytest.warns(UserWarning, match="k_best is larger than"):
        st.fit(X_train, y_train_missing_labels)
    assert st.termination_condition_ == "all_labeled"

k_best 超过未标注样本数时警告，首轮全标注，终止条件为 all_labeled。

---

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_prefitted_throws_error()（第151-160行）

def test_prefitted_throws_error():
    knn = KNeighborsClassifier()
    knn.fit(X_train, y_train)
    st = SelfTrainingClassifier(knn)
    with pytest.raises(
        NotFittedError,
        match="This SelfTrainingClassifier instance is not fitted yet",
    ):
        st.predict(X_train)

传入已拟合估计器，clone 会复制其拟合状态，但 SelfTrainingClassifier 自身未拟合，调用 predict 抛 NotFittedError。

---

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_strings_dtype()（第206-214行）

def test_strings_dtype():
    clf = SelfTrainingClassifier(KNeighborsClassifier())
    X, y = make_blobs(n_samples=30, random_state=0, cluster_std=0.1)
    labels_multiclass = ["one", "two", "three"]
    y_strings = np.take(labels_multiclass, y)
    with pytest.raises(ValueError, match="dtype"):
        clf.fit(X, y_strings)

字符串 dtype（kind in ["U", "S"]）被显式拒绝，引导用户使用 object dtype + -1 表示未标注。

---

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_verbose() 与 test_verbose_k_best()（第217-253行）

验证 verbose=True 时打印迭代信息，verbose=False 静默。

---

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_estimator_meta_estimator()（第286-315行）

def test_estimator_meta_estimator():
    estimator = StackingClassifier(
        estimators=[("clf_1", LogisticRegression()), ("clf_2", LogisticRegression())],
        final_estimator=LogisticRegression(),
        cv=2,
    )
    assert hasattr(estimator, "predict_proba")
    clf = SelfTrainingClassifier(estimator=estimator)
    clf.fit(X_train, y_train_missing_labels)
    clf.predict_proba(X_test)

    estimator = StackingClassifier(
        estimators=[("svc_1", SVC()), ("svc_2", SVC())],
        final_estimator=SVC(),
        cv=2,
    )
    assert not hasattr(estimator, "predict_proba")
    clf = SelfTrainingClassifier(estimator=estimator)
    with pytest.raises(AttributeError):
        clf.fit(X_train, y_train_missing_labels)

验证“拟合后才暴露 predict_proba”的元估计器：StackingClassifier 含 LogisticRegression 时有 predict_proba，自训练能工作；含 SVC(probability=False) 时无 predict_proba，fit 时抛 AttributeError。

---

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_self_training_estimator_attribute_error()（第318-355行）

def test_self_training_estimator_attribute_error():
    estimator = SVC(gamma="scale")
    self_training = SelfTrainingClassifier(estimator)
    with pytest.raises(AttributeError, match="has no attribute 'predict_proba'"):
        self_training.fit(X_train, y_train_missing_labels)

    self_training = SelfTrainingClassifier(estimator=DecisionTreeClassifier())
    outer_msg = "This 'SelfTrainingClassifier' has no attribute 'decision_function'"
    inner_msg = "'DecisionTreeClassifier' object has no attribute 'decision_function'"
    with pytest.raises(AttributeError, match=outer_msg) as exec_info:
        self_training.fit(X_train, y_train_missing_labels).decision_function(X_train)
    assert isinstance(exec_info.value.__cause__, AttributeError)
    assert inner_msg in str(exec_info.value.__cause__)

两种错误场景：

1. fit 内部调用 predict_proba，底层无此方法 → 直接 AttributeError

2. decision_function 由 available_if 条件暴露，底层无此方法 → 包装后的 AttributeError，__cause__ 链保留原始错误信息

---

17.14.7 元数据路由测试

源码路径：sklearn/semi_supervised/tests/test_self_training.py - test_routing_passed_metadata_not_supported()（第358-387行）

@pytest.mark.filterwarnings("ignore:y contains no unlabeled samples:UserWarning")
@pytest.mark.parametrize(
    "method", ["decision_function", "predict_log_proba", "predict_proba", "predict"]
)
def test_routing_passed_metadata_not_supported(method):
    """Test that the right error message is raised when metadata is passed while
    not supported when `enable_metadata_routing=False`."""
    est = SelfTrainingClassifier(estimator=SimpleEstimator())
    with pytest.raises(
        ValueError, match="is only supported if enable_metadata_routing=True"
    ):
        est.fit([[1], [1]], [1, 1], sample_weight=[1], prop="a")

    est = SelfTrainingClassifier(estimator=SimpleEstimator())
    with pytest.raises(
        ValueError, match="is only supported if enable_metadata_routing=True"
    ):
        # make sure that the estimator thinks it is already fitted
        est.fitted_params_ = True
        getattr(est, method)([[1]], sample_weight=[1], prop="a")

验证未启用全局路由配置时，向 fit/predict/... 传递额外参数（如 sample_weight）抛出明确错误。

---

17.15 设计中的取舍

17.15.1 为什么 LabelPropagation 和 LabelSpreading 不共用一个类，而是通过 _variant 区分？

硬钳制与软钳制在数学本质上不同：前者是约束优化问题（固定已标注节点），后者是正则化优化问题（软约束）。合并为单一类会导致参数空间混淆（alpha 对 LabelPropagation 无意义）、默认 max_iter 差异大（1000 vs 30）、收敛性质不同。通过继承共享 90% 通用逻辑（核调度、迭代骨架、预测接口），仅在 _build_graph 和钳制分支差异化，符合“组合优于继承”的变体设计原则。

---

17.15.2 为什么 SelfTrainingClassifier 要求 estimator 仅有 fit 而非 fit + predict_proba？

为了支持“拟合后才暴露 predict_proba”的元估计器（如 StackingClassifier）。如果在 __init__ 校验阶段就要求 predict_proba，这类元估计器会被错误拒绝。scikit-learn 选择在运行时（fit 内部调用 predict_proba）才检查，配合 available_if 实现动态接口暴露，体现了“鸭子类型”与“后期绑定”的 Python 动态哲学。

---

17.15.3 为什么 LabelSpreading 的闭式解用 (I - alpha * S)^(-1) Y 而 LabelPropagation 用分块矩阵求解？

LabelSpreading 的软钳制迭代 \(F^{(t)} = \alpha S F^{(t-1)} + (1-\alpha)Y\) 是线性定常系统，直接求解线性方程组得闭式解。LabelPropagation 的硬钳制将已标注节点固定，相当于在方程组中移除已标注行列，只对未标注子系统求解，因此需要分块矩阵 \(T_{uu}, T_{ul}\)。两种数学推导路径不同，但最终都验证了迭代算法的正确性。

---

17.16 动手练习

17.16.1 练习 1：对比硬钳制与软钳制的迭代传播差异

1. 阅读 sklearn/semi_supervised/_label_propagation.py 第 198-293 行的 fit() 方法。

2. 回答问题：

   • 硬钳制分支中 np.where(unlabeled, self.label_distributions_, y_static) 完成了什么操作？

   • 软钳制分支中 alpha * label_distributions_ + y_static 的数学含义是什么？

   • 如果 alpha 趋近于 0，LabelSpreading 的行为会退化成什么？

17.16.2 练习 2：推导闭式解并验证迭代一致性

1. 阅读 test_label_propagation.py 第 63-112 行的两个闭式解测试。

2. 回答问题：

   • LabelSpreading 的闭式解为什么是 (I - alpha * S)⁻¹ Y？

   • LabelPropagation 的闭式解中，Tuu 和 Tul 子矩阵分别代表什么？

   • 为什么两种闭式解都要在最终做行归一化？

17.16.3 练习 3：追踪自训练的 k_best 选择逻辑

1. 阅读 sklearn/semi_supervised/_self_training.py 第 220-270 行。

2. 回答问题：

   • np.argpartition(-max_proba, n_to_select)[:n_to_select] 返回的是索引还是布尔掩码？

   • selected_full = np.nonzero(~has_label)[0][selected] 完成了什么映射？

   • 当 n_to_select == max_proba.shape[0] 时，为什么退化为 np.ones_like(max_proba, dtype=bool)？

---

17.17 本章小结

这一章中我们学习了 scikit-learn 半监督学习模块的核心实现。首先我们了解了模块入口的最小暴露设计，其次深入剖析了 BaseLabelPropagation 基类如何通过模板方法模式统一核函数调度、标签分布初始化、迭代传播主循环与收敛判据，接着对比了 LabelPropagation 的硬钳制（行归一化概率转移矩阵）与 LabelSpreading 的软钳制（归一化图拉普拉斯 + alpha 加权融合）两种图构建与传播策略的数学本质与工程实现差异，然后通过测试代码中的闭式解验证理解了迭代算法的数学正确性保障，最后详细解析了 SelfTrainingClassifier 作为元估计器的伪标注迭代流水线、条件方法暴露机制与元数据路由实现。

本章核心概念速查表：

| 概念 | 解释 |

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

| BaseLabelPropagation._get_kernel() | 核函数调度中心：RBF 生成 O(N²) 稠密亲和矩阵，KNN 生成 O(kN) 稀疏图，callable 直接调用用户函数 |

| LabelPropagation._build_graph() | 行归一化构造概率转移矩阵 P = D⁻¹W，稀疏/稠密分支分别使用 diags 缩放与广播除法 |

| LabelSpreading._build_graph() | 归一化图拉普拉斯变换 S = -L_norm = D⁻¹ᐟ² W D⁻¹ᐟ²，对角线置零防止自传播 |

| BaseLabelPropagation.fit() | 迭代传播主循环：硬钳制强制重置标注样本，软钳制 alpha 加权融合，收敛判据为分布变化和 < tol |

| BaseLabelPropagation.predict_proba() | 跨图传播生成新样本概率，KNN 分支取近邻分布直接求和，稠密分支 safe_sparse_dot 后行归一化 |

| SelfTrainingClassifier.fit() | 迭代伪标注：每轮拟合已标注子集，对未标注子集 predict_proba 后按 threshold 或 k_best 选择高置信样本 |

| termination_condition_ | 三种终止原因：max_iter 耗尽、no_change 无新增样本、all_labeled 全量标注完成 |

| labeled_iter_ | 记录每个样本获得标签的迭代轮次：初始标注样本为 0，从未标注为 -1 |

| available_if | 条件方法暴露装饰器，按底层估计器能力动态挂载 predict/decision_function/score 等方法 |

| get_metadata_routing() | MetadataRouter + MethodMapping 声明 7 组 callee→caller 映射，支持 sample_weight 等元数据路由 |

| 闭式解验证 | label_propagation_closed_form 用 Tuu/Tul 子矩阵求解，label_spreading_closed_form 用 (I-αS)⁻¹Y 验证迭代结果 |

| test_self_training.py __main__ | 模块级测试数据准备：train_test_split 切分 iris，50 个标注样本 + 字符串标签映射，供后续所有测试复用 |

下一章中，我们将学习核近似与核岭回归，探索随机傅里叶特征、Nyström 方法如何将非线性核映射到显式特征空间，兼顾表达力与计算效率。

17.18 架构与数据流图

graph TD A[__init__] --> B[_label_propagation] B --> C[_label_propagation]

sequenceDiagram participant U as 调用者 participant E as __init__ participant C as _label_propagation U->>E: 调用入口 E->>C: 传递参数 C-->>U: 返回结果

graph LR I[输入] --> P[参数校验] P --> T[核心处理] T --> O[输出]

graph TD L1[用户 API 层] --> L2[算法/服务层] L2 --> L3[数据结构层] L3 --> L4[运行时与依赖层]

上述图分别展示模块依赖、调用时序、数据流和架构分层。

第 18 章 —— 核近似与核岭回归 —— 延展“线性模型的非线性边界”

18.1 学习目标

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

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

• 理解核近似方法（随机傅里叶特征、Nyström、张量草图、显式特征映射）的数学原理与适用场景

• 掌握 RBFSampler、SkewedChi2Sampler、AdditiveChi2Sampler、PolynomialCountSketch、Nystroem 等近似器的源码实现细节

• 了解核岭回归 的对偶求解机制、预计算核支持、Cholesky 分解与奇异矩阵回退逻辑

• 熟悉稀疏/稠密矩阵统一处理、数据类型传播、Array API 多后端兼容等工程化设计模式

• 能阅读并编写核近似与核岭回归的单元测试，验证数学正确性、数值稳定性与跨后端一致性

18.2 生活类比

想象核近似是一家高维特征压缩工厂：原始数据 = 复杂的高维设计图纸（难以直接加工）；核函数 = 设计规范（定义图纸间的相似度）；随机傅里叶特征 = 随机投影仪（用随机光栅将图纸投影到低维屏幕）；Nyström 方法 = 关键图纸抽样（挑选代表性图纸作为基准，其他图纸只需与基准比对）；张量草图 = 哈希压缩机（将高阶特征交互通过哈希桶压缩，再用 FFT 快速重组）；显式特征映射 = 标准化展开图（按数学公式将每个特征展开为正弦/余弦/平方根分量）；核岭回归 = 对偶空间求解器（在样本相似度空间而非特征空间求解，支持预计算相似度矩阵）。就像工程师用投影、抽样、哈希等手段将复杂设计简化为可加工的标准件，核近似将隐式高维核映射显式化，兼顾表达力与计算效率。

18.3 源码地图

sklearn/kernel_approximation.py

├── PolynomialCountSketch

│ ├── init # 参数初始化

│ ├── fit # 随机哈希初始化

│ ├── transform # FFT 加速张量草图特征映射

│ └── sklearn_tags # 支持稀疏输入

├── RBFSampler

│ ├── init # 参数初始化

│ ├── fit # 自适应 gamma 与随机权重/偏移采样

│ ├── transform # 就地操作链：点积→加偏移→cos→归一化

│ └── sklearn_tags # 稀疏输入、dtype 保持

├── SkewedChi2Sampler

│ ├── init # 参数初始化

│ ├── fit # sech 分布采样权重、均匀采样偏移

│ ├── transform # 对数域变换、点积、余弦映射

│ └── sklearn_tags # 正值输入、dtype 保持

├── AdditiveChi2Sampler

│ ├── init # 参数初始化

│ ├── fit # 仅参数校验

│ ├── transform # 稠密/稀疏分派到 _transform_dense/_transform_sparse

│ ├── _transform_dense # 稠密矩阵显式周期采样映射

│ ├── _transform_sparse # 稀疏矩阵显式周期采样映射

│ ├── get_feature_names_out # 语义化特征命名

│ └── sklearn_tags # 无需拟合、正值输入、稀疏支持

├── Nystroem

│ ├── init # 参数初始化

│ ├── fit # 子采样、核矩阵 SVD、归一化矩阵构造

│ ├── transform # 新样本嵌入

│ ├── _get_kernel_params # 灵活内核参数管理

│ └── sklearn_tags # Array API 支持、稀疏输入、dtype 保持

sklearn/kernel_ridge.py

├── KernelRidge

│ ├── init # 参数初始化

│ ├── _get_kernel # 统一内核计算入口

│ ├── fit # Cholesky 求解对偶系数、奇异回退、样本权重、多输出

│ ├── predict # 核矩阵乘法预测

│ └── sklearn_tags # 稀疏输入、预计算核 pairwise 标记

sklearn/tests/test_kernel_approximation.py

├── test_polynomial_count_sketch # 无偏性与精度验证

├── test_polynomial_count_sketch_dense_sparse # 稠密/稀疏一致性

├── test_additive_chi2_sampler # 精确核匹配与稀疏等价

├── test_additive_chi2_sampler_sample_steps # 参数组合校验

├── test_additive_chi2_sampler_wrong_sample_steps # 异常处理

├── test_skewed_chi2_sampler # 对数域验证与边界检查

├── test_additive_chi2_sampler_exceptions # 负值输入异常

├── test_rbf_sampler # RBF 近似精度

├── test_rbf_sampler_fitted_attributes_dtype # dtype 传播

├── test_rbf_sampler_dtype_equivalence # float32/64 等价性

├── test_rbf_sampler_gamma_scale # 自适应 gamma 校验

├── test_skewed_chi2_sampler_fitted_attributes_dtype # dtype 传播

├── test_skewed_chi2_sampler_dtype_equivalence # float32/64 等价性

├── test_input_validation # 列表/稀疏输入兼容性

├── test_nystroem_approximation # 基础功能与可调用内核

├── test_nystroem_approximation_array_api # 多后端一致性

├── test_nystroem_default_parameters # 默认参数行为

├── test_nystroem_singular_kernel # 奇异核矩阵处理

├── test_nystroem_poly_kernel_params # 多项式核参数透传

├── test_nystroem_callable # 可调用内核与参数冲突检查

├── test_nystroem_precomputed_kernel # 预计算核路径

├── test_nystroem_component_indices # 基向量索引记录

├── test_get_feature_names_out # 通用特征名规范

├── test_additivechi2sampler_get_feature_names_out # 加性卡方特征名语义

├── _linear_kernel # 线性核辅助函数

└── main # 全局代码

sklearn/tests/test_kernel_ridge.py

├── test_kernel_ridge # 线性核下与 Ridge 系数对齐

├── test_kernel_ridge_sparse # 稀疏输入一致性

├── test_kernel_ridge_singular_kernel # 奇异核回退 lstsq

├── test_kernel_ridge_precomputed # 预计算核等价性

├── test_kernel_ridge_precomputed_kernel_unchanged # 预计算核不修改输入

├── test_kernel_ridge_sample_weights # 样本权重原始/预计算等价

├── test_kernel_ridge_multi_output # 多目标回归等价性

└── main # 全局代码

18.4 PolynomialCountSketch —— 利用 FFT 加速的张量草图多项式核近似

18.4.1 核心概念与类型定义

PolynomialCountSketch 实现了张量草图，它通过 Count Sketch 哈希将高维多项式特征映射压缩到低维空间，并利用快速傅里叶变换（FFT）高效实现多项式核的特征映射近似。其目标核函数为 \(K(X, Y) = (\gamma \langle X, Y \rangle + c_0)^d\)。我们先看构造函数签名：

源码路径：sklearn/kernel_approximation.py - PolynomialCountSketch.__init__（第39-55行）

def __init__(
    self, *, gamma=1.0, degree=2, coef0=0, n_components=100, random_state=None
):
    self.gamma = gamma
    self.degree = degree
    self.coef0 = coef0
    self.n_components = n_components
    self.random_state = random_state

参数图纸解读：

• gamma：核函数缩放因子，控制内积权重。

• degree：多项式次数，决定特征交互阶数（如 degree=2 即二阶交互）。

• coef0：常数项，引入偏置，相当于在输入特征中追加一个常数维度。

• n_components：输出特征维度，通常需大于输入特征数以保证近似质量。

• random_state：控制随机哈希初始化的可复现性。

该类继承自 ClassNamePrefixFeaturesOutMixin、TransformerMixin 和 BaseEstimator，标准实现了 fit/transform 接口。

18.4.2 随机哈希初始化 (fit)

源码路径：sklearn/kernel_approximation.py - PolynomialCountSketch.fit（第87-118行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
    """Fit the model with X.

    Initializes the internal variables. The method needs no information
    about the distribution of data, so we only care about n_features in X.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Training data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    y : array-like of shape (n_samples,) or (n_samples, n_outputs), \
            default=None
        Target values (None for unsupervised transformations).

    Returns
    -------
    self : object
        Returns the instance itself.
    """
    X = validate_data(self, X, accept_sparse="csc")
    random_state = check_random_state(self.random_state)

    n_features = X.shape[1]
    if self.coef0 != 0:
        n_features += 1

    self.indexHash_ = random_state.randint(
        0, high=self.n_components, size=(self.degree, n_features)
    )

    self.bitHash_ = random_state.choice(a=[-1, 1], size=(self.degree, n_features))
    self._n_features_out = self.n_components
    return self

逐行解析：

1. validate_data(self, X, accept_sparse="csc")：验证输入数据，接受 CSC 稀疏矩阵格式（便于按列访问哈希索引）。

2. check_random_state(self.random_state)：获取确定性的随机数生成器。

3. n_features = X.shape[1]：获取输入特征数。

4. if self.coef0 != 0: n_features += 1：若有常数项，逻辑上视为额外增加一维特征（后续 transform 中会显式拼接）。

5. self.indexHash_ = random_state.randint(...)：生成形状为 (degree, n_features) 的哈希桶索引，取值范围 [0, n_components)。这是 Count Sketch 的核心：将每个原始特征在每个多项式次数下映射到某个输出维度。

6. self.bitHash_ = random_state.choice(a=[-1, 1], ...)：生成同形状的随机符号（+1/-1），构成 2-wise 独立哈希函数族，保证无偏估计。

7. self._n_features_out = self.n_components：设置输出维度属性。

这段代码仅初始化随机哈希参数，不依赖数据分布，体现了“无状态特征映射”的设计：哈希函数固定后，映射即为确定性线性变换。

18.4.3 稀疏/稠密统一的特征映射 (transform)

源码路径：sklearn/kernel_approximation.py - PolynomialCountSketch.transform（第120-170行）

def transform(self, X):
    """Generate the feature map approximation for X.

    Parameters
    ----------
    X : {array-like}, shape (n_samples, n_features)
        New data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    Returns
    -------
    X_new : array-like, shape (n_samples, n_components)
        Returns the instance itself.
    """

    check_is_fitted(self)
    X = validate_data(self, X, accept_sparse="csc", reset=False)

    X_gamma = np.sqrt(self.gamma) * X

    if sp.issparse(X_gamma) and self.coef0 != 0:
        X_gamma = sp.hstack(
            [X_gamma, np.sqrt(self.coef0) * np.ones((X_gamma.shape[0], 1))],
            format="csc",
        )

    elif not sp.issparse(X_gamma) and self.coef0 != 0:
        X_gamma = np.hstack(
            [X_gamma, np.sqrt(self.coef0) * np.ones((X_gamma.shape[0], 1))]
        )

    if X_gamma.shape[1] != self.indexHash_.shape[1]:
        raise ValueError(
            "Number of features of test samples does not"
            " match that of training samples."
        )

    count_sketches = np.zeros((X_gamma.shape[0], self.degree, self.n_components))

    if sp.issparse(X_gamma):
        for j in range(X_gamma.shape[1]):
            for d in range(self.degree):
                iHashIndex = self.indexHash_[d, j]
                iHashBit = self.bitHash_[d, j]
                count_sketches[:, d, iHashIndex] += (
                    (iHashBit * X_gamma[:, [j]]).toarray().ravel()
                )

    else:
        for j in range(X_gamma.shape[1]):
            for d in range(self.degree):
                iHashIndex = self.indexHash_[d, j]
                iHashBit = self.bitHash_[d, j]
                count_sketches[:, d, iHashIndex] += iHashBit * X_gamma[:, j]

    # For each same, compute a count sketch of phi(x) using the polynomial
    # multiplication (via FFT) of p count sketches of x.
    count_sketches_fft = fft(count_sketches, axis=2, overwrite_x=True)
    count_sketches_fft_prod = np.prod(count_sketches_fft, axis=1)
    data_sketch = np.real(ifft(count_sketches_fft_prod, overwrite_x=True))

    return data_sketch

逐行解析：

1. check_is_fitted(self)：确保已调用 fit。

2. X = validate_data(...)：验证输入，reset=False 禁止重置 n_features_in_。

3. X_gamma = np.sqrt(self.gamma) * X：输入缩放，将 \(\gamma\) 吸入特征向量，后续内积自动带入 \(\gamma\) 因子。

4. 稀疏/稠密分支拼接常数项：若 coef0 != 0，拼接一列 \(\sqrt{c_0}\)。稀疏用 sp.hstack 保持 CSC 格式；稠密用 np.hstack。

5. 特征数一致性检查，防止维度错位。

6. count_sketches = np.zeros((n_samples, degree, n_components))：分配三维累加器，第二维 degree 对应多项式各阶次的 Count Sketch。

7. 稀疏分支：按列遍历（CSC 高效），X_gamma[:, [j]] 保持二维切片，toarray().ravel() 转为一维加到对应哈希桶。

8. 稠密分支：直接向量化加法 iHashBit * X_gamma[:, j] 累加到 count_sketches[:, d, iHashIndex]。

9. FFT 加速多项式乘法核心：

   • fft(count_sketches, axis=2, overwrite_x=True)：对最后一维（哈希桶维度）做 FFT，利用卷积定理将时域卷积转为频域逐元素相乘。

   • np.prod(..., axis=1)：沿 degree 维相乘，等价于多项式展开中各阶次特征的张量积累加。

   • np.real(ifft(..., overwrite_x=True))：逆变换回时域，取实部（理论上虚部为数值噪声）。

10. 返回形状 (n_samples, n_components) 的近似特征映射。

代码总结：这段代码实现了张量草图的核心计算流程：输入缩放 → Count Sketch 哈希累加（稀疏/稠密统一） → FFT 频域多项式乘法 → IFFT 得到显式低维嵌入。overwrite_x=True 避免了中间数组拷贝，体现了工程上的内存优化。

18.4.4 标签系统与全局代码

源码路径：sklearn/kernel_approximation.py - PolynomialCountSketch.__sklearn_tags__（第172-176行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.input_tags.sparse = True
    return tags

声明支持稀疏输入（CSC 格式），使元估计器（如 Pipeline）能正确路由稀疏数据。

源码路径：sklearn/kernel_approximation.py - __main__（第1-30行）

模块级导入与常量定义，无可执行逻辑。

---

18.5 RBFSampler —— 基于随机傅里叶特征的 RBF 核近似

18.5.1 核心概念与类型定义

RBFSampler 基于 Bochner 定理 与 Random Kitchen Sinks 技术：RBF 核 \(K(x, y) = \exp(-\gamma \|x-y\|^2)\) 的傅里叶变换是高斯分布。从 \(\mathcal{N}(0, 2\gamma I)\) 采样随机权重 \(w\)，从 \([0, 2\pi)\) 采样随机偏移 \(b\)，显式特征映射为 \(z(x) = \sqrt{2/D} \cos(x w^T + b)\)，使 \(z(x)^T z(y) \approx K(x, y)\)。

源码路径：sklearn/kernel_approximation.py - RBFSampler.__init__（第181-197行）

def __init__(self, *, gamma=1.0, n_components=100, random_state=None):
    self.gamma = gamma
    self.n_components = n_components
    self.random_state = random_state

参数中 gamma 支持 'scale' 字符串（v1.2 新增），表示自适应计算 \(\gamma = 1 / (n_{\text{features}} \cdot \text{Var}(X))\)。

18.5.2 自适应 gamma 与随机投影采样 (fit)

源码路径：sklearn/kernel_approximation.py - RBFSampler.fit（第227-262行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
    """Fit the model with X.

    Samples random projection according to n_features.

    Parameters
    ----------
    X : {array-like, sparse matrix}, shape (n_samples, n_features)
        Training data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    y : array-like, shape (n_samples,) or (n_samples, n_outputs), \
            default=None
        Target values (None for unsupervised transformations).

    Returns
    -------
    self : object
        Returns the instance itself.
    """
    X = validate_data(self, X, accept_sparse="csr")
    random_state = check_random_state(self.random_state)
    n_features = X.shape[1]
    sparse = sp.issparse(X)
    if self.gamma == "scale":
        # var = E[X^2] - E[X]^2 if sparse
        X_var = (X.multiply(X)).mean() - (X.mean()) ** 2 if sparse else X.var()
        self._gamma = 1.0 / (n_features * X_var) if X_var != 0 else 1.0
    else:
        self._gamma = self.gamma
    self.random_weights_ = (2.0 * self._gamma) ** 0.5 * random_state.normal(
        size=(n_features, self.n_components)
    )

    self.random_offset_ = random_state.uniform(0, 2 * np.pi, size=self.n_components)

    if X.dtype == np.float32:
        # Setting the data type of the fitted attribute will ensure the
        # output data type during `transform`.
        self.random_weights_ = self.random_weights_.astype(X.dtype, copy=False)
        self.random_offset_ = self.random_offset_.astype(X.dtype, copy=False)
    self._n_features_out = self.n_components
    return self

逐行解析：

1. accept_sparse="csr"：接受 CSR 稀疏矩阵（便于行向量点积）。

2. sparse = sp.issparse(X)：判断稀疏性。

3. 自适应 gamma 计算：

   • 稀疏矩阵方差：(X.multiply(X)).mean() - (X.mean())**2，利用 \(Var(X) = E[X^2] - E[X]^2\) 避免显式稠密化。

   • 稠密矩阵直接用 X.var()。

   • _gamma = 1.0 / (n_features * X_var)，若方差为 0 退回 1.0。

4. random_weights_：形状 (n_features, n_components)，从 \(\mathcal{N}(0, 1)\) 采样并缩放 \(\sqrt{2\gamma}\)，预乘 \(\sqrt{2\gamma}\) 使得 transform 中只需做 \(X W\) 而无需额外缩放。

5. random_offset_：形状 (n_components,)，均匀采样 \([0, 2\pi)\)。

6. 数据类型保持：若输入为 float32，将拟合属性转为 float32（copy=False 避免拷贝），保证 transform 输出 dtype 与输入一致。这是 scikit-learn 统一的 dtype 传播 模式。

这段代码完成了随机傅里叶特征的参数采样与 dtype 固化，为 transform 的高效就地操作铺路。

18.5.3 高效前向投影 (transform)

源码路径：sklearn/kernel_approximation.py - RBFSampler.transform（第264-278行）

def transform(self, X):
    """Apply the approximate feature map to X.

    Parameters
    ----------
    X : {array-like, sparse matrix}, shape (n_samples, n_features)
        New data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    Returns
    -------
    X_new : array-like, shape (n_samples, n_components)
        Returns the instance itself.
    """
    check_is_fitted(self)

    X = validate_data(self, X, accept_sparse="csr", reset=False)
    projection = safe_sparse_dot(X, self.random_weights_)
    projection += self.random_offset_
    np.cos(projection, projection)
    projection *= (2.0 / self.n_components) ** 0.5
    return projection

逐行解析：

1. safe_sparse_dot(X, self.random_weights_)：稀疏安全点积，自动处理 CSR 稀疏与稠密权重矩阵的乘法，返回稠密投影向量。

2. projection += self.random_offset_：就地加偏移（广播机制）。

3. np.cos(projection, projection)：就地余弦，out=projection 复用内存，避免分配新数组。

4. projection *= (2.0 / self.n_components) ** 0.5：就地归一化缩放 \(\sqrt{2/D}\)。

5. 返回最终特征映射。

代码总结：这是一个极简且高性能的就地操作链：点积 → 加偏移 → 就地 cos → 就地缩放。三次就地操作（+=、np.cos(..., out=...)、*=）将内存占用压缩到单个 (n_samples, n_components) 缓冲区，体现了 scikit-learn 在数值计算层的工程极致。

18.5.4 标签系统

源码路径：sklearn/kernel_approximation.py - RBFSampler.__sklearn_tags__（第280-285行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.input_tags.sparse = True
    tags.transformer_tags.preserves_dtype = ["float64", "float32"]
    return tags

声明：支持稀疏输入；保持 float64/float32 dtype 不变。这使得元估计器可感知其 dtype 保持能力，避免不必要的类型转换。

---

18.6 SkewedChi2Sampler —— 偏态卡方核的随机傅里叶近似

18.6.1 核心概念

目标核为偏态乘法卡方核 \(K(x, y) = \prod_i \frac{2 (x_i+c)(y_i+c)}{x_i+y_i+2c}\)，定义域 \(x_i > -c\)。利用 \(\text{sech}\) 分布随机权重 \(w \sim \frac{1}{\pi} \log \tan(\frac{\pi}{2} u)\) 与均匀偏移 \(b\)，在对数域完成随机傅里叶映射。

18.6.2 数值稳健的拟合 (fit)

源码路径：sklearn/kernel_approximation.py - SkewedChi2Sampler.fit（第335-360行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
    """Fit the model with X.

    Samples random projection according to n_features.

    Parameters
    ----------
    X : array-like, shape (n_samples, n_features)
        Training data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    y : array-like, shape (n_samples,) or (n_samples, n_outputs), \
            default=None
        Target values (None for unsupervised transformations).

    Returns
    -------
    self : object
        Returns the instance itself.
    """
    X = validate_data(self, X)
    random_state = check_random_state(self.random_state)
    n_features = X.shape[1]
    uniform = random_state.uniform(size=(n_features, self.n_components))
    # transform by inverse CDF of sech
    self.random_weights_ = 1.0 / np.pi * np.log(np.tan(np.pi / 2.0 * uniform))
    self.random_offset_ = random_state.uniform(0, 2 * np.pi, size=self.n_components)

    if X.dtype == np.float32:
        # Setting the data type of the fitted attribute will ensure the
        # output data type during `transform`.
        self.random_weights_ = self.random_weights_.astype(X.dtype, copy=False)
        self.random_offset_ = self.random_offset_.astype(X.dtype, copy=False)

    self._n_features_out = self.n_components
    return self

逐行解析：

1. validate_data(self, X)：不接受稀疏（该核需逐元素对数变换，稀疏零值会导致 \(\log(0)\) 问题）。

2. uniform = random_state.uniform(...)：采样均匀分布 \(u \in [0,1)\)。

3. sech 分布逆变换采样：1.0 / np.pi * np.log(np.tan(np.pi / 2.0 * uniform))，这是 \(\text{sech}\) 分布的逆累积分布函数（ICDF），生成随机权重。

4. random_offset_：均匀采样偏移。

5. dtype 传播：同 RBFSampler，输入 float32 导致属性转为 float32。

18.6.3 对数域变换与投影 (transform)

源码路径：sklearn/kernel_approximation.py - SkewedChi2Sampler.transform（第362-385行）

def transform(self, X):
    """Apply the approximate feature map to X.

    Parameters
    ----------
    X : array-like, shape (n_samples, n_features)
        New data, where `n_samples` is the number of samples
        and `n_features` is the number of features. All values of X must be
        strictly greater than "-skewedness".

    Returns
    -------
    X_new : array-like, shape (n_samples, n_components)
        Returns the instance itself.
    """
    check_is_fitted(self)
    X = validate_data(
        self, X, copy=True, dtype=[np.float64, np.float32], reset=False
    )
    if (X <= -self.skewedness).any():
        raise ValueError("X may not contain entries smaller than -skewedness.")

    X += self.skewedness
    np.log(X, X)
    projection = safe_sparse_dot(X, self.random_weights_)
    projection += self.random_offset_
    np.cos(projection, projection)
    projection *= np.sqrt(2.0) / np.sqrt(self.n_components)
    return projection

逐行解析：

1. validate_data(..., copy=True, dtype=[np.float64, np.float32])：强制拷贝并限制 dtype，因为后续要做就地对数变换（np.log(X, X)），不能修改用户原始数据。

2. 边界检查：X <= -skewedness 触发 ValueError，保证定义域合法。

3. X += self.skewedness：就地平移 \(X \leftarrow X + c\)。

4. np.log(X, X)：就地取对数 \(X \leftarrow \log(X + c)\)。

5. safe_sparse_dot：点积（此处 X 为稠密，但统一接口）。

6. 后续与 RBFSampler 一致：加偏移 → 就地 cos → 归一化 \(\sqrt{2/D}\)。

代码总结：该实现展示了对数域数值稳健性设计：强制拷贝保护输入、就地对数变换避免中间数组、统一的点积/三角/缩放管线。

18.6.4 标签系统

源码路径：sklearn/kernel_approximation.py - SkewedChi2Sampler.__sklearn_tags__（第387-392行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.transformer_tags.preserves_dtype = ["float64", "float32"]
    return tags

声明 dtype 保持，不声明 sparse=True（因对数变换不支持稀疏零值），positive_only 隐含在验证逻辑中。

---

18.7 AdditiveChi2Sampler —— 加性卡方核的显式周期采样映射

18.7.1 无状态设计与参数验证

AdditiveChi2Sampler 近似加性卡方核 \(K(x, y) = \sum_i \frac{2 x_i y_i}{x_i + y_i}\)。采用确定性显式特征展开，无需随机采样，因此无状态（requires_fit=False），fit 仅做参数校验。

源码路径：sklearn/kernel_approximation.py - AdditiveChi2Sampler.__init__（第397-412行）

def __init__(self, *, sample_steps=2, sample_interval=None):
    self.sample_steps = sample_steps
    self.sample_interval = sample_interval

• sample_steps：采样阶数，每个特征展开为 \(2 \times \text{sample\_steps} - 1\) 个分量。

• sample_interval：采样间隔，sample_steps 在 {1,2,3} 时有默认值（0.8, 0.5, 0.4），否则必填。

源码路径：sklearn/kernel_approximation.py - AdditiveChi2Sampler.fit（第414-435行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
    """Only validates estimator's parameters.

    This method allows to: (i) validate the estimator's parameters and
    (ii) be consistent with the scikit-learn transformer API.

    Parameters
    ----------
    X : array-like, shape (n_samples, n_features)
        Training data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    y : array-like, shape (n_samples,) or (n_samples, n_outputs), \
            default=None
        Target values (None for unsupervised transformations).

    Returns
    -------
    self : object
        Returns the transformer.
    """
    X = validate_data(self, X, accept_sparse="csr", ensure_non_negative=True)

    if self.sample_interval is None and self.sample_steps not in (1, 2, 3):
        raise ValueError(
            "If sample_steps is not in [1, 2, 3],"
            " you need to provide sample_interval"
        )

    return self

逐行解析：

1. accept_sparse="csr", ensure_non_negative=True：接受 CSR 稀疏，强制非负校验（卡方核定义域）。

2. 参数校验：若 sample_interval 未提供且 sample_steps 不在 {1,2,3}，抛出 ValueError。

3. 返回 self，不存储任何拟合属性，体现无状态设计。

18.7.2 确定性特征展开 (transform)

源码路径：sklearn/kernel_approximation.py - AdditiveChi2Sampler.transform（第437-472行）

def transform(self, X):
    """Apply approximate feature map to X.

    Parameters
    ----------
    X : {array-like, sparse matrix}, shape (n_samples, n_features)
        Training data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    Returns
    -------
    X_new : {ndarray, sparse matrix}, \
           shape = (n_samples, n_features * (2*sample_steps - 1))
        Whether the return value is an array or sparse matrix depends on
        the type of the input X.
    """
    X = validate_data(
        self, X, accept_sparse="csr", reset=False, ensure_non_negative=True
    )
    sparse = sp.issparse(X)

    if self.sample_interval is None:
        # See figure 2 c) of "Efficient additive kernels via explicit feature maps"
        # <http://www.robots.ox.ac.uk/~vedaldi/assets/pubs/vedaldi11efficient.pdf>
        # A. Vedaldi and A. Zisserman, Pattern Analysis and Machine Intelligence,
        # 2011
        if self.sample_steps == 1:
            sample_interval = 0.8
        elif self.sample_steps == 2:
            sample_interval = 0.5
        elif self.sample_steps == 3:
            sample_interval = 0.4
        else:
            raise ValueError(
                "If sample_steps is not in [1, 2, 3],"
                " you need to provide sample_interval"
            )
    else:
        sample_interval = self.sample_interval

    # zeroth component
    # 1/cosh = sech
    # cosh(0) = 1.0
    transf = self._transform_sparse if sparse else self._transform_dense
    return transf(X, self.sample_steps, sample_interval)

逐行解析：

1. 验证输入，reset=False，非负校验。

2. sparse = sp.issparse(X)：分派稠密/稀疏实现。

3. 默认 sample_interval 查表：引用文献 Figure 2(c) 的经验值。

4. 根据稀疏性选择 _transform_sparse 或 _transform_dense。

18.7.3 稠密/稀疏双实现

源码路径：sklearn/kernel_approximation.py - AdditiveChi2Sampler._transform_dense（第500-525行）

@staticmethod
def _transform_dense(X, sample_steps, sample_interval):
    non_zero = X != 0.0
    X_nz = X[non_zero]

    X_step = np.zeros_like(X)
    X_step[non_zero] = np.sqrt(X_nz * sample_interval)

    X_new = [X_step]

    log_step_nz = sample_interval * np.log(X_nz)
    step_nz = 2 * X_nz * sample_interval

    for j in range(1, sample_steps):
        factor_nz = np.sqrt(step_nz / np.cosh(np.pi * j * sample_interval))

        X_step = np.zeros_like(X)
        X_step[non_zero] = factor_nz * np.cos(j * log_step_nz)
        X_new.append(X_step)

        X_step = np.zeros_like(X)
        X_step[non_zero] = factor_nz * np.sin(j * log_step_nz)
        X_new.append(X_step)

    return np.hstack(X_new)

数学公式对应：

• 第 0 项：\(\sqrt{x \cdot \Delta}\)

• 第 \(j\) 项：\(\sqrt{\frac{2x\Delta}{\cosh(\pi j \Delta)}} \cos(j \Delta \log x)\) 与 \(\sin\) 对

• 仅对非零元素计算（non_zero 掩码），零值保持 0，稀疏友好。

源码路径：sklearn/kernel_approximation.py - AdditiveChi2Sampler._transform_sparse（第527-560行）

@staticmethod
def _transform_sparse(X, sample_steps, sample_interval):
    indices = X.indices.copy()
    indptr = X.indptr.copy()

    data_step = np.sqrt(X.data * sample_interval)
    X_step = sp.csr_matrix(
        (data_step, indices, indptr), shape=X.shape, dtype=X.dtype, copy=False
    )
    X_new = [X_step]

    log_step_nz = sample_interval * np.log(X.data)
    step_nz = 2 * X.data * sample_interval

    for j in range(1, sample_steps):
        factor_nz = np.sqrt(step_nz / np.cosh(np.pi * j * sample_interval))

        data_step = factor_nz * np.cos(j * log_step_nz)
        X_step = sp.csr_matrix(
            (data_step, indices, indptr), shape=X.shape, dtype=X.dtype, copy=False
        )
        X_new.append(X_step)

        data_step = factor_nz * np.sin(j * log_step_nz)
        X_step = sp.csr_matrix(
            (data_step, indices, indptr), shape=X.shape, dtype=X.dtype, copy=False
        )
        X_new.append(X_step)

    return sp.hstack(X_new)

稀疏实现精髓：

• 复用 indices、indptr，仅变换 data 数组，保持稀疏结构不变。

• copy=False 避免拷贝索引数组。

• sp.hstack 水平拼接所有分量，输出仍为 CSR 稀疏矩阵，实现稀疏输入稀疏输出的零拷贝语义。

18.7.4 语义化特征命名 (get_feature_names_out)

源码路径：skernel_approximation.py - AdditiveChi2Sampler.get_feature_names_out（第474-498行）

def get_feature_names_out(self, input_features=None):
    """Get output feature names for transformation.

    Parameters
    ----------
    input_features : array-like of str or None, default=None
        Only used to validate feature names with the names seen in :meth:`fit`.

    Returns
    -------
    feature_names_out : ndarray of str objects
        Transformed feature names.
    """
    # Note that passing attributes="n_features_in_" forces check_is_fitted
    # to check if the attribute is present. Otherwise it will pass on this
    # stateless estimator (requires_fit=False)
    check_is_fitted(self, attributes="n_features_in_")
    input_features = _check_feature_names_in(
        self, input_features, generate_names=True
    )
    est_name = self.__class__.__name__.lower()

    names_list = [f"{est_name}_{name}_sqrt" for name in input_features]

    for j in range(1, self.sample_steps):
        cos_names = [f"{est_name}_{name}_cos{j}" for name in input_features]
        sin_names = [f"{est_name}_{name}_sin{j}" for name in input_features]
        names_list.extend(cos_names + sin_names)

    return np.asarray(names_list, dtype=object)

命名规范：additivechi2sampler_<原始名>_sqrt、_cos{j}、_sin{j}。顺序严格对应 _transform_dense/_sparse 的拼接顺序，便于模型解释与特征重要性分析。

18.7.5 标签系统

源码路径：sklearn/kernel_approximation.py - AdditiveChi2Sampler.__sklearn_tags__（第562-568行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.requires_fit = False
    tags.input_tags.positive_only = True
    tags.input_tags.sparse = True
    return tags

声明：无需拟合（requires_fit=False）、仅接受正值输入（positive_only=True）、支持稀疏输入。这是无状态 Transformer 的标准标签配置。

---

18.8 Nystroem —— 基于子采样的 Nyström 低秩核近似

18.8.1 核心流程 (fit)

Nyström 方法通过抽样 \(m\) 个基向量，计算基向量间核矩阵 \(K_{bb}\) 的 SVD 分解 \(U \Sigma V^T\)，构造归一化矩阵 \(U \Sigma^{-1/2} V^T\)，新样本嵌入为 \(Z = K_{xb} \cdot \text{normalization}^T\)。

源码路径：sklearn/kernel_approximation.py - Nystroem.fit（第623-685行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
    """Fit estimator to data.

    Samples a subset of training points, computes kernel
    on these and computes normalization matrix.

    Parameters
    ----------
    X : array-like, shape (n_samples, n_features)
        Training data, where `n_samples` is the number of samples
        and `n_features` is the number of features.

    y : array-like, shape (n_samples,) or (n_samples, n_outputs), \
            default=None
        Target values (None for unsupervised transformations).

    Returns
    -------
    self : object
        Returns the instance itself.
    """
    xp, _, device = get_namespace_and_device(X)
    X = validate_data(self, X, accept_sparse="csr")
    rnd = check_random_state(self.random_state)
    n_samples = X.shape[0]

    # get basis vectors
    if self.n_components > n_samples:
        # XXX should we just bail?
        n_components = n_samples
        warnings.warn(
            "n_components > n_samples. This is not possible.\n"
            "n_components was set to n_samples, which results"
            " in inefficient evaluation of the full kernel."
        )

    else:
        n_components = self.n_components
    n_components = min(n_samples, n_components)
    inds = rnd.permutation(n_samples)
    basis_inds = xp.asarray(inds[:n_components], dtype=xp.int64, device=device)
    if sp.issparse(X):
        basis = X[basis_inds]
    else:
        basis = _safe_indexing(X, basis_inds, axis=0)

    basis_kernel = pairwise_kernels(
        basis,
        metric=self.kernel,
        filter_params=True,
        n_jobs=self.n_jobs,
        **self._get_kernel_params(),
    )

    # sqrt of kernel matrix on basis vectors
    _, _, dtype = _find_floating_dtype_allow_sparse(basis_kernel, Y=None, xp=xp)
    basis_kernel = xp.asarray(basis_kernel, dtype=dtype, device=device)
    U, S, V = xp.linalg.svd(basis_kernel)
    S = xp.clip(S, 1e-12, None)
    self.normalization_ = U / xp.sqrt(S) @ V
    self.components_ = basis
    self.component_indices_ = basis_inds
    self._n_features_out = n_components
    return self

逐行解析：

1. xp, _, device = get_namespace_and_device(X)：Array API 兼容层入口，获取数组命名空间（NumPy/CuPy/JAX...）与设备。

2. validate_data(..., accept_sparse="csr")：接受 CSR 稀疏。

3. n_components > n_samples 保护：截断并警告。

4. inds = rnd.permutation(n_samples)：随机排列索引。

5. basis_inds = xp.asarray(...)：将索引转为当前后端数组（支持 GPU 设备）。

6. 稀疏/稠密分支索引：稀疏用 X[basis_inds]（CSR 支持花式索引）；稠密用 _safe_indexing 统一处理。

7. pairwise_kernels(...)：计算基向量间核矩阵，支持 n_jobs 并行。

8. _find_floating_dtype_allow_sparse：推断浮点 dtype（兼容稀疏）。

9. xp.asarray(basis_kernel, dtype=dtype, device=device)：确保核矩阵在正确后端/设备上。

10. U, S, V = xp.linalg.svd(basis_kernel)：后端无关的 SVD 调用。

11. S = xp.clip(S, 1e-12, None)：裁剪奇异值，防止除以极小值导致数值爆炸（对应测试 test_nystroem_singular_kernel）。

12. self.normalization_ = U / xp.sqrt(S) @ V：构造归一化矩阵 \(U \Sigma^{-1/2} V^T\)（注意：\(V\) 已为 \(V^T\) 形式，xp.linalg.svd 返回 \(V^H\)）。

13. 存储 components_、component_indices_。

代码总结：这是 Array API 多后端兼容的典范：从索引采样、核计算、SVD 到矩阵乘法，全程使用 xp 命名空间与 device，零硬编码 NumPy 调用，实现了“一次编写，多后端运行”。

18.8.2 高效嵌入新样本 (transform)

源码路径：sklearn/kernel_approximation.py - Nystroem.transform（第687-708行）

def transform(self, X):
    """Apply feature map to X.

    Computes an approximate feature map using the kernel
    between some training points and X.

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Data to transform.

    Returns
    -------
    X_transformed : ndarray of shape (n_samples, n_components)
        Transformed data.
    """
    check_is_fitted(self)

    xp, _, device = get_namespace_and_device(X)
    X = validate_data(self, X, accept_sparse="csr", reset=False)

    kernel_params = self._get_kernel_params()
    embedded = pairwise_kernels(
        X,
        self.components_,
        metric=self.kernel,
        filter_params=True,
        n_jobs=self.n_jobs,
        **kernel_params,
    )
    dtype = _find_matching_floating_dtype(embedded, xp=xp)
    embedded = xp.asarray(embedded, dtype=dtype, device=device)
    return embedded @ self.normalization_.T

逐行解析：

1. 获取后端 xp 与 device。

2. pairwise_kernels(X, self.components_, ...)：计算新样本与基向量的核矩阵 \(K_{xb}\)。

3. _find_matching_floating_dtype + xp.asarray：统一 dtype/device。

4. embedded @ self.normalization_.T：矩阵乘法完成嵌入 \(Z = K_{xb} N^T\)。

18.8.3 灵活内核参数管理 (_get_kernel_params)

源码路径：sklearn/kernel_approximation.py - Nystroem._get_kernel_params（第710-728行）

def _get_kernel_params(self):
    params = self.kernel_params
    if params is None:
        params = {}
    if not callable(self.kernel) and self.kernel != "precomputed":
        for param in KERNEL_PARAMS[self.kernel]:
            if getattr(self, param) is not None:
                params[param] = getattr(self, param)
    else:
        if (
            self.gamma is not None
            or self.coef0 is not None
            or self.degree is not None
        ):
            raise ValueError(
                "Don't pass gamma, coef0 or degree to "
                "Nystroem if using a callable "
                "or precomputed kernel"
            )

    return params

设计意图：

• 字符串内核：自动收集 gamma/degree/coef0 非空参数透传给 pairwise_kernels。

• 可调用内核 / precomputed：禁止传入上述参数，避免歧义（用户应通过 kernel_params 字典显式传递）。

18.8.4 标签系统

源码路径：sklearn/kernel_approximation.py - Nystroem.__sklearn_tags__（第730-736行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.array_api_support = True
    tags.input_tags.sparse = True
    tags.transformer_tags.preserves_dtype = ["float64", "float32"]
    return tags

声明：array_api_support=True（核心标记）、稀疏支持、dtype 保持。这是全模块唯一设置 array_api_support=True 的类，标志着其对跨后端计算的完整支持。

---

18.9 KernelRidge —— 核岭回归的对偶求解与预计算核支持

18.9.1 模型公式与统一内核入口

核岭回归目标：\(\min_\alpha \| K\alpha - y \|^2 + \alpha \alpha^T K \alpha\)，闭式解 \((K + \alpha I) \alpha = y\) → dual_coef_ = \alpha。预测：\(\hat{y} = K_{test, train} \alpha\)。

源码路径：sklearn/kernel_ridge.py - KernelRidge._get_kernel（第98-104行）

def _get_kernel(self, X, Y=None):
    if callable(self.kernel):
        params = self.kernel_params or {}
    else:
        params = {"gamma": self.gamma, "degree": self.degree, "coef0": self.coef0}
    return pairwise_kernels(X, Y, metric=self.kernel, filter_params=True, **params)

统一分派：可调用内核用 kernel_params 字典；字符串内核自动打包 gamma/degree/coef0。

18.9.2 标签系统

源码路径：sklearn/kernel_ridge.py - KernelRidge.__sklearn_tags__（第106-110行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.input_tags.sparse = True
    tags.input_tags.pairwise = self.kernel == "precomputed"
    return tags

动态标签：pairwise=True 仅当 kernel='precomputed' 时设为 True。这意味着：

• 普通内核：X 为特征矩阵，fit/predict 接收 (n_samples, n_features)。

• 预计算核：X 为核矩阵，fit 接收 (n_samples, n_samples)，predict 接收 (n_samples, n_samples_fitted)。元估计器据此调整数据路由。

18.9.3 Cholesky 求解与奇异回退 (fit)

源码路径：sklearn/kernel_ridge.py - KernelRidge.fit（第116-150行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y, sample_weight=None):
    """Fit Kernel Ridge regression model.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Training data. If kernel == "precomputed" this is instead
        a precomputed kernel matrix, of shape (n_samples, n_samples).

    y : array-like of shape (n_samples,) or (n_samples, n_targets)
        Target values.

    sample_weight : float or array-like of shape (n_samples,), default=None
        Individual weights for each sample, ignored if None is passed.

    Returns
    -------
    self : object
        Returns the instance itself.
    """
    # Convert data
    X, y = validate_data(
        self, X, y, accept_sparse=("csr", "csc"), multi_output=True, y_numeric=True
    )
    if sample_weight is not None and not isinstance(sample_weight, float):
        sample_weight = _check_sample_weight(sample_weight, X)

    K = self._get_kernel(X)
    alpha = np.atleast_1d(self.alpha)

    ravel = False
    if len(y.shape) == 1:
        y = y.reshape(-1, 1)
        ravel = True

    copy = self.kernel == "precomputed"
    self.dual_coef_ = _solve_cholesky_kernel(K, y, alpha, sample_weight, copy)
    if ravel:
        self.dual_coef_ = self.dual_coef_.ravel()

    self.X_fit_ = X

    return self

逐行解析：

1. validate_data(..., accept_sparse=("csr", "csc"), multi_output=True, y_numeric=True)：同时接受 CSR/CSC，支持多输出回归，强制 y 为数值类型。

2. 样本权重标准化：_check_sample_weight 处理数组/标量统一。

3. K = self._get_kernel(X)：计算核矩阵（或直接使用预计算核）。

4. alpha = np.atleast_1d(self.alpha)：支持标量或数组（多目标不同正则化）。

5. ravel 标记：记录 y 是否为一维，拟合后恢复 dual_coef_ 形状。

6. 关键 copy 逻辑：copy = self.kernel == "precomputed"。

   • 预计算核：copy=True，必须复制核矩阵，因为 _solve_cholesky_kernel 会原地修改对角线（加 \(\alpha\)），不能破坏用户传入的核矩阵。

   • 普通内核：copy=False，避免多余拷贝，_solve_cholesky_kernel 可原地操作临时核矩阵。

7. self.dual_coef_ = _solve_cholesky_kernel(...)：委托给 sklearn.linear_model._ridge._solve_cholesky_kernel，内部含 Cholesky 分解失败回退至 lstsq 的逻辑（测试 test_kernel_ridge_singular_kernel 验证）。

8. self.X_fit_ = X：保存训练数据（或预计算核矩阵）供 predict 使用。

18.9.4 预测阶段 (predict)

源码路径：sklearn/kernel_ridge.py - KernelRidge.predict（第152-165行）

def predict(self, X):
    """Predict using the kernel ridge model.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Samples. If kernel == "precomputed" this is instead a
        precomputed kernel matrix, shape = [n_samples,
        n_samples_fitted], where n_samples_fitted is the number of
        samples used in the fitting for this estimator.

    Returns
    -------
    C : ndarray of shape (n_samples,) or (n_samples, n_targets)
        Returns predicted values.
    """
    check_is_fitted(self)
    X = validate_data(self, X, accept_sparse=("csr", "csc"), reset=False)
    K = self._get_kernel(X, self.X_fit_)
    return np.dot(K, self.dual_coef_)

逐行解析：

1. 计算测试样本与训练样本（或预计算核矩阵行）的核矩阵 \(K_{test, train}\)。

2. np.dot(K, self.dual_coef_)：矩阵乘法得到预测值，自动处理单输出/多输出（dual_coef_ 形状为 (n_samples,) 或 (n_samples, n_targets)）。

---

18.10 核近似测试套件 —— 从数学正确性到跨后端一致性的全方位验证

18.10.1 无偏性与精度验证

测试 test_polynomial_count_sketch（第30-65行）对比精确多项式核与近似核：

error = kernel - kernel_approx
assert np.abs(np.mean(error)) <= 0.05  # 无偏性：均值误差 ≤ 0.05
np.abs(error, out=error)
assert np.max(error) <= 0.1  # 最大绝对误差 ≤ 0.1
assert np.mean(error) <= 0.05

测试 test_rbf_sampler 同理，容差更严（均值 ≤ 0.01）。

测试 test_skewed_chi2_sampler 在对数域验证相对误差。

测试 test_additive_chi2_sampler 使用解析解精确匹配（assert_array_almost_equal(kernel, kernel_approx, 1)）。

18.10.2 稠密/稀疏数值一致性

test_polynomial_count_sketch_dense_sparse、test_additive_chi2_sampler 等使用 @pytest.mark.parametrize("csr_container", CSR_CONTAINERS) 参数化，确保 CSR、COO、LIL 等多种稀疏容器与稠密结果位级一致（assert_allclose/assert_array_equal）。

18.10.3 数据类型传播与等价性

test_rbf_sampler_fitted_attributes_dtype / test_rbf_sampler_dtype_equivalence：

• 验证 float32 输入导致 random_weights_/random_offset_ 为 float32。

• 验证 float32 与 float64 相同随机种子下结果数值等价（assert_allclose）。

SkewedChi2Sampler 同理。

18.10.4 Array API 多后端一致性 (test_nystroem_approximation_array_api)

源码路径：sklearn/tests/test_kernel_approximation.py - test_nystroem_approximation_array_api（第340-395行）

@pytest.mark.parametrize(
    "array_namespace, device, dtype_name", yield_namespace_device_dtype_combinations()
)
@pytest.mark.parametrize(
    "kernel", list(kernel_metrics()) + [_linear_kernel, "precomputed"]
)
@pytest.mark.parametrize("n_components", [2, 100])
def test_nystroem_approximation_array_api(
    array_namespace, device, dtype_name, kernel, n_components
):
    xp = _array_api_for_tests(array_namespace, device)
    rnd = np.random.RandomState(0)
    n_samples = 10
    # Ensure full-rank linear kernel to limit the impact of device-specific
    # rounding discrepancies.
    n_features = 2 * n_samples
    X_np = rnd.uniform(size=(n_samples, n_features)).astype(dtype_name)
    if kernel == "precomputed":
        X_np = rbf_kernel(X_np[:n_components])

    X_xp = xp.asarray(X_np, device=device)

    nystroem = Nystroem(n_components=n_components, kernel=kernel, random_state=0)
    X_np_transformed = nystroem.fit_transform(X_np)

    with config_context(array_api_dispatch=True):
        X_xp_transformed = nystroem.fit_transform(X_xp)
        X_xp_transformed_np = _convert_to_numpy(X_xp_transformed, xp=xp)

        for attribute_name in ["components_", "normalization_"]:
            xp_attr, _, device_attr = get_namespace_and_device(
                getattr(nystroem, attribute_name)
            )
            assert xp_attr is xp
            assert device_attr == array_device(X_xp)

    atol = _atol_for_type(dtype_name)
    assert_allclose(X_np_transformed, X_xp_transformed_np, atol=atol)

关键点解析：

1. yield_namespace_device_dtype_combinations()：生成 (NumPy/CuPy/JAX, CPU/GPU, float32/float64) 组合。

2. n_features = 2 * n_samples：构造满秩线性核矩阵，避免奇异值截断导致的跨后端数值差异放大。

3. config_context(array_api_dispatch=True)：启用 Array API 分派路径。

4. 属性后端/设备校验：components_、normalization_ 必须位于正确后端与设备。

5. _atol_for_type(dtype_name)：根据 dtype 设置容差（float32 宽容度更大）。

18.10.5 边界条件与错误处理

• test_nystroem_singular_kernel：重复样本导致奇异核矩阵，SVD 裁剪保证有限输出。

• test_additive_chi2_sampler_exceptions / test_skewed_chi2_sampler：负值输入触发 ValueError。

• test_nystroem_callable / test_nystroem_precomputed_kernel：内核参数冲突检查（传入 gamma 等报错）。

• test_nystroem_precomputed_kernel：预计算核路径验证。

18.10.6 特征名语义正确性

• test_get_feature_names_out：所有近似器输出特征名符合 {estimator_lower}{index} 规范。

• test_additivechi2sampler_get_feature_names_out：验证 _sqrt、_cos{j}、_sin{j} 后缀顺序。

18.10.7 KernelRidge 回归专项测试

• test_kernel_ridge：线性核下与 Ridge(fit_intercept=False) 系数级对齐。

• test_kernel_ridge_sparse：稀疏输入与稠密结果一致。

• test_kernel_ridge_precomputed / ..._kernel_unchanged：预计算核路径等价且不修改输入（copy=True 生效）。

• test_kernel_ridge_sample_weights：样本权重在原始/预计算核下等价。

• test_kernel_ridge_multi_output：多目标回归与独立单目标拟合等价。

---

18.11 设计中的取舍

18.11.1 为什么 PolynomialCountSketch 使用 FFT 而非直接计算张量积？

直接计算 \(d\) 次多项式核的显式特征映射维度为 \(O(n_{\text{features}}^d)\)，指数级爆炸。Count Sketch 将其哈希压缩到固定 \(n_{\text{components}}\) 维，FFT 利用卷积定理在 \(O(n_{\text{components}} \log n_{\text{components}})\) 时间内完成频域多项式乘法，以可控的近似误差换取指数级的计算与存储降低。

18.11.2 为什么 RBFSampler/SkewedChi2Sampler 采用就地操作链？

就地操作（+=、np.cos(out=...)、*=）将内存占用从 \(O(3 \times n_{\text{samples}} \times n_{\text{components}})\) 压缩到 \(O(1 \times n_{\text{samples}} \times n_{\text{components}})\)，避免中间数组分配与 GC 压力。在大规模数据（如 \(10^5\) 样本 \(\times\) \(10^4\) 组件）下，这决定了能否在内存中完成计算。

18.11.3 为什么 AdditiveChi2Sampler 采用无状态设计（requires_fit=False）？

加性卡方核的显式特征映射有解析式确定性公式，无需从数据分布学习参数（随机采样仅用于 RBF/偏态卡方等无解析映射的核）。无状态设计使其可直接用于 transform，且天然支持流式/增量处理。但为了 API 一致性，建议调用 fit_transform 触发参数校验。

18.11.4 为什么 Nystroem 是唯一设置 array_api_support=True 的近似器？

Nyström 方法的核心计算（核矩阵构造、SVD、矩阵乘法）均为标准线性代数算子，完美契合 Array API 标准接口。而 PolynomialCountSketch 依赖 scipy.fft（非标准 API）、RBFSampler/SkewedChi2Sampler 依赖 safe_sparse_dot 与特定随机采样分布，当前生态尚无统一的稀疏/FFT/随机数 Array API 实现，故暂不开启。

18.11.5 为什么 KernelRidge 在 kernel='precomputed' 时强制 copy=True？

_solve_cholesky_kernel 需在核矩阵对角线原地加 \(\alpha\) 以求解 \((K + \alpha I) \alpha = y\)。若不复制，用户传入的预计算核矩阵会被破坏，违反“输入不可变”约定。普通内核由内部临时生成核矩阵，原地修改安全且节省内存。

---

18.12 动手练习

1. 阅读核近似器核心实现

   • 阅读 sklearn/kernel_approximation.py 中以下方法的实现：

     1. PolynomialCountSketch.transform（第120-170行）：理解 FFT 如何加速张量草图计算

     2. RBFSampler.fit（第227-262行）：理解自适应 gamma 计算与随机权重/偏移采样

     3. Nystroem.fit（第623-685行）：理解子采样、核矩阵 SVD 与归一化矩阵构造

   • 回答问题：

     • PolynomialCountSketch 中 count_sketches 的形状含义是什么？为什么要在 axis=2 上做 FFT？

     • RBFSampler 中 gamma='scale' 时如何计算方差？稀疏矩阵为何用 (X.multiply(X)).mean() - (X.mean())**2？

     • Nystroem.fit 中 normalization_ = U / sqrt(S) @ V 的数学含义是什么？为什么要裁剪奇异值？

2. 分析核岭回归求解路径

   • 阅读 sklearn/kernel_ridge.py 与 sklearn/linear_model/_ridge.py（_solve_cholesky_kernel）的相关代码：

     1. KernelRidge.fit（第116-150行）：数据验证、核矩阵计算、求解委托

     2. KernelRidge.predict（第152-165行）：预测阶段核矩阵计算

   • 回答问题：

     • KernelRidge 为何在 kernel='precomputed' 时设置 copy=True？普通内核为何 copy=False？

     • _solve_cholesky_kernel 如何处理 Cholesky 分解失败？这对应测试 test_kernel_ridge_singular_kernel 的什么场景？

     • 多输出回归时 dual_coef_ 的形状如何变化？ravel 标记何时生效？

3. 设计核近似器扩展测试

   • 参考 sklearn/tests/test_kernel_approximation.py 中的测试模式，为 AdditiveChi2Sampler 设计以下测试用例（仅描述测试思路，不写代码）：

     1. 验证 sample_steps=1,2,3 时默认 sample_interval 与文献一致性

     2. 验证 get_feature_names_out 输出顺序：先 _sqrt，再按 j=1..sample_steps-1 依次 _cos{j}、_sin{j}

     3. 验证稀疏输入下零值特征不产生额外分量（non_zero 掩码逻辑）

     4. 验证 sample_interval 显式传入时覆盖默认值，且 sample_steps=4 不再报错

   • 回答问题：

     • 如何构造已知解析解的加性卡方核数据来验证近似精度？

     • 测试稠密/稀疏数值一致性时，为何使用 CSR_CONTAINERS 参数化？

4. 探究 Array API 多后端一致性测试

   • 阅读 test_nystroem_approximation_array_api（第340-395行）的实现：

     1. 理解 yield_namespace_device_dtype_combinations 如何生成测试组合

     2. 理解 config_context(array_api_dispatch=True) 如何启用分派路径

     3. 理解 _atol_for_type 如何根据 dtype 设置容差

   • 回答问题：

     • 为何测试中生成 n_features = 2 * n_samples 的数据？这与线性核满秩有何关系？

     • get_namespace_and_device 在验证 components_ 与 normalization_ 属性时检查什么？

     • 若要在该测试中增加 PolynomialCountSketch 的 Array API 支持验证，需要修改哪些源码？

---

18.13 本章小结

这一章中我们学习了核近似与核岭回归的核心源码实现。首先，我们剖析了 PolynomialCountSketch 如何利用 Count Sketch 哈希与 FFT 加速实现多项式核的张量草图近似；其次，我们深入理解了 RBFSampler 基于 Bochner 定理的随机傅里叶特征、自适应 gamma 与就地操作链的高性能设计；接着，我们探讨了 SkewedChi2Sampler 在对数域的数值稳健变换与 sech 分布采样；然后，我们解读了 AdditiveChi2Sampler 的确定性显式周期采样映射、稠密/稀疏双实现与语义化特征命名；之后，我们详细分析了 Nystroem 的 Nyström 低秩近似流程、Array API 多后端兼容架构与灵活内核参数管理；最后，我们揭示了 KernelRidge 的对偶求解机制、预计算核优化、Cholesky 分解奇异回退与样本权重/多输出支持。同时我们也梳理了测试套件如何从无偏性精度、稠密稀疏一致性、dtype 传播、Array API 跨后端一致性、边界异常处理到特征名语义等多维度保障正确性。

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

| 概念 | 解释 |

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

| PolynomialCountSketch | 张量草图+FFT加速多项式核近似，支持稀疏输入 |

| RBFSampler | 随机傅里叶特征近似RBF核，自适应gamma、dtype保持、就地操作链 |

| SkewedChi2Sampler | 偏态卡方核随机傅里叶近似，sech分布采样、对数域变换、正值约束 |

| AdditiveChi2Sampler | 加性卡方核显式周期采样映射，无状态设计、稠密/稀疏双实现、语义化特征名 |

| Nystroem | Nyström低秩核近似，子采样+SVD归一化、Array API多后端支持、灵活内核参数 |

| KernelRidge | 核岭回归对偶求解，Cholesky分解+奇异回退lstsq、预计算核优化、样本权重与多输出支持 |

| 稀疏/稠密统一处理 | safe_sparse_dot、accept_sparse、CSR/CSC容器参数化测试保证数值一致性 |

| 数据类型传播 | float32输入导致拟合属性与输出dtype一致，避免隐式升级 |

| Array API 兼容 | Nystroem通过get_namespace_and_device/xp.linalg实现跨NumPy/CuPy/JAX统一 |

| 数学正确性验证 | 无偏性误差界、精确核匹配、相对误差容差、奇异矩阵鲁棒性 |

| 边界条件与错误处理 | 负值输入检查、参数冲突校验、样本数不足警告、预计算核不修改输入 |

下一章中，我们将学习判别分析与保序回归，走近“统计模型的经典传承”，解析 LDA/QDA 的判别函数与保序回归的 PAVA 算法和插值预测，展现统计学习在 scikit-learn 中的优雅实现。

18.14 架构与数据流图

graph TD A[kernel_approximation] --> B[kernel_approximation] B --> C[kernel_approximation]

sequenceDiagram participant U as 调用者 participant E as kernel_approximation participant C as kernel_approximation U->>E: 调用入口 E->>C: 传递参数 C-->>U: 返回结果

graph LR I[输入] --> P[参数校验] P --> T[核心处理] T --> O[输出]

graph TD L1[用户 API 层] --> L2[算法/服务层] L2 --> L3[数据结构层] L3 --> L4[运行时与依赖层]

上述图分别展示模块依赖、调用时序、数据流和架构分层。

第 19 章 —— 判别分析与保序回归 —— 走近“统计模型的经典传承”

19.1 学习目标

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

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

• 理解线性判别分析（LDA）与二次判别分析（QDA）的统计原理及判别函数推导

• 掌握LDA三大求解器（SVD、最小二乘、特征值分解）的适用场景与数值实现差异

• 掌握QDA的协方差估计策略（SVD/特征分解）、正则化机制（reg_param/shrinkage）及秩缺失处理

• 理解保序回归的PAVA算法原理、Cython加速实现与单调性自动检测（Spearman相关+Fisher变换）

• 熟悉LDA降维投影、QDA二次决策边界、IsotonicRegression插值预测与序列化机制

• 能阅读并对比判别分析与保序回归的测试策略：概率校验、正交性验证、正则化必要性、边界条件覆盖

想象你是一位城市交通规划师，需要为三类车辆设计导航系统：LDA 是共享地图的线性导航，所有车辆（类别）共用同一张路况图（共享协方差矩阵），规划直线路径（线性决策边界），提供三种路线规划算法——SVD 适合高维稀疏路网（不显式建图）、最小二乘适合任意路况图（支持收缩/自定义估计器）但不支持降维、特征分解适合需规划主干道（支持降维）的场景。QDA 是专属地图的曲线导航，每车辆（类别）持有私有地图（类别协方差），规划曲线路径（二次决策边界），但地图可能不全（秩缺失），需正则化填补——SVD 模式靠 reg_param 抬高海拔（奇异值平方加常数）却无法凭空造路（n_samples <= n_features 根本限制），Eigen 模式靠 shrinkage 向标准地图靠拢（Ledoit-Wolf/固定收缩）可修复断路。IsotonicRegression 是单调攀登的“台阶修路工”：PAVA 算法像推土机，将乱序高程（y）推平成不降/不升台阶（单调拟合），合并相邻违规段（Pool Adjacent Violators），再铺设线性插值桥梁（interp1d）供新车通行，自动识别上坡/下坡（Spearman+Fisher 置信区间）。就像导航系统需权衡地图精度与计算速度，判别分析在共享/独立协方差、线性/二次边界、三大求解器间权衡；保序回归在单调约束与数据拟合间平衡，Cython PAVA 保证 O(n) 极速铺路。

19.2 源码地图

graph TD subgraph discriminant_analysis["sklearn/discriminant_analysis.py"] direction TB utils["工具函数"] utils --> _cov["_cov(): 协方差估计(经验/Ledoit-Wolf/自定义)"] utils --> _class_means["_class_means(): 类均值计算(Array API/NumPy双路径)"] utils --> _class_cov["_class_cov(): 加权类内协方差矩阵"] mixin["DiscriminantAnalysisPredictionMixin"] mixin --> decision_function["decision_function(): 统一决策函数接口(二分类返回log likelihood ratio)"] mixin --> predict["predict(): argmax预测"] mixin --> predict_proba["predict_proba(): softmax/expit概率估计"] mixin --> predict_log_proba["predict_log_proba(): log_softmax对数概率"] LDA["LinearDiscriminantAnalysis"] LDA --> init_LDA["__init__(): 参数初始化(solver/shrinkage/priors等)"] LDA --> svd["_solve_svd(): SVD求解器(两阶段SVD、不显式计算协方差、Array API兼容)"] LDA --> lstsq["_solve_lstsq(): 最小二乘求解器(直接求解线性系统、不支持transform)"] LDA --> eigen["_solve_eigen(): 特征值分解求解器(广义特征值问题Sb v=λ Sw v、支持降维)"] LDA --> fit_LDA["fit(): 统一拟合入口(参数验证、先验处理、求解器分发)"] LDA --> transform["transform(): 判别子空间投影(svd: (X-xbar)@scalings; eigen: X@scalings)"] LDA --> predict_proba_LDA["predict_proba(): 二分类_expit/多分类_softmax"] LDA --> predict_log_proba_LDA["predict_log_proba(): 对数概率(数值稳定性处理)"] LDA --> get_feature_names["get_feature_names_out(): 特征名输出(lineardiscriminantanalysis0...)"] LDA --> tags["__sklearn_tags__(): Array API支持标记"] LDA --> decision_LDA["decision_function(): 线性决策函数(覆盖父类文档)"] QDA["QuadraticDiscriminantAnalysis"] QDA --> init_QDA["__init__(): 参数初始化(solver/shrinkage/reg_param等)"] QDA --> svd_QDA["_solve_svd(): 逐类SVD(奇异值平方为方差、reg_param正则化)"] QDA --> eigen_QDA["_solve_eigen(): 逐类特征分解(shrinkage/自定义协方差估计器)"] QDA --> fit_QDA["fit(): 逐类拟合(秩缺失检测、svd强制n>d、eigen可正则化修复)"] QDA --> _decision["_decision_function(): 二次判别函数(马氏距离+log先验)"] QDA --> decision_QDA["decision_function(): 二次决策函数(覆盖父类文档)"] end subgraph isotonic["sklearn/isotonic.py"] direction TB check_inc["check_increasing(): Spearman相关系数+Fisher变换置信区间自动判断单调方向"] iso_reg["isotonic_regression(): 统一入口(SciPy>=1.12用优化器,旧版回退Cython PAVA)"] IR["IsotonicRegression"] IR --> init_IR["__init__(): 参数(y_min/y_max/increasing/out_of_bounds)"] IR --> check_shape["_check_input_data_shape(): 输入形状验证(仅1D或单特征2D)"] IR --> build_y["_build_y(): 核心拟合流程(排序、去重_make_unique、PAVA、trim_duplicates)"] IR --> build_f["_build_f(): interp1d线性插值函数构建(out_of_bounds处理)"] IR --> fit_IR["fit(): 验证输入、调用_build_y/_build_f、存储阈值用于pickle"] IR --> _transform["_transform(): 统一预测/变换逻辑(clip/nan/raise越界处理)"] IR --> predict["predict()/transform(): 公共预测接口"] IR --> get_feature_names_IR["get_feature_names_out(): 输出特征名"] IR --> serialize["__getstate__/__setstate__(): 序列化支持(f_不可pickle,存阈值重建)"] IR --> tags_IR["__sklearn_tags__(): 1D数组输入标记"] IR --> make_unique["_make_unique(): 合并重复X值(加权平均y、按精度容差去重)"] end subgraph cython["sklearn/_isotonic.pyx"] direction TB pava["_inplace_contiguous_isotonic_regression(): PAVA核心(O(n)回溯合并块、原地修改)"] make_unique_cy["_make_unique(): Cython加速去重(floating类型模板、resolution容差)"] end subgraph tests_disc["sklearn/tests/test_discriminant_analysis.py"] direction TB test_lda_predict["test_lda_predict(): 9种solver/shrinkage组合验证fit/predict/proba一致性"] test_lda_proba["test_lda_predict_proba(): 理论后验概率验证(ESL书P.127公式、atol=1e-2)"] test_lda_priors["test_lda_priors(): 先验处理(负值报错、列表支持、非归一化警告重归一化)"] test_lda_coefs["test_lda_coefs(): 三求解器系数一致性(decimal=1)"] test_lda_ortho["test_lda_orthogonality(): 降维正交性验证(风筝形均值、类内协方差为单位阵)"] test_qda_reg["test_qda_regularization(): 秩缺失触发LinAlgError、reg_param/shrinkage修复、svd根本限制"] test_cov["test_covariance(): _cov工具函数验证(经验/auto对称性)"] test_qda_coefs["test_qda_coefs(): SVD/Eigen旋转缩放一致性"] test_qda_priors["test_qda_priors(): QDA先验概率影响预测分布"] test_cov_eq["协方差估计器等价性测试: shrinkage=0.5 vs ShrunkCovariance、auto vs StandardizedLedoitWolf"] end subgraph tests_iso["sklearn/tests/test_isotonic.py"] direction TB test_perm["test_permutation_invariance(): 打乱顺序拟合预测不变(回归测试样本权重排序Bug)"] test_ties_sec["test_isotonic_regression_ties_secondary_(): 对标R isotone包secondary ties方法"] test_ties_diff["test_isotonic_regression_with_ties_in_differently_sized_groups(): 不同组大小并列处理"] test_oob["test_isotonic_regression_oob_raise/clip/nan(): 三种越界模式验证"] test_fast["test_fast_predict(): trim_duplicates优化预测一致性(1000样本随机数据)"] test_dtype["test_isotonic_dtype(): int32/64/float32/64输入输出dtype匹配"] test_make_unique_tol["test_make_unique_tolerance(): 浮点精度容差去重验证(float64保留1e-14,float32合并)"] test_inf_slope["test_isotonic_non_regression_inf_slope(): 极小数值不产生inf斜率"] test_feat_names["test_get_feature_names_out(): 特征名输出isotonicregression0"] test_check_inc["test_check_increasing(): Spearman+Fisher CI跨零报警、小样本/极值处理"] test_shape["test_input_shape_validation(): 1D/2D单特征输入兼容、多特征报错"] end

19.3 LDA/QDA 核心框架 —— 统计判别模型的统一骨架

判别分析模块的设计核心在于通过类继承体系实现接口复用与功能分层。LinearDiscriminantAnalysis 继承自 LinearClassifierMixin、TransformerMixin 和 BaseEstimator，这意味着它不仅具备线性分类器的决策函数，还支持降维变换和流水线集成。QuadraticDiscriminantAnalysis 则继承自 DiscriminantAnalysisPredictionMixin、ClassifierMixin 和 BaseEstimator，专注于二次判别分类。两者均实现了 sklearn_tags 以支持 Array API 兼容性，体现了 scikit-learn 向多数组后端演进的架构决策。参数验证采用 _parameter_constraints 字典声明式约束，配合 _fit_context 装饰器管理拟合上下文，实现了编译期般的严格检查与嵌套验证跳过的运行时灵活性。核心协方差估计工具函数 _cov、_class_means、_class_cov 构成了统计计算的基础设施：_cov 支持经验协方差、Ledoit-Wolf 自动收缩、自定义估计器三种模式；_class_means 计算各类别均值，支持 Array API 与 NumPy 双路径实现；_class_cov 计算加权类内协方差矩阵，按先验概率加权求和。

19.3.1 类型定义详解

19.3.1.1 核心参数约束

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis._parameter_constraints (第187-210行)
_parameter_constraints: dict = {
    "solver": [StrOptions({"svd", "lsqr", "eigen"})],
    "shrinkage": [StrOptions({"auto"}), Interval(Real, 0, 1, closed="both"), None],
    "n_components": [Interval(Integral, 1, None, closed="left"), None],
    "priors": ["array-like", None],
    "store_covariance": ["boolean"],
    "tol": [Interval(Real, 0, None, closed="left")],
    "covariance_estimator": [HasMethods("fit"), None],
}

逐行注释：

1. solver 限制为三种字符串选项："svd"、"lsqr"、"eigen"，对应三大求解器

2. shrinkage 支持字符串 'auto'、0-1 间的浮点数或 None，控制协方差收缩强度

3. n_components 必须为正整数，限制降维输出维度上界

4. priors 接受数组或 None，指定类先验概率

5. store_covariance 布尔值，控制是否显式存储协方差矩阵（仅 SVD 求解器受影响）

6. tol 非负实数，SVD 求解器的奇异值显著性阈值

7. covariance_estimator 需具有 fit 方法的自定义估计器对象，与 shrinkage 互斥

代码总结：这段代码定义了 LDA 的参数约束字典，使用声明式验证替代手工 if-else 检查。StrOptions、Interval、HasMethods 等约束类在 _fit_context 装饰器作用下，会在 fit 调用前自动验证参数合法性，保证了接口的健壮性与可维护性。

19.3.1.2 协方差估计工具函数

# 第 19 章 —— sklearn/discriminant_analysis.py - _cov (第43-80行)
def _cov(X, shrinkage=None, covariance_estimator=None):
    """Estimate covariance matrix (using optional covariance_estimator)."""
    if covariance_estimator is None:
        shrinkage = "empirical" if shrinkage is None else shrinkage
        if isinstance(shrinkage, str):
            if shrinkage == "auto":
                sc = StandardScaler()  # standardize features
                X = sc.fit_transform(X)
                s = ledoit_wolf(X)[0]
                # rescale
                s = sc.scale_[:, np.newaxis] * s * sc.scale_[np.newaxis, :]
            elif shrinkage == "empirical":
                s = empirical_covariance(X)
        elif isinstance(shrinkage, Real):
            s = shrunk_covariance(empirical_covariance(X), shrinkage)
    else:
        if shrinkage is not None and shrinkage != 0:
            raise ValueError(
                "covariance_estimator and shrinkage parameters "
                "are not None. Only one of the two can be set."
            )
        covariance_estimator.fit(X)
        if not hasattr(covariance_estimator, "covariance_"):
            raise ValueError(
                "%s does not have a covariance_ attribute"
                % covariance_estimator.__class__.__name__
            )
        s = covariance_estimator.covariance_
    return s

逐行注释：

1. 若未提供自定义估计器，根据 shrinkage 参数选择估计策略：None 默认为 'empirical'

2. shrinkage='auto'：标准化特征后使用 Ledoit-Wolf 估计，再反标准化还原尺度

3. shrinkage='empirical'：直接计算经验协方差矩阵

4. shrinkage 为浮点数：在经验协方差基础上应用固定收缩强度

5. 若提供 covariance_estimator：校验与 shrinkage 互斥，调用其 fit 方法并提取 covariance_ 属性

6. 返回估计的协方差矩阵 s

代码总结：_cov 函数实现了三种协方差估计模式的统一调度，体现了策略模式。'auto' 模式通过标准化消除特征尺度影响，利用 Ledoit-Wolf 理论最优收缩，再反标准化还原，兼顾了数值稳定性与尺度不变性。自定义估计器接口扩展了灵活性，使上层求解器无需关心具体估计细节。

# 第 19 章 —— sklearn/discriminant_analysis.py - _class_means (第82-105行)
def _class_means(X, y):
    """Compute class means."""
    xp, is_array_api_compliant = get_namespace(X)
    classes, y = xp.unique_inverse(y)
    means = xp.zeros((classes.shape[0], X.shape[1]), device=device(X), dtype=X.dtype)

    if is_array_api_compliant:
        for i in range(classes.shape[0]):
            means[i, :] = xp.mean(X[y == i], axis=0)
    else:
        cnt = np.bincount(y)
        np.add.at(means, y, X)
        means /= cnt[:, None]
    return means

逐行注释：

1. 获取命名空间 xp 与 Array API 兼容性标志，支持 NumPy/CuPy/PyTorch 等后端

2. xp.unique_inverse(y) 返回唯一类别标签 classes 与整数编码后的 y

3. 预分配均值数组 means，形状 (n_classes, n_features)，保持设备与 dtype 一致

4. Array API 路径：逐类使用 xp.mean 计算均值，通用但可能较慢

5. NumPy 路径：np.bincount(y) 统计各类样本数，np.add.at(means, y, X) 向量化累加样本，最后按计数归一化

6. 返回类均值矩阵

代码总结：_class_means 计算各类别均值，展示了 Array API 兼容层的双路径实现。NumPy 路径利用 np.bincount 与 np.add.at 实现向量化加速，避免 Python 循环；Array API 路径保证跨后端一致行为。这种设计体现了 scikit-learn "性能优先、兼容兜底" 的工程哲学。

# 第 19 章 —— sklearn/discriminant_analysis.py - _class_cov (第107-140行)
def _class_cov(X, y, priors, shrinkage=None, covariance_estimator=None):
    """Compute weighted within-class covariance matrix."""
    classes = np.unique(y)
    cov = np.zeros(shape=(X.shape[1], X.shape[1]))
    for idx, group in enumerate(classes):
        Xg = X[y == group, :]
        cov += priors[idx] * np.atleast_2d(_cov(Xg, shrinkage, covariance_estimator))
    return cov

逐行注释：

1. 获取唯一类别标签 classes

2. 初始化零协方差矩阵 cov，形状 (n_features, n_features)

3. 遍历各类别：提取该类样本 Xg，调用 _cov 估计类协方差

4. np.atleast_2d 确保单样本类也能生成二维协方差矩阵

5. 按先验概率 priors[idx] 加权求和累积到 cov

6. 返回加权类内协方差矩阵

代码总结：_class_cov 实现了 LDA 共享协方差假设的数学基础——加权类内协方差矩阵。通过遍历类别、调用统一的 _cov 接口、按先验加权，实现了协方差估计策略与加权逻辑的解耦。该矩阵是 LDA 所有求解器的核心输入（SVD 除外，它隐式处理）。

---

19.4 LDA 三大求解器 —— 线性判别的数值优化三重奏

LDA 提供三种求解器，对应不同的数学推导路径与适用场景。_solve_lstsq 基于最小二乘直接求解线性系统，支持任意协方差估计器但不支持降维；_solve_eigen 求解广义特征值问题，同时支持分类与降维；_solve_svd 通过两阶段 SVD 避免显式计算协方差矩阵，适合高维数据并实验性支持 Array API。

flowchart TD subgraph LDA_Solvers["LDA 三大求解器对比"] direction TB SVD["_solve_svd()\n两阶段SVD\n不显式计算协方差\n支持Array API\n适合高维数据\n支持降维"] LSQR["_solve_lstsq()\n最小二乘求解线性系统\n支持任意协方差估计器\n不支持降维\n仅用于分类"] EIGEN["_solve_eigen()\n广义特征值分解\nSb v = λ Sw v\n同时支持分类与降维\n计算explained_variance_ratio"] end SVD -.->|共享输入| LSQR SVD -.->|共享输入| EIGEN LSQR -.->|共享输入| EIGEN

19.4.1 逐行解析关键函数

19.4.1.1 最小二乘求解器

# 第 19 章 —— sklearn/discriminant_analysis.py - _solve_lstsq (第142-185行)
def _solve_lstsq(self, X, y, shrinkage, covariance_estimator):
    """Least squares solver."""
    self.means_ = _class_means(X, y)
    self.covariance_ = _class_cov(
        X, y, self.priors_, shrinkage, covariance_estimator
    )
    self.coef_ = linalg.lstsq(self.covariance_, self.means_.T)[0].T
    self.intercept_ = -0.5 * np.diag(np.dot(self.means_, self.coef_.T)) + np.log(
        self.priors_
    )

逐行注释：

1. 计算类均值 means_，形状 (n_classes, n_features)

2. 计算加权类内协方差 covariance_，形状 (n_features, n_features)

3. 求解线性系统 covariance_ @ coef_.T = means_.T，得到系数矩阵 coef_，形状 (n_classes, n_features)

4. 计算截距项：-0.5 * diag(means_ @ coef_.T) + log(priors_)，对应判别函数常数项

代码总结：这段代码实现了基于最小二乘的 LDA 求解器。它直接利用协方差矩阵与类均值求解线性判别系数，数学上等价于最大化后验概率的线性判别函数。由于未计算特征向量，不支持 transform 降维，但支持任意协方差估计器（收缩、自定义），适合仅需分类且协方差需正则化的场景。

19.4.1.2 特征值分解求解器

# 第 19 章 —— sklearn/discriminant_analysis.py - _solve_eigen (第187-233行)
def _solve_eigen(self, X, y, shrinkage, covariance_estimator):
    """Eigenvalue solver."""
    self.means_ = _class_means(X, y)
    self.covariance_ = _class_cov(
        X, y, self.priors_, shrinkage, covariance_estimator
    )

    Sw = self.covariance_  # within scatter
    St = _cov(X, shrinkage, covariance_estimator)  # total scatter
    Sb = St - Sw  # between scatter

    evals, evecs = linalg.eigh(Sb, Sw)
    self.explained_variance_ratio_ = np.sort(evals / np.sum(evals))[::-1][
        : self._max_components
    ]
    evecs = evecs[:, np.argsort(evals)[::-1]]  # sort eigenvectors

    self.scalings_ = evecs
    self.coef_ = np.dot(self.means_, evecs).dot(evecs.T)
    self.intercept_ = -0.5 * np.diag(np.dot(self.means_, self.coef_.T)) + np.log(
        self.priors_
    )

逐行注释：

1. 计算类均值 means_ 与类内协方差 Sw

2. 计算总散度矩阵 St，类间散度 Sb = St - Sw

3. 求解广义特征值问题 Sb @ v = λ * Sw @ v，得到特征值 evals 与特征向量 evecs

4. 计算解释方差比：特征值归一化后降序排列，取前 _max_components 个

5. 特征向量按特征值降序排列，存入 scalings_

6. 计算判别系数 coef_ = means_ @ evecs @ evecs.T 与截距项

代码总结：这段代码实现了基于 Fisher 判别准则的特征值分解求解器。通过求解广义特征值问题 Sb v = λ Sw v，找到最大化类间散度与类内散度比率的投影方向。scalings_ 存储判别方向（特征向量），explained_variance_ratio_ 量化各判别成分捕获的类间方差比例。该求解器同时支持分类与降维，是 transform 方法的数学基础。

19.4.1.3 SVD 求解器

# 第 19 章 —— sklearn/discriminant_analysis.py - _solve_svd (第235-300行)
def _solve_svd(self, X, y):
    """SVD solver."""
    xp, is_array_api_compliant = get_namespace(X)

    if is_array_api_compliant:
        svd = xp.linalg.svd
    else:
        svd = scipy.linalg.svd

    n_samples, _ = X.shape
    n_classes = self.classes_.shape[0]

    self.means_ = _class_means(X, y)
    if self.store_covariance:
        self.covariance_ = _class_cov(X, y, self.priors_)

    Xc = []
    for idx, group in enumerate(self.classes_):
        Xg = X[y == group]
        Xc.append(Xg - self.means_[idx, :])

    self.xbar_ = self.priors_ @ self.means_

    Xc = xp.concat(Xc, axis=0)

    # 1) within (univariate) scaling by with classes std-dev
    std = xp.std(Xc, axis=0)
    # avoid division by zero in normalization
    std[std == 0] = 1.0
    fac = xp.asarray(1.0 / (n_samples - n_classes), dtype=X.dtype, device=device(X))

    # 2) Within variance scaling
    X = xp.sqrt(fac) * (Xc / std)
    # SVD of centered (within)scaled data
    _, S, Vt = svd(X, full_matrices=False)

    rank = xp.sum(xp.astype(S > self.tol, xp.int32))
    # Scaling of within covariance is: V' 1/S
    scalings = (Vt[:rank, :] / std).T / S[:rank]
    fac = 1.0 if n_classes == 1 else 1.0 / (n_classes - 1)

    # 3) Between variance scaling
    # Scale weighted centers
    X = (
        (xp.sqrt((n_samples * self.priors_) * fac)) * (self.means_ - self.xbar_).T
    ).T @ scalings
    # Centers are living in a space with n_classes-1 dim (maximum)
    # Use SVD to find projection in the space spanned by the
    # (n_classes) centers
    _, S, Vt = svd(X, full_matrices=False)

    if self._max_components == 0:
        self.explained_variance_ratio_ = xp.empty((0,), dtype=S.dtype)
    else:
        self.explained_variance_ratio_ = (S**2 / xp.sum(S**2))[
            : self._max_components
        ]

    rank = xp.sum(xp.astype(S > self.tol * S[0], xp.int32))
    self.scalings_ = scalings @ Vt.T[:, :rank]
    coef = (self.means_ - self.xbar_) @ self.scalings_
    self.intercept_ = -0.5 * xp.sum(coef**2, axis=1) + xp.log(self.priors_)
    self.coef_ = coef @ self.scalings_.T
    self.intercept_ -= self.xbar_ @ self.coef_.T

逐行注释：

1. 获取命名空间 xp，根据是否 Array API 兼容选择 svd 实现

2. 计算各类中心化数据 Xc = [Xg - mean_k]，整体均值 xbar_ = priors_ @ means_

3. 拼接所有中心化数据 Xc，计算特征标准差 std，零标准差置 1 避免除零

4. 计算缩放因子 fac = 1/(n_samples - n_classes)，对 Xc 做类内标准化并加权

5. 对标准化后的 X 做 SVD，得到奇异值 S 和右奇异向量 Vt

6. 根据容差 tol 确定秩 rank，构建类内缩放矩阵 scalings = (Vt[:rank]/std).T / S[:rank]

7. 计算类间缩放因子 fac = 1/(n_classes-1)，对加权类中心差 (means_ - xbar_) 施加 scalings

8. 对变换后的类中心矩阵再次 SVD，得到判别方向的奇异值 S 和向量 Vt

9. 计算解释方差比 S^2 / sum(S^2)，确定最终秩

10. 组合两阶段缩放：scalings_ = scalings @ Vt.T[:, :rank]

11. 计算判别系数 coef_ 与截距 intercept_，修正整体均值偏移

代码总结：这段代码实现了 LDA 的 SVD 求解器，核心创新在于两阶段 SVD 避免显式计算协方差矩阵。第一阶段对类内标准化数据 SVD 得到白化变换 scalings，第二阶段对变换后的类中心 SVD 得到判别方向。这种设计使得计算复杂度主要取决于 n_samples 而非 n_features，极其适合高维数据 (n_features > n_samples)。实验性支持 Array API 使其能在 CuPy、PyTorch 等后端运行。

---

19.5 QDA 求解器与正则化 —— 二次判别的几何与稳定性

QDA 为每个类别单独估计协方差矩阵，形成二次决策边界。_solve_svd 对每个类做 SVD，奇异值平方为方差，支持 reg_param 正则化；_solve_eigen 对类协方差做特征分解，支持 shrinkage 与自定义估计器。秩缺失检测通过 tol 阈值判断，SVD 模式强制要求 n_samples > n_features，Eigen 模式可通过正则化修复。

flowchart TD subgraph QDA_Solvers["QDA 两大求解器对比"] direction TB SVD_QDA["_solve_svd()\n逐类SVD\n奇异值平方为方差\nreg_param正则化\n要求 n_samples > n_features\n秩缺失无法修复"] EIGEN_QDA["_solve_eigen()\n逐类特征分解\n支持shrinkage/自定义协方差估计器\n正则化可修复秩缺失\n支持 n_samples <= n_features"] end SVD_QDA -.->|逐类独立拟合| EIGEN_QDA

19.5.1 逐行解析关键函数

19.5.1.1 QDA SVD 求解器

# 第 19 章 —— sklearn/discriminant_analysis.py - QuadraticDiscriminantAnalysis._solve_svd (第515-540行)
def _solve_svd(self, X):
    """SVD solver for Quadratic Discriminant Analysis."""
    n_samples, n_features = X.shape

    mean = X.mean(0)
    Xc = X - mean
    # Xc = U * S * V.T
    _, S, Vt = np.linalg.svd(Xc, full_matrices=False)
    scaling = (S**2) / (n_samples - 1)  # scalings are squared singular values
    scaling = ((1 - self.reg_param) * scaling) + self.reg_param
    rotation = Vt.T

    cov = None
    if self.store_covariance:
        # cov = V * (S^2 / (n-1)) * V.T
        cov = scaling * Vt.T @ Vt

    return scaling, rotation, cov

逐行注释：

1. 计算类内中心化数据 Xc = X - mean

2. 对 Xc 做 SVD，得到奇异值 S 和右奇异向量 Vt

3. 奇异值平方除以 (n-1) 得到方差估计 scaling（即特征值）

4. 应用正则化：scaling = (1-reg_param)*scaling + reg_param，将极小奇异值拉升防止奇异

5. 旋转矩阵 rotation = Vt.T 即主轴方向

6. 若 store_covariance，重构协方差矩阵 cov = V @ diag(scaling) @ V.T

7. 返回缩放向量、旋转矩阵、协方差矩阵

代码总结：这段代码实现了 QDA 的逐类 SVD 求解器。每个类别独立计算协方差的谱分解，scaling 存储主轴方差，rotation 存储主轴方向。reg_param 通过向奇异值平方添加常数实现岭类正则化，防止协方差矩阵奇异。但 SVD 本身要求 n_samples > n_features，否则奇异值数量不足导致根本性秩缺失，正则化也无法修复。

19.5.1.2 QDA 特征值分解求解器

# 第 19 章 —— sklearn/discriminant_analysis.py - QuadraticDiscriminantAnalysis._solve_eigen (第500-513行)
def _solve_eigen(self, X):
    """Eigenvalue solver."""
    n_samples, n_features = X.shape

    cov = _cov(X, self.shrinkage, self.covariance_estimator)
    scaling, rotation = linalg.eigh(cov)  # scalings are eigenvalues
    rotation = rotation[:, np.argsort(scaling)[::-1]]  # sort eigenvectors
    scaling = scaling[np.argsort(scaling)[::-1]]  # sort eigenvalues
    return scaling, rotation, cov

逐行注释：

1. 调用 _cov 估计协方差矩阵，支持 shrinkage（'auto' 或浮点数）与自定义估计器

2. 对协方差矩阵做特征分解 eigh，特征值即 scaling，特征向量即 rotation

3. 按特征值降序排列特征向量与特征值

4. 返回缩放向量（特征值）、旋转矩阵（特征向量）、协方差矩阵

代码总结：这段代码实现了 QDA 的逐类特征值分解求解器。通过 _cov 统一接口支持收缩估计与自定义估计器，shrinkage='auto' 触发 Ledoit-Wolf 自动收缩，有效缓解小样本高维下的协方差估计不稳。特征分解天然得到主轴方差（特征值）与方向（特征向量），支持通过正则化修复秩缺失。

19.5.1.3 QDA 统一拟合入口与秩缺失处理

# 第 19 章 —— sklearn/discriminant_analysis.py - QuadraticDiscriminantAnalysis.fit (第542-610行)
def fit(self, X, y):
    """Fit the model according to the given training data and parameters."""
    X, y = validate_data(self, X, y)
    check_classification_targets(y)
    self.classes_ = np.unique(y)
    n_samples, n_features = X.shape
    n_classes = len(self.classes_)
    if n_classes < 2:
        raise ValueError(
            "The number of classes has to be greater than one. Got "
            f"{n_classes} class."
        )
    if self.priors is None:
        _, cnts = np.unique(y, return_counts=True)
        self.priors_ = cnts / float(n_samples)
    else:
        self.priors_ = np.array(self.priors)

    if self.solver == "svd":
        if self.shrinkage is not None:
            raise NotImplementedError("shrinkage not supported with 'svd' solver.")
        if self.covariance_estimator is not None:
            raise ValueError(
                "covariance_estimator is not supported with solver='svd'. "
                "Try solver='eigen' instead."
            )
        specific_solver = self._solve_svd
    elif self.solver == "eigen":
        specific_solver = self._solve_eigen

    means = []
    cov = []
    scalings = []
    rotations = []
    for class_idx, class_label in enumerate(self.classes_):
        X_class = X[y == class_label, :]
        if len(X_class) == 1:
            raise ValueError(
                "y has only 1 sample in class %s, covariance is ill defined."
                % str(self.classes_[class_idx])
            )

        mean_class = X_class.mean(0)
        means.append(mean_class)

        scaling_class, rotation_class, cov_class = specific_solver(X_class)

        rank = np.sum(scaling_class > self.tol)
        if rank < n_features:
            n_samples_class = X_class.shape[0]
            if self.solver == "svd" and n_samples_class <= n_features:
                raise linalg.LinAlgError(
                    f"The covariance matrix of class {class_label} is not full "
                    f"rank. When using `solver='svd'` the number of samples in "
                    f"each class should be more than the number of features, but "
                    f"class {class_label} has {n_samples_class} samples and "
                    f"{n_features} features. Try using `solver='eigen'` and "
                    f"setting the parameter `shrinkage` for regularization."
                )
            else:
                msg_param = "shrinkage" if self.solver == "eigen" else "reg_param"
                raise linalg.LinAlgError(
                    f"The covariance matrix of class {class_label} is not full "
                    f"rank. Increase the value of `{msg_param}` to reduce the "
                    f"collinearity.",
                )

        cov.append(cov_class)
        scalings.append(scaling_class)
        rotations.append(rotation_class)

    if self.store_covariance:
        self.covariance_ = cov
    self.means_ = np.asarray(means)
    self.scalings_ = scalings
    self.rotations_ = rotations
    return self

逐行注释：

1. 验证输入数据，检查分类目标，获取类别标签

2. 计算或验证先验概率 priors_

3. 根据 solver 选择具体求解器，校验参数互斥性（svd 不支持 shrinkage/estimator）

4. 逐类循环：提取类样本 X_class，检查单样本类报错

5. 计算类均值，调用具体求解器获取 scaling、rotation、cov

6. 秩缺失检测：统计 scaling > tol 的数量，若小于 n_features 则秩缺失

7. SVD 模式下 n_samples <= n_features：直接报错，建议改用 eigen+shrinkage

8. Eigen 模式或 SVD 其他情况：提示增加 shrinkage 或 reg_param 正则化

9. 收集各类结果，存储 means_、scalings_、rotations_、可选 covariance_

代码总结：这段代码实现了 QDA 的统一拟合流程，核心在于逐类拟合与秩缺失的差异化处理。SVD 求解器受限于奇异值数量不超过样本数，当 n_samples <= n_features 时根本性秩缺失无法通过 reg_param 修复，必须改用 Eigen 求解器配合 shrinkage。Eigen 求解器通过 _cov 的收缩机制（Ledoit-Wolf 或固定收缩）可使协方差矩阵满秩，修复秩缺失。这种设计体现了数值线性代数约束对算法选择的指导意义。

---

19.6 判别函数与预测逻辑 —— 从后验概率到决策边界

判别分析的预测逻辑统一封装在 DiscriminantAnalysisPredictionMixin 中。LDA 的决策函数为线性形式 X @ coef_.T + intercept_，二分类特殊处理返回 log likelihood ratio；QDA 的决策函数基于马氏距离计算二次型。两者共享 predict、predict_proba、predict_log_proba 实现。

flowchart TD subgraph Prediction["判别分析预测流程"] direction TB mixin["DiscriminantAnalysisPredictionMixin"] mixin --> decision["decision_function()\n调用子类_decision_function\n二分类返回log likelihood ratio"] mixin --> predict["predict()\nargmax决策"] mixin --> proba["predict_proba()\nexp(predict_log_proba)"] mixin --> log_proba["predict_log_proba()\nLogSoftmax数值稳定"] LDA_dec["LinearDiscriminantAnalysis._decision_function\n线性: X @ coef_.T + intercept_"] QDA_dec["QuadraticDiscriminantAnalysis._decision_function\n二次: 马氏距离 + log|Σ| + log先验"] decision -.->|委托| LDA_dec decision -.->|委托| QDA_dec end

19.6.1 逐行解析关键函数

19.6.1.1 预测混入类核心方法

# 第 19 章 —— sklearn/discriminant_analysis.py - DiscriminantAnalysisPredictionMixin.decision_function (第142-165行)
def decision_function(self, X):
    """Apply decision function to an array of samples."""
    y_scores = self._decision_function(X)
    if len(self.classes_) == 2:
        return y_scores[:, 1] - y_scores[:, 0]
    return y_scores

逐行注释：

1. 调用子类实现的 _decision_function 获取各类得分 y_scores，形状 (n_samples, n_classes)

2. 二分类时返回差值 y_scores[:, 1] - y_scores[:, 0]，即 log likelihood ratio

3. 多分类直接返回各类得分矩阵

# 第 19 章 —— sklearn/discriminant_analysis.py - DiscriminantAnalysisPredictionMixin.predict_log_proba (第180-195行)
def predict_log_proba(self, X):
    """Estimate log class probabilities."""
    scores = self._decision_function(X)
    log_likelihood = scores - scores.max(axis=1)[:, np.newaxis]
    return log_likelihood - np.log(
        np.exp(log_likelihood).sum(axis=1)[:, np.newaxis]
    )

逐行注释：

1. 获取判别得分 scores

2. 数值稳定性处理：减去每行最大值 log_likelihood = scores - max(scores)

3. LogSoftmax：log_likelihood - log(sum(exp(log_likelihood)))

代码总结：这段代码实现了判别分析的统一预测接口。decision_function 将子类的具体判别函数（LDA 线性、QDA 二次）适配为统一接口，二分类返回标量 log likelihood ratio。predict_log_proba 使用 LogSoftmax 技巧将任意尺度的判别得分转换为对数概率，减去行最大值避免数值溢出。

19.6.1.2 LDA 线性决策函数

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis._decision_function (第340-360行)
def _decision_function(self, X):
    """Apply decision function to an array of samples."""
    check_is_fitted(self)
    X = validate_data(self, X, reset=False)
    return np.dot(X, self.coef_.T) + self.intercept_

逐行注释：

1. 检查拟合状态，验证输入数据

2. 返回线性决策函数：X @ coef_.T + intercept_，形状 (n_samples, n_classes)

代码总结：这段代码实现了 LDA 的核心线性判别函数。coef_ 和 intercept_ 在 fit 阶段由求解器计算，_decision_function 直接应用线性变换。该方法被基类 decision_function 调用，二分类时返回 log likelihood ratio。

19.6.1.3 LDA 概率预测

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis.predict_proba (第360-375行)
def predict_proba(self, X):
    """Estimate probability."""
    check_is_fitted(self)
    xp, _ = get_namespace(X)
    decision = self.decision_function(X)
    if size(self.classes_) == 2:
        proba = _expit(decision, xp=xp)
        return xp.stack([1 - proba, proba], axis=1)
    else:
        return softmax(decision)

逐行注释：

1. 检查拟合状态，获取命名空间 xp

2. 计算决策函数值 decision

3. 二分类：使用 _expit (sigmoid) 计算正类概率，堆叠为 [1-p, p]

4. 多分类：使用 softmax 归一化

代码总结：这段代码实现了 LDA 的概率预测，体现了 Array API 兼容设计。二分类使用 sigmoid (_expit)，多分类使用 softmax，均通过 xp 命名空间适配不同数组后端。

19.6.1.4 LDA 对数概率预测

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis.predict_log_proba (第377-395行)
def predict_log_proba(self, X):
    """Estimate log probability."""
    xp, _ = get_namespace(X)
    prediction = self.predict_proba(X)

    smallest_normal = xp.finfo(prediction.dtype).smallest_normal
    prediction[prediction == 0.0] += smallest_normal
    return xp.log(prediction)

逐行注释：

1. 获取命名空间 xp，调用 predict_proba 获取概率

2. 数值稳定性处理：将 0 概率替换为该 dtype 的最小正规数 smallest_normal

3. 返回对数概率

代码总结：这段代码实现了 LDA 的对数概率预测，通过添加 smallest_normal 避免 log(0) 产生 -inf。相比基类的 LogSoftmax 实现，LDA 复用 predict_proba 结果，保证了二分类 sigmoid 与多分类 softmax 的数值一致性。

19.6.1.5 LDA 特征名称输出

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis.get_feature_names_out (第410-430行)
def get_feature_names_out(self, input_features=None):
    """Get output feature names for transformation."""
    check_is_fitted(self, "scalings_")
    class_name = self.__class__.__name__.lower()
    return np.asarray(
        [f"{class_name}{i}" for i in range(self._n_features_out)], dtype=object
    )

逐行注释：

1. 检查 scalings_ 存在性（已拟合）

2. 类名小写 lineardiscriminantanalysis

3. 生成特征名数组：lineardiscriminantanalysis0, lineardiscriminantanalysis1, ...

代码总结：这段代码实现了 ClassNamePrefixFeaturesOutMixin 接口，为降维输出生成语义化特征名。前缀使用类名小写，后缀为组件索引，便于下游流水线阶段引用。

19.6.1.6 LDA Array API 支持标记

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis.__sklearn_tags__ (第432-440行)
def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.array_api_support = True
    return tags

逐行注释：

1. 调用父类方法获取标签字典

2. 设置 array_api_support = True，标记支持 Array API

3. 返回更新后的标签

代码总结：这段代码声明了 LDA 支持 Array API 标准，使其可在 CuPy、PyTorch 等兼容数组库上运行。这是 scikit-learn 向多数组后端演进的关键标记。

19.6.1.7 LDA 决策函数文档覆盖

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis.decision_function (第397-410行)
def decision_function(self, X):
    """Apply decision function to an array of samples.

    The decision function is equal (up to a constant factor) to the
    log-posterior of the model, i.e. `log p(y = k | x)`. In a binary
    classification setting this instead corresponds to the difference
    `log p(y = 1 | x) - log p(y = 0 | x)`. See :ref:`lda_qda_math`.

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Array of samples (test vectors).

    Returns
    -------
    y_scores : ndarray of shape (n_samples,) or (n_samples, n_classes)
        Decision function values related to each class, per sample.
        In the two-class case, the shape is `(n_samples,)`, giving the
        log likelihood ratio of the positive class.
    """
    # Only overrides for the docstring.
    return super().decision_function(X)

逐行注释：

1. 仅覆盖父类文档字符串，提供 LDA 专用的数学解释

2. 调用父类 decision_function 实现（即 DiscriminantAnalysisPredictionMixin.decision_function）

代码总结：这段代码通过覆盖 docstring 为 LDA 提供了专门的数学文档，阐明决策函数与对数后验的关系，二分类时对应 log likelihood ratio。实际逻辑复用父类实现，体现了文档与实现分离的设计。

19.6.1.8 QDA 初始化

# 第 19 章 —— sklearn/discriminant_analysis.py - QuadraticDiscriminantAnalysis.__init__ (第615-640行)
def __init__(
    self,
    *,
    solver="svd",
    shrinkage=None,
    priors=None,
    reg_param=0.0,
    store_covariance=False,
    tol=1.0e-4,
    covariance_estimator=None,
):
    self.solver = solver
    self.shrinkage = shrinkage
    self.priors = priors
    self.reg_param = reg_param
    self.store_covariance = store_covariance
    self.tol = tol
    self.covariance_estimator = covariance_estimator

逐行注释：

1. 关键字参数初始化，强制使用命名参数

2. solver 默认 'svd'，可选 'eigen'

3. shrinkage 支持 'auto'、float、None，仅 eigen 支持

4. reg_param 默认 0.0，仅 svd 支持，正则化强度

5. store_covariance 控制是否存储类协方差

6. tol 秩缺失判断阈值

7. covariance_estimator 自定义协方差估计器，仅 eigen 支持

代码总结：这段代码定义了 QDA 的构造函数，参数设计体现了两种求解器的能力差异：svd 独享 reg_param，eigen 独享 shrinkage 和 covariance_estimator。_parameter_constraints 会在 fit 前验证这些互斥关系。

19.6.1.9 QDA 二次判别函数

# 第 19 章 —— sklearn/discriminant_analysis.py - QuadraticDiscriminantAnalysis._decision_function (第612-630行)
def _decision_function(self, X):
    # return log posterior, see eq (4.12) p. 110 of the ESL.
    check_is_fitted(self)

    X = validate_data(self, X, reset=False)
    norm2 = []
    for i in range(len(self.classes_)):
        R = self.rotations_[i]
        S = self.scalings_[i]
        Xm = X - self.means_[i]
        X2 = np.dot(Xm, R * (S ** (-0.5)))
        norm2.append(np.sum(X2**2, axis=1))
    norm2 = np.array(norm2).T  # shape = [len(X), n_classes]
    u = np.asarray([np.sum(np.log(s)) for s in self.scalings_])
    return -0.5 * (norm2 + u) + np.log(self.priors_)

逐行注释：

1. 验证输入数据，逐类循环计算

2. 获取第 i 类的旋转矩阵 R 和缩放向量 S

3. 中心化：Xm = X - means_[i]

4. 白化变换：X2 = Xm @ R * S^(-0.5)，即投影到主轴并按标准差缩放

5. 计算马氏距离平方：norm2 = sum(X2^2, axis=1)

6. 收集所有类的 norm2，转置为 (n_samples, n_classes)

7. 计算协方差对数行列式项：u = sum(log(s)) 对应 log|Σ|

8. 返回对数后验：-0.5 * (马氏距离^2 + log|Σ|) + log(prior)

代码总结：这段代码实现了 QDA 的核心判别函数，数学对应高斯判别分析的对数后验概率。通过 rotations_ 和 scalings_ 实现马氏距离计算：先旋转到主轴坐标系，再按标准差缩放，最后计算欧氏距离平方。u 项对应协方差行列式的对数，惩罚方差大的类。最终加上对数先验得到对数后验，predict 取 argmax 即贝叶斯最优决策。

19.6.1.10 QDA 决策函数文档覆盖

# 第 19 章 —— sklearn/discriminant_analysis.py - QuadraticDiscriminantAnalysis.decision_function (第632-650行)
def decision_function(self, X):
    """Apply decision function to an array of samples.

    The decision function is equal (up to a constant factor) to the
    log-posterior of the model, i.e. `log p(y = k | x)`. In a binary
    classification setting this instead corresponds to the difference
    `log p(y = 1 | x) - log p(y = 0 | x)`. See :ref:`lda_qda_math`.

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Array of samples (test vectors).

    Returns
    -------
    C : ndarray of shape (n_samples,) or (n_samples, n_classes)
        Decision function values related to each class, per sample.
        In the two-class case, the shape is `(n_samples,)`, giving the
        log likelihood ratio of the positive class.
    """
    # Only overrides for the docstring.
    return super().decision_function(X)

逐行注释：

1. 仅覆盖父类文档字符串，提供 QDA 专用的数学解释

2. 调用父类 decision_function 实现（即 DiscriminantAnalysisPredictionMixin.decision_function）

代码总结：这段代码通过覆盖 docstring 为 QDA 提供了专门的数学文档，阐明决策函数与对数后验的关系，二分类时对应 log likelihood ratio。实际逻辑复用父类实现，体现了文档与实现分离的设计。

---

19.7 LDA 降维与特征输出 —— 判别子空间投影

LDA 的 transform 方法将数据投影到判别子空间，仅 svd 和 eigen 求解器支持。SVD 模式以整体均值为中心投影，Eigen 模式直接投影到特征向量。explained_variance_ratio_ 量化各判别成分捕获的类间方差比例。

flowchart TD subgraph LDA_Transform["LDA 降维投影流程"] direction TB check["检查 solver != 'lsqr'"] svd_proj["SVD: (X - xbar_) @ scalings_\n以整体均值为中心"] eigen_proj["Eigen: X @ scalings_\n直接投影到特征向量"] slice["截取前 _max_components 列"] end check --> svd_proj check --> eigen_proj svd_proj --> slice eigen_proj --> slice

19.7.1 逐行解析关键函数

19.7.1.1 降维变换

# 第 19 章 —— sklearn/discriminant_analysis.py - LinearDiscriminantAnalysis.transform (第377-395行)
def transform(self, X):
    """Project data to maximize class separation."""
    if self.solver == "lsqr":
        raise NotImplementedError(
            "transform not implemented for 'lsqr' solver (use 'svd' or 'eigen')."
        )
    check_is_fitted(self)
    X = validate_data(self, X, reset=False)

    if self.solver == "svd":
        X_new = (X - self.xbar_) @ self.scalings_
    elif self.solver == "eigen":
        X_new = X @ self.scalings_

    return X_new[:, : self._max_components]

逐行注释：

1. lsqr 求解器不支持降维，直接报错

2. 检查拟合状态，验证输入数据

3. SVD 模式：(X - xbar_) @ scalings_，以整体均值 xbar_ 为中心投影

4. Eigen 模式：X @ scalings_，直接投影到特征向量（特征向量已为正交基）

5. 截取前 _max_components 个成分

代码总结：这段代码实现了 LDA 的降维投影。两种求解器的投影中心不同：SVD 基于整体均值 xbar_（白化后的坐标系），Eigen 基于原始坐标系原点（特征向量正交基）。这种差异源于两阶段 SVD 的数学推导：第一阶段白化以 xbar_ 为中心，第二阶段在白化空间投影。explained_variance_ratio_ 在两种模式下分别由奇异值平方和特征值归一化得到，反映判别方向的类间分离能力。

---

19.8 PAVA 算法核心 —— 保序回归的 Cython 极速实现

保序回归的核心是 Pool Adjacent Violators Algorithm (PAVA)，在 sklearn/_isotonic.pyx 中用 Cython 实现，达到 O(n) 线性时间复杂度。_make_unique 预处理合并重复 X 值，按浮点精度容差去重。

flowchart TD subgraph PAVA["PAVA 算法流程"] direction TB init["初始化: target[i]=i, 每个点自成块\nw[i]=权重, y[i]=目标值"] loop["主循环 i=0..n-1\nk=target[i]+1 (下一块起始)"] check["y[i] < y[k]?\n是: i=k 继续\n否: 进入违规合并"] merge["违规合并:\n累积 sum_wy, sum_w\nk=target[k]+1 跳跃\n直到 k==n 或 prev_y < y[k]"] update["更新块起始 i:\ny[i]=sum_wy/sum_w\nw[i]=sum_w\ntarget[i]=k-1, target[k-1]=i"] backtrack["回溯: i=target[i-1]\n保证单次遍历 O(n)"] recon["重构解:\n遍历target链表\n块内所有值设为块均值 y[i]"] end init --> loop loop --> check check -- 是 --> loop check -- 否 --> merge merge --> update update --> backtrack backtrack --> loop loop -- 结束 --> recon

19.8.1 逐行解析关键函数

19.8.1.1 PAVA 核心实现

# 第 19 章 —— sklearn/_isotonic.pyx - _inplace_contiguous_isotonic_regression (第15-60行)
def _inplace_contiguous_isotonic_regression(floating[::1] y, floating[::1] w):
    cdef:
        Py_ssize_t n = y.shape[0], i, k
        floating prev_y, sum_wy, sum_w
        Py_ssize_t[::1] target = np.arange(n, dtype=np.intp)

    # target describes a list of blocks.  At any time, if [i..j] (inclusive) is
    # an active block, then target[i] := j and target[j] := i.

    # For "active" indices (block starts):
    # w[i] := sum{w_orig[j], j=[i..target[i]]}
    # y[i] := sum{y_orig[j]*w_orig[j], j=[i..target[i]]} / w[i]

    with nogil:
        i = 0
        while i < n:
            k = target[i] + 1
            if k == n:
                break
            if y[i] < y[k]:
                i = k
                continue
            sum_wy = w[i] * y[i]
            sum_w = w[i]
            while True:
                # We are within a decreasing subsequence.
                prev_y = y[k]
                sum_wy += w[k] * y[k]
                sum_w += w[k]
                k = target[k] + 1
                if k == n or prev_y < y[k]:
                    # Non-singleton decreasing subsequence is finished,
                    # update first entry.
                    y[i] = sum_wy / sum_w
                    w[i] = sum_w
                    target[i] = k - 1
                    target[k - 1] = i
                    if i > 0:
                        # Backtrack if we can.  This makes the algorithm
                        # single-pass and ensures O(n) complexity.
                        i = target[i - 1]
                    # Otherwise, restart from the same point.
                    break
        # Reconstruct the solution.
        i = 0
        while i < n:
            k = target[i] + 1
            y[i + 1 : k] = y[i]
            i = k