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

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

第 ② 行创建 RandomState 实例而非依赖全局种子，保证函数自身可复现。

第 ③ 行使用 COO（Coordinate）格式而非直接 CSR——COO 构造最简单，只需 (data, (row, col)) 三元组。

第 ④ 行非零值使用 rng.randn（正态分布）而非 rng.rand（均匀分布），模拟真实高维稀疏特征的"重尾"分布。

第 ⑤、⑥ 行通过 rng.randint 生成行/列索引，注意 randint(upper) 返回 [0, upper) 区间。

第 ⑦ 行指定矩阵形状。

第 ⑧ 行通过 toarray() 与 tocsr() 同时返回稠密与 CSR 版本，受主流程的 opts.dense 标志控制选择哪个。

这段代码用 COO 格式便捷地构造指定非零元数量的稀疏矩阵，同步提供稠密与 CSR 双版本，让后续基准测试能在相同数据分布下对比不同存储格式的性能。

86.7.5 主流程：参数解析、估计器构造、计时循环与格式化输出

源码路径：benchmarks/bench_random_projections.py - __main__（89-174 行）

op.add_option(  # ① --n-times：重复次数
    "--n-times",
    dest="n_times",
    default=5,
    type=int,
    help="Benchmark results are average over n_times experiments",
)

op.add_option(  # ② --n-components 与 --density 支持 'auto'
    "--n-components",
    dest="n_components",
    default="auto",
    help="Size of the random subspace. ('auto' or int > 0)",
)

# 第 86 章 —— ...

(opts, args) = op.parse_args()
opts.n_components = type_auto_or_int(opts.n_components)  # ③ 转换
opts.density = type_auto_or_float(opts.density)
selected_transformers = opts.selected_transformers.split(",")  # ④ 字符串转列表

n_nonzeros = int(opts.ratio_nonzeros * opts.n_features)  # ⑤ 计算非零元数

if opts.n_components == "auto":  # ⑥ JL 引理自动维度
    print(
        "n_components \t= %s (auto)"
        % johnson_lindenstrauss_min_dim(n_samples=opts.n_samples, eps=opts.eps)
    )

transformers = {}  # ⑦ 估计器字典

gaussian_matrix_params = {  # ⑧ Gaussian 投影参数
    "n_components": opts.n_components,
    "random_state": opts.random_seed,
}
transformers["GaussianRandomProjection"] = GaussianRandomProjection(**gaussian_matrix_params)

sparse_matrix_params = {  # ⑨ Sparse 投影参数（额外 density 与 eps）
    "n_components": opts.n_components,
    "random_state": opts.random_seed,
    "density": opts.density,
    "eps": opts.eps,
}
transformers["SparseRandomProjection"] = SparseRandomProjection(**sparse_matrix_params)

X_dense, X_sparse = make_sparse_random_data(  # ⑩ 生成数据对
    opts.n_samples, opts.n_features, n_nonzeros, random_state=opts.random_seed
)
X = X_dense if opts.dense else X_sparse  # ⑪ 根据 opts.dense 切换

for name in selected_transformers:  # ⑫ 遍历选中的投影器
    print("Perform benchmarks for %s..." % name)

    for iteration in range(opts.n_times):  # ⑬ n_times 次重复
        print("\titer %s..." % iteration, end="")
        time_to_fit, time_to_transform = bench_scikit_transformer(
            X_dense, transformers[name]  # ⑭ 始终传入稠密数据
        )
        time_fit[name].append(time_to_fit)
        time_transform[name].append(time_to_transform)
        print("done")

for name in sorted(selected_transformers):
    print_row(name, np.mean(time_fit[name]), np.mean(time_transform[name]))  # ⑮ 排序输出

逐行解析这段代码：

第 ① 行使用 optparse 定义 --n-times 选项，默认 5 次重复。

第 ② 行 --n-components 与 --density 都是字符串类型，因为它们可能取值 'auto'。

第 ③ 行在 parse_args 后通过 type_auto_or_int / type_auto_or_float 把字符串转为合适的 Python 类型。

第 ④ 行用逗号分割 --transformers 字符串为列表。

第 ⑤ 行根据 ratio_nonzeros 与 n_features 计算非零元总数。

第 ⑥ 行当 n_components == "auto" 时，调用 johnson_lindenstrauss_min_dim 计算 JL 引理给出的理论最小维度。

第 ⑦ 行初始化估计器字典。

第 ⑧ 行 Gaussian 投影只需 n_components 与 random_state。

第 ⑨ 行 Sparse 投影额外需要 density（非零比例）与 eps（JL 精度参数）。

第 ⑩ 行生成稠密与稀疏数据对。

第 ⑪ 行根据命令行 --dense 标志选择输入类型。注意只有 make_sparse_random_data 会切换数据格式，而 bench_scikit_transformer 总是接收稠密数据——这是因为我们想测试算法对稠密输入的性能，而不是稀疏输入处理（稀疏输入的测试由其他模块负责）。

第 ⑫ 行遍历命令行指定的所有投影器。

第 ⑬ 行每个投影器重复 n_times 次。

第 ⑭ 行调用 bench_scikit_transformer 返回二元组。

第 ⑮ 行调用 print_row 打印格式化输出，按 np.mean 取多次重复均值。

这段代码通过 optparse 配置实验参数、johnson_lindenstrauss_min_dim 给出理论下界、双投影器字典支持命令行选择、sklearn.clone 隔离状态、fit/transform 分离计时，构建了完整的"理论 + 实验"双轨评测框架。

86.7.6 模块依赖关系图

graph TD A["CLI --n-components auto"] --> B["type_auto_or_int convert"] B --> C["johnson_lindenstrauss_min_dim"] C --> D["GaussianRandomProjection"] B --> D C --> E["SparseRandomProjection"] B --> E F["make_sparse_random_data"] --> G["X_dense"] F --> H["X_sparse CSR"] G --> I["bench_scikit_transformer"] H --> I I --> J["fit timing"] I --> K["transform timing"] J --> L["compute_time delta"] K --> L["compute_time delta"] L --> M["np.mean aggregate"] M --> N["print_row format"] K --> N["print_row format"]

86.8 设计中的取舍

问：为什么不用 time.perf_counter 而用 datetime.now()？ 答：time.perf_counter() 返回单调时钟的纳秒级浮点数，确实精度更高，但 datetime.now() 返回的 datetime 对象支持 timedelta 运算，让 compute_time 函数可以复用 delta.seconds 与 delta.microseconds。这是用一点点精度换取代码可读性的取舍——在毫秒级的基准测试里，两者的差异可以忽略。

问：HashingVectorizer 为什么不监控 transform 阶段？ 答：基准测试聚焦端到端性能，因为业务上 fit_transform 几乎总是串联调用。分离 fit 与 transform 会让脚本复杂化，且 HashingVectorizer 的 fit 是空操作（无状态），无分离必要。

问：为什么 bench_feature_expansions.py 把 PolynomialFeatures 实例化放在循环外？ 答：因为我们要测量 transform 阶段的计算成本，而非构造与 fit 成本。把对象创建移到循环外是隔离测量目标的标准做法。如果把构造放进来，结果会包含约 10 毫秒的固定初始化开销，淹没维度/密度的变化信号。

问：为什么随机投影默认传入稠密数据而非稀疏数据？ 答：因为我们想比较的是 GaussianRandomProjection 与 SparseRandomProjection 两者之间的差异，而不是它们处理稀疏输入的能力。固定输入格式让两者的对比聚焦在投影矩阵生成与矩阵乘法阶段，与数据存储格式解耦。

问：为什么 bench_sample_without_replacement.py 用 optparse 而非 argparse？ 答：这是历史原因——该脚本早于 argparse 成为标准库（Python 2.7 时代）。两者功能等价，但 optparse 的 add_option 接口更接近传统 Unix 工具风格（如 getopt），对维护者而言迁移成本不高。

86.9 动手练习

1. 对比向量化器内存占用模式

   • 阅读 bench_text_vectorizers.py 中 memory_usage 的调用方式，回答：

     1. 为什么 HashingVectorizer 在大 n-gram 配置下内存优势最明显？

     2. 若将 n_repeat 从 3 增到 10，time 的标准差会如何变化？原因是什么？

     3. 如何修改脚本以同时统计 fit 与 transform 的分阶段耗时？

2. 分析多项式特征扩展的稀疏优势消失点

   • 阅读 bench_feature_expansions.py 的绘图逻辑，回答：

     1. 在 density=1.0 时，CSR 与 Dense 曲线为何几乎重合？

     2. 若将 degree 改为 3，稀疏矩阵的非零元增长规律会如何变化？对计时曲线有何影响？

     3. 脚本中为何在每个 trial 开始前不重新实例化 PolynomialFeatures？是否会影响结果？

3. 剖析无放回采样四种 sklearn 策略的适用边界

   • 阅读 bench_sample_without_replacement.py 中 sampling_algorithm 字典构建，回答：

     1. tracking_selection 与 reservoir_sampling 在 ratio → 1 时性能差异的根本原因？

     2. method='auto' 的内部派发逻辑依据是什么？（提示：查看 sklearn 源码 sklearn/utils/_random.pyx）

     3. 为何 numpy-permutation 在小 ratio 时极慢？其时间复杂度是多少？

4. 随机投影基准中 fit 与 transform 的性能分离

   • 阅读 bench_random_projections.py 的 bench_scikit_transformer 与主流程，回答：

     1. GaussianRandomProjection.fit 的主要开销在哪里？为何随 n_components 线性增长？

     2. SparseRandomProjection 的 density 参数如何影响 transform 阶段的稀疏矩阵乘法速度？

     3. 若在主循环外预先 fit 好投影器，再在循环内仅 transform，结果会有何不同？这对应什么实际场景？

86.10 本章小结

这一章中我们学习了特征工程与采样领域的独立基准脚本。首先，我们了解了文本向量化器基准如何通过 run_vectorizer 闭包实现"状态隔离的端到端计时"，并用 itertools.product 构建 3×4 笛卡尔积实验网格；其次，我们剖析了 bench_feature_expansions.py 如何通过 3 trials × 3 densities × 6 dimensionalities 的三重循环构建稀疏优势消失的完整曲面；接着，我们深入了 bench_sample_without_replacement.py 中 6 种算法（4 种 sklearn 策略 + Python 内置 + NumPy permutation）的 lambda 适配层与 gc.collect 严谨计时模式；最后，我们见证了 bench_random_projections.py 如何通过 JL 引理给出最小投影维度的理论下界，并分离 fit（投影矩阵生成）与 transform（矩阵乘法）两阶段进行精细化对比。

下表汇总了本章涉及的核心概念：

| 概念 | 解释 |

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

| bench_text_vectorizers.py | 三足鼎立：Count/Tfidf/Hashing × word/char n-gram，timeit+memory_profiler 双指标，DataFrame 多级索引聚合 |

| run_vectorizer 闭包 | 每次实例化新向量化器，隔离状态，聚焦端到端 fit_transform 性能 |

| bench_feature_expansions.py | 维度爆炸实验：稀疏 CSR vs 稠密 Dense × 6 维度 × 3 密度，PolynomialFeatures(degree=2) 计时曲面 |

| bench_sample_without_replacement.py | 公平抽签对决：6 种算法（Python/NumPy/sklearn 4 策略）× 采样比例，optparse 配置，gc.collect 严谨计时 |

| sample_without_replacement 4 策略 | auto/ tracking_selection(集合标记)/ reservoir_sampling(水池)/ pool(全排列切片) |

| bench_random_projections.py | JL 引理免训练降维：Gaussian(密集高斯) vs Sparse(稀疏 Achlioptas) × fit/transform 分离计时 × auto 维度计算 |

| johnson_lindenstrauss_min_dim | 给定样本数与 eps，理论计算最小投影维度，保证距离近似保持 |

| make_sparse_random_data | COO 格式构造指定非零元数的稀疏矩阵，同步返回稠密/CSR 对照组 |

| bench_scikit_transformer | clone 深拷贝估计器隔离状态，分别计时 fit(含投影矩阵生成) 与 transform(纯矩阵乘法) |

| type_auto_or_float/int | 命令行参数解析工具，支持 'auto' 关键字与数值的双模式自动转换 |

| compute_time | 微秒级计时工具函数，将 datetime.timedelta 转换为浮点秒数 |

| print_row | 格式化输出工具，对齐打印 Transformer 名称、fit 耗时与 transform 耗时 |

感谢你读到了这里。基准测试的精髓不在于追求极致性能，而在于用严谨的方法揭示不同算法在不同条件下的性能边界，从而为工程决策提供数据支撑。

第 87 章 —— 异常检测基准脚本 —— 寻找"离群之马"的效率

87.1 学习目标

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

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

• 理解异常检测基准测试的设计模式与评估指标

• 掌握 IsolationForest、LocalOutlierFactor、IsotonicRegression 在基准测试中的配置与性能分析方法

• 了解合成数据生成器如何探测算法的最佳/最差时间复杂度场景

• 熟悉 argparse 子命令模式实现基准运行与结果可视化的分离

• 学会使用 joblib.parallel_config 控制预测阶段的并行度以分离训练/推理性能

87.2 生活类比

想象异常检测基准测试是一场侦探选拔赛，其中六个数据集就是六个犯罪现场：KDD 网络入侵现场藏匿着恶意连接，森林覆盖现场混杂着珍稀物种与异常地块，航天飞机与航天飞船的遥测数据中潜伏着机械隐患，而合成的高斯簇与均匀噪声则像训练场里的标准关卡。各位侦探的办案手法各具特色，其中 IsolationForest 像一个随机砍伐树林的猎人在特征空间里随机切出隔离墙，正常的树木深藏在密林深处需要很多刀才能砍倒，而异常点就像孤立在外的灌木几刀就能孤立出来。LocalOutlierFactor 则像是社区巡逻队长，他会逐一审视每个居民的邻居密度，那些明显比邻居孤僻的个体立刻会被标记为异类。IsotonicRegression 则是位单调校准师，他不找异常，而是把混乱的预测分数重新排列成一条单调上升的曲线，让置信度随着分数的提升严格递增。为了让选拔更加公平，合成数据生成器就是模拟训练场的造物主。普通的对数扰动数据像日常出警，逻辑回归型数据像分类概率校准演练，而那精心设计的病理数据则像一个充满陷阱的迷宫，能让老练的 PAVA 算法也陷入 O(n²) 的泥潭。参数扫描就是让侦探在不同团队规模（n_jobs）、案件量（n_samples）、线索维度（n_features）下比拼办案速度，毕竟神探也不能在海量线索中卡壳。最后的评估环节像极了破案率与误报率的权衡：ROC 曲线下方面积越大，越说明这位侦探既能抓住真凶又不冤枉好人。而那 1000 次空循环的缓存清理，则像出警前清理警械箱，确保每次计时都从干净状态起步。

87.3 源码地图

benchmarks/bench_isolation_forest.py
├── print_outlier_ratio()           # 打印目标分布与异常比例
├── __main__                       # 核心流程：数据获取→预处理→训练→ROC评估
│   ├── fetch_kddcup99/fetch_openml/fetch_covtype  # 多数据源加载
│   ├── LabelBinarizer              # 分类特征独热编码
│   ├── IsolationForest(n_jobs=-1)  # 全核并行训练
│   ├── decision_function + roc_curve/auc  # 异常评分与ROC计算
│   └── matplotlib 可视化对比       # 多条ROC曲线叠加展示

benchmarks/bench_isolation_forest_predict.py
├── get_data()                     # 合成高斯簇+均匀噪声数据生成器
├── bench()                        # 基准运行：四维参数网格搜索
│   ├── 固定训练集/变化测试集/特征/污染率/n_jobs
│   ├── 1000次空循环清理缓存        # 减少测量噪声
│   ├── parallel_config控制预测并行度 # 分离训练/推理性能
│   └── defaultdict累积结果→DataFrame保存CSV
├── plot()                         # 可视化对比：seaborn双分支折线图
│   ├── 读取两分支CSV合并对比
│   └── n_samples_test vs predict_time 按n_jobs分组
└── argparse子命令分派             # bench/plot 模式分离

benchmarks/bench_lof.py
├── __main__                       # 全量数据转导式评估流程
│   ├── 复用相同数据集加载/预处理逻辑
│   ├── LocalOutlierFactor(n_neighbors=20)  # 就地异常检测
│   ├── negative_outlier_factor_ 属性评分  # 无需预测阶段
│   └── 训练集上直接计算ROC AUC      # 体现转导式学习特性

benchmarks/bench_isotonic.py
├── generate_perturbed_logarithm_dataset()  # 对数趋势+噪声
├── generate_logistic_dataset()             # 逻辑回归型单调关系
├── generate_pathological_dataset()         # 触发O(n²)病理数据
├── bench_isotonic_regression()             # 单次迭代高精度计时
│   ├── gc.collect() 消除GC抖动
│   └── timeit.default_timer 纳秒级计时
└── __main__                        # 指数级规模扫描+可选对数坐标绘图

87.4 隔离森林基准 —— ROC 曲线背后的"狩猎游戏"

这一节我们先剖析隔离森林基准，它是隔离森林在六个经典异常检测数据集上的"选秀现场"。我们将沿着"数据获取 → 预处理 → 训练 → ROC 评估"的链条，逐行解读源码。

在隔离森林里，"隔离"是核心隐喻：异常点天然具有"易于被随机切分孤立"的特性，而正常点则密集地抱团在一起，需要更多次的随机切分才能将它们彼此分开。隔离森林正是基于这一直觉构建随机树集成。在基准测试中，我们需要验证它在多个真实数据集上的异常检测能力，并用 ROC 曲线直观呈现。

87.4.1 辅助函数与主流程概览

下图展示了隔离森林基准的整体执行流程，从数据加载、预处理、模型训练到 ROC 曲线绘制的完整链路：

flowchart TB subgraph 数据加载["数据加载阶段"] A["fetch_kddcup99<br/>fetch_openml<br/>fetch_covtype"] --> B["数据集子集选择<br/>http/smtp/SA/SF/shuttle/forestcover"] end subgraph 预处理["特征工程阶段"] C["LabelBinarizer<br/>独热编码分类特征"] --> D["数据类型转换<br/>X.astype float"] D --> E["数据切分<br/>50%训练 / 50%测试"] end subgraph 训练评估["训练与评估阶段"] F["IsolationForest<br/>n_jobs=-1, random_state=1"] --> G["model.fit X_train<br/>训练阶段计时"] G --> H["-decision_function X_test<br/>异常评分"] H --> I["roc_curve y_test, scoring<br/>计算FPR/TPR"] I --> J["auc fpr, tpr<br/>计算AUC面积"] end subgraph 可视化["可视化阶段"] J --> K["多条ROC曲线叠加<br/>ax_roc.plot"] K --> L["plt.show<br/>展示结果"] end A --> C B --> C E --> F

先看辅助函数 print_outlier_ratio，它帮助我们一眼了解目标分布：

源码路径：benchmarks/bench_isolation_forest.py - print_outlier_ratio()（30-42行）

def print_outlier_ratio(y):
    """
    Helper function to show the distinct value count of element in the target.
    Useful indicator for the datasets used in bench_isolation_forest.py.
    """
    # 使用np.unique统计目标数组y中所有唯一值及其出现次数
    # 返回的uniq是唯一值数组，cnt是对应每个唯一值的计数数组
    uniq, cnt = np.unique(y, return_counts=True)
    print("----- Target count values: ")  # 打印表头分隔符
    # 遍历每个唯一标签及其计数，逐行打印分布情况
    for u, c in zip(uniq, cnt):
        print("------ %s -> %d occurrences" % (str(u), c))
    # 计算异常样本占比：取计数数组中的最小值（即最少类别的样本数）除以总样本数
    # 这个比例反映了数据集中异常样本的稀缺程度，是评估异常检测难度的关键指标
    print("----- Outlier ratio: %.5f" % (np.min(cnt) / len(y)))

这段代码定义了 print_outlier_ratio 函数，专门打印目标变量的分布情况。它通过 np.unique 找出所有唯一标签并计数，然后逐行打印，最后用最小类别的计数除以总样本数得到异常比例。这个比例是评估异常检测难度的关键指标——若比例极低（如 smtp 数据集），随机切分测试集可能完全不含异常，此时 roc_curve 会发出警告。

接下来看主流程的核心循环：

源码路径：benchmarks/bench_isolation_forest.py - __main__（44-160行）

# 第 87 章 —— 设置全局随机种子为1，确保数据打乱和模型训练的可复现性
# 第 87 章 —— 这个值在后续所有使用random_state参数的地方都会被引用
random_state = 1

# 第 87 章 —— 创建ROC曲线绘制画布，尺寸为8x5英寸，用于后续叠加显示六个数据集的ROC曲线
fig_roc, ax_roc = plt.subplots(1, 1, figsize=(8, 5))

# 第 87 章 —— 设置是否绘制决策函数直方图的开关，默认关闭以减少输出噪音
with_decision_function_histograms = False

# 第 87 章 —— 定义六个候选数据集名称，对应KDD Cup 99的子集和OpenML/UCSD数据集
datasets = ["http", "smtp", "SA", "SF", "shuttle", "forestcover"]

# 第 87 章 —— 遍历每个数据集进行完整的训练和评估流程
for dat in datasets:
    print("====== %s ======" % dat)  # 打印当前处理的数据集名称作为分隔标识
    print("--- Fetching data...")

    # 对于KDD Cup 99的四个子集（http/smtp/SA/SF），使用subset参数指定协议类型
    # percent10=True使用10%抽样版本以加速下载和处理
    # random_state=random_state确保每次运行都获得相同的数据子集
    if dat in ["http", "smtp", "SF", "SA"]:
        dataset = fetch_kddcup99(
            subset=dat, shuffle=True, percent10=True, random_state=random_state
        )
        X = dataset.data  # 特征矩阵
        y = dataset.target  # 目标标签

    # shuttle数据集从OpenML下载，需要转换为整数类型标签
    if dat == "shuttle":
        dataset = fetch_openml("shuttle", as_frame=False)  # 从OpenML获取shuttle数据
        X = dataset.data
        y = dataset.target.astype(np.int64)  # 转换为64位整数
        X, y = sh(X, y, random_state=random_state)  # 打乱数据顺序以消除原始顺序偏差
        s = y != 4  # 创建布尔掩码，标记非类别4的样本
        X = X[s, :]  # 移除类别4的样本
        y = y[s]
        y = (y != 1).astype(int)  # 将类别1设为正常（0），其余设为异常（1）
        print("----- ")

    # forestcover数据集包含森林覆盖类型，类别2为正常，类别4为异常
    if dat == "forestcover":
        dataset = fetch_covtype(shuffle=True, random_state=random_state)  # 下载森林覆盖数据
        X = dataset.data
        y = dataset.target
        s = (y == 2) + (y == 4)  # 创建布尔掩码，只保留类别2和4
        X = X[s, :]  # 应用掩码过滤样本
        y = y[s]
        y = (y != 2).astype(int)  # 类别2为正常（0），类别4为异常（1）
        print_outlier_ratio(y)  # 打印过滤后的异常比例

    print("--- Vectorizing data...")

    # SF数据集的第二列（服务类型）为字符串，需要独热编码处理
    if dat == "SF":
        lb = LabelBinarizer()  # 创建独热编码器
        x1 = lb.fit_transform(X[:, 1].astype(str))  # 对服务类型列进行拟合和转换
        X = np.c_[X[:, :1], x1, X[:, 2:]]  # 替换原列：保留第一列，插入编码后的列，追加其余列
        y = (y != b"normal.").astype(int)  # 将字节字符串标签转换为0/1整数
        print_outlier_ratio(y)

    # SA数据集有三列字符串特征需要编码：协议类型、服务类型、标志位
    if dat == "SA":
        lb = LabelBinarizer()
        x1 = lb.fit_transform(X[:, 1].astype(str))  # 编码协议类型
        x2 = lb.fit_transform(X[:, 2].astype(str))  # 编码服务类型
        x3 = lb.fit_transform(X[:, 3].astype(str))  # 编码标志位
        # 将三列编码后的特征与原始数值特征拼接，形成完整的特征矩阵
        X = np.c_[X[:, :1], x1, x2, x3, X[:, 4:]]
        y = (y != b"normal.").astype(int)
        print_outlier_ratio(y)

    # http和smtp数据集的标签处理相对简单，直接转换即可
    if dat in ("http", "smtp"):
        y = (y != b"normal.").astype(int)  # 将"normal."标签转为0，其余转为1
        print_outlier_ratio(y)

    # 获取数据集规模信息，用于后续的数据切分
    n_samples, n_features = X.shape
    n_samples_train = n_samples // 2  # 按50%比例切分训练集

    X = X.astype(float)  # 统一转换为浮点类型，便于后续数值计算
    X_train = X[:n_samples_train, :]  # 前一半作为训练集
    X_test = X[n_samples_train:, :]  # 后一半作为测试集
    y_train = y[:n_samples_train]
    y_test = y[n_samples_train:]

这段代码定义了主流程中的数据加载与预处理阶段。六个数据集各有定制化处理：KDD Cup 99 子集直接通过 fetch_kddcup99 加载；shuttle 数据集需要移除类别 4 并将类别 1 作为正常标签；forestcover 则精选类别 2 与 4 并以类别 2 为正常。每个分支都使用 random_state=1 控制随机性以保证可复现。LabelBinarizer 将字符串型分类特征转为 0/1 向量，避免模型将离散值误读为有序关系。数据按时间或顺序等量切分（前一半训练、后一半测试），这是一种简单但能避免信息泄露的切分方式。

    print("--- Fitting the IsolationForest estimator...")
    # 创建隔离森林模型实例，使用全部CPU核心进行训练
    # random_state=random_state确保每次运行产生相同的随机树结构
    model = IsolationForest(n_jobs=-1, random_state=random_state)
    tstart = time()  # 记录训练开始时间
    model.fit(X_train)  # 在训练集上拟合模型
    fit_time = time() - tstart  # 计算训练耗时（秒）
    tstart = time()  # 记录预测开始时间

    # decision_function返回"正常程度"分数，越正表示越正常
    # 取负后变为"异常程度"，越负表示越异常，符合roc_curve的要求（分数越高越可能为正类）
    scoring = -model.decision_function(X_test)

    print("--- Preparing the plot elements...")
    # 当开关打开时，绘制决策函数在各数据集上的分布直方图
    if with_decision_function_histograms:
        fig, ax = plt.subplots(3, sharex=True, sharey=True)  # 创建3行子图，共享坐标轴
        bins = np.linspace(-0.5, 0.5, 200)  # 定义直方图的区间边界
        ax[0].hist(scoring, bins, color="black")  # 绘制全部样本的直方图
        ax[0].set_title("Decision function for %s dataset" % dat)
        ax[1].hist(scoring[y_test == 0], bins, color="b", label="normal data")  # 正常样本分布
        ax[1].legend(loc="lower right")
        ax[2].hist(scoring[y_test == 1], bins, color="r", label="outliers")  # 异常样本分布
        ax[2].legend(loc="lower right")

    # 计算ROC曲线和AUC值
    predict_time = time() - tstart  # 计算预测耗时
    fpr, tpr, thresholds = roc_curve(y_test, scoring)  # 计算假阳性率、真阳性率和阈值
    auc_score = auc(fpr, tpr)  # 计算ROC曲线下面积（AUC）
    # 格式化标签字符串，包含数据集名、AUC值、训练时间和预测时间
    label = "%s (AUC: %0.3f, train_time= %0.2fs, test_time= %0.2fs)" % (
        dat,
        auc_score,
        fit_time,
        predict_time,
    )
    print(label)  # 在控制台输出性能指标
    ax_roc.plot(fpr, tpr, lw=1, label=label)  # 将当前数据集的ROC曲线叠加到主图

# 第 87 章 —— 设置坐标轴范围和标签，完成图表装饰
ax_roc.set_xlim([-0.05, 1.05])
ax_roc.set_ylim([-0.05, 1.05])
ax_roc.set_xlabel("False Positive Rate")  # X轴标签：假阳性率
ax_roc.set_ylabel("True Positive Rate")  # Y轴标签：真阳性率
ax_roc.set_title("Receiver operating characteristic (ROC) curves")  # 图表标题
ax_roc.legend(loc="lower right")  # 图例放在右下角
fig_roc.tight_layout()  # 调整布局以减少留白
plt.show()  # 显示最终图表

这段代码定义了主流程中的训练、评估与可视化阶段。注意 n_jobs=-1 让隔离森林训练阶段使用全部 CPU 核；预测阶段虽然也调用了 decision_function，但因为是单进程 Python 代码，n_jobs 不生效。scoring = -model.decision_function 是关键技巧：隔离森林的 decision_function 返回的是"正常程度"，越正越正常，取负后变为"异常程度"，越负越异常，符合 roc_curve 要求的"得分越高越可能为正类"。

87.4.2 隔离森林基准的设计取舍

隔离森林基准采用训练/测试分离的评估范式，主要基于以下考量：作为归纳式学习算法，隔离森林构建的隔离树结构可以独立应用于新数据，因此可以通过分离训练集与测试集来衡量模型的泛化能力。这种设计模拟了真实场景中的部署需求——模型在已知数据上训练后，需要对未知数据进行预测。此外，50/50 的固定切分比例简化了实验设计，避免了交叉验证带来的复杂度，同时确保每个数据集都有足够大的测试集来可靠地估计 ROC 曲线。

87.5 隔离森林预测专项 —— 参数扫描下的"推理体检"

如果说上一节关注的是隔离森林的"狩猎准确率"，那么这一节就专注于它的"办案速度"。预测基准脚本的核心思想是分离训练与预测的并行度：训练时固定用全部 CPU 核确保模型快速构建，预测时遍历不同 n_jobs 值来量化并行推理的加速比。

87.5.1 合成数据生成与基准执行

下图展示了隔离森林预测基准的整体架构，包括数据生成、四维参数扫描和并行控制的核心流程：

flowchart TB subgraph 数据生成["数据生成模块 get_data"] A["n_samples_train=1000"] --> B["生成两个高斯簇<br/>X+2 和 X-2"] C["n_samples_test=1k/10k/50k"] --> D["同样构造两个高斯簇"] E["contamination=0.01/0.1/0.5"] --> F["均匀分布噪声作为异常点<br/>替换部分测试样本"] end subgraph 参数扫描["四维参数网格扫描"] G["n_features=10/100/1000"] --> H["IsolationForest训练<br/>n_jobs=-1全核并行"] H --> I["1000次空循环<br/>清理CPU缓存预热"] I --> J["parallel_config<br/>threading, n_jobs=1/2/3/4"] J --> K["decision_function计时<br/>记录predict_time"] end subgraph 结果输出["结果汇总"] K --> L["defaultdict累积结果"] L --> M["DataFrame.to_csv<br/>保存为CSV文件"] end subgraph 可视化对比["plot函数跨分支对比"] M --> N["读取PR分支CSV"] M --> O["读取main分支CSV"] N --> P["pd.merge合并<br/>按n_samples_test和n_jobs"] O --> P P --> Q["seaborn双子图<br/>左右并排对比"] end subgraph 子命令分派["argparse子命令"] R["bench子命令<br/>运行参数扫描"] --> S["bench_results/branch参数"] T["plot子命令<br/>生成对比图"] --> U["pr_name/main_name/image_path参数"] end

先看合成数据生成器，它模拟了"两个正常簇 + 均匀分布噪声"的经典场景：

源码路径：benchmarks/bench_isolation_forest_predict.py - get_data()（40-65行）

def get_data(
    n_samples_train, n_samples_test, n_features, contamination=0.1, random_state=0
):
    """
    合成异常检测基准数据的生成函数。
    构造两个高斯簇作为正常样本，再注入均匀分布的离群点作为异常样本。
    """
    # 创建随机数生成器，使用固定种子确保每次调用生成相同的数据
    rng = np.random.RandomState(random_state)

    # 生成标准差为0.3的基础扰动数据（n_samples_train行，n_features列）
    # 这个较小的标准差确保数据点聚集在原点附近
    X = 0.3 * rng.randn(n_samples_train, n_features)
    # 将数据分为两部分，一部分中心在+2，一部分中心在-2，形成两个分离的高斯簇
    # 这模拟了真实场景中正常数据的聚类结构
    X_train = np.r_[X + 2, X - 2]

    # 同样方式生成测试集的基础扰动数据
    X = 0.3 * rng.randn(n_samples_test, n_features)
    X_test = np.r_[X + 2, X - 2]

    # 根据污染率计算离群点数量
    n_outliers = int(np.floor(contamination * n_samples_test))
    # 在[-4, 4]区间内均匀采样，生成覆盖范围比正常簇更广的离群点
    # 正常簇的中心在±2，标准差0.3，99%数据落在[1.1, 2.9]和[-2.9, -1.1]区间内
    # 而离群点在[-4, 4]均匀分布，会出现在正常簇范围之外
    X_outliers = rng.uniform(low=-4, high=4, size=(n_outliers, n_features))

    # 随机选择测试集中的部分样本索引，用离群点替换它们
    # replace=False确保每个索引只被选择一次
    outlier_idx = rng.choice(np.arange(0, n_samples_test), n_outliers, replace=False)
    X_test[outlier_idx, :] = X_outliers  # 用离群点覆盖选中的正常样本

    return X_train, X_test  # 返回训练集和测试集

这段代码定义了 get_data 函数，构造了典型的双簇异常检测场景。两个高斯簇（中心分别为 +2 和 -2）作为正常数据，均匀分布的离群点作为异常值，且异常点散布在比正常簇更广的 [-4, 4] 区间内。这种设计的妙处在于：正常数据有清晰的聚类结构，异常数据则散布在边界以外，让隔离森林的"随机切分隔离"机制能高效工作。

接下来看 bench 函数的四维参数网格搜索：

源码路径：benchmarks/bench_isolation_forest_predict.py - bench()（120-165行）

def bench(args):
    """
    执行隔离森林预测阶段的基准测试。
    通过四维参数网格搜索，量化不同n_jobs配置下的预测性能。
    """
    results_dir = Path(args.bench_results)  # 结果保存目录
    branch = args.branch  # 分支名称，用于CSV文件命名
    random_state = 1  # 固定随机种子确保可复现

    # 使用defaultdict自动初始化列表，避免手动检查键是否存在
    results = defaultdict(list)

    n_samples_train = 1000  # 训练集固定为1000样本

    # 外层循环：遍历测试集规模（1k/10k/50k），评估算法在不同数据量下的扩展性
    for n_samples_test in [1000, 10000, 50000]:
        # 中层循环：遍历特征维度（10/100/1000），评估高维数据的计算开销
        for n_features in [10, 100, 1000]:
            # 内层循环：遍历异常比例（1%/10%/50%），评估污染率对性能的影响
            for contamination in [0.01, 0.1, 0.5]:
                # 最内层循环：遍历并行作业数（1/2/3/4），量化并行推理的加速比
                for n_jobs in [1, 2, 3, 4]:
                    # 根据当前参数组合生成对应的训练集和测试集
                    X_train, X_test = get_data(
                        n_samples_train, n_samples_test, n_features, contamination, random_state
                    )

                    print("--- Fitting the IsolationForest estimator...")
                    # 训练阶段固定使用全部CPU核心（n_jobs=-1）
                    # 这是因为训练是O(n_estimators × log(n))的并行友好任务
                    model = IsolationForest(n_jobs=-1, random_state=random_state)
                    tstart = time()  # 记录训练开始时间
                    model.fit(X_train)  # 在训练集上拟合模型
                    fit_time = time() - tstart  # 计算训练耗时

                    # 缓存预热：执行1000次空循环，让CPU缓存、分支预测器等进入稳定状态
                    # 这样可以避免首次预测因冷启动而显著慢于稳态测量值
                    for _ in range(1000):
                        1 + 1  # 纯计算操作，不产生实际效果

                    # 使用joblib的parallel_config上下文管理器控制预测阶段的并行度
                    # "threading"指定使用线程级并行，n_jobs参数动态调整线程数
                    # 这样可以在不重建模型的情况下测试不同并行配置的性能
                    with parallel_config("threading", n_jobs=n_jobs):
                        tstart = time()  # 记录预测开始时间
                        model.decision_function(X_test)  # 执行异常评分预测
                        predict_time = time() - tstart  # 计算预测耗时

                    # 将所有结果累积到defaultdict中，便于后续转换为DataFrame
                    results["predict_time"].append(predict_time)
                    results["fit_time"].append(fit_time)
                    results["n_samples_train"].append(n_samples_train)
                    results["n_samples_test"].append(n_samples_test)
                    results["n_features"].append(n_features)
                    results["contamination"].append(contamination)
                    results["n_jobs"].append(n_jobs)

    # 将累积的结果转换为DataFrame，并保存为CSV文件
    df = pd.DataFrame(results)
    df.to_csv(results_dir / f"{branch}.csv", index=False)

这段代码定义了 bench 函数。它构建了 3×3×3×4 = 108 种参数组合的全因子实验。关键工程细节有三处：

1. 训练固定 n_jobs=-1 而预测遍历 n_jobs：训练阶段是 O(n_estimators × log(n)) 的并行构建，全核最快；预测阶段是 O(n_samples × tree_depth) 的并行推理，可控变量便于量化加速比。

2. for _ in range(1000): 1 + 1 空循环：这是隔离基准测试中的"缓存预热"技巧，让 CPU 缓存、分支预测器等恢复到稳定状态，避免首次预测因冷启动而显著慢于稳态。

3. parallel_config("threading", n_jobs=n_jobs)：隔离森林预测底层使用 Cython + OpenMP，可通过 joblib 的 parallel_config 在不重新构建模型的情况下动态调整线程数。

最后看 plot 函数的双分支对比绘图：

源码路径：benchmarks/bench_isolation_forest_predict.py - plot()（67-118行）

def plot(args):
    """
    绘制跨分支的性能对比图。
    读取PR分支和main分支的CSV结果，以n_samples_test为X轴绘制双分支折线图。
    """
    import matplotlib.pyplot as plt
    import seaborn as sns

    bench_results = Path(args.bench_results)  # 结果目录路径
    pr_name = args.pr_name  # PR分支名称
    main_name = args.main_name  # main分支名称
    image_path = args.image_path  # 图片保存路径

    results_path = Path(bench_results)  # 解析为Path对象
    pr_path = results_path / f"{pr_name}.csv"  # PR分支CSV路径
    main_path = results_path / f"{main_name}.csv"  # main分支CSV路径
    image_path = results_path / image_path  # 图片保存完整路径

    # 读取两个CSV文件，并为每条记录添加分支标识列
    df_pr = pd.read_csv(pr_path).assign(branch=pr_name)
    df_main = pd.read_csv(main_path).assign(branch=main_name)

    # 使用pd.merge按共同列合并两个DataFrame
    # suffixes参数用于处理重名列，添加后缀区分来源
    merged_data = pd.merge(
        df_pr,
        df_main,
        on=["n_samples_test", "n_jobs"],  # 以测试样本数和并行作业数为键进行合并
        suffixes=("_pr", "_main"),  # PR分支的列加_pr后缀，main分支加_main后缀
    )

    # 设置seaborn的绘图风格
    sns.set(style="whitegrid", context="notebook", font_scale=1.5)

    # 创建双子图布局，共享X轴和Y轴
    fig, axes = plt.subplots(1, 2, figsize=(18, 6), sharex=True, sharey=True)

    # 调试输出：打印n_jobs的唯一值
    print(merged_data["n_jobs"].unique())

    # 左图：PR分支的预测时间
    ax = axes[0]
    # 使用seaborn的lineplot绘制折线图，按n_jobs分组着色
    sns.lineplot(
        data=merged_data,
        x="n_samples_test",  # X轴：测试样本数
        y="predict_time_pr",  # Y轴：PR分支的预测时间
        hue="n_jobs",  # 按并行作业数分组（不同颜色）
        style="n_jobs",  # 按并行作业数区分线型
        markers="o",  # 使用圆形标记
        ax=ax,
        legend="full",  # 显示完整图例
    )
    ax.set_title(f"Predict Time vs. n_samples_test - {pr_name} branch")  # 设置子图标题
    ax.set_ylabel("Predict Time (Seconds)")  # Y轴标签
    ax.set_xlabel("n_samples_test")  # X轴标签

    # 右图：main分支的预测时间
    ax = axes[1]
    sns.lineplot(
        data=merged_data,
        x="n_samples_test",
        y="predict_time_main",
        hue="n_jobs",
        style="n_jobs",
        markers="X",  # 使用X形标记
        dashes=True,  # 使用虚线
        ax=ax,
        legend=None,  # 右图不显示图例（避免重复）
    )
    ax.set_title(f"Predict Time vs. n_samples_test - {main_name} branch")
    ax.set_ylabel("Predict Time")
    ax.set_xlabel("n_samples_test")

    # 调整布局并保存图片
    plt.tight_layout()
    fig.savefig(image_path, bbox_inches="tight")  # bbox_inches="tight"减少多余留白
    print(f"Saved image to {image_path}")

这段代码定义了 plot 函数，专用于跨分支对比。它读取两个 CSV（PR 与 main）并以 n_samples_test 与 n_jobs 为键合并，得到每个组合在两个分支上的预测时间。通过左右双子图的并排呈现，seaborn 自动按 n_jobs 分组上色，让读者一眼看出"测试样本量增大时，加速比如何变化"以及"两个分支是否有性能回归"。

整个脚本的入口使用 argparse 子命令模式：

源码路径：benchmarks/bench_isolation_forest_predict.py - __main__（167-195行）

if __name__ == "__main__":
    # 创建顶层参数解析器
    parser = argparse.ArgumentParser()

    # 创建子命令解析器，允许同一个脚本处理不同的子命令
    subparsers = parser.add_subparsers()

    # 定义bench子命令：用于运行参数扫描基准测试
    bench_parser = subparsers.add_parser("bench")
    bench_parser.add_argument("bench_results")  # 结果保存目录（位置参数）
    bench_parser.add_argument("branch")  # 分支名称（位置参数）
    bench_parser.set_defaults(func=bench)  # 将bench函数绑定为bench子命令的处理器

    # 定义plot子命令：用于绘制跨分支对比图
    plot_parser = subparsers.add_parser("plot")
    plot_parser.add_argument("bench_results")  # 结果目录
    plot_parser.add_argument("pr_name")  # PR分支名称
    plot_parser.add_argument("main_name")  # main分支名称
    plot_parser.add_argument("image_path")  # 图片保存路径
    plot_parser.set_defaults(func=plot)  # 将plot函数绑定为plot子命令的处理器

    # 解析命令行参数并调用对应的处理函数
    args = parser.parse_args()
    args.func(args)  # 根据子命令类型分派执行bench或plot函数

这段代码定义了入口点的子命令分派逻辑。bench 与 plot 两条子命令相互独立：前者负责采集数据生成 CSV，后者负责对比 CSV 生成图片。这种解耦设计是 CI/CD 性能回归检测的关键：先在 PR 分支运行 bench，再在 main 分支运行 bench，最后用 plot 把两次结果拼成一目了然的对比图。

87.5.2 隔离森林预测基准的设计取舍

在隔离森林预测基准中，训练阶段固定使用 n_jobs=-1 而预测阶段遍历 n_jobs 的设计背后有明确的工程考量。训练阶段涉及构建 N 棵随机树，每棵树的构建过程相互独立，非常适合并行化，因此使用全部 CPU 核心可以让模型快速构建，减少实验启动的开销。相比之下，预测阶段是批量推理任务，可控变量正是评估"并行推理加速比"的关键。通过 parallel_config("threading", n_jobs=n_jobs) 上下文管理器，可以在不重建模型的情况下动态切换线程数，这种设计让性能对比更加纯粹——每次测量只受单一因素影响。

87.6 LOF 基准 —— 局部离群因子的"异常雷达"

与隔离森林的"训练/测试分离"不同，LOF（局部离群因子）基准采用了全量数据转导式评估的范式。这是因为 LOF 的核心思想是"局部密度对比"——每个点的异常得分都依赖于全体数据点的邻域关系，不存在一个"独立于数据"的可复用模型。

87.6.1 转导式评估流程概览

下图展示了 LOF 基准的转导式评估流程，与隔离森林的归纳式评估形成鲜明对比：

flowchart LR subgraph 数据加载["数据加载"] A["fetch_kddcup99<br/>fetch_openml<br/>fetch_covtype"] --> B["数据预处理<br/>独热编码+类型转换"] end subgraph 评估["转导式评估（无测试集）"] B --> C["LocalOutlierFactor<br/>n_neighbors=20"] C --> D["model.fit X<br/>在全量数据上训练"] D --> E["negative_outlier_factor_<br/>直接获取异常分数"] end subgraph ROC计算["ROC评估"] E --> F["-negative_outlier_factor_<br/>取负值使越小越异常"] F --> G["roc_curve y, scoring<br/>在全量标签上计算"] G --> H["auc fpr, tpr<br/>得到AUC分数"] end subgraph 可视化["可视化"] H --> I["plt.plot多条ROC曲线<br/>叠加展示"] end

源码路径：benchmarks/bench_lof.py - __main__（1-100行）

# 第 87 章 —— 设置全局随机种子为2，用于控制SA数据集中异常点的随机选择
# 第 87 章 —— 与隔离森林基准使用不同的种子，增加了测试场景的多样性
random_state = 2

# 第 87 章 —— 定义与隔离森林基准相同的六个数据集
datasets = ["http", "smtp", "SA", "SF", "shuttle", "forestcover"]

# 第 87 章 —— 创建图表用于绘制所有数据集的ROC曲线
plt.figure()

# 第 87 章 —— 遍历每个数据集进行转导式评估
for dataset_name in datasets:
    print("loading data")  # 打印状态信息

    # 数据加载阶段：与bench_isolation_forest.py高度相似的加载逻辑
    if dataset_name in ["http", "smtp", "SA", "SF"]:
        # 从KDD Cup 99下载指定子集，percent10=True使用10%抽样版
        # random_state确保异常选择的可复现性
        dataset = fetch_kddcup99(
            subset=dataset_name, percent10=True, random_state=random_state
        )
        X = dataset.data
        y = dataset.target

    # shuttle数据集的处理：移除类别4，将类别1设为正常
    if dataset_name == "shuttle":
        dataset = fetch_openml("shuttle", as_frame=False)
        X = dataset.data
        y = dataset.target.astype(np.int64)
        s = y != 4  # 创建掩码移除类别4
        X = X[s, :]
        y = y[s]
        y = (y != 1).astype(int)  # 类别1为正常，其余为异常

    # forestcover数据集的处理：只保留类别2和4，以类别2为正常
    if dataset_name == "forestcover":
        dataset = fetch_covtype()  # 注意：此处未使用random_state参数
        X = dataset.data
        y = dataset.target
        s = (y == 2) + (y == 4)  # 创建掩码只保留类别2和4
        X = X[s, :]
        y = y[s]
        y = (y != 2).astype(int)  # 类别2为正常，类别4为异常

    print("vectorizing data")

    # 独热编码处理字符串分类特征（与隔离森林基准一致）
    if dataset_name == "SF":
        lb = LabelBinarizer()
        x1 = lb.fit_transform(X[:, 1].astype(str))
        X = np.c_[X[:, :1], x1, X[:, 2:]]
        y = (y != b"normal.").astype(int)

    if dataset_name == "SA":
        lb = LabelBinarizer()
        x1 = lb.fit_transform(X[:, 1].astype(str))
        x2 = lb.fit_transform(X[:, 2].astype(str))
        x3 = lb.fit_transform(X[:, 3].astype(str))
        X = np.c_[X[:, :1], x1, x2, x3, X[:, 4:]]
        y = (y != b"normal.").astype(int)

    if dataset_name == "http" or dataset_name == "smtp":
        y = (y != b"normal.").astype(int)

    X = X.astype(float)  # 统一转换为浮点类型

    print("LocalOutlierFactor processing...")
    # 创建LOF模型，使用默认的20个邻居计算局部密度
    model = LocalOutlierFactor(n_neighbors=20)
    tstart = time()  # 记录开始时间
    model.fit(X)  # 在全量数据上拟合
    fit_time = time() - tstart  # 计算耗时

    # 获取异常分数：negative_outlier_factor_是LOF训练后立即可用的属性
    # 取负值使越负的分数表示越异常，符合roc_curve的要求
    scoring = -model.negative_outlier_factor_

    # 在全量数据上计算ROC曲线和AUC值
    fpr, tpr, thresholds = roc_curve(y, scoring)
    AUC = auc(fpr, tpr)

    # 绘制当前数据集的ROC曲线，包含AUC值和训练时间
    plt.plot(
        fpr,
        tpr,
        lw=1,
        label="ROC for %s (area = %0.3f, train-time: %0.2fs)"
        % (dataset_name, AUC, fit_time),
    )

# 第 87 章 —— 设置坐标轴范围和标签
plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("Receiver operating characteristic")
plt.legend(loc="lower right")
plt.show()

这段代码定义了 LOF 基准的核心循环。与隔离森林基准相比，它有三个本质区别。下表对比了两种基准在评估范式上的差异：

| 维度 | IsolationForest 基准 | LOF 基准 |

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

| 数据切分 | 50/50 训练/测试 | 全量数据 |

| 评估指标 | 测试集 ROC AUC | 训练集 ROC AUC |

| 评分来源 | decision_function(X_test) | negative_outlier_factor_ 属性 |

| 打乱数据 | shuffle=True | 无shuffle |

| 随机性来源 | random_state控制整体随机 | 仅控制SA数据集异常选择 |

没有 shuffle 操作也耐人寻味——因为 LOF 在全量数据上评估，没有"测试集"概念，自然不需要打乱顺序。随机性仅来自 SA 数据集内部异常点的随机选择，因此 random_state=2 就足以保证结果可复现。negative_outlier_factor_ 是 LOF 训练后立即可用的属性，无需再次调用 predict，这是 LOF 与隔离森林在"评分接口"上的根本差异。

整体而言，LOF 基准的代码结构比隔离森林基准更简洁：少了训练/测试切分逻辑，少了 decision_function 调用，少了 with_decision_function_histograms 开关。这反映了 LOF 算法本身的设计哲学——它是一个"就地异常检测器"（in-place detector），没有清晰的"训练-预测"边界。

87.6.2 LOF 基准的设计取舍

LOF 选择转导式评估而非归纳式评估，本质上是由其算法原理决定的。LOF 计算每个点与其 k 个最近邻的平均局部密度之比，这个比值依赖于所有数据点的相对位置关系。如果将数据切分为训练集和测试集，测试集中每个点的邻居可能不完整，导致异常分数不可靠。因此，LOF 必须使用全量数据进行评估，这使得它更像是一种"数据探索工具"而非"预测模型"。这种设计权衡换来的是算法的简洁性——不需要显式的训练/预测分离，也不需要处理冷启动问题。

87.7 保序回归基准 —— PAVA 算法的"复杂度陷阱"

这一节我们剖析 bench_isotonic.py，它专注于保序回归（Isotonic Regression）的 PAVA（Pool Adjacent Violators Algorithm）算法在不同数据模式下的复杂度验证。

87.7.1 数据生成器与基准流程

下图展示了保序回归基准的架构，包括三类合成数据生成器、纳秒级计时和复杂度可视化：

flowchart TB subgraph 数据生成["三类合成数据生成器"] A["perturbed_logarithm<br/>对数趋势+均匀扰动"] --> A1["50*log1+range<br/>单调递增对数"] A1 --> A2["±50随机噪声<br/>模拟一般回归场景"] B["logistic<br/>逻辑回归型单调"] --> B1["sort normal<br/>排序正态分布"] B1 --> B2["expit sigmoid<br/>概率阈值生成0/1"] B2 --> B3["非严格单调序列<br/>适合概率校准"] C["pathological<br/>病理数据"] --> C1["递增+V字+递增<br/>三段拼接"] C1 --> C2["触发PAVA最坏情况<br/>O(n²)池合并循环"] end subgraph 计时["高精度计时模块"] D["bench_isotonic_regression Y"] --> D1["gc.collect<br/>主动垃圾回收"] D1 --> D2["default_timer纳秒级计时"] D2 --> D3["isotonic_regression Y<br/>单次调用执行"] D3 --> D4["返回elapsed时间"] end subgraph 主循环["指数级规模扫描"] E["exponent: log_min→log_max"] --> E1["n = 10**exponent<br/>指数增长规模"] E1 --> E2["生成对应规模数据"] E2 --> E3["迭代N次取平均<br/>减少测量噪声"] E3 --> E4["保存n和mean_time"] end subgraph 可视化["可选可视化输出"] E4 --> F["plt.plot散点连线"] F --> G["plt.loglog双对数坐标"] G --> H["斜率1=O(n<br/>斜率2=O(n²"] end

先看三个合成数据生成器，它们分别代表了三种复杂度场景：

源码路径：benchmarks/bench_isotonic.py - generate_perturbed_logarithm_dataset()（25-27行）

def generate_perturbed_logarithm_dataset(size):
    """
    生成"对数趋势 + 均匀扰动"的合成数据。
    模拟一般的有序回归场景，数据整体呈现单调趋势但带有随机噪声。
    """
    # np.arange(size)生成0到size-1的整数序列
    # 50.0 * np.log(1 + np.arange(size))创造单调递增的对数趋势
    # 随着size增大，对数项的增量逐渐减小，体现边际递减效应
    # np.random.randint(-50, 50, size=size)叠加±50的均匀随机噪声
    # 最终返回的是有噪点的单调递增序列
    return np.random.randint(-50, 50, size=size) + 50.0 * np.log(1 + np.arange(size))

这段代码定义了 generate_perturbed_logarithm_dataset，生成"对数趋势 + 均匀扰动"的数据。50.0 * np.log(1 + np.arange(size)) 创造了一个单调递增的对数趋势，np.random.randint(-50, 50, size=size) 叠加了 ±50 的随机噪声。这种数据模拟的是一般的有序回归场景——大多数现实数据都遵循某种单调趋势但带有噪声。

源码路径：benchmarks/bench_isotonic.py - generate_logistic_dataset()（29-31行）

def generate_logistic_dataset(size):
    """
    生成逻辑回归型的单调序列数据。
    模拟分类概率校准场景，数据呈现S型单调递增趋势。
    """
    # np.random.normal生成size个标准正态分布随机数
    # np.sort将随机数排序，确保输入单调递增
    # 这是必要的，因为保序回归要求输入x必须单调
    X = np.sort(np.random.normal(size=size))

    # expit是scipy的sigmoid函数，将正态分布的X映射到(0,1)区间
    # np.random.random生成(0,1)均匀分布的随机数
    # 比较操作生成布尔数组，再转换为整数0/1
    # 结果是一个"非严格单调"序列：局部可能有波动但整体递增
    # 这模拟了概率校准场景：预测概率与真实概率的关系
    return np.random.random(size=size) < expit(X)

这段代码定义了 generate_logistic_dataset，模拟分类概率校准场景。expit(X) 是 sigmoid 函数，将正态分布的 X 映射到 (0, 1) 区间作为概率阈值；随机数小于该阈值则为 1，否则为 0。结果是一个"非严格"单调序列——局部波动但整体递增，适合验证保序回归对概率校准的拟合能力。

源码路径：benchmarks/bench_isotonic.py - generate_pathological_dataset()（33-37行）

def generate_pathological_dataset(size):
    """
    生成触发PAVA算法O(n²)最坏情况的病理数据。
    这种数据模式会让原始PAVA实现陷入密集的池合并循环。
    """
    # np.arange(size)生成0到size-1的递增序列
    # np.arange(-(size-1), size)生成-(size-1)到size-1的序列，形成"V字形"
    # np.arange(-(size-1), 1)生成-(size-1)到0的递增序列
    # 三段拼接形成"递增→V字→递增"的锯齿状模式
    # 这种模式迫使PAVA算法在每个violator点都触发跨多个块的池合并
    # Triggers O(n^2) complexity on the original implementation.
    return np.r_[
        np.arange(size), np.arange(-(size - 1), size), np.arange(-(size - 1), 1)
    ]

这段代码定义了 generate_pathological_dataset，是触发 PAVA 最坏情况的关键。它由三段拼接而成：np.arange(size) 是 0 到 size-1 的递增段；np.arange(-(size-1), size) 是 -(size-1) 到 size-1 的"V 字形"段；np.arange(-(size-1), 1) 是 -(size-1) 到 0 的递增段。这种"递增-V字-递增"的组合让原始 PAVA 实现陷入 O(n²) 的池合并循环，正是 PAVA 算法复杂度分析的经典病理数据。

接下来看高精度计时函数：

源码路径：benchmarks/bench_isotonic.py - bench_isotonic_regression()（45-53行）

def bench_isotonic_regression(Y):
    """
    对输入数据执行单次保序回归迭代，报告总耗时（秒）。
    使用纳秒级计时器和GC控制确保测量稳定性。
    """
    # 主动触发垃圾回收，消除CPython引用计数GC的不确定性
    # CPython的GC在累积到一定阈值时会触发回收，导致某次迭代突然变慢
    # 主动调用可让基准测试更稳定
    gc.collect()

    # 使用timeit的纳秒级计时器，在不同平台选择最精确的时钟源
    # Linux上使用clock_gettime，Windows上使用QueryPerformanceCounter
    tstart = default_timer()

    # 执行单次保序回归调用
    isotonic_regression(Y)

    # 返回经过的精确时间
    return default_timer() - tstart

这段代码定义了 bench_isotonic_regression，是单次 PAVA 调用的高精度计时封装。gc.collect() 是关键工程细节——CPython 的引用计数 GC 在累积到一定阈值时会触发回收，导致某次迭代突然变慢几十倍，主动调用可让基准测试更稳定。timeit.default_timer 在不同平台选择最精确的计时器（Linux 上是 clock_gettime，Windows 上是 QueryPerformanceCounter）。

最后看主流程的指数级规模扫描：

源码路径：benchmarks/bench_isotonic.py - __main__（55-95行）

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Isotonic Regression benchmark tool")
    parser.add_argument("--seed", type=int, help="RNG seed")  # 随机数生成器种子
    parser.add_argument(
        "--iterations",
        type=int,
        required=True,
        help="Number of iterations to average timings over for each problem size",
    )
    parser.add_argument(
        "--log_min_problem_size",
        type=int,
        required=True,
        help="Base 10 logarithm of the minimum problem size",
    )
    parser.add_argument(
        "--log_max_problem_size",
        type=int,
        required=True,
        help="Base 10 logarithm of the maximum problem size",
    )
    parser.add_argument(
        "--show_plot", action="store_true", help="Plot timing output with matplotlib"
    )
    parser.add_argument(
        "--dataset",
        choices=DATASET_GENERATORS.keys(),  # 限制为预定义的三个生成器之一
        required=True
    )

    args = parser.parse_args()

    np.random.seed(args.seed)  # 设置NumPy随机种子确保可复现性

    timings = []  # 初始化结果列表，存储(规模, 平均耗时)元组

    # 指数级规模扫描：遍历从log_min到log_max的所有指数
    # 例如：log_min=2, log_max=5 意味着扫描 10², 10³, 10⁴, 10⁵
    for exponent in range(args.log_min_problem_size, args.log_max_problem_size):
        n = 10**exponent  # 计算当前规模：10的exponent次方
        Y = DATASET_GENERATORS[args.dataset](n)  # 调用对应的数据生成器

        # 对同一规模执行多次迭代，取平均值以减少测量噪声
        # 单次计时可能受系统调度、缓存状态等影响，平均值更稳定
        time_per_iteration = [
            bench_isotonic_regression(Y) for i in range(args.iterations)
        ]
        timing = (n, np.mean(time_per_iteration))  # 打包为元组
        timings.append(timing)  # 累积到结果列表

        # 如果不是绘图模式，直接打印规模和时间
        if not args.show_plot:
            print(n, np.mean(time_per_iteration))

    if args.show_plot:
        # 使用zip(*timings)解压元组列表，分别获取x值和y值
        # *zip(*timings)将[(n1,t1), (n2,t2), ...]解压为(n1,n2,...), (t1,t2,...)
        # plt.plot接收解包后的参数绘制散点连线图
        plt.plot(*zip(*timings))
        plt.title("Average time taken running isotonic regression")
        plt.xlabel("Number of observations")
        plt.ylabel("Time (s)")
        plt.axis("tight")
        # plt.loglog将X轴和Y轴都设置为对数刻度
        # 在双对数坐标下，O(n)复杂度表现为斜率1的直线
        # O(n²)复杂度表现为斜率2的直线，便于观察算法复杂度
        plt.loglog()
        plt.show()

这段代码定义了入口点与主循环。指数级规模扫描（10**exponent）是验证算法复杂度的标准做法：在双对数坐标下，线性复杂度 O(n) 表现为斜率 1 的直线，平方复杂度 O(n²) 表现为斜率 2 的直线。这让"病理数据是否真的触发了 O(n²) 行为"一目了然。

--dataset 参数使用 argparse 的 choices 限定三个生成器名之一，避免拼写错误。--show_plot 是 action="store_true" 形式——仅在命令行传入该标志时才为真，默认 False。

87.7.2 保序回归基准的设计取舍

保序回归基准选择纳秒级计时器（timeit.default_timer）而非普通的 time.time，主要基于算法特性的考量。保序回归在单次调用上的执行时间非常短（微秒到毫秒级），尤其是对于小规模数据，time.time 的毫秒级精度可能引入较大误差。而隔离森林基准对整个 fit + predict 流程计时，每次运行时间在毫秒到秒级，time.time 的精度已经足够。此外，病理数据生成器中三段拼接的设计是有意为之——它强迫 PAVA 算法在每个 violator 点都触发跨多个块的池合并，将最坏情况的复杂度暴露出来，这正是基准测试想要验证的核心问题。

87.8 隔离森林基准的设计取舍

隔离森林基准采用训练/测试分离的评估范式，主要基于以下考量：作为归纳式学习算法，隔离森林构建的隔离树结构可以独立应用于新数据，因此可以通过分离训练集与测试集来衡量模型的泛化能力。这种设计模拟了真实场景中的部署需求——模型在已知数据上训练后，需要对未知数据进行预测。相比之下，LOF 是一种转导式算法，每个点的异常得分都依赖于全体数据的邻域关系，不存在"独立于训练数据"的预测阶段，因此必须使用全量数据进行评估。

在隔离森林预测基准中，训练阶段固定使用 n_jobs=-1 而预测阶段遍历 n_jobs 的设计背后有明确的工程考量。训练阶段涉及构建 N 棵随机树，每棵树的构建过程相互独立，非常适合并行化，因此使用全部 CPU 核心可以让模型快速构建，减少实验启动的开销，让对比重点集中在预测阶段。预测阶段是批量推理任务，可控变量正是评估"并行推理加速比"的关键。通过 parallel_config("threading", n_jobs=n_jobs) 上下文管理器，可以在不重建模型的情况下动态切换线程数，这种设计让性能对比更加纯粹——每次测量只受单一因素影响。

保序回归基准选择纳秒级计时器（timeit.default_timer）而非普通的 time.time，主要基于算法特性的考量。保序回归在单次调用上的执行时间非常短（微秒到毫秒级），尤其是对于小规模数据，time.time 的毫秒级精度可能引入较大误差。主动调用 gc.collect() 消除 CPython 引用计数 GC 的不确定性，避免某次迭代因 GC 触发而突然变慢。而隔离森林基准对整个 fit + predict 流程计时，每次运行时间在毫秒到秒级，time.time 的毫秒级精度已经足够。

关于 1000 次空循环的预热技巧，这是社区中流传的经验式方法，目标是让 CPU 缓存、分支预测器、内存分配器等进入稳定状态。相比于正式的 timeit 模块（它会测量多次取最优值），这里的预热更轻量——1000 次空循环足够触发 CPU 微指令流水线预热，又不会让代码变得冗长。更标准的做法是使用 timeit 的 repeat 或 number 参数，但会让这个基准脚本多出几十行代码。

87.9 动手练习

1. 对比 IsolationForest 与 LOF 的评估范式差异

   阅读 benchmarks/bench_isolation_forest.py 与 benchmarks/bench_lof.py 的主流程。首先思考为什么 IsolationForest 要拆分训练/测试集，而 LOF 使用全量数据。这是因为隔离森林是归纳式算法，构建的模型可以独立应用于新数据；而 LOF 是转导式算法，每个点的异常得分都依赖于全体数据。其次，LOF 的 negative_outlier_factor_ 属性与 IsolationForest 的 decision_function 有何本质区别？前者是就地计算的结果属性，后者是模型方法。最后，如果要在 LOF 基准中引入测试集，需要修改哪些核心逻辑？

2. 分析隔离森林预测基准的参数扫描设计

   阅读 benchmarks/bench_isolation_forest_predict.py 的 bench 函数。首先思考为什么训练阶段固定 n_jobs=-1 而预测阶段遍历 [1,2,3,4]？这是为了分离关注点，让训练时间固定而聚焦于预测并行度的量化。其次，get_data 生成的数据分布特点是什么？为何选择两个高斯簇+均匀噪声？这种设计模拟了真实异常检测场景——正常数据聚集成簇，异常点散布在边界外。最后，1000次空循环 for _ in range(1000): 1+1 的作用是什么？是否有更标准的替代方案？

3. 探究保序回归的病理数据构造原理

   阅读 benchmarks/bench_isotonic.py 的 generate_pathological_dataset。首先思考该函数生成的序列模式为何能触发原始 PAVA 算法的 O(n²) 最坏情况？关键在于三段拼接形成的"递增-V字-递增"锯齿状模式。其次，对比三种数据生成器，它们分别模拟了哪些真实应用场景？perturbed_logarithm 模拟一般有序回归，logistic 模拟概率校准，pathological 触发最坏情况。最后，如果将 isotonic_regression 替换为 IsotonicRegression().fit，基准结果会有何变化？

87.10 本章小结

这一章我们深入剖析了 scikit-learn 异常检测专题的四个独立基准脚本，揭示了异常检测算法在性能评估与复杂度验证上的工程智慧。首先介绍了 IsolationForest 在六个经典数据集上的 ROC 评估流程，它展现了"训练/测试分离 + 并行加速 + AUC 量化"的归纳式异常检测基准范式；其次分析了 IsolationForest 预测专项，它通过四维参数网格与 parallel_config 上下文管理展示了"分离训练/推理并行度"的精细实验设计；然后剖析了 LOF 基准，它以转导式学习范式在全量数据上计算 negative_outlier_factor_；最后探讨了保序回归基准，它用三类合成数据生成器（一般/校准/病理）精准探测 PAVA 算法的复杂度边界。

| 概念 | 解释 |

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

| IsolationForest基准 | 随机分割树集成，训练/测试分离，ROC/AUC评估异常检测能力 |

| IsolationForest预测基准 | 四维参数网格扫描，分离训练/推理计时，支持跨分支性能对比 |

| LocalOutlierFactor基准 | 转导式学习，全量数据拟合评估，邻域密度比值作为异常分数 |

| IsotonicRegression基准 | PAVA算法复杂度验证，三类合成数据探测最佳/平均/最坏情况 |

| 合成数据生成器 | perturbed_logarithm(一般)/logistic(校准)/pathological(O²陷阱) |

| 并行控制 | joblib.parallel_config 精确控制预测阶段 n_jobs，消除训练干扰 |

| 缓存清理技巧 | 1000次空循环预热，减少内存分配/缓存命中对计时的干扰 |

| argparse子命令模式 | bench/plot 分离，CI/CD友好，支持历史版本性能回归分析 |

下一章中，我们将学习 NMF、收敛性与 t-SNE 专题，深入探讨非负矩阵分解的自定义求解器、逻辑回归收敛性、t-SNE 大规模实验与多项式核近似等计算密集型任务的基准脚本。

第 88 章 —— 独立基准脚本：NMF、收敛性与 t-SNE 专题 —— 挑战"计算密集巅峰"

88.1 学习目标

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

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

• 理解 NMF 投影梯度求解器的核心数学原理与实现细节

• 掌握逻辑回归不同求解器收敛性的基准测试方法与评估指标

• 了解 SAGA 求解器在多种数据集、惩罚类型与数据精度下的性能对比实验设计

• 理解 t-SNE 大规模基准中嵌入质量评估（nn_accuracy）与 PCA 预处理的作用

• 掌握多项式核近似中 PolynomialCountSketch 与 Nystroem 的精度-效率权衡分析方法

• 能够阅读并修改基准测试脚本以适配新的算法或数据集

88.2 生活类比

想象这些基准测试脚本是一场"算法奥林匹克"的裁判系统。在每一场赛事中，裁判们都需要用公正、可量化、可复现的方式评判算法的真实水平。

NMF 基准就像一场举重比赛的三种姿势对决：坐标下降（CD）像标准抓举，运动员稳扎稳打，每一步都精确可控；投影梯度（PG）像挺举，爆发力强但技术难度高，需要精妙的线搜索才能找到最佳发力点；乘法更新（MU）像力量举，纯粹的乘性迭代，简洁但收敛行为独特。裁判 _PGNMF 甚至复原了已退役的"经典姿势"（已从 scikit-learn 0.19 版本移除的 Projected Gradient 求解器），让它重新站上擂台，与现役选手同台较量。

逻辑回归收敛性基准则是一场马拉松全程计时赛：七位选手（liblinear、SAG、Newton-CG、LBFGS、SGD 等）同场竞技，每隔几公里（迭代轮次）就有裁判记录一次配速（损失值）与心率（准确率）。最后，四位裁判分别绘制"配速-时间"曲线，从目标函数下降速度、训练集准确率提升、测试集准确率演变、目标函数对数残差四个维度，全方位剖析哪位选手最稳健、最适合长跑。

SAGA 多数据集基准则是十项全能积分榜：同一位运动员（SAGA 求解器）要在 rcv1、digits、iris 三个赛场、两种规则（l1/l2）、两种装备（float32/float64）下全能比拼。多位裁判（joblib.Parallel）同时执裁，最后汇总成 JSON 总成绩单，让跨数据集、跨精度、跨惩罚的比较一目了然。

t-SNE 基准犹如一场地图绘制精度大赛：参赛者（sklearn TSNE vs bhtsne）需要将 784 维的高维地形投影到 2D 平面，裁判用"最近邻保持率"（nn_accuracy）检查邻里关系有无扭曲——原空间中 A 的最近邻是 B，嵌入后是否还是？得分越高，说明地形保真度越好。PCA 预处理就像先用地形等高线仪把 784 维压缩到 50 维，再让选手绘制地图，大大降低了比赛难度。

多项式核近似基准则是压缩算法压缩率测试：PolynomialCountSketch 像"流式压缩"，边读边压，几乎没有训练开销；Nystroem 像"字典压缩"，需先扫描全文档建立码本，训练阶段相当耗时。考核指标双重——压缩后分类精度（准确率）与压缩速度（可扩展性），看谁在保持精度的同时能处理更大的数据。

88.3 源码地图

benchmarks/bench_plot_nmf.py
├── _norm 函数 (78-82行)                      # 点积实现的欧氏范数
├── _nls_subproblem 函数 (85-136行)            # 非负最小二乘子问题求解器
│   ├── 计算 WtX, WtW                          # 预计算常量矩阵
│   ├── 投影梯度下降主循环 (max_iter)            # 外层迭代
│   ├── 线搜索内循环 (20次)                     # 步长调整 (gamma, beta, sigma)
│   ├── 投影步骤 (Hn *= Hn > 0)                  # 非负约束
│   └── 充分下降条件检查                         # Armijo 规则变体
├── _fit_projected_gradient 函数 (176-215行)     # W/H 交替优化主流程
│   ├── 计算梯度 gradW, gradH                    # 使用 safe_sparse_dot
│   ├── 交替调用 _nls_subproblem 更新 W/H        # 转置技巧复用子问题求解器
│   ├── 动态调整容差 tolW/tolH                   # 加速收敛
│   └── 最终修正 W                               # 处理边界情况
├── _PGNMF 类 (138-215行)                       # 自定义投影梯度 NMF 求解器
│   ├── __init__                                 # 初始化参数 (solver='pg', nls_max_iter等)
│   ├── fit                                      # 调用 fit_transform
│   ├── transform                                # 固定 H 求 W (调用 _nls_subproblem)
│   ├── inverse_transform                        # W × H 重构
│   ├── fit_transform                            # 同时更新 W/H (调用 _fit_projected_gradient)
│   └── _fit_transform                           # 核心拟合逻辑 (初始化、收敛判断、警告)
├── bench_one 函数 (242-274行)                    # 单次实验执行与缓存
│   ├── 拷贝初始 W0/H0                           # 避免污染
│   ├── 实例化分类器并 fit_transform             # 核心计算
│   ├── 计算 beta_divergence (beta=2)            # Frobenius 范数损失
│   └── 返回 loss, duration                      # 供 run_bench 聚合
├── run_bench 函数 (276-322行)                    # 基准测试主控流程
│   ├── 遍历求解器、初始化方式、迭代次数           # 网格搜索参数空间
│   ├── 调用 bench_one (joblib 缓存)              # 单次实验执行
│   ├── 收集 loss/time/init/method                # 结果记录
│   └── 绘图调用 plot_results                     # 可视化输出
├── plot_results 函数 (324-352行)                 # 分面绘制 loss-time 曲线
├── load_20news / load_faces                     # 数据加载器 (稀疏/稠密)
├── build_clfs                                   # 构建三类求解器对比组 (CD, PG, MU)
└── __main__                                     # 入口: 设置超参、加载数据、运行基准、绘图

benchmarks/bench_rcv1_logreg_convergence.py
├── get_loss 函数 (18-26行)                      # 逻辑损失 + L2 正则化计算
├── bench_one 函数 (29-52行)                      # 单次训练与评估 (joblib 缓存)
│   ├── 设置 max_iter/random_state                # 兼容不同求解器参数名
│   ├── 计时 fit                                  # 训练耗时
│   ├── 提取 coef_/intercept_                     # 统一不同模型属性
│   ├── 调用 get_loss 计算训练目标                # 目标函数值
│   └── 计算 train/test score                     # 准确率评估
├── bench 函数 (55-84行)                          # 多求解器、多迭代次数批量实验
│   ├── 遍历 clfs 列表                            # 每个求解器一组迭代范围
│   ├── gc.collect()                              # 内存清理减少噪声
│   ├── 聚合 losses/scores/durations              # 记录曲线数据
│   └── 打印进度                                  # 监控实验状态
├── 绘图函数族                                    # 四种可视化视角
│   ├── plot_train_losses (87-94)                 # 目标函数随时间下降曲线
│   ├── plot_train_scores (96-104)                # 训练集准确率随时间
│   ├── plot_test_scores (106-114)                # 测试集准确率随时间
│   └── plot_dloss (116-128)                      # 对数目标差距 (log(p - p*))
├── get_max_squared_sum (130-132)                 # 计算 max ||x_i||^2 用于步长
├── 数据准备区 (RCV1 加载、标签二值化、切分)       # 23149 样本作测试集
├── 求解器配置 clfs 列表                          # 7 个基线求解器 + 可选 lightning
│   ├── LR-liblinear (primal/dual)                # 坐标下降
│   ├── LR-SAG                                    # 随机平均梯度
│   ├── LR-newton-cg / LR-lbfgs                   # 二阶/拟牛顿
│   └── SGD                                       # 随机梯度下降
└── __main__ (134-220)                            # 执行 bench -> 绘制 4 图 -> show

benchmarks/bench_saga.py
├── fit_single 函数 (22-78行)                      # 单求解器/单数据集/单精度完整实验
│   ├── 参数映射 (penalty -> alpha/beta/l1_ratio)   # 统一数学定义
│   ├── 实例化模型 (Lightning 或 sklearn)          # 统一接口适配
│   ├── 多轮 max_iter 循环 (step=2)                 # 记录轨迹
│   ├── X_train.max()                              # 预热缓存
│   ├── time.clock() 计时 fit                      # CPU 时间
│   ├── 计算目标函数值 (log_loss + 正则项)          # 统一数学定义
│   ├── 计算准确率                                 # 评估指标
│   └── 返回 (lr, times, train_scores, test_scores, accuracies)
├── _predict_proba (180-188)                       # Lightning 多类概率补全 (softmax)
├── exp 函数 (80-124)                              # 实验调度器
│   ├── 加载数据集 (rcv1/digits/iris/20news)        # 单/多目标标签转换
│   ├── 截取 n_samples                             # 控制规模
│   ├── Parallel(delayed(fit_single))              # 并行跑所有 solver x dtype
│   └── 序列化 JSON 到 bench_saga.json              # 持久化中间结果
├── plot 函数 (126-178)                            # 从 JSON 读取并绘制 4 子图
│   ├── pandas DataFrame 聚合                      # 整理长表格
│   ├── 4 子图: train_obj/test_obj/acc/iter_obj     # 多维度视角
│   ├── 样式映射 (color/linestyle/alpha)           # 区分 solver/dtype
│   └── 保存 PNG                                   # 批量输出图表
└── __main__ (180-202)                             # 网格: penalty x n_samples -> exp -> plot

benchmarks/bench_tsne_mnist.py
├── load_data 函数 (cached) (54-68)                # fetch_openml MNIST + 归一化 + 缓存
├── nn_accuracy 函数 (37-45)                       # 最近邻保持率 (核心质量指标)
│   ├── NearestNeighbors(n_jobs=-1)                # 并行最近邻搜索
│   ├── 比较原始空间 vs 嵌入空间邻居索引            # 一致性检查
│   └── 返回平均匹配率                              # 0-1 分数
├── tsne_fit_transform (48-51)                     # 统一接口包装 (返回 embedding + n_iter)
├── sanitize (53-54)                               # 文件名安全化
├── argparse 参数解析 (70-90)                      # order/perplexity/bhtsne/all/profile
├── PCA 预处理 (可选) (92-98)                       # 降至 50 维加速
├── 方法注册列表 methods (100-120)                  # sklearn TSNE + 可选 bhtsne
│   ├── sklearn TSNE (init='pca', n_iter=1000)     # 精确/巴恩斯-赫特
│   └── bhtsne (run_bh_tsne, 无 PCA、无 n_iter 报告) # 参考实现
├── 数据规模扫描 data_size (122)                    # [100, 500, 1k, 5k, 10k, 70k]
├── 主循环 (124-150)                               # 遍历规模 x 方法
│   ├── 计时 method(X_train)                       # fit_transform 耗时
│   ├── 计算 nn_accuracy                           # 质量评估
│   ├── 记录 JSON (method, duration, n_samples)     # 持久化日志
│   └── 保存 .npy (原始数据、标签、嵌入结果)        # 供可视化脚本使用
└── 实时打印进度                                    # 监控长时间实验

benchmarks/plot_tsne_mnist.py
├── argparse (labels/embedding 路径)               # 指定可视化文件
├── 加载 .npy                                      # 嵌入坐标 + 标签
├── 按类别散点图 (alpha=0.2)                        # 10 类数字可视化
└── plt.legend/show                                # 交互式查看

benchmarks/bench_plot_polynomial_kernel_approximation.py
├── 数据准备: load_digits + train_test_split (70/30)
├── 基准模型评估                                    # 两条虚线基准线
│   ├── LinearSVC (原始空间)                        # lsvm_score
│   └── SVC(kernel='poly')                          # ksvm_score (精度上界)
├── PolynomialCountSketch 精度扫描 (68-152行)      # n_components 20..400
│   ├── n_runs=5 平均                              # 抗随机性
│   ├── Pipeline(PS + LinearSVC)                   # 组合估计器
│   └── 记录 ps_svm_scores                         # 准确率曲线
├── Nystroem 精度扫描 (155-190行)                   # 同协议对比
│   ├── kernel='poly', gamma=1, degree=2, coef0=0
│   └── 记录 ny_svm_scores
├── 精度对比图 (fig1)                               # 4 条曲线 + 虚线基准
├── 可扩展性实验 (大规模随机数据 10k x 100)        # out_dims 500..6000
│   ├── time(ps.fit_transform)                     # 仅变换时间 (含 fit)
│   ├── time(ny.fit_transform)                     # Nystroem 含昂贵 fit 阶段
│   └── 记录 ps_svm_times / ny_svm_times
├── 可扩展性对比图 (fig2)                           # 时间随 n_components 增长
│   ├── PS: O(n(d + k log k)) 预期线性增长
│   └── Nystroem: O(n(dk + k^2)) 预期超线性
└── plt.show()                                     # 显示双图

88.4 NMF 基准与自定义投影梯度求解器 —— 矩阵分解的"深度解剖"

NMF（Non-Negative Matrix Factorization，非负矩阵分解）是 scikit-learn 中一个重要算法。它的目标是将一个非负矩阵 X 分解为两个非负矩阵 W 和 H 的乘积，即 X ≈ W @ H。在文本挖掘、图像处理、推荐系统等领域有广泛应用。本节聚焦的 bench_plot_nmf.py 是一个特殊的基准脚本——它不仅测试了 scikit-learn 中现存的两种求解器（Coordinate Descent 与 Multiplicative Update），还复活了一个曾经存在、后来被移除的求解器：Projected Gradient。

88.4.1 为什么需要自定义 _PGNMF 类？

scikit-learn 自 0.19 版本起移除了 Projected Gradient 求解器，但该算法在文献中（如 C.-J. Lin 2007 年的论文）仍然是非常经典的参考实现。_PGNMF 类正是为了在受控的基准环境中重现该算法，与当前的两大求解器在稀疏/稠密数据、不同初始化方式下进行"三方对决"。

源码路径：benchmarks/bench_plot_nmf.py - _PGNMF.__init__()（138-152行）

class _PGNMF(NMF):
    """Non-Negative Matrix Factorization (NMF) with projected gradient solver.

    This class is private and for comparison purpose only.
    It may change or disappear without notice.
    """
    # ① 继承自 sklearn NMF，复用参数管理与初始化逻辑
    def __init__(
        self,
        n_components=None,
        solver="pg",                   # ② 固定为投影梯度
        init=None,
        tol=1e-4,
        max_iter=200,
        random_state=None,
        alpha=0.0,
        l1_ratio=0.0,
        nls_max_iter=10,               # ③ 子问题最大迭代次数
    ):
        # ④ 调用父类构造，统一参数入口
        super().__init__(
            n_components=n_components,
            init=init,
            solver=solver,
            tol=tol,
            max_iter=max_iter,
            random_state=random_state,
            alpha_W=alpha,             # ⑤ W 的正则化系数
            alpha_H=alpha,             # ⑥ H 的正则化系数
            l1_ratio=l1_ratio,
        )
        # ⑦ 自定义子求解器的最大迭代次数
        self.nls_max_iter = nls_max_iter

这段代码定义了 _PGNMF 类，它继承自 scikit-learn 内置的 NMF 类，通过继承复用父类的所有参数管理、初始化、__repr__ 等基础设施。构造函数中重点关注三个关键设计：

• 第 2 行的 solver="pg"：这是该类唯一标识，将求解器固定为投影梯度。

• 第 5、6 行的 alpha_W 与 alpha_H：将原本一个 alpha 参数复制给 W 与 H 两个矩阵，允许后续扩展。

• 第 7 行的 nls_max_iter：这是 PG 求解器特有的超参数，控制内部非负最小二乘（Non-negative Least Squares, NLS）子问题的最大迭代次数。

通过继承机制，_PGNMF 既复用了 scikit-learn 的标准接口（fit、transform、fit_transform、inverse_transform），又定制了核心求解逻辑。

88.4.2 NLS 子问题：投影梯度下降的核心引擎

NMF 的整体优化是一个非凸问题，但固定 W 求 H 或固定 H 求 W 后，问题转化为凸的非负最小二乘子问题。_nls_subproblem 正是求解这一子问题的工作马。

源码路径：benchmarks/bench_plot_nmf.py - _nls_subproblem()（85-136行）

def _nls_subproblem(
    X, W, H, tol, max_iter, alpha=0.0, l1_ratio=0.0, sigma=0.01, beta=0.1
):
    # ① 预计算两个常量矩阵，避免内循环重复计算
    WtX = safe_sparse_dot(W.T, X)
    WtW = np.dot(W.T, W)

    # ② 初始化线搜索步长 gamma（论文中称为 alpha）
    gamma = 1
    # ④ 外层循环开始，n_iter 是实际执行的迭代次数
    for n_iter in range(1, max_iter + 1):
        # ③ 计算当前梯度：∇f(H) = W^T W H - W^T X
        grad = np.dot(WtW, H) - WtX
        # ④ 加入 L1/L2 正则化梯度
        if alpha > 0 and l1_ratio == 1.0:
            grad += alpha
        elif alpha > 0:
            grad += alpha * (l1_ratio + (1 - l1_ratio) * H)

        # ⑤ 检查投影梯度范数收敛
        # 投影梯度 = grad * mask(grad < 0 OR H > 0)
        if _norm(grad * np.logical_or(grad < 0, H > 0)) < tol:
            break

        Hp = H

        # ⑥ 线搜索内循环，最多 20 次
        for inner_iter in range(20):
            # ⑦ 梯度下降步
            Hn = H - gamma * grad
            # ⑧ 投影到非负空间：负值截断为 0
            Hn *= Hn > 0
            d = Hn - H                                  # ⑨ 步进方向
            gradd = np.dot(grad.ravel(), d.ravel())     # ⑩ 梯度内积
            dQd = np.dot(np.dot(WtW, d).ravel(), d.ravel())  # ⑪ 二次型
            # ⑫ 充分下降条件（Armijo 规则的变体）
            suff_decr = (1 - sigma) * gradd + 0.5 * dQd < 0
            if inner_iter == 0:
                # ⑬ 首轮检测：若不满足，标记为需要减小 gamma
                decr_gamma = not suff_decr

            if decr_gamma:
                if suff_decr:
                    H = Hn
                    break
                else:
                    gamma *= beta    # ⑭ 步长缩小
            elif not suff_decr or (Hp == Hn).all():
                H = Hp
                break
            else:
                gamma /= beta        # ⑮ 步长放大
                Hp = Hn

    # ⑯ 达到最大迭代时发出警告
    if n_iter == max_iter:
        warnings.warn("Iteration limit reached in nls subproblem.", ConvergenceWarning)

    return H, grad, n_iter

这段代码定义了 NMF 优化中最核心的子问题求解器——投影梯度下降法求解非负最小二乘。整个流程可分为四大部分：

第一部分（第 1-2 行）：常量矩阵预计算。WtX 与 WtW 在整个循环中不变，提前计算能避免每次迭代重新做矩阵乘法，显著加速迭代。

第二部分（第 3-7 行）：梯度计算与收敛判断。第 5 行的 grad * np.logical_or(grad < 0, H > 0) 极为精妙——这是投影梯度的数学定义。在非负约束下，若某点的梯度方向试图把它推离可行域（grad < 0），或该点已在边界（H > 0），则保留梯度分量；否则该分量为零。投影梯度范数小于容差时，认为达到最优。

第三部分（第 8-21 行）：Armijo 线搜索。这是工程实现的核心难点。每次外层迭代进入内层线搜索（最多 20 次），尝试找到一个合适的步长 gamma：

• 第 12 行的 suff_decr 是充分下降条件：（1 - sigma) * gradd + 0.5 * dQd < 0，对应经典 Armijo 规则。其中 sigma=0.01 控制"充分性"的严格程度，gradd 是梯度与步进方向的内积，dQd 是步进方向上的二次型。

• 第 13-15 行：首次进入内循环时记录"是否需要缩小 gamma"，若不满足下降条件则持续缩小 gamma *= beta（默认 beta=0.1）；若首轮已满足，则后续可以尝试放大 gamma /= beta 以加速收敛。

第四部分（第 22-25 行）：边界处理。若达到最大迭代次数仍不收敛，则发出警告。

88.4.3 W/H 交替优化与动态容差调整

外层 _fit_projected_gradient 负责在 W 与 H 之间交替优化。

源码路径：benchmarks/bench_plot_nmf.py - _fit_projected_gradient()（176-215行）

def _fit_projected_gradient(X, W, H, tol, max_iter, nls_max_iter, alpha, l1_ratio):
    # ① 计算 W 和 H 的梯度（safe_sparse_dot 支持稀疏输入）
    gradW = np.dot(W, np.dot(H, H.T)) - safe_sparse_dot(X, H.T, dense_output=True)
    gradH = np.dot(np.dot(W.T, W), H) - safe_sparse_dot(W.T, X, dense_output=True)

    # ② 初始化梯度范数作为参考尺度
    init_grad = squared_norm(gradW) + squared_norm(gradH.T)
    # ③ max(0.001, tol) 保证至少迭代一次
    tolW = max(0.001, tol) * np.sqrt(init_grad)
    tolH = tolW

    for n_iter in range(1, max_iter + 1):
        # ④ 投影梯度范数判断
        proj_grad_W = squared_norm(gradW * np.logical_or(gradW < 0, W > 0))
        proj_grad_H = squared_norm(gradH * np.logical_or(gradH < 0, H > 0))

        if (proj_grad_W + proj_grad_H) / init_grad < tol**2:
            break

        # ⑤ 通过转置技巧复用 _nls_subproblem
        # 原问题: min ||X - W H^T||^2，固定 H 求 W
        # 转置后: min ||X^T - H^T W^T||^2，固定 H^T 求 W^T
        Wt, gradWt, iterW = _nls_subproblem(
            X.T, H.T, W.T, tolW, nls_max_iter, alpha=alpha, l1_ratio=l1_ratio
        )
        W, gradW = Wt.T, gradWt.T           # ⑥ 转置回来

        # ⑦ 若子问题 1 步收敛，放宽 W 的容差以加速
        if iterW == 1:
            tolW = 0.1 * tolW

        # ⑧ 固定 W 求 H（不需转置）
        H, gradH, iterH = _nls_subproblem(
            X, W, H, tolH, nls_max_iter, alpha=alpha, l1_ratio=l1_ratio
        )
        # ⑨ 同样的动态容差启发式
        if iterH == 1:
            tolH = 0.1 * tolH

    # ⑩ 修正边界负零值
    H[H == 0] = 0

    # ⑪ 若未收敛，最后再精炼一次 W
    if n_iter == max_iter:
        Wt, _, _ = _nls_subproblem(
            X.T, H.T, W.T, tolW, nls_max_iter, alpha=alpha, l1_ratio=l1_ratio
        )
        W = Wt.T

    return W, H, n_iter

这段代码展示了交替最小化的核心思想：固定 H 求 W，再固定 W 求 H，周而复始。其中第 5-6 行的"转置技巧"非常精巧——_nls_subproblem 默认求解 min ||X - W H||，但更新 W 时问题变为 min ||X - W H^T||，通过将 X、H、W 同时转置，可以复用同一求解器，最后再把结果转置回来。

第 7、9 行的 tolW = 0.1 * tolW 是一个巧妙的工程优化：当子问题在 1 步内就收敛（即梯度已经很小），说明该子问题比较容易，可以放宽其容差以加速整体迭代。这是基于"难度感知"的动态容差调整启发式。

第 11 行：若达到最大迭代仍未收敛，最后再调用一次子问题精炼 W，保证返回的 W 与 H 尽可能匹配。

88.4.4 基准测试的主控流程

run_bench 是基准测试的总指挥，它遍历"求解器 × 初始化方式 × 迭代次数"三维网格，调用 bench_one 执行单次实验。

源码路径：benchmarks/bench_plot_nmf.py - run_bench()（276-322行）

def run_bench(X, clfs, plot_name, n_components, tol, alpha, l1_ratio):
    start = time()
    results = []
    # ① 遍历每组 (求解器名, 类型, 迭代范围, 参数)
    for name, clf_type, iter_range, clf_params in clfs:
        print("Training %s:" % name)
        # ② 三种初始化方式：nndsvd, nndsvdar, random
        for rs, init in enumerate(("nndsvd", "nndsvdar", "random")):
            print("    %s %s: " % (init, " " * (8 - len(init))), end="")
            # ③ 用同一随机种子生成初始 W, H
            W, H = _initialize_nmf(X, n_components, init, 1e-6, rs)

            # ④ 遍历预设的迭代范围
            for max_iter in iter_range:
                # ⑤ 注入超参数
                clf_params["alpha"] = alpha
                clf_params["l1_ratio"] = l1_ratio
                clf_params["max_iter"] = max_iter
                clf_params["tol"] = tol
                clf_params["random_state"] = rs
                clf_params["init"] = "custom"
                clf_params["n_components"] = n_components

                # ⑥ 调用 bench_one (joblib 缓存)
                this_loss, duration = bench_one(
                    name, X, W, H, X.shape, clf_type, clf_params, init, n_components, rs
                )

                init_name = "init='%s'" % init
                results.append((name, this_loss, duration, init_name))
                # ⑦ 用点号表示进度
                print(".", end="")
                sys.stdout.flush()
            print(" ")

    # ⑧ 用 pandas DataFrame 组织结果
    results_df = pandas.DataFrame(results, columns="method loss time init".split())
    print("Total time = %0.3f sec\n" % (time() - start))

    # ⑨ 调用 plot_results 绘制分面曲线
    plot_results(results_df, plot_name)
    return results_df

run_bench 通过三层嵌套循环构造完整的参数网格。最关键的工程实践是第 6 行：bench_one 被 joblib.Memory 装饰器包装，所有相同参数的实验结果会自动缓存——重复运行时不必重算，大幅加速调试。

下面这张图展示了 NMF 基准测试的整体数据流：

graph TD A[run_bench 启动] --> B[遍历求解器] B --> C[遍历初始化方式] C --> D[遍历迭代次数] D --> E[_initialize_nmf 生成 W0, H0] E --> F[bench_one 执行 fit_transform] F --> G[计算 beta_divergence loss] G --> H[joblib 缓存结果] H --> I[追加到 results 列表] I --> J{迭代次数遍历完?} J -->|否| D J -->|是| K{初始化方式遍历完?} K -->|否| C K -->|是| L[pandas DataFrame 聚合] L --> M[plot_results 分面绘图]

88.4.5 单次实验与缓存策略

bench_one 是最小执行单元，通过 joblib 缓存避免重复计算。

源码路径：benchmarks/bench_plot_nmf.py - bench_one()（242-274行）

@ignore_warnings(category=ConvergenceWarning)
# 第 88 章 —— use joblib to cache the results.
# 第 88 章 —— X_shape is specified in arguments for avoiding hashing X
@mem.cache(ignore=["X", "W0", "H0"])
def bench_one(
    name, X, W0, H0, X_shape, clf_type, clf_params, init, n_components, random_state
):
    # ① 拷贝初始矩阵避免污染
    W = W0.copy()
    H = H0.copy()

    # ② 实例化估计器
    clf = clf_type(**clf_params)
    # ③ 计时
    st = time()
    W = clf.fit_transform(X, W=W, H=H)
    end = time()
    # ④ 提取最终 H（fit_transform 已更新）
    H = clf.components_

    # ⑤ 计算 beta-divergence (beta=2 即 Frobenius 范数平方)
    this_loss = _beta_divergence(X, W, H, 2.0, True)
    duration = end - st
    return this_loss, duration

@mem.cache(ignore=["X", "W0", "H0"]) 装饰器是关键：它告诉 joblib 不要哈希 X、W0、H0 这三个大数据（因为哈希代价太高），而是通过 X_shape 和 random_state 等小参数来标识缓存键。如果这些参数都一致，joblib 会从磁盘读取之前保存的结果——这在反复调试可视化代码时极为有用。

88.5 逻辑回归收敛性基准 —— 求解器的"马拉松赛跑"

与 NMF 基准的"自定义求解器复活"不同，逻辑回归基准测试的是 scikit-learn 中现有的七种求解器在大规模真实数据（RCV1）上的收敛行为。

88.5.1 统一的目标函数度量

不同求解器内部优化的目标可能略有差异（如对截距的处理、L2 正则化的归一化方式）。get_loss 提供了一个统一的标尺。

源码路径：benchmarks/bench_rcv1_logreg_convergence.py - get_loss()（18-26行）

def get_loss(w, intercept, myX, myy, C):
    n_samples = myX.shape[0]
    w = w.ravel()
    # ① 逻辑损失：log(1 + exp(-y * (w·x + b)))
    p = np.mean(np.log(1.0 + np.exp(-myy * (myX.dot(w) + intercept))))
    print("%f + %f" % (p, w.dot(w) / 2.0 / C / n_samples))
    # ② L2 正则化项
    p += w.dot(w) / 2.0 / C / n_samples
    return p

get_loss 实现了带 L2 正则化的逻辑回归目标函数。需要注意第 1 行使用 np.log(1 + exp(...)) 而非 log(1 + exp(-y*score)) 的 max(0,...) 形式——前者数学上等价但数值稳定性略差，适合做"事后度量"而非训练时优化。第 2 行用 1/C/n_samples 归一化，这与 LogisticRegression(C=...) 的内部实现完全对齐。

88.5.2 单次训练与评估

源码路径：benchmarks/bench_rcv1_logreg_convergence.py - bench_one()（29-52行）

@m.cache()
def bench_one(name, clf_type, clf_params, n_iter):
    clf = clf_type(**clf_params)
    try:
        # ① 兼容 max_iter 参数（部分求解器如 liblinear 使用 n_iter）
        clf.set_params(max_iter=n_iter, random_state=42)
    except Exception:
        clf.set_params(n_iter=n_iter, random_state=42)

    # ② 计时训练过程
    st = time.time()
    clf.fit(X, y)
    end = time.time()

    # ③ 还原正则化系数 C
    try:
        C = 1.0 / clf.alpha / n_samples
    except Exception:
        C = clf.C

    # ④ 兼容不同模型对截距的命名
    try:
        intercept = clf.intercept_
    except Exception:
        intercept = 0.0

    # ⑤ 计算统一目标值与准确率
    train_loss = get_loss(clf.coef_, intercept, X, y, C)
    train_score = clf.score(X, y)
    test_score = clf.score(X_test, y_test)
    duration = end - st

    return train_loss, train_score, test_score, duration

这里有几个工程亮点：

• 第 3 行的 try/except：SGDClassifier 使用 max_iter，而 liblinear 等可能使用 n_iter，通过异常捕获做兼容。

• 第 7-9 行：通过 clf.alpha（SGD 风格）或 clf.C（LogisticRegression 风格）反推正则化系数。

• 第 12-14 行：截距（intercept）在不同模型中可能为 0，统一处理。

88.5.3 四视角可视化

基准脚本最有价值的设计是四张图组成的多维度评估体系：

源码路径：benchmarks/bench_rcv1_logreg_convergence.py - 绘图函数族（87-128行）

def plot_train_losses(clfs):
    # ① 训练损失随时间下降曲线
    plt.figure()
    for name, _, _, train_losses, _, _, durations in clfs:
        plt.plot(durations, train_losses, "-o", label=name)
        plt.legend(loc=0); plt.xlabel("seconds"); plt.ylabel("train loss")

def plot_train_scores(clfs):
    # ② 训练集准确率随时间变化
    plt.figure()
    for name, _, _, _, train_scores, _, durations in clfs:
        plt.plot(durations, train_scores, "-o", label=name)
        plt.legend(loc=0); plt.xlabel("seconds"); plt.ylabel("train score")
        plt.ylim((0.92, 0.96))

def plot_test_scores(clfs):
    # ③ 测试集准确率随时间变化
    plt.figure()
    for name, _, _, _, _, test_scores, durations in clfs:
        plt.plot(durations, test_scores, "-o", label=name)
        plt.legend(loc=0); plt.xlabel("seconds"); plt.ylabel("test score")
        plt.ylim((0.92, 0.96))

def plot_dloss(clfs):
    # ④ 对数目标差距：log|p - p*|，揭示收敛阶数
    plt.figure()
    pobj_final = []
    for name, _, _, train_losses, _, _, durations in clfs:
        pobj_final.append(train_losses[-1])

    indices = np.argsort(pobj_final)
    pobj_best = pobj_final[indices[0]]   # 取所有求解器中的最优目标值

    for name, _, _, train_losses, _, _, durations in clfs:
        log_pobj = np.log(abs(np.array(train_losses) - pobj_best)) / np.log(10)
        plt.plot(durations, log_pobj, "-o", label=name)
        plt.legend(loc=0); plt.xlabel("seconds"); plt.ylabel("log(best - train_loss)")

四张图各有侧重：

• plot_train_losses：直接看目标函数下降速度，能区分"线性收敛"（CD、SGD）与"超线性收敛"（Newton-CG、LBFGS）。

• plot_train_scores / plot_test_scores：从准确率角度评估泛化能力。固定 y 轴范围（0.92-0.96）让差异更显著。

• plot_dloss：对数尺度下的目标残差，最关键的一张图——如果某求解器在双对数图上呈现直线且斜率更陡，意味着它的收敛阶数更高（如 Newton-CG 的二次收敛性在此会非常明显）。

下面这张图揭示了逻辑回归基准测试中的数据流转：

graph LR A[RCV1 数据加载] --> B[标签二值化 CCAT vs rest] B --> C[切分训练/测试集] C --> D[7 个求解器配置] D --> E[bench_one: 训练+计时] E --> F[get_loss 计算目标] E --> G[score 计算准确率] F --> H[plot_train_losses] F --> I[plot_dloss] G --> J[plot_train_scores] G --> K[plot_test_scores]

88.5.4 求解器配置的精妙之处

求解器列表是基准脚本的核心，每个求解器的迭代范围都做了精心调整，让所有曲线在同一张图上"刚好能跑完"：

源码路径：benchmarks/bench_rcv1_logreg_convergence.py - clfs 列表（155-200行附近）

# 第 88 章 —— max_iter range
sgd_iter_range = list(range(1, 121, 10))         # SGD: 1..121，步长 10
newton_iter_range = list(range(1, 25, 3))         # Newton-CG: 1..25，步长 3
lbfgs_iter_range = list(range(1, 242, 12))        # LBFGS: 1..242，步长 12
liblinear_iter_range = list(range(1, 37, 3))      # liblinear (primal): 1..37
liblinear_dual_iter_range = list(range(1, 85, 6)) # liblinear (dual): 1..85
sag_iter_range = list(range(1, 37, 3))            # SAG: 1..37

迭代范围的精心设计反映了各求解器的计算成本：SGD 每轮极快但需多轮迭代，Newton-CG 每轮昂贵但收敛迅速，因此前者迭代范围设为 121、后者仅 25。

88.6 SAGA 求解器多数据集基准 —— 目标函数的"全景透视"

如果说 NMF 基准是"对比三种 NMF 求解器"，逻辑回归基准是"对比多种 LR 求解器"，那 SAGA 基准则是"深入一种求解器在不同战场上的全维表现"。

88.6.1 fit_single：单实验的原子单元

源码路径：benchmarks/bench_saga.py - fit_single()（22-78行）

def fit_single(
    solver, X, y, penalty="l2", single_target=True, C=1, max_iter=10,
    skip_slow=False, dtype=np.float64,
):
    if skip_slow and solver == "lightning" and penalty == "l1":
        print("skip_slowping l1 logistic regression with solver lightning.")
        return

    print("Solving %s logistic regression with penalty %s, solver %s."
          % ("binary" if single_target else "multinomial", penalty, solver))

    if solver == "lightning":
        from lightning.classification import SAGAClassifier

    # ① 多类策略映射
    # sklearn 原生 saga 支持 multinomial；其他求解器 + 单目标用 OvR
    if single_target or solver not in ["sag", "saga"]:
        multi_class = "ovr"
    else:
        multi_class = "multinomial"

    # ② 数据类型转换
    X = X.astype(dtype)
    y = y.astype(dtype)

    X_train, X_test, y_train, y_test = train_test_split(
        X, y, random_state=42, stratify=y
    )
    n_samples = X_train.shape[0]
    n_classes = np.unique(y_train).shape[0]
    test_scores = [1]; train_scores = [1]; accuracies = [1 / n_classes]
    times = [0]

    # ③ penalty 参数映射
    if penalty == "l2":
        l1_ratio = 0
        alpha = 1.0 / (C * n_samples)  # sklearn 风格
        beta = 0
        lightning_penalty = None
    else:  # l1
        l1_ratio = 1
        alpha = 0.0
        beta = 1.0 / (C * n_samples)  # lightning 风格
        lightning_penalty = "l1"

    # ④ 多次迭代记录轨迹
    for this_max_iter in range(1, max_iter + 1, 2):
        print("[%s, %s, %s] Max iter: %s" % (...))
        if solver == "lightning":
            # ⑤ Lightning 接口
            lr = SAGAClassifier(loss="log", alpha=alpha, beta=beta,
                                penalty=lightning_penalty, tol=-1,
                                max_iter=this_max_iter)
        else:
            # ⑥ sklearn 接口
            lr = LogisticRegression(solver=solver, C=C, l1_ratio=l1_ratio,
                                    fit_intercept=False, tol=0,
                                    max_iter=this_max_iter, random_state=42)
            if multi_class == "ovr":
                lr = OneVsRestClassifier(lr)

        # ⑦ 预热 CPU 缓存，保证计时公平
        X_train.max()
        t0 = time.clock()
        lr.fit(X_train, y_train)
        train_time = time.clock() - t0

        # ⑧ 计算训练/测试集目标与准确率
        scores = []
        for X, y in [(X_train, y_train), (X_test, y_test)]:
            try:
                y_pred = lr.predict_proba(X)
            except NotImplementedError:
                # Lightning 在多分类下不实现 predict_proba
                y_pred = _predict_proba(lr, X)
            # ⑨ 提取 coef（兼容 OneVsRestClassifier）
            if isinstance(lr, OneVsRestClassifier):
                coef = np.concatenate([est.coef_ for est in lr.estimators_])
            else:
                coef = lr.coef_
            # ⑩ 统一目标函数：log_loss + 0.5*alpha*||w||^2 + beta*||w||_1
            score = log_loss(y, y_pred, normalize=False) / n_samples
            score += 0.5 * alpha * np.sum(coef**2) + beta * np.sum(np.abs(coef))
            scores.append(score)
        train_score, test_score = tuple(scores)

        # ⑪ 准确率
        y_pred = lr.predict(X_test)
        accuracy = np.sum(y_pred == y_test) / y_test.shape[0]
        test_scores.append(test_score); train_scores.append(train_score)
        accuracies.append(accuracy); times.append(train_time)
    return lr, times, train_scores, test_scores, accuracies

这段代码有三个工程精髓：

第一，参数空间映射。第 12-17 行通过一个统一的 penalty 字符串映射到三套参数体系：

• scikit-learn 用 l1_ratio 与 alpha = 1/(C*n)

• lightning 用 alpha 与 beta，参数语义不同但目标函数一致

第二，统一目标函数。第 38 行 score = log_loss + 0.5*alpha*||w||^2 + beta*||w||_1 是真正"统一标尺"——无论底层用哪个库、哪个求解器，目标函数都是相同的，确保公平比较。

第三，CPU 缓存预热。第 30 行 X_train.max() 是一个看似无关但至关重要的细节：在第一次 fit 之前先访问一遍数据，确保数据被加载到 CPU 缓存，避免后续的迭代把数据加载时间算进训练耗时。

88.6.2 exp 函数：并行实验调度器

源码路径：benchmarks/bench_saga.py - exp()（80-124行）

def exp(solvers, penalty, single_target, n_samples=30000, max_iter=20,
        dataset="rcv1", n_jobs=1, skip_slow=False):
    dtypes_mapping = {"float64": np.float64, "float32": np.float32}

    if dataset == "rcv1":
        rcv1 = fetch_rcv1()
        lbin = LabelBinarizer()
        lbin.fit(rcv1.target_names)
        X = rcv1.data; y = rcv1.target
        y = lbin.inverse_transform(y)  # 转成类别标签
        le = LabelEncoder()
        y = le.fit_transform(y)
        if single_target:
            y_n = y.copy()
            y_n[y > 16] = 1; y_n[y <= 16] = 0
            y = y_n
    elif dataset == "digits":
        X, y = load_digits(return_X_y=True)
        if single_target:
            y_n = y.copy(); y_n[y < 5] = 1; y_n[y >= 5] = 0
            y = y_n
    elif dataset == "iris":
        iris = load_iris()
        X, y = iris.data, iris.target
    elif dataset == "20newspaper":
        ng = fetch_20newsgroups_vectorized()
        X = ng.data; y = ng.target
        if single_target:
            y_n = y.copy(); y_n[y > 4] = 1; y_n[y <= 16] = 0
            y = y_n

    X = X[:n_samples]
    y = y[:n_samples]

    # ① 并行执行所有 solver x dtype 组合
    out = Parallel(n_jobs=n_jobs, mmap_mode=None)(
        delayed(fit_single)(
            solver, X, y,
            penalty=penalty, single_target=single_target, dtype=dtype,
            C=1, max_iter=max_iter, skip_slow=skip_slow,
        )
        for solver in solvers
        for dtype in dtypes_mapping.values()
    )

    # ② 收集结果，序列化到 JSON
    res = []
    idx = 0
    for dtype_name in dtypes_mapping.keys():
        for solver in solvers:
            if not (skip_slow and solver == "lightning" and penalty == "l1"):
                lr, times, train_scores, test_scores, accuracies = out[idx]
                this_res = dict(solver=solver, penalty=penalty, dtype=dtype_name,
                                single_target=single_target,
                                times=times, train_scores=train_scores,
                                test_scores=test_scores, accuracies=accuracies)
                res.append(this_res)
            idx += 1

    # ③ 持久化中间结果（断点续跑关键）
    with open("bench_saga.json", "w+") as f:
        json.dump(res, f)

exp 函数是实验调度器，最关键的特性是第 31 行的 Parallel(delayed(fit_single))(...)：它将"solver × dtype"的所有组合并行执行。n_jobs=1 是默认值，但在多核机器上可以设为 -1 使用全部核心。

第 49 行的 with open("bench_saga.json", "w+") 持久化中间结果是非常好的工程实践——如果中途中断，已经完成的实验无需重跑。

88.6.3 Lightning 的 predict_proba 补全

源码路径：benchmarks/bench_saga.py - _predict_proba()（180-188行）

def _predict_proba(lr, X):
    """Predict proba for lightning for n_classes >=3."""
    # ① 计算 logit：X · coef^T + intercept
    pred = safe_sparse_dot(X, lr.coef_.T)
    if hasattr(lr, "intercept_"):
        pred += lr.intercept_
    # ② softmax 转概率
    return softmax(pred)

这段代码补全了 Lightning 在多分类场景下缺失的 predict_proba 方法。它先计算线性得分，再通过 softmax 函数转换为概率分布——这是逻辑回归的标准数学定义。

88.6.4 plot 函数：四子图全景透视

源码路径：benchmarks/bench_saga.py - plot()（126-178行）

def plot(outname=None):
    import pandas as pd
    with open("bench_saga.json", "r") as f:
        f = json.load(f)
    res = pd.DataFrame(f)
    res.set_index(["single_target"], inplace=True)
    grouped = res.groupby(level=["single_target"])

    # ① 颜色与样式映射
    colors = {"saga": "C0", "liblinear": "C1", "lightning": "C2"}
    linestyles = {"float32": "--", "float64": "-"}
    alpha = {"float64": 0.5, "float32": 1}

    for idx, group in grouped:
        single_target = idx
        # ② 4 子图布局
        fig, axes = plt.subplots(figsize=(12, 4), ncols=4)

        # 子图 1：训练目标随时间
        ax = axes[0]
        for scores, times, solver, dtype in zip(
            group["train_scores"], group["times"], group["solver"], group["dtype"]
        ):
            ax.plot(times, scores, label="%s - %s" % (solver, dtype),
                    color=colors[solver], alpha=alpha[dtype],
                    marker=".", linestyle=linestyles[dtype])
        ax.set_xlabel("Time (s)"); ax.set_ylabel("Training objective (relative to min)")
        ax.set_yscale("log")

        # 子图 2：测试目标随时间（结构与子图 1 类似）
        ax = axes[1]
        # ... 相同的循环结构 ...

        # 子图 3：测试准确率随时间
        ax = axes[2]
        for accuracy, times, solver, dtype in zip(...):
            ax.plot(times, accuracy, label="%s - %s" % (solver, dtype), ...)
        ax.set_xlabel("Time (s)"); ax.set_ylabel("Test accuracy")
        ax.legend()

        # 子图 4：目标随迭代次数
        ax = axes[3]
        for scores, times, solver, dtype in zip(...):
            ax.plot(np.arange(len(scores)), scores, label="%s - %s" % ..., ...)
        ax.set_yscale("log"); ax.set_xlabel("# iterations")
        ax.set_ylabel("Objective function"); ax.legend()

        plt.suptitle(name)
        plt.savefig(outname)

四子图的设计让 SAGA 求解器的性能特征无所遁形：

• 子图 1（训练目标-时间）：谁优化得最快？

• 子图 2（测试目标-时间）：谁泛化得最稳？

• 子图 3（测试准确率-时间）：谁的预测最准？

• 子图 4（目标-迭代次数）：谁的迭代效率最高（每轮的进步）？

颜色（colors）和线型（linestyles）的统一映射让多 solver × 多 dtype 的对比清晰可读。

88.7 t-SNE 大规模基准 —— 嵌入质量的"保真度检验"

t-SNE 是高维数据可视化的明星算法，但它的时间复杂度是 O(n²)，难以扩展到大数据集。本节基准测试回答两个核心问题：（1）sklearn 的 t-SNE 与参考实现 bhtsne 性能差距多大？（2）PCA 预处理是否值得？

88.7.1 数据加载与内存映射

源码路径：benchmarks/bench_tsne_mnist.py - load_data()（54-68行）

@memory.cache
def load_data(dtype=np.float32, order="C", shuffle=True, seed=0):
    """Load the data, then cache and memmap the train/test split"""
    print("Loading dataset...")
    data = fetch_openml("mnist_784", as_frame=True)

    # ① check_array 做类型与内存顺序转换
    X = check_array(data["data"], dtype=dtype, order=order)
    y = data["target"]

    if shuffle:
        X, y = _shuffle(X, y, random_state=seed)

    # ② 归一化到 [0, 1]
    X /= 255
    return X, y

@memory.cache 装饰器让首次下载后永久缓存到磁盘；mmap_mode="r" 让 joblib 在缓存时使用内存映射，避免将整个数据集一次性加载到内存。

88.7.2 嵌入质量评估：nn_accuracy

源码路径：benchmarks/bench_tsne_mnist.py - nn_accuracy()（37-45行）

def nn_accuracy(X, X_embedded, k=1):
    """Accuracy of the first nearest neighbor"""
    # ① 并行搜索原始空间 1-NN
    knn = NearestNeighbors(n_neighbors=1, n_jobs=-1)
    _, neighbors_X = knn.fit(X).kneighbors()
    # ② 并行搜索嵌入空间 1-NN
    _, neighbors_X_embedded = knn.fit(X_embedded).kneighbors()
    # ③ 索引一致率 = 保真度指标
    return np.mean(neighbors_X == neighbors_X_embedded)

nn_accuracy 是 t-SNE 基准的灵魂指标。它的逻辑是：在原始 50 维空间中，每个点的最近邻是谁；在嵌入 2 维空间中，每个点的最近邻又是谁；两者一致的比率就是保真度。这个指标直接量化了 t-SNE 是否"扭曲"了局部邻域结构。

88.7.3 主实验流程

源码路径：benchmarks/bench_tsne_mnist.py - __main__()（70-150行）

if __name__ == "__main__":
    # ① 命令行参数解析
    parser = argparse.ArgumentParser("Benchmark for t-SNE")
    parser.add_argument("--order", type=str, default="C", help="Order of the input data")
    parser.add_argument("--perplexity", type=float, default=30)
    parser.add_argument("--bhtsne", action="store_true", ...)
    parser.add_argument("--all", action="store_true",
                        help="if set, run with the whole MNIST. Note it will take up to 1 hour.")
    parser.add_argument("--profile", action="store_true",
                        help="if set, run with a memory profiler.")
    parser.add_argument("--verbose", type=int, default=0)
    parser.add_argument("--pca-components", type=int, default=50,
                        help="Number of principal components for preprocessing.")
    args = parser.parse_args()

    print("Used number of threads: {}".format(_openmp_effective_n_threads()))
    X, y = load_data(order=args.order)

    # ② PCA 预处理（可选）
    if args.pca_components > 0:
        t0 = time()
        X = PCA(n_components=args.pca_components).fit_transform(X)
        print("PCA preprocessing down to {} dimensions took {:0.3f}s".format(
            args.pca_components, time() - t0))

    methods = []
    # ③ 注册 sklearn TSNE
    tsne = TSNE(n_components=2, init="pca", perplexity=args.perplexity,
                verbose=args.verbose, n_iter=1000)
    methods.append(("sklearn TSNE", lambda data: tsne_fit_transform(tsne, data)))

    # ④ 注册 bhtsne（可选）
    if args.bhtsne:
        try:
            from bhtsne.bhtsne import run_bh_tsne
        except ImportError as e:
            raise ImportError(...) from e

        def bhtsne(X):
            """Wrapper for the reference lvdmaaten/bhtsne implementation."""
            # PCA 预处理已在主流程完成
            n_iter = -1  # TODO find a way to report the number of iterations
            return (
                run_bh_tsne(X, use_pca=False, perplexity=args.perplexity,
                            verbose=args.verbose > 0),
                n_iter,
            )
        methods.append(("lvdmaaten/bhtsne", bhtsne))

    if args.profile:
        try:
            from memory_profiler import profile
            methods = [(n, profile(m)) for n, m in methods]
        except ImportError as e:
            raise ImportError(...) from e

    # ⑤ 数据规模扫描
    data_size = [100, 500, 1000, 5000, 10000]
    if args.all:
        data_size.append(70000)  # 加上全量数据

    results = []
    basename = os.path.basename(os.path.splitext(__file__)[0])
    log_filename = os.path.join(LOG_DIR, basename + ".json")

    # ⑥ 主循环：遍历规模 × 遍历方法
    for n in data_size:
        X_train = X[:n]
        y_train = y[:n]
        n = X_train.shape[0]
        for name, method in methods:
            print("Fitting {} on {} samples...".format(name, n))
            t0 = time()
            # ⑦ 持久化原始数据与标签
            np.save(os.path.join(LOG_DIR, "mnist_{}_{}.npy".format("original", n)), X_train)
            np.save(os.path.join(LOG_DIR, "mnist_{}_{}.npy".format("original_labels", n)), y_train)
            # ⑧ 执行 t-SNE
            X_embedded, n_iter = method(X_train)
            duration = time() - t0
            # ⑨ 计算保真度
            precision_5 = nn_accuracy(X_train, X_embedded)
            print("Fitting {} on {} samples took {:.3f}s in {:d} iterations, "
                  "nn accuracy: {:0.3f}".format(name, n, duration, n_iter, precision_5))
            # ⑩ 记录结果
            results.append(dict(method=name, duration=duration, n_samples=n))
            with open(log_filename, "w", encoding="utf-8") as f:
                json.dump(results, f)
            # ⑪ 持久化嵌入结果
            method_name = sanitize(name)
            np.save(op.join(LOG_DIR, "mnist_{}_{}.npy".format(method_name, n)), X_embedded)

这段代码展示了工程级基准脚本的最佳实践：

第 8-13 行（PCA 预处理）：默认将 784 维降至 50 维，显著加速 t-SNE。这是工程经验值。

第 35 行（数据规模扫描）：默认 [100, 500, 1000, 5000, 10000]，加 --all 后追加 70000。这种"小到大"的扫描模式能清晰展示算法的扩展性曲线。

第 49-50 行（实时持久化）：每次实验结束立刻写入 JSON 与 .npy。中途中断后，下次运行可以从断点继续（虽然代码没有自动断点续跑逻辑，但持久化数据不会丢失）。

第 45 行的 n_iter = -1：bhtsne 的 C++ 二进制接口无法报告迭代次数，这是一个工程局限。

88.7.4 可视化脚本 plot_tsne_mnist.py

源码路径：benchmarks/plot_tsne_mnist.py - __main__()（1-30行）

import argparse
import os.path as op
import matplotlib.pyplot as plt
import numpy as np

LOG_DIR = "mnist_tsne_output"

if __name__ == "__main__":
    parser = argparse.ArgumentParser("Plot benchmark results for t-SNE")
    parser.add_argument("--labels", type=str,
                        default=op.join(LOG_DIR, "mnist_original_labels_10000.npy"),
                        help="1D integer numpy array for labels")
    parser.add_argument("--embedding", type=str,
                        default=op.join(LOG_DIR, "mnist_sklearn_TSNE_10000.npy"),
                        help="2D float numpy array for embedded data")
    args = parser.parse_args()

    # ① 加载嵌入与标签
    X = np.load(args.embedding)
    y = np.load(args.labels)

    # ② 按类别散点图，10 个数字 10 种颜色
    for i in np.unique(y):
        mask = y == i
        plt.scatter(X[mask, 0], X[mask, 1], alpha=0.2, label=int(i))
    plt.legend(loc="best")
    plt.show()

这段代码简洁明了——它把基准测试产出的 .npy 嵌入结果与标签加载后，按数字类别分色散点展示。alpha=0.2 让密集区域透明度叠加，更清晰展示聚类结构。

88.8 多项式核近似基准 —— CountSketch 与 Nystroem 的"近似比拼"

核方法（如 SVM with polynomial kernel）能捕捉非线性关系，但计算核矩阵的开销是 O(n²)，难以扩展。显式特征映射（Explicit Feature Map Approximation）通过将核函数近似为低维特征空间中的点积，把核 SVM 转化为线性 SVM，从而用 O(n) 的复杂度处理大规模数据。本节基准测试两种代表性方法：PolynomialCountSketch（数据无关）与 Nystroem（数据相关）。

88.8.1 数据准备与基准线

源码路径：benchmarks/bench_plot_polynomial_kernel_approximation.py - 数据加载与基准（33-48行）

# 第 88 章 —— ① 加载 digits 数据集
X, y = load_digits()["data"], load_digits()["target"]
# 第 88 章 —— ② 70/30 划分
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.7)

# 第 88 章 —— ③ 设置 n_components 扫描范围
out_dims = range(20, 400, 20)

# 第 88 章 —— ④ 基准 1：原始空间上的线性 SVM（下界）
lsvm = LinearSVC().fit(X_train, y_train)
lsvm_score = 100 * lsvm.score(X_test, y_test)

# 第 88 章 —— ⑤ 基准 2：多项式核 SVM（精度上界）
ksvm = SVC(kernel="poly", degree=2, gamma=1.0).fit(X_train, y_train)
ksvm_score = 100 * ksvm.score(X_test, y_test)

两条基准线的设计非常关键：

• 线性 SVM：告诉我们"不做核近似，直接用线性模型"的效果。

• 核 SVM：告诉我们"完整核方法的精度上界"。

后续的近似方法（如 PS + LinearSVC）应该位于这两条线之间，越接近核 SVM 越好。

88.8.2 PolynomialCountSketch 精度扫描

源码路径：benchmarks/bench_plot_polynomial_kernel_approximation.py - PS 精度测试（68-152行）

ps_svm_scores = []
n_runs = 5    # ① 多次运行取平均抵消随机性

# 第 88 章 —— ② 遍历 n_components
for k in out_dims:
    score_avg = 0
    for _ in range(n_runs):
        # ③ Pipeline 组合：PS 近似 + 线性 SVM
        ps_svm = Pipeline([
            ("PS", PolynomialCountSketch(degree=2, n_components=k)),
            ("SVM", LinearSVC()),
        ])
        score_avg += ps_svm.fit(X_train, y_train).score(X_test, y_test)
    ps_svm_scores.append(100 * score_avg / n_runs)

n_runs=5 是抗随机性的关键——PolynomialCountSketch 内部使用随机哈希，多次运行取平均能消除单次实验的偶然性。Pipeline 巧妙地把"特征映射 + 线性 SVM"组合成一个整体，自动完成 fit/transform 的衔接。

88.8.3 Nystroem 精度扫描

源码路径：benchmarks/bench_plot_polynomial_kernel_approximation.py - Nystroem 精度测试（155-190行）

ny_svm_scores = []
n_runs = 5

for k in out_dims:
    score_avg = 0
    for _ in range(n_runs):
        # Nystroem 参数：多项式核，degree=2, gamma=1, coef0=0
        ny_svm = Pipeline([
            ("NY", Nystroem(kernel="poly", gamma=1.0, degree=2, coef0=0, n_components=k)),
            ("SVM", LinearSVC()),
        ])
        score_avg += ny_svm.fit(X_train, y_train).score(X_test, y_test)
    ny_svm_scores.append(100 * score_avg / n_runs)

Nystroem 与 PS 的协议完全一致——同样跑 5 次取平均——让两者对比绝对公平。

88.8.4 精度对比图

源码路径：benchmarks/bench_plot_polynomial_kernel_approximation.py - 绘图（70-152行尾部）

fig, ax = plt.subplots(figsize=(6, 4))
ax.set_title("Accuracy results")
# 第 88 章 —— ① PS 曲线（橙色）
ax.plot(out_dims, ps_svm_scores, label="PolynomialCountSketch + linear SVM", c="orange")
# 第 88 章 —— ② Nystroem 曲线（蓝色）
ax.plot(out_dims, ny_svm_scores, label="Nystroem + linear SVM", c="blue")
# 第 88 章 —— ③ 线性 SVM 基准（黑色虚线，下界）
ax.plot([out_dims[0], out_dims[-1]], [lsvm_score, lsvm_score],
        label="Linear SVM", c="black", dashes=[2, 2])
# 第 88 章 —— ④ 核 SVM 基准（红色虚线，上界）
ax.plot([out_dims[0], out_dims[-1]], [ksvm_score, ksvm_score],
        label="Poly-kernel SVM", c="red", dashes=[2, 2])
ax.legend()
ax.set_xlabel("N_components for PolynomialCountSketch and Nystroem")
ax.set_ylabel("Accuracy (%)")
fig.tight_layout()

四条曲线在一张图上清晰呈现：PS 与 Nystroem 应在两条虚线之间上升，越接近核 SVM 虚线越好。

88.8.5 可扩展性实验

源码路径：benchmarks/bench_plot_polynomial_kernel_approximation.py - 扩展性测试（155-190行）

# 第 88 章 —— ① 生成大规模随机数据：10000 样本 × 100 特征
fakeData = np.random.randn(10000, 100)
fakeDataY = np.random.randint(0, high=10, size=(10000))

# 第 88 章 —— ② 更大的 n_components 范围
out_dims = range(500, 6000, 500)

# 第 88 章 —— ③ 测量 PolynomialCountSketch 的 fit_transform 时间
ps_svm_times = []
for k in out_dims:
    ps = PolynomialCountSketch(degree=2, n_components=k)
    start = time()
    ps.fit_transform(fakeData, None)
    ps_svm_times.append(time() - start)

# 第 88 章 —— ④ 测量 Nystroem 的 fit_transform 时间
ny_svm_times = []
for k in out_dims:
    ny = Nystroem(kernel="poly", gamma=1.0, degree=2, coef0=0, n_components=k)
    start = time()
    ny.fit_transform(fakeData, None)
    ny_svm_times.append(time() - start)

# 第 88 章 —— ⑤ 绘图
fig, ax = plt.subplots(figsize=(6, 4))
ax.set_title("Scalability results")
ax.plot(out_dims, ps_svm_times, label="PolynomialCountSketch", c="orange")
ax.plot(out_dims, ny_svm_times, label="Nystroem", c="blue")
ax.legend()
ax.set_xlabel("N_components for PolynomialCountSketch and Nystroem")
ax.set_ylabel("fit_transform time (s/10.000 samples)")
fig.tight_layout()
plt.show()

扩展性实验的设计非常考究：

• 数据规模：10000 × 100。10000 样本让差异可见，100 特征避免维度灾难。

• n_components 范围：500 到 6000。这是 PS 优势开始显现的范围。

• 计时对象：fit_transform。包含 PS 的初始化与变换、Nystroem 的训练与变换。

理论复杂度对比：

| 方法 | 复杂度 | 训练开销 | 数据相关性 |

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

| PolynomialCountSketch | O(n(d + k log k)) | 几乎无（仅随机变量初始化） | 数据无关 |

| Nystroem | O(n(dk + k²)) | 较高（需要数据子集 SVD） | 数据相关 |

在 out_dims=500..6000 区间，PS 的 O(n(d + k log k)) 接近线性增长，而 Nystroem 的 O(n(dk + k²)) 呈超线性。当 k 较小时，Nystroem 的常数项可能更小，因此初期可能更快；但随着 k 增大，PS 的优势将逐渐显现。

下面这张图揭示了多项式核近似基准的双实验范式：

graph TD A[digits 数据 70/30 切分] --> B[两条基准线] B --> B1[LinearSVC 原始空间] B --> B2[SVC poly 核空间] A --> C[精度实验] C --> C1[PS: 5 次平均准确率] C --> C2[Nystroem: 5 次平均准确率] C --> D[精度对比图 4 曲线] A --> E[扩展性实验] E --> E1[fakeData 10k x 100] E1 --> F[PS: fit_transform 计时] E1 --> G[Nystroem: fit_transform 计时] F --> H[扩展性对比图] G --> H

88.9 设计中的取舍

• 为什么 NMF 基准不复用 sklearn 内置的 Projected Gradient，而是自定义 _PGNMF 类？

  scikit-learn 自 0.19 版本起移除了 Projected Gradient 求解器（参见 _PGNMF 类注释第 23-25 行 "This class is private and for comparison purpose only. It may change or disappear without notice"）。移除的原因是该算法的实现复杂度高且收敛速度并不显著优于 CD 与 MU，但作为学术对比的经典参考，它仍有保留价值。_PGNMF 通过继承 NMF 基类复用参数管理与初始化基础设施，仅重写求解逻辑，避免了重复造轮子。基准测试的本质是"在受控环境中观察三种求解器的真实差异"，自定义 PG 实现正是为了这一目的。

• 逻辑回归基准的 plot_dloss 取所有求解器最优值作为参考，而不是各自的最优值，这样设计的原因是什么？

  plot_dloss 的核心目的是揭示各求解器的收敛阶数——线性收敛与超线性收敛在对数残差图上有截然不同的斜率。如果每个求解器用各自的最优值作为参考，每条曲线都从 0 开始，无法看出"谁先达到更小的残差"。统一参考线后，曲线越陡峭下降意味着该求解器越接近全局最优。多求解器共享同一参考线还能让它们在终点处的差异清晰可见——某求解器曲线"提前进入平台"，说明它最终停留在了一个比最优解差的局部极小。

• SAGA 基准为何使用 time.clock() 而非 time.time()？为什么需要 X_train.max() 预热？

  time.clock()（在 Python 3.3+ 已被 time.process_time() 取代）测量的是进程占用的 CPU 时间，而非挂钟时间。SAGA 基准对比的是不同求解器的计算效率，因此 CPU 时间更准确（不受系统其他进程干扰）。X_train.max() 预热是为了让数据进入 CPU 缓存，避免后续 fit 第一次迭代时把数据加载时间算进计时中。这种细节是高质量基准测试的标志——只有排除所有干扰因素，才能保证对比的公平性。

• t-SNE 基准为何将 PCA 预处理作为默认开启项？50 维是经验值吗？

  PCA 预处理的依据是 t-SNE 的局部保持特性：它只关心高维空间的局部邻域结构，对全局距离不敏感。784 维的 MNIST 数据中，50 个主成分已经能解释绝大部分方差（通常 >90%）。将维度从 784 降到 50 后，t-SNE 的 O(n² d) 计算量减少 15 倍以上，邻域搜索（K-D 树）的速度也显著加快。这是工程经验值——50 是 sklearn 文档明确推荐的默认值。基准脚本通过 --pca-components 参数允许用户覆盖这一选择。

• 多项式核近似基准为何设计成"双实验范式"（精度图 + 扩展性图）？

  单张图无法同时回答"精度如何"与"速度如何"两个问题。精度图在小维度范围（20-400）扫描，n_runs=5 取平均，主要观察"精度提升曲线"。扩展性图在大维度范围（500-6000）扫描，每次仅计时一次，主要观察"计算成本曲线"。两者结合才能完整评估 PolynomialCountSketch 与 Nystroem 的精度-效率权衡——精度图告诉你"近似效果如何"，扩展性图告诉你"训练成本如何"。这种"双视角"的实验设计是机器学习基准测试的通用范式。

88.10 动手练习

1. 深入 NMF 投影梯度求解器

   阅读 benchmarks/bench_plot_nmf.py 中 _nls_subproblem (85-136行) 与 _fit_projected_gradient (176-215行)：

   1. 解释 grad * np.logical_or(grad < 0, H > 0) 这一投影梯度范数的几何含义。

   2. 线搜索中 suff_decr = (1 - sigma) * gradd + 0.5 * dQd < 0 对应哪种经典线搜索条件？参数 sigma=0.01, beta=0.1 的典型作用是什么？

   3. _fit_projected_gradient 为何交替更新 W 和 H？tolW/tolH 动态缩小 0.1 倍的启发式依据是什么？

2. 复现与扩展逻辑回归收敛性基准

   基于 benchmarks/bench_rcv1_logreg_convergence.py：

   1. 修改脚本加入 solver='saga' 与 solver='liblinear' (dual=False) 在 l1 正则下的对比，观察目标函数下降曲线差异。

   2. 将 plot_dloss 中的 pobj_best 改为所有求解器最终目标值的最小值，而非当前最小值，分析对曲线形状的影响。

   3. 尝试在 bench_one 中记录 clf.n_iter_ (实际迭代数) 与设定 max_iter 的偏差，绘制"实际迭代数 vs 设定迭代数"散点图。

3. SAGA 基准的多维度分析

   阅读 benchmarks/bench_saga.py 中 fit_single 与 exp：

   1. 为什么 LogisticRegression 加上 OneVsRestClassifier 时 multi_class='ovr'，而原生 saga 支持 multinomial？这对目标函数计算有何影响？

   2. _predict_proba 为何仅用于 Lightning？实现 softmax(safe_sparse_dot(X, coef_.T) + intercept_) 的数学依据是什么？

   3. 实验设计中 n_samples=[100000, 300000, 500000, 800000, None] 为何包含 None？它对应 RCV1 全量数据多少样本？

4. t-SNE 嵌入质量与性能权衡

   分析 benchmarks/bench_tsne_mnist.py 与 plot_tsne_mnist.py：

   1. nn_accuracy 使用 k=1 最近邻，若改为 k=5 或 k=10 会如何改变质量评估的敏感度？

   2. PCA 预处理维度 pca_components=50 是经验值，设计实验对比 30/50/100/无 PCA 对最终 nn_accuracy 与运行时的影响。

   3. bhtsne 实现为何无法报告 n_iter？其 C++ 核心 bh_tsne 二进制接口限制是什么？

5. 多项式核近似：理论复杂度与实测对比

   结合 benchmarks/bench_plot_polynomial_kernel_approximation.py 两个实验段：

   1. 精度实验中 n_runs=5 取平均，PolynomialCountSketch 的方差来源是什么？Nystroem 是否也有随机性？

   2. 扩展性实验中 fakeData = np.random.randn(10000, 100) 为何选此规模？若改为稀疏矩阵，两者相对性能会如何反转？

   3. 理论复杂度 PS: O(n(d + k log k)) vs Nystroem: O(n(dk + k²))，在 out_dims=500..6000 下实测斜率比是否符合预期？为何 Nystroem 曲线在小 k 时可能更快？

88.11 本章小结

这一章我们深入剖析了五个独立基准脚本，从自定义求解器复活（_PGNMF）到大规模 t-SNE 评估，从逻辑回归收敛性的四视角观察到 SAGA 多数据集全景透视，再到多项式核近似的双实验范式。每一个脚本都展示了"算法奥林匹克"裁判系统的不同侧面。

首先，我们学习了 NMF 自定义投影梯度求解器（_PGNMF），它通过继承 scikit-learn 内置 NMF 类复用基础设施，仅重写求解逻辑，让已移除的算法重新站上擂台。其次，我们掌握了逻辑回归收敛性基准的四视角可视化——train_loss、train_score、test_score、log|p-p*|——全方位剖析求解器的优化行为与泛化能力。接着，我们了解了 SAGA 多数据集基准的统一目标函数度量，无论底层是 scikit-learn 还是 lightning，score 计算始终一致。然后，我们理解了 t-SNE 大规模基准中的核心质量指标 nn_accuracy（最近邻保持率）与 PCA 预处理的加速价值。最后，我们掌握了多项式核近似的精度-效率权衡——PolynomialCountSketch 的数据无关特性让 O(n(d+k log k)) 接近线性增长，而 Nystroem 的数据相关性导致 O(n(dk + k²)) 呈超线性。

| 概念 | 解释 |

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

| _PGNMF | 复现已移除的投影梯度 NMF 求解器，用于对比 CD/MU 求解器在稀疏/稠密数据上的收敛行为 |

| _nls_subproblem | 非负最小二乘子问题求解器，采用投影梯度下降 + Armijo 线搜索，含投影步骤与充分下降条件 |

| run_bench / bench_one | 基准测试主控流程：网格搜索求解器×初始化×迭代数，joblib 缓存单次实验，聚合 loss/time 绘图 |

| get_loss (logreg) | 逻辑损失 + L2 正则项：log(1+exp(-y(w·x+b))) + ||w||²/(2Cn)，统一度量不同求解器目标值 |

| 收敛性四视图 | train_loss/time、train_score/time、test_score/time、log(p-p*)/time —— 全方位剖析求解器动态 |

| fit_single (SAGA) | 单实验原子单元：统一 penalty→alpha/beta 映射、适配 Lightning/sklearn 接口、记录迭代轨迹 |

| exp + Parallel | 实验调度器：跨数据集/精度/求解器并行执行，结果持久化 JSON 供绘图复用 |

| nn_accuracy | t-SNE 嵌入质量核心指标：原始空间与嵌入空间 1-NN 索引一致率，量化局部结构保真度 |

| PCA 预处理 + 规模扫描 | 先降维再 t-SNE，从 100 到 70k 样本递增，对比 sklearn 与 bhtsne 运行时与迭代数 |

| PolynomialCountSketch vs Nystroem | 两类多项式核显式映射：PS 数据无关 O(n(d+klog k))、Nystroem 数据相关 O(n(dk+k²))，精度-速度权衡 |

| 双图评估范式 | 精度图(小维度扫描+多次平均) + 扩展性图(大维度扫描+计时) —— 覆盖实用与极限两种场景 |

下一章中，我们将学习并行成对距离基准 —— 释放"多核算力"，通过 bench_plot_parallel_pairwise.py 剖析欧氏距离与 RBF 核计算在单核与多核下的加速比，展示 sklearn 并行基础设施的性能收益。

第 89 章 —— 并行距离与核计算 —— 多核时代的"加速引擎"

89.1 学习目标

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

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

• 理解 pairwise_distances 与 pairwise_kernels 的并行化实现机制

• 掌握使用 n_jobs 参数控制并行度（单核 vs 全核）的性能基准设计方法

• 了解 joblib.Parallel 在成对距离计算与核计算中的任务分发策略

• 能分析内存带宽密集型（欧氏距离）与计算密集型（RBF核）任务的并行扩展性差异

• 能编写基准脚本量化并行加速比随数据规模的变化趋势

89.2 生活类比

我们可以把并行成对距离计算想象成一场"全员核对身份证"的马拉松赛事。每位参赛选手对应一个样本，核对的不是真实的证件，而是高维特征向量之间的两两差异。当选手们两两核对身份证号码的差异时，欧氏距离的计算几乎全在搬运和比对数据——核验本身很简单，难点在于把成千上万条记录读出来再算差值，因此它属于内存带宽密集型操作；而 RBF 核的核对则需要在此基础上进一步计算复杂的指纹匹配分数，也就是 exp(-γ * ||x-y||²) 的指数运算，CPU 必须深度参与大量的浮点计算，因此它属于计算密集型任务。如果我们设置 n_jobs=1，就如同指派一个裁判独自处理所有选手的核对工作，过程严格串行，速度较慢但结果可预测；如果设置 n_jobs=-1，则相当于调动全体裁判分组上场，每个核负责一部分比较任务，理论速度会成倍提升。然而现实中的运动会终究受限于身份证发放速度，对应到硬件层面，欧氏距离的并行加速比往往受限于内存子载速率而非 CPU 的运算能力；RBF 核则因计算密度高，通常能获得更接近线性的并行加速收益。这场"多核运动会"贯穿两个核心维度：当样本规模从 1000 递增到 5000 时，n_jobs=1 与 n_jobs=-1 的耗时曲线如何分叉，正是观察"发放速度瓶颈"与"裁判人手瓶颈"差异的最佳窗口：发放速度瓶颈类似于身份证读卡机的吞吐上限，裁判人手瓶颈则像签核台的 CPU 计算周期。两种瓶颈决定了赛事最终的提速空间，也决定了我们能在这场多核运动会中收获多少回报。

89.3 源码地图

graph TD A[脚本入口<br/>bench_plot_parallel_pairwise.py] --> B[plot func 调用<br/>euclidean_distances] A --> C[plot func 调用<br/>rbf_kernels] B --> D[plot func<br/>14-34行] C --> D D --> E[check_random_state 0<br/>固定随机种子] D --> F[sample_sizes 循环<br/>1000~5000 步长 1000] D --> G[euclidean_distances 包装<br/>36-38行] D --> H[rbf_kernels 包装<br/>40-42行] G --> I[pairwise_distances<br/>metric euclidean] H --> J[pairwise_kernels<br/>metric rbf gamma 0.1] D --> K[time.time 计时] D --> L[matplotlib.pyplot 绘图] A --> M[plt.show 显示]

源码地图本身即是理解脚本结构的最佳架构图。从 check_random_state(0) 出发，经 sample_sizes 循环，到 matplotlib.pyplot 收尾，构成了一个完整的"数据生成 → 计时执行 → 结果可视化"流水线。

89.4 并行距离与核计算 —— 多核时代的"加速引擎"

本章围绕 pairwise_distances 与 pairwise_kernels 在多核环境下的并行加速机制展开，覆盖三个核心要点：单核与全核的对比基线设计、控制变量法的基准流程、以及两类任务在并行扩展性上的本质差异。

为什么要对比 n_jobs=1 与 n_jobs=-1？n_jobs=1 代表单线程执行，是我们评估并行收益的性能基线。n_jobs=-1 则代表调用全部可用 CPU 核心，用于展示 pairwise_distances 与 pairwise_kernels 在多核环境下的加速效果。通过两者的耗时之比（即加速比），我们能够量化这两个函数在不同数据规模下的并行可扩展性（scalability）。这个比值如果接近 CPU 核心数，意味着近似线性扩展；如果远低于核心数，则说明存在内存带宽、任务调度或锁竞争等瓶颈。基准脚本通过将这两个参数封装在 plot 函数的统一循环中，使得每一档样本量都能直接产出一对可比数据点。

核心基准流程是如何设计的？基准脚本采用控制变量法：固定特征维度为 300（中等规模），仅逐步增加样本量（从 1000 到 5000，步长 1000）。这种设计思路是隔离"数据规模增长"这一关键变量对并行效率的影响。使用 check_random_state(0) 固定随机种子，确保所有实验基于相同的数据分布，结果具备可复现性。脚本分别测试两类代表性任务：欧氏距离（内存带宽敏感型）和 RBF 核（计算密集型），最终以折线图形式呈现加速比随样本量的变化趋势。这种"固定维度、变化样本量"的设计是 scikit-learn 基准测试的经典范式，便于横向对比不同 metric 的扩展曲线。

pairwise_distances 与 pairwise_kernels 的并行化机制有何异同？两者均通过 n_jobs 参数控制并行度，底层依赖 joblib.Parallel 实现任务分发。然而，两者的性能瓶颈截然不同。欧氏距离计算涉及大量数据搬运（O(N²) 次内存访问），并行加速比通常受限于内存带宽——当所有核都在"等待数据搬运"时，再多的核也无法提升速度。RBF 核计算包含 exp(-γ * ||x-y||²) 的指数运算，计算密度高，每个数据点需要更多 CPU 周期处理，因此通常能获得更接近线性的并行加速比。这正是脚本挑选这两类任务作为基准对象的原因：它们恰好代表了并行计算的两种典型瓶颈场景。

graph TD A[plot func 启动] --> B[生成随机数据 X] B --> C{选择并行度} C -->|n_jobs=1| D[单核顺序执行] C -->|n_jobs=-1| E[joblib.Parallel 任务分发] E --> F[工作进程池] F --> G[分块计算距离/核矩阵] G --> H[归约结果] D --> I[记录单核耗时] H --> J[记录多核耗时] I --> K[绘制对比折线图] J --> K

89.4.1 核心基准驱动函数 plot()

源码路径：benchmarks/bench_plot_parallel_pairwise.py - plot()（14-34行）

def plot(func):                                               # ① 接受一个待测试的成对计算函数作为参数
    random_state = check_random_state(0)                      # ② 调用 sklearn 工具方法，传入固定种子 0，确保可复现
    one_core = []                                             # ③ 初始化空列表，用于存储 n_jobs=1 下的耗时记录
    multi_core = []                                           # ④ 初始化空列表，用于存储 n_jobs=-1 下的耗时记录
    sample_sizes = range(1000, 6000, 1000)                     # ⑤ 生成样本量梯度，从 1000 到 5000，步长为 1000

    for n_samples in sample_sizes:                             # ⑥ 遍历每一档样本量，进入主测试循环
        X = random_state.rand(n_samples, 300)                 # ⑦ 调用随机状态对象，生成形状为 (n_samples, 300) 的均匀分布矩阵

        start = time.time()                                   # ⑧ 调用 time.time 获取当前墙钟时间，作为单核测试起点
        func(X, n_jobs=1)                                     # ⑨ 调用传入的待测试函数，强制 n_jobs=1 单线程串行执行
        one_core.append(time.time() - start)                  # ⑩ 计算耗时差值并追加到单核耗时列表

        start = time.time()                                   # ⑪ 再次记录当前墙钟时间，作为多核测试起点
        func(X, n_jobs=-1)                                    # ⑫ 调用传入的待测试函数，启用 n_jobs=-1 自动使用全部 CPU 核心
        multi_core.append(time.time() - start)                # ⑬ 计算耗时差值并追加到多核耗时列表

    plt.figure("scikit-learn parallel %s benchmark results" % func.__name__)  # ⑭ 创建 matplotlib 画布，标题动态拼接函数名
    plt.plot(sample_sizes, one_core, label="one core")        # ⑮ 绘制单核耗时随样本量变化的折线，标签为 one core
    plt.plot(sample_sizes, multi_core, label="multi core")    # ⑯ 绘制多核耗时随样本量变化的折线，标签为 multi core
    plt.xlabel("n_samples")                                   # ⑰ 设置 X 轴标签为 n_samples
    plt.ylabel("Time (s)")                                    # ⑱ 设置 Y 轴标签为 Time s，单位为秒
    plt.title("Parallel %s" % func.__name__)                  # ⑲ 设置图表标题，动态拼接函数名
    plt.legend()                                              # ⑳ 调用 legend 显示图例，区分 one core 与 multi core 两条曲线

这段代码定义了基准驱动函数 plot(func)，它接受一个成对计算函数（euclidean_distances 或 rbf_kernels），通过固定随机种子生成 5 个不同规模的随机数据集（样本量 1000~5000），分别测试单核（n_jobs=1）与多核（n_jobs=-1）模式下的墙钟耗时，最终用 matplotlib 绘制对比折线图。函数内部采用"先单核、后多核"的固定顺序计时，避免并行分支的预热效应污染单核基线。

graph LR A[plot func] --> B[check_random_state 0] B --> C[初始化 one_core 与 multi_core 列表] C --> D[循环 sample_sizes] D --> E[生成随机数据 X] E --> F[n_jobs 1 单核计时] F --> G[n_jobs -1 多核计时] G --> H{遍历完所有样本量?} H -->|否| D H -->|是| I[matplotlib 绘制折线图] I --> J[设置标签与图例]

89.4.2 欧氏距离包装函数 euclidean_distances()

源码路径：benchmarks/bench_plot_parallel_pairwise.py - euclidean_distances()（36-38行）

def euclidean_distances(X, n_jobs):                           # ① 定义薄包装函数，参数为特征矩阵与并行度
    return pairwise_distances(X, metric="euclidean", n_jobs=n_jobs)  # ② 透传 n_jobs 至 sklearn 底层 pairwise_distances，显式指定欧氏距离度量

这段代码定义了一个薄包装函数 euclidean_distances，它直接调用 sklearn.metrics.pairwise.pairwise_distances 并显式指定 metric="euclidean"。n_jobs 参数从外层 plot 函数透传至底层，触发 joblib.Parallel 的任务分发逻辑。函数体仅一行，体现了基准脚本"只测基础设施、不改业务逻辑"的纯净设计。

graph LR A[euclidean_distances 包装] --> B[调用 pairwise_distances] B --> C[metric euclidean] B --> D[n_jobs 透传] C --> E[joblib.Parallel 分发] D --> E E --> F[返回距离矩阵]

89.4.3 RBF 核包装函数 rbf_kernels()

源码路径：benchmarks/bench_plot_parallel_pairwise.py - rbf_kernels()（40-42行）

def rbf_kernels(X, n_jobs):                                  # ① 定义薄包装函数，参数为特征矩阵与并行度
    return pairwise_kernels(X, metric="rbf", n_jobs=n_jobs, gamma=0.1)  # ② 透传 n_jobs 至底层 pairwise_kernels，指定 RBF 度量与带宽参数 0.1

这段代码定义了 rbf_kernels 包装函数，调用 pairwise_kernels 并设置 metric="rbf"、gamma=0.1（径向基函数带宽参数的常用默认值）。它代表了计算密集型的并行任务，与欧氏距离形成对比。两者的函数结构完全对称，唯一差异在于 metric 与 gamma 参数的取值，这种对称设计便于横向对比加速曲线。

graph LR A[rbf_kernels 包装] --> B[调用 pairwise_kernels] B --> C[metric rbf] B --> D[gamma 0.1] B --> E[n_jobs 透传] C --> F[joblib.Parallel 分发] D --> F E --> F F --> G[exp 指数运算] G --> H[返回核矩阵]

89.4.4 脚本主流程（隐式 __main__）

源码路径：benchmarks/bench_plot_parallel_pairwise.py - 主流程（44-46行）

plot(euclidean_distances)                                     # ① 首次调用 plot 函数，将欧氏距离包装函数传入作为被测对象
plot(rbf_kernels)                                             # ② 再次调用 plot 函数，将 RBF 核包装函数传入作为被测对象
plt.show()                                                    # ③ 调用 matplotlib 的 show 方法，阻塞弹出两幅独立的折线图

这段代码以隐式 __main__ 的形式在模块底部执行主流程：依次调用 plot 函数测试两类任务，最终通过 plt.show() 弹出两幅独立的折线图。这种设计允许用户直观对比欧氏距离与 RBF 核在不同样本量下的并行加速效果。模块级调用使得脚本作为独立入口直接可执行，无需 if __name__ == "__main__" 保护。

graph TD A[模块加载完成] --> B[plot 调用 euclidean_distances] A --> C[plot 调用 rbf_kernels] B --> D[生成图1 欧氏距离对比] C --> E[生成图2 RBF核对比] D --> F[plt.show 阻塞显示] E --> F

89.5 设计中的取舍

为什么选择欧氏距离与 RBF 核作为基准对象？ 脚本刻意挑选了欧氏距离与 RBF 核这两种典型代表进行对比：前者代表"内存带宽密集型"任务，后者代表"计算密集型"任务。选择这两种是因为它们在 scikit-learn 中应用最广泛，且具有截然不同的瓶颈特征，能够完整揭示并行基础设施在不同负载下的表现差异。如果对所有 metric 都进行基准测试，会产生大量冗余数据且难以聚焦核心结论，这正是"代表性 vs 全覆盖"的经典权衡。

为什么使用 time.time() 而非 time.perf_counter() 计时？ 前者返回挂钟时间，会受系统时钟调整影响；后者是单调时钟，更适合测量短时间间隔。但对于本脚本中秒级甚至更长的成对距离计算，time.time() 的精度已经足够，且能直观反映用户感知的"等待时间"，属于"实现简洁性 vs 计时精度"的权衡。

为什么选择维度 300 与样本量 1000~5000？ 300 是一个中等规模维度，介于低维（如 10）与高维（如 10000）之间。在这个维度下，欧氏距离的内存访问开销与 RBF 核的计算开销都足够显著，能够在 1000~5000 的样本量区间内产生可观测的耗时差异。若维度太低，两类任务都会过于迅速，难以区分并行效率；若维度太高，则会过早触达内存带宽瓶颈，无法完整观察扩展曲线，这是"观测灵敏度 vs 现实代表性"的取舍。

graph LR A[基准脚本设计取舍] --> B[代表性 vs 全覆盖] A --> C[实现简洁性 vs 计时精度] A --> D[观测灵敏度 vs 现实代表性] B --> B1[仅选欧氏与RBF] B --> B2[避免冗余数据] C --> C1[time.time 挂钟时间] C --> C2[time.perf_counter 单调时钟] D --> D1[维度 300 中等规模] D --> D2[样本量 1000~5000]

89.6 动手练习

1. 阅读并行基准脚本核心逻辑：阅读 benchmarks/bench_plot_parallel_pairwise.py 全文，回答三个问题：plot 函数中为何选择 range(1000, 6000, 1000) 作为样本量梯度？如果改为 range(100, 1100, 100) 会观测到什么现象？euclidean_distances 与 rbf_kernels 分别调用了 sklearn.metrics.pairwise 中的哪两个核心函数？它们的 n_jobs 参数是如何透传到底层的？脚本使用 time.time() 而非 time.perf_counter() 计时，对于并行基准测试可能带来什么影响？

2. 动手修改基准脚本探究扩展性：修改 benchmarks/bench_plot_parallel_pairwise.py 进行三项实验：新增 manhattan_distances 基准（调用 pairwise_distances(metric='manhattan')），对比其与欧氏距离的并行加速比差异；固定样本量为 5000，增加特征维度（如 100, 300, 1000, 3000），观测高维下并行效率变化，解释内存带宽压力如何影响加速比；将 n_jobs=-1 改为具体核心数（如 2, 4, 8），绘制加速比随核心数变化的曲线，判断是否存在超线性加速或饱和点。

3. 深入底层并行实现：阅读 sklearn/metrics/pairwise.py 与 sklearn/metrics/_pairwise_fast.pyx 源码，回答三个问题：pairwise_distances 与 pairwise_kernels 如何根据 metric 选择并行策略？哪些 metric 支持并行，哪些不支持？_pairwise_fast.pyx 中的 OpenMP 并行实现（如 parallel_pairwise）是如何划分工作负载的？是否存在负载不均问题？对比 joblib.Parallel（进程/线程池）与 Cython nogil + OpenMP（线程级并行）在成对计算中的优劣势，scikit-learn 采用哪种方式？

graph TD A[动手练习三层递进] --> B[第一层 脚本阅读] A --> C[第二层 修改实验] A --> D[第三层 源码深挖] B --> B1[理解 plot 函数设计] B --> B2[理解 n_jobs 透传] C --> C1[新增 manhattan 基准] C --> C2[高维扩展性测试] C --> C3[核心数梯度测试] D --> D1[pairwise.py 并行策略] D --> D2[_pairwise_fast.pyx OpenMP] D --> D3[joblib vs nogil 对比]

89.7 本章小结

本章我们围绕 pairwise_distances 与 pairwise_kernels 在多核环境下的并行加速机制展开，重点学习了脚本通过 plot 函数固定维度、变化样本量的基准设计思路；euclidean_distances 与 rbf_kernels 两个包装函数如何透传 n_jobs 参数触发 joblib.Parallel 的任务分发；欧氏距离作为内存带宽密集型任务与 RBF 核作为计算密集型任务在并行扩展性上的差异；从任务分发到结果归约的完整调度链；以及为何选择这两种代表性 metric 与 time.time() 计时方式的工程合理性。

下面用一张表格汇总本章涉及的关键概念与解释：

| 概念 | 解释 |

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

| bench_plot_parallel_pairwise.py | 并行成对距离/核计算性能基准脚本，量化多核加速收益 |

| plot(func) | 基准驱动函数，固定维度变化样本量，对比单核/多核墙钟时间并绘图 |

| euclidean_distances | 封装 pairwise_distances(metric='euclidean')，内存带宽敏感型并行任务 |

| rbf_kernels | 封装 pairwise_kernels(metric='rbf')，计算密集型并行任务，加速比通常更高 |

| n_jobs=1 vs -1 | 单线程基线 vs 全核并行，直观展示 joblib.Parallel 调度效果 |

| check_random_state(0) | 固定随机种子，确保实验数据分布一致，结果可复现 |

| 加速比随样本量变化 | 大样本量下多核收益显著，小样本量调度开销可能抵消并行优势 |

下一章中，我们将学习构建系统基础架构，理解 scikit-learn 如何打造一条稳定、可复现的"工程流水线"，从依赖管理到环境锁定的完整机制。

graph LR A[本章核心要点] --> B[基准设计思路] A --> C[n_jobs 透传机制] A --> D[扩展性差异分析] A --> E[调度链理解] A --> F[工程合理性] B --> B1[固定维度变化样本量] C --> C1[joblib.Parallel 任务分发] D --> D1[欧氏内存带宽敏感] D --> D2[RBF 计算密集] E --> E1[分发到归约全链路] F --> F1[代表性 metric 选取] F --> F2[time.time 计时策略]

第 90 章 —— 项目总览与环境搭建 —— 走进"scikit-learn 基准测试世界"

90.1 学习目标

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

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

• 理解 scikit-learn 基准测试套件的整体定位与作用

• 明确基准测试运行所需的核心依赖环境与最低版本约束

• 掌握项目源码结构、开发指南及社区协作流程的导航入口

• 建立从源码阅读到性能评测的全局认知框架

90.2 生活类比

想象 scikit-learn 基准测试套件是一座"性能奥林匹克体育场"。README.md 就是这座体育场的导览图与入场规则手册——上面印着徽章、赛事日程和观众须知，是观众与新运动员进入这座体育场的必经之路。核心依赖则好比比赛认证的标准器材：Python/NumPy/SciPy 是跑道、计时器与起跑枪，决定了基础赛道是否平整；joblib 与 threadpoolctl 则是裁判组与计分系统，确保多组选手同时比赛时不会相互干扰，并且能精确统计成绩。版本约束就像参赛资格审核标准，低于规定版本的"器材"根本无法入场参赛。安装指南则充当观众入场指引，pip 和 conda 就像体育场东西两侧的闸机，持有任一票据均可顺利入场。测试命令好比正式比赛的发令枪——一旦在终端敲下 pytest sklearn，全场数百个测试用例便同时起跑，争夺"代码健康"的最高荣誉。贡献指南则是运动员注册手册，详细规定动作姿势、号码牌样式与犯规罚则，避免新提交的 PR 因为格式问题被直接罚下场。就像奥运会需要统一的场地标准、器材认证和竞赛规则才能让各国选手公平较量，基准测试也需要固定的环境版本、统一的度量标准和可复现的运行流程，才能横向对比各算法"选手"的真实实力。

90.3 源码地图

本章节我们以 README.md 为切入点展开分析，它是整个基准测试套件的"导航地图"。根据细纲定义，本章仅解析 README.md 文件第1-50行的内容，对应于文件开头的徽章墙与版本常量定义部分。下面以纯文本树状结构展示该区域的代码骨架：

README.md
├── 项目元数据与徽章（第1-30行）
│   ├── 文件模式声明（第1行）
│   ├── 徽章引用声明（第3行）
│   ├── Azure 徽章定义（第5-6行）
│   ├── CircleCI 徽章定义（第8-9行）
│   ├── Codecov 徽章定义（第11-12行）
│   ├── Nightly wheels 徽章定义（第14-15行）
│   ├── Ruff 徽章定义（第17-18行）
│   ├── PythonVersion 徽章定义（第20-21行）
│   ├── PyPI 徽章定义（第23-24行）
│   ├── DOI 徽章定义（第26-27行）
│   └── Benchmark 徽章定义（第29-30行）
└── 版本常量定义（第32-42行）
    ├── PythonMinVersion（第32行）
    ├── NumPyMinVersion（第33行）
    ├── SciPyMinVersion（第34行）
    ├── JoblibMinVersion（第35行）
    ├── ThreadpoolctlMinVersion（第36行）
    ├── MatplotlibMinVersion（第37行）
    ├── Scikit-ImageMinVersion（第38行）
    ├── PandasMinVersion（第39行）
    ├── SeabornMinVersion（第40行）
    ├── PytestMinVersion（第41行）
    └── PlotlyMinVersion（第42行）

90.4 README 元信息：体育场的"入场大屏"

在正式进入技术细节之前，我们先把目光投向 README.md 文件最顶部的"徽章墙"——这些徽章是 scikit-learn 项目健康度的实时仪表盘。接下来的解析将严格限定在 README.md 第1-50行范围内，即徽章墙与紧随其后的版本常量定义区。作为整座"性能奥林匹克体育场"的正门入口，README 充当了项目的 __main__——它是开发者与用户接触 scikit-learn 基准测试世界的第一站，所有后续的安装、运行、贡献行为都从此处发散。

源码路径：README.md - 项目元数据（第1-50行）

.. -*- mode: rst -*-                                # [第1行] Emacs/Vim 文件局部变量，声明文档采用 reStructuredText 模式
                                                     # [第2行] 空行
|Azure| |Codecov| |CircleCI| |Nightly wheels| |Ruff| |PythonVersion| |PyPI| |DOI| |Benchmark|   # [第3行] 一次性引用 9 个待定义的徽章变量
                                                     # [第4行] 空行
.. |Azure| image:: https://dev.azure.com/scikit-learn/scikit-learn/_apis/build/status/scikit-learn.scikit-learn?branchName=main   # [第5行] Azure CI 构建状态徽章
   :target: https://dev.azure.com/scikit-learn/scikit-learn/_build/latest?definitionId=1&branchName=main                       # [第6行] 点击跳转至 Azure 构建详情
                                                     # [第7行] 空行
.. |CircleCI| image:: https://circleci.com/gh/scikit-learn/scikit-learn/tree/main.svg?style=shield   # [第8行] CircleCI 备用构建状态徽章
   :target: https://circleci.com/gh/scikit-learn/scikit-learn                                         # [第9行] 点击跳转至 CircleCI 项目页
                                                     # [10行] 空行
.. |Codecov| image:: https://codecov.io/gh/scikit-learn/scikit-learn/branch/main/graph/badge.svg?token=Pk8G9gg3y9   # [11行] Codecov 测试覆盖率徽章
   :target: https://codecov.io/gh/scikit-learn/scikit-learn                                             # [12行] 点击跳转至 Codecov 项目页
                                                     # [13行] 空行
.. |Nightly wheels| image:: https://github.com/scikit-learn/scikit-learn/actions/workflows/wheels.yml/badge.svg?event=schedule   # [14行] 每晚构建 wheel 包状态徽章
   :target: https://github.com/scikit-learn/scikit-learn/actions?query=workflow%3A%22Wheel+builder%22+event%3Aschedule        # [15行] 点击跳转至 GitHub Actions 工作流页
                                                     # [16行] 空行
.. |Ruff| image:: https://img.shields.io/badge/code%20style-ruff-000000.svg   # [17行] Ruff 代码风格徽章
   :target: https://github.com/astral-sh/ruff                                # [18行] 点击跳转至 Ruff 项目主页
                                                     # [19行] 空行
.. |PythonVersion| image:: https://img.shields.io/pypi/pyversions/scikit-learn.svg   # [20行] Python 版本兼容范围徽章
   :target: https://pypi.org/project/scikit-learn/                              # [21行] 点击跳转至 PyPI 项目页
                                                     # [22行] 空行
.. |PyPI| image:: https://img.shields.io/pypi/v/scikit-learn   # [23行] PyPI 最新发布版本徽章
   :target: https://pypi.org/project/scikit-learn              # [24行] 点击跳转至 PyPI 项目页
                                                     # [25行] 空行
.. |DOI| image:: https://zenodo.org/badge/21369/scikit-learn/scikit-learn.svg   # [26行] Zenodo DOI 学术引用徽章
   :target: https://zenodo.org/badge/latestdoi/21369/scikit-learn/scikit-learn  # [27行] 点击跳转至 Zenodo 引用页
                                                     # [28行] 空行
.. |Benchmark| image:: https://img.shields.io/badge/Benchmarked%20by-asv-blue   # [29行] asv 基准测试入口徽章
   :target: https://scikit-learn.org/scikit-learn-benchmarks                     # [30行] 点击跳转至 scikit-learn-benchmarks 网站

这一段 RST 源码在第1行通过 Emacs/Vim 文件局部变量声明文档采用 reStructuredText 模式，使编辑器能够启用对应的语法高亮与缩进规则；第3行则以一行简洁的语法同时引用了 9 个后续会定义的 RST 替换变量名，这种"声明一次、引用一次"的写法是 RST 文档系统管理可复用元素的惯用手法。第5-6行开始进入真正的徽章定义区：每一对连续的行（第5-6、第8-9、第11-12……第29-30行）都构成一个标准的 RST 替换变量声明，其中 .. |Name| image:: 指令负责告诉渲染引擎"当遇到 |Name| 时，插入这张图片"，而紧跟其后的 :target: 缩进指令则指定了用户点击图片后浏览器应跳转到的目标 URL，二者搭配即可在文档中渲染出一枚带超链接的徽章。从语义内容上，徽章涵盖了 CI 构建状态（Azure、CircleCI）、测试覆盖率（Codecov）、每晚构建产物（Nightly wheels）、代码风格工具（Ruff）、Python 版本兼容范围（PythonVersion）、PyPI 发布版本（PyPI）、学术引用 DOI（DOI）以及基准测试入口（Benchmark）共 9 个维度，构成了一张浓缩的"项目健康仪表盘"。

下表汇总了 9 枚徽章的具体含义，便于读者对照查看。

| 徽章名称 | 含义 | 跳转目标 |

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

| Azure | 主分支 Azure Pipelines 构建状态 | Azure DevOps 构建详情页 |

| CircleCI | CircleCI 备用构建状态 | CircleCI 项目页 |

| Codecov | 代码覆盖率统计 | Codecov 项目页 |

| Nightly wheels | 每晚自动构建的 wheel 包状态 | GitHub Actions 工作流页 |

| Ruff | 代码风格检查工具标识 | Ruff 项目主页 |

| PythonVersion | 支持的 Python 版本范围 | PyPI 项目页 |

| PyPI | 最新发布版本号 | PyPI 项目页 |

| DOI | Zenodo 数字对象标识符 | Zenodo 引用页 |

| Benchmark | asv 性能基准测试入口 | scikit-learn-benchmarks 网站 |