Sklearn-源码解析-书-v1-0-三十五-
Sklearn 源码解析（书）v1.0（三十五）
下图展示了分解基准类在 ASV 框架下的完整执行流程，从参数网格展开到计时与评分的每一步都清晰可见：
80.6 流形学习基准 —— 高维嵌入的"计算马拉松"
流形学习是非线性降维的代表，t-SNE 是其中最著名的算法。本节我们看 TSNEBenchmark 如何通过 timeout 护栏与 method 参数网格，量化 exact 与 barnes_hut 两种 t-SNE 实现的性能差异。
80.6.1 TSNEBenchmark：t-SNE 算法的性能极限挑战
t-SNE 的核心挑战是其 O(n²) 的计算复杂度，这使得它在中等规模数据上就需要数分钟。我们看 TSNEBenchmark 如何通过 timeout 与数据规模控制应对这一挑战：
源码路径：asv_benchmarks/benchmarks/manifold.py - TSNEBenchmark
class TSNEBenchmark(Estimator, Benchmark):
"""
Benchmarks for t-SNE.
"""
# 单一参数：近似方法
param_names = ["method"]
params = (["exact", "barnes_hut"],)
# 关键：3 分钟超时护栏
timeout = 180000
def setup_cache(self):
super().setup_cache()
# 根据方法动态调整样本数
def make_data(self, params):
(method,) = params
# exact 方法只能处理小数据
n_samples = 500 if method == "exact" else None
return _digits_dataset(n_samples=n_samples)
# 实例化 TSNE
def make_estimator(self, params):
(method,) = params
estimator = TSNE(
random_state=0,
method=method, # 'exact' 或 'barnes_hut'
)
return estimator
# 通过访问内部属性 kl_divergence_ 获取 KL 散度
def make_scorers(self):
self.train_scorer = lambda _, __: self.estimator.kl_divergence_
self.test_scorer = lambda _, __: self.estimator.kl_divergence_
这段代码定义/展示了 TSNEBenchmark 的三项核心设计。第一，timeout=180000 设置 3 分钟超时护栏，防止 exact 方法在大数据上无法完成测试。第二，make_data 中根据 method 动态调整 n_samples：exact 方法只能处理 500 样本（O(n²) 限制），barnes_hut 则可以使用完整 digits 数据集（1797 样本）。第三，make_scorers 通过访问 self.estimator.kl_divergence_ 这个 scikit-learn 内部属性，把 t-SNE 的最终优化目标（KL 散度）作为质量评估指标——这是 ASV 框架通过访问 sklearn 内部状态实现自定义评分的一个典型示例。
80.6.2 流形学习基准的特殊性
80.7 设计中的取舍
为什么 KMeansBenchmark 同时测试 init 和 algorithm 两个维度？这是因为 init 参数考察的是初始化策略对聚类质量的影响——k-means++ 通过智能选择初始质心让收敛更快、质量更稳，而 random 则是基准线。algorithm 参数考察的是底层求解器的实现差异——lloyd 是经典迭代，elkan 利用三角不等式避免不必要的距离计算。在高维或大规模数据上，elkan 通常能显著加速；但在稀疏数据上，由于三角不等式的剪枝效果有限，lloyd 反而可能更快。同时测试两个维度，才能完整勾勒出 KMeans 的性能画像。
为什么 TSNEBenchmark 设置 timeout = 180000 而其他基准如 CrossValidationBenchmark 仅 20000？这是因为 timeout 反映的是算法的预期耗时量级。CrossValidationBenchmark 用的是 5 万样本的合成数据，16 折交叉验证在现代机器上大约几十秒即可完成，因此 20 秒的护栏足够触发"超时报警"。但 TSNEBenchmark 即使是 barnes_hut 方法在 1797 样本上也需要数十秒，exact 方法在 500 样本上同样需要数十秒；3 分钟的护栏给 t-SNE 留出充分空间，同时仍能在意外死循环时及时终止测试。timeout 不是越短越好，而是要在"防止无限等待"与"允许合理耗时"之间取得平衡。
为什么 DictionaryLearningBenchmark 的 fit_algorithm 不包含 'lgd'？原大纲中曾提到 'lgd'（Lagrangian Gradient Descent），但实际代码中 DictionaryLearning 仅支持 'lars' 与 'cd' 两种算法。这是因为 lgd 已被弃用并从 scikit-learn 中移除，参数网格必须与实际可用的 API 保持一致。这也体现了基准代码的演进特性——它需要随着 scikit-learn 主代码库的算法更新而同步调整。
为什么 PCABenchmark 固定 n_components=32 而不是更高维度？这是因为 MNIST 数据集本身具有显著的低维流形结构（手写数字主要分布在低维子空间中），32 维足以捕获绝大部分方差信息；同时，固定 n_components 能避免降维维度本身成为性能差异的干扰变量——三种 SVD 求解器（full/arpack/randomized）的核心差异在于"求解策略"而非"目标维度"，固定维度让比较聚焦于算法本身的实现效率。此外，32 这个数值也确保了 arpack 与 randomized 求解器有足够的截断迭代空间（n_components 太小会导致 arpack 的截断迭代失去意义）。
80.8 动手练习
-
聚类基准：KMeans 与 MiniBatchKMeans 的算法与参数影响
-
阅读
asv_benchmarks/benchmarks/cluster.py中 KMeansBenchmark 和 MiniBatchKMeansBenchmark。 -
回答问题：
-
KMeansBenchmark 中
init参数（'k-means++' vs 'random'）对聚类收敛速度和惯性的影响是什么？ -
为什么 KMeansBenchmark 使用
lloyd和elkan两种算法进行对比？elkan 算法在什么情况下更优？ -
MiniBatchKMeansBenchmark 引入了
batch_size参数，请解释其对训练速度和聚类质量的权衡。代码中默认 batch_size 是多少？ -
两个基准都使用了
train_scorer/test_scorer跟踪inertia，这说明了无监督基准在评估聚类质量时的共同策略是什么？
-
-
-
维度降低基准：PCA 与字典学习的求解器与算法对比
-
阅读
asv_benchmarks/benchmarks/decomposition.py中 PCABenchmark、DictionaryLearningBenchmark 和 MiniBatchDictionaryLearningBenchmark。 -
回答问题：
-
PCABenchmark 中
svd_solver参数（full, arpack, randomized）在 MNIST 数据上的性能差异主要体现在哪些方面？ -
DictionaryLearningBenchmark 与 MiniBatchDictionaryLearningBenchmark 在数据规模和训练策略上有什么关键区别？
-
为什么 DictionaryLearningBenchmark 的
fit_algorithm参数包括 'lars', 'cd'？这两种算法在稀疏编码中的适用场景是什么？ -
所有分解基准类均使用
make_pca_scorers或make_dict_learning_scorers生成评分器，这说明了它们在评估降维质量时的共同目标是什么？
-
-
-
流形学习与监督基准对比：t-SNE 的特殊性与线性模型的共性
-
阅读
asv_benchmarks/benchmarks/manifold.py中 TSNEBenchmark 和asv_benchmarks/benchmarks/linear_model.py中 LogisticRegressionBenchmark。 -
回答问题：
-
TSNEBenchmark 为什么设置了
timeout = 180000？这与其他基准（如 CrossValidationBenchmark 的 timeout=20000）有什么不同？ -
TSNEBenchmark 只测试
method参数（exact vs barnes_hut），为什么不包含n_jobs或数据规模作为可调参数？ -
在
make_scorers方法中，如何访问 TSNE 模型的 KL 散度？这说明了 ASV 框架在评估非标准指标时的扩展方式是什么？ -
对比 TSNEBenchmark（继承 Benchmark 和 Estimator）与 LogisticRegressionBenchmark（继承 Predictor, Estimator, Benchmark），二者在执行流程上的根本区别是什么？
-
-
80.9 本章小结
这一章中我们系统学习了聚类、维度降低与流形学习三大无监督学习基准的测试设计原则与具体实现。首先，我们了解了 KMeansBenchmark 与 MiniBatchKMeansBenchmark 如何通过 representation、algorithm、init、batch_size 等参数网格考察聚类算法的速度-质量权衡；其次，我们剖析了 PCABenchmark 与 DictionaryLearningBenchmark 如何在 MNIST、Olivetti 等数据集上测试不同 SVD 求解器与字典学习算法的训练效率；接着，我们深入了 TSNEBenchmark 如何通过 timeout 护栏与数据规模动态调整应对 t-SNE 的 O(n²) 复杂度挑战；最后，我们梳理了 ASV 框架下 make_data、make_estimator、make_scorers、setup_cache 等钩子的标准模式及其在无监督基准中的差异化应用。
下表总结了本章涉及的核心基准类及其关键特性：
| 概念 | 解释 |
|------|------|
| KMeansBenchmark | 对比 representation（dense/sparse）、algorithm（lloyd/elkan）和 init（random/k-means++）三维参数，固定 n_clusters=20，n_init=1，max_iter 动态调整（sparse:30, dense:100），使用 neg_mean_inertia 评估聚类质量 |
| MiniBatchKMeansBenchmark | 测试 batch_size 对 MiniBatchKMeans 性能的影响，参数网格为 representation × init，batch_size=1000，禁用早停与标签计算以保证计时公平 |
| PCABenchmark | 测试 svd_solver（full/arpack/randomized）在 MNIST 数据上的 PCA 性能，固定 n_components=32，调用 make_pca_scorers 评估重构误差 |
| DictionaryLearningBenchmark | 测试 fit_algorithm（lars/cd）和 n_jobs 对字典学习训练效率的影响，固定字典大小=15，使用 Olivetti 人脸数据集 |
| MiniBatchDictionaryLearningBenchmark | 测试 fit_algorithm（lars/cd）和 n_jobs 对 MiniBatch 字典学习的影响，增加 batch_size=3 参数，其他设置同 DictionaryLearningBenchmark |
| TSNEBenchmark | 对比 exact 和 barnes_hut 方法在 digits 数据上的 t-SNE 性能，exact 用 500 样本，barnes_hut 用完整 1797 样本，timeout=180000 防超时，通过 kl_divergence_ 评估 KL 散度 |
第 81 章 —— 独立基准脚本：梯度提升专题 —— 挑战"集成学习天花板"
81.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解梯度提升决策树（GBDT）基准测试的核心流程与设计模式
-
掌握 HistGradientBoosting 与 LightGBM/XGBoost/CatBoost 的参数等价映射机制
-
了解多线程环境下 OpenMP 线程控制与并行扩展性评测方法
-
理解真实异构数据集（Adult）与合成纯分类特征数据集的基准测试差异
-
掌握大规模数据集（HIGGS）下的缓存、子采样与特征交互约束配置技巧
-
能编写自定义基准脚本对比多个梯度提升库的训练/推理性能
81.2 生活类比
想象这是一场"梯度提升库奥林匹克大赛"。scikit-learn 的 HistGradientBoosting、LightGBM、XGBoost、CatBoost 四位选手各怀绝技，要在速度与精度两个维度上一决高下，而我们要为他们设计公平且富有挑战性的比赛。具体而言，bench_hist_gradient_boosting.py 充当标准赛道的角色，用合成数据和统一规则考察各库在标准分类/回归任务上的基础速度与精度，如同百米跑道上比试爆发力；bench_hist_gradient_boosting_threading.py 则是接力赛/多核挑战赛，通过调节线程数（即"接力棒人数"），考察哪位选手能更高效地利用多核 CPU，评估并行扩展能力；bench_hist_gradient_boosting_higgsboson.py 堪称马拉松/重量级比赛，2.6GB 的 HIGGS 真实大数据登场，考察选手的内存管理、IO 效率与大规模并行能力，谁能在这场耐力赛中保持稳定，谁就是真正的高手；bench_hist_gradient_boosting_adult.py 是复杂地形越野赛，使用真实异构表格数据（数值+类别混合），考察选手对原生分类特征的处理能力与预处理集成度；最后，bench_hist_gradient_boosting_categorical_only.py 是极限障碍赛，纯离散序数特征登场，考察选手在无连续数值、全部分箱场景下的鲁棒性。而贯穿所有比赛的"通用翻译官"get_equivalent_estimator，将 sklearn 的通用参数"实时同传"成 LightGBM/XGBoost/CatBoost 的专属方言，确保比赛规则一致、参数等价，让比赛结果真正反映各库的性能差异。
81.3 源码地图
benchmarks/bench_hist_gradient_boosting.py
├── 核心流程
│ ├── one_run() # 单次基准运行：训练、评分、记录耗时
│ ├── get_estimator_and_data() # 生成合成数据与获取对应 Estimator 类
│ └── __main__ # 参数解析、循环调用 one_run、绘图展示
├── 模型转换
│ └── get_equivalent_estimator() (调用) # 来自 sklearn.ensemble._hist_gradient_boosting.utils，跨库参数映射
benchmarks/bench_hist_gradient_boosting_higgsboson.py
├── 核心流程
│ ├── load_data() # 下载、缓存、解析 HIGGS 数据集 (joblib.Memory)
│ ├── fit() # 统一训练计时封装
│ ├── predict() # 统一预测计时与指标计算 (ROC AUC, Accuracy)
│ └── __main__ # 参数配置、交互约束设置、多库对比流程
├── 模型转换
│ └── get_equivalent_estimator() (调用) # sklearn 模型参数转 LightGBM/XGBoost/CatBoost
benchmarks/bench_hist_gradient_boosting_threading.py
├── 核心流程
│ ├── one_run() # 线程级基准：threadpool_limits 控制 OpenMP 线程数
│ ├── get_estimator_and_data() # 生成数据与基础 sklearn Estimator
│ └── __main__ # 线程列表生成、循环测试、可视化扩展性曲线
├── 模型转换
│ └── get_equivalent_estimator() (调用) # 并设置各库线程参数 (num_threads, nthread, thread_count)
benchmarks/bench_hist_gradient_boosting_adult.py
├── 核心流程
│ ├── fit() # 通用训练封装
│ ├── predict() # 通用预测封装
│ └── __main__ # fetch_openml、ColumnTransformer+OrdinalEncoder、原始生分类特征支持对比
├── 模型转换
│ └── get_equivalent_estimator() (调用) # 设置 LightGBM max_cat_to_onehot=1 禁用 One-Hot
benchmarks/bench_hist_gradient_boosting_categorical_only.py
├── 核心流程
│ ├── fit() # 通用训练封装
│ ├── predict() # 简化预测：仅计时，不计算指标
│ └── __main__ # make_classification + KBinsDiscretizer 生成纯分类数据、全特征 categorical=True
├── 模型转换
│ └── get_equivalent_estimator() (调用) # LightGBM categorical_feature=list(range(n_features))
81.4 梯度提升库大对决 —— HistGB 与 LightGBM/XGBoost/CatBoost 的"擂台赛"
bench_hist_gradient_boosting.py 是整个梯度提升基准测试的"标准赛道"，它以合成数据为考题，统一调用 one_run 函数对四个主流梯度提升库进行训练/推理耗时与精度的全面对比。
one_run(n_samples) 是基准测试的原子操作单元，在固定数据规模和线程配置下，完成"训练→计时→评分→计时"的标准化流程。脚本通过命令行参数控制是否同时对比 LightGBM、XGBoost、CatBoost，最多支持四个选手同场竞技。
源码路径：benchmarks/bench_hist_gradient_boosting.py - one_run()（74-145行）
def one_run(n_samples):
X_train = X_train_[:n_samples] # ① 截取指定规模的训练子集
X_test = X_test_[:n_samples] # ② 截取对应测试子集
y_train = y_train_[:n_samples] # ③ 同步训练标签
y_test = y_test_[:n_samples] # ④ 同步测试标签
if sample_weight is not None: # ⑤ 若启用样本权重
sample_weight_train = sample_weight_train_[:n_samples] # ⑥ 同步样本权重
else:
sample_weight_train = None
assert X_train.shape[0] == n_samples # ⑦ 断言：保证切片大小一致
assert X_test.shape[0] == n_samples
print("Data size: %d samples train, %d samples test." % (n_samples, n_samples))
print("Fitting a sklearn model...")
tic = time() # ⑧ 启动训练计时器
est = Estimator( # ⑨ 实例化 HistGB 估计器
learning_rate=lr,
max_iter=n_trees,
max_bins=max_bins,
max_leaf_nodes=n_leaf_nodes,
early_stopping=False, # ⑩ 关闭早停，保证固定轮数
random_state=0, # ⑪ 固定随机种子确保可复现
verbose=0,
)
loss = args.loss
if args.problem == "classification":
if loss == "default":
loss = "log_loss" # ⑫ 分类默认对数损失
else:
if loss == "default":
loss = "squared_error" # ⑬ 回归默认平方损失
est.set_params(loss=loss) # ⑭ 设置损失函数
est.fit(X_train, y_train, sample_weight=sample_weight_train) # ⑮ 执行训练
sklearn_fit_duration = time() - tic # ⑯ 记录训练耗时
tic = time()
sklearn_score = est.score(X_test, y_test) # ⑰ 计算测试集评分
sklearn_score_duration = time() - tic # ⑱ 记录评分耗时
print("score: {:.4f}".format(sklearn_score))
print("fit duration: {:.3f}s,".format(sklearn_fit_duration))
print("score duration: {:.3f}s,".format(sklearn_score_duration))
lightgbm_score = None # ⑲ 初始化 LightGBM 指标占位
lightgbm_fit_duration = None
lightgbm_score_duration = None
if args.lightgbm: # ⑳ 可选：对比 LightGBM
print("Fitting a LightGBM model...")
lightgbm_est = get_equivalent_estimator( # ㉑ 翻译为等价的 LightGBM 模型
est, lib="lightgbm", n_classes=args.n_classes
)
tic = time()
lightgbm_est.fit(X_train, y_train, sample_weight=sample_weight_train)
lightgbm_fit_duration = time() - tic
tic = time()
lightgbm_score = lightgbm_est.score(X_test, y_test)
lightgbm_score_duration = time() - tic
print("score: {:.4f}".format(lightgbm_score))
print("fit duration: {:.3f}s,".format(lightgbm_fit_duration))
print("score duration: {:.3f}s,".format(lightgbm_score_duration))
xgb_score = None # ㉒ 初始化 XGBoost 指标占位
xgb_fit_duration = None
xgb_score_duration = None
if args.xgboost: # ㉓ 可选：对比 XGBoost
print("Fitting an XGBoost model...")
xgb_est = get_equivalent_estimator(est, lib="xgboost", n_classes=args.n_classes)
tic = time()
xgb_est.fit(X_train, y_train, sample_weight=sample_weight_train)
xgb_fit_duration = time() - tic
tic = time()
xgb_score = xgb_est.score(X_test, y_test)
xgb_score_duration = time() - tic
print("score: {:.4f}".format(xgb_score))
print("fit duration: {:.3f}s,".format(xgb_fit_duration))
print("score duration: {:.3f}s,".format(xgb_score_duration))
cat_score = None # ㉔ 初始化 CatBoost 指标占位
cat_fit_duration = None
cat_score_duration = None
if args.catboost: # ㉕ 可选：对比 CatBoost
print("Fitting a CatBoost model...")
cat_est = get_equivalent_estimator(
est, lib="catboost", n_classes=args.n_classes
)
tic = time()
cat_est.fit(X_train, y_train, sample_weight=sample_weight_train)
cat_fit_duration = time() - tic
tic = time()
cat_score = cat_est.score(X_test, y_test)
cat_score_duration = time() - tic
print("score: {:.4f}".format(cat_score))
print("fit duration: {:.3f}s,".format(cat_fit_duration))
print("score duration: {:.3f}s,".format(cat_score_duration))
return ( # ㉖ 返回 12 项指标供主循环汇总
sklearn_score,
sklearn_fit_duration,
sklearn_score_duration,
lightgbm_score,
lightgbm_fit_duration,
lightgbm_score_duration,
xgb_score,
xgb_fit_duration,
xgb_score_duration,
cat_score,
cat_fit_duration,
cat_score_duration,
)
这段代码定义了 one_run 函数，在给定数据规模下完成 sklearn、LightGBM、XGBoost、CatBoost 四个选手的完整对比。它通过 get_equivalent_estimator 实现跨库参数映射，让对比结果真正反映各库性能差异。one_run 每次返回 12 个关键指标（每个库 3 个：评分、训练耗时、推理耗时），主循环汇总这些指标后调用 matplotlib 绘制三联图：评分曲线、训练耗时曲线、推理耗时曲线，横轴为样本数（对数刻度）。注意每个未启用的库（如 args.lightgbm=False）对应指标均初始化为 None，避免在主循环汇总时因缺失值报错——这是 Python 中常见的"占位+可选赋值"模式。
下图展示了 one_run 内部"数据切片→sklearn→LightGBM→XGBoost→CatBoost"四阶段串行对比的控制流：
注意四个库的训练与评分阶段在 one_run 内是严格串行执行的（同一函数体、按代码顺序运行），并非并行。这种设计保证了计时的纯净性——每次只有单一库在占用 CPU，但代价是无法捕获并发场景下的资源竞争。脚本最终通过 for n_samples in n_samples_list 外层循环扫描 1000、10000、100000、500000、1000000、5000000、10000000 这七个数据规模点（实际取不超过 args.n_samples_max 的部分），记录各库在每个规模下的耗时与精度，最终用 plt.subplots(3, sharex=True) 绘制三联对比图。
get_estimator_and_data() 负责"装备准备"，根据任务类型生成数据并返回对应的 sklearn 估计器类。
源码路径：benchmarks/bench_hist_gradient_boosting.py - get_estimator_and_data()（54-72行）
def get_estimator_and_data():
if args.problem == "classification": # ① 分类任务
X, y = make_classification( # ② 生成合成分类数据
args.n_samples_max * 2, # ③ 生成两倍样本，后续切分
n_features=args.n_features,
n_classes=args.n_classes,
n_clusters_per_class=1, # ④ 每类一个高斯团簇
n_informative=args.n_classes, # ⑤ 信息特征数等于类别数
random_state=0,
)
return X, y, HistGradientBoostingClassifier # ⑥ 返回分类器类
elif args.problem == "regression": # ⑦ 回归任务
X, y = make_regression(
args.n_samples_max * 2, n_features=args.n_features, random_state=0
)
return X, y, HistGradientBoostingRegressor # ⑧ 返回回归器类
这段代码定义了数据准备函数，根据问题类型选择合成数据生成器，返回 (X, y, EstimatorClass) 三元组。n_clusters_per_class=1 与 n_informative=args.n_classes 是关键参数：前者保证数据可分，后者让信息特征数与类别数匹配，确保每个类别都有足够的判别信号，模型既能学到模式又不至于过于简单。
get_estimator_and_data 的"分类 vs 回归"分支结构如下图所示：
主流程 __main__ 部分（147-230行）串联起整个基准测试：参数解析 → 数据生成 → 缺失值注入 → 样本切分 → 循环调用 one_run → 结果可视化。其核心逻辑是"用 for n_samples in n_samples_list 循环扫描 1000、10000、100000、500000、1000000、5000000、10000000 这七个数据规模点，记录各库在每个规模下的耗时与精度，最终用 plt.subplots(3, sharex=True) 绘制三联对比图。
fit() 与 predict() 函数在 bench_hist_gradient_boosting_higgsboson.py 中以通用计时封装的形式出现：
源码路径：benchmarks/bench_hist_gradient_boosting_higgsboson.py - fit()（63-70行）
def fit(est, data_train, target_train, libname):
print(f"Fitting a {libname} model...") # ① 输出库标识
tic = time() # ② 启动计时
est.fit(data_train, target_train) # ③ 执行训练
toc = time() # ④ 结束计时
print(f"fitted in {toc - tic:.3f}s") # ⑤ 输出耗时
源码路径：benchmarks/bench_hist_gradient_boosting_higgsboson.py - predict()（72-81行）
def predict(est, data_test, target_test):
if args.no_predict: # ① 可选：跳过预测阶段
return
tic = time() # ② 启动预测计时
predicted_test = est.predict(data_test) # ③ 硬预测（类别）
predicted_proba_test = est.predict_proba(data_test) # ④ 软预测（概率）
toc = time()
roc_auc = roc_auc_score(target_test, predicted_proba_test[:, 1]) # ⑤ 计算 ROC AUC
acc = accuracy_score(target_test, predicted_test) # ⑥ 计算准确率
print(f"predicted in {toc - tic:.3f}s, ROC AUC: {roc_auc:.4f}, ACC: {acc:.4f}")
这段代码定义了两个计时辅助函数，将"计时"和"指标计算"封装成独立单元，让主流程的对比代码更清晰。fit 只做训练计时，predict 同时完成预测计时和 ROC/AUC/ACC 指标计算，遵循"单一职责"原则。fit 函数的设计极其精简——只接收估计器、训练数据与库名三个核心参数，避免了对每个库的不同 fit 签名做特殊处理；predict 则在计算两个指标时分别使用 predict（硬预测，输出类别标签）和 predict_proba（软预测，输出概率分布）的第二列作为 ROC AUC 的输入——这是二分类问题的标准做法。
fit 与 predict 在 HIGGS 脚本中的协作流程如下图所示：
81.5 模型转换机制 —— get_equivalent_estimator 的"跨库翻译官"
get_equivalent_estimator 是连接 sklearn 与各梯度提升库的"翻译官"，它位于 sklearn.ensemble._hist_gradient_boosting.utils 模块中，负责将 HistGradientBoosting 的核心参数映射为目标库（LightGBM/XGBoost/CatBoost）的等价参数。在所有五个基准脚本中，这一函数都是核心调用点。
以下是各脚本中调用 get_equivalent_estimator 的位置与上下文：
源码路径：benchmarks/bench_hist_gradient_boosting.py - get_equivalent_estimator()（96, 110, 124, 138行）
# 第 81 章 —— 在 one_run() 函数中，针对不同库的三次调用
lightgbm_est = get_equivalent_estimator(est, lib="lightgbm", n_classes=args.n_classes) # ① 翻译为 LightGBM
xgb_est = get_equivalent_estimator(est, lib="xgboost", n_classes=args.n_classes) # ② 翻译为 XGBoost
cat_est = get_equivalent_estimator(est, lib="catboost", n_classes=args.n_classes) # ③ 翻译为 CatBoost
源码路径：benchmarks/bench_hist_gradient_boosting_higgsboson.py - get_equivalent_estimator()（95, 103, 111, 119行）
# 第 81 章 —— 在 __main__ 中，针对不同库的三次调用
est = get_equivalent_estimator(est, lib="lightgbm", n_classes=n_classes) # ① 翻译为 LightGBM
est = get_equivalent_estimator(est, lib="xgboost", n_classes=n_classes) # ② 翻译为 XGBoost
est = get_equivalent_estimator(est, lib="catboost", n_classes=n_classes) # ③ 翻译为 CatBoost
源码路径：benchmarks/bench_hist_gradient_boosting_threading.py - get_equivalent_estimator()（128, 152, 167, 182行）
# 第 81 章 —— 在 one_run() 中，调用后需手动设置各库线程参数
lightgbm_est = get_equivalent_estimator(est, lib="lightgbm", n_classes=args.n_classes)
lightgbm_est.set_params(num_threads=n_threads) # ⑦ 设置 LightGBM 线程数
xgb_est = get_equivalent_estimator(est, lib="xgboost", n_classes=args.n_classes)
xgb_est.set_params(nthread=n_threads) # ⑧ 设置 XGBoost 线程数
cat_est = get_equivalent_estimator(est, lib="catboost", n_classes=args.n_classes)
cat_est.set_params(thread_count=n_threads) # ⑨ 设置 CatBoost 线程数
源码路径：benchmarks/bench_hist_gradient_boosting_adult.py - get_equivalent_estimator()（89行）
# 第 81 章 —— 在 __main__ 中，针对 LightGBM 的调用
est = get_equivalent_estimator(est, lib="lightgbm", n_classes=n_classes)
est.set_params(max_cat_to_onehot=1) # dont use OHE # ⑦ 禁用 One-Hot 编码
源码路径：benchmarks/bench_hist_gradient_boosting_categorical_only.py - get_equivalent_estimator()（87行）
# 第 81 章 —— 在 __main__ 中，针对 LightGBM 的调用
est = get_equivalent_estimator(est, lib="lightgbm", n_classes=2)
est.set_params(max_cat_to_onehot=1) # dont use OHE # ⑧ 禁用 One-Hot 编码
get_equivalent_estimator 接收三个关键参数：est（sklearn HistGB 估计器实例）、lib（目标库名）、n_classes（类别数）。它读取 est 的参数（如 learning_rate、max_iter、max_leaf_nodes、max_bins、loss 等），然后返回一个已实例化的目标库估计器对象，参数值与原 sklearn 估计器保持一致。这种"翻译官"模式确保了对比的公平性——所有库都在等价配置下竞技，性能差异真正反映算法实现而非参数设置。
get_equivalent_estimator 在整个基准测试体系中的位置与调用层次如下图所示：
81.6 线程扩展性基准 —— 多核时代的"并行加速赛"
bench_hist_gradient_boosting_threading.py 是多核时代的"并行加速赛"脚本。它使用 threadpool_limits 强制控制 OpenMP 线程数，遍历 1、2、4、...、max_threads 等线程配置，对比四个梯度提升库在 fit 与 score 阶段的并行扩展效率。
get_estimator_and_data() 是数据与基础估计器准备函数，与标准赛脚本保持一致但参数有调整：
源码路径：benchmarks/bench_hist_gradient_boosting_threading.py - get_estimator_and_data()（72-90行）
def get_estimator_and_data():
if args.problem == "classification": # ① 分类任务
X, y = make_classification(
args.n_samples * 2, # ② 注意：使用 args.n_samples 而非 n_samples_max
n_features=args.n_features,
n_classes=args.n_classes,
n_clusters_per_class=1,
n_informative=args.n_features // 2, # ③ 半数特征为信息特征（区别于标准赛）
random_state=0,
)
return X, y, HistGradientBoostingClassifier
elif args.problem == "regression": # ④ 回归任务
X, y = make_regression(
args.n_samples_max * 2, n_features=args.n_features, random_state=0
)
return X, y, HistGradientBoostingRegressor
注意第 ③ 行 n_informative=args.n_features // 2 与标准赛不同：线程测试脚本采用"半数信息特征"配置，确保合成数据有足够难度（更高的特征交互复杂度），让并行收益更明显，避免在简单任务上线程数翻倍没有可观测的加速比。
get_estimator_and_data 在线程测试脚本中的分支流程如下图所示：
one_run(n_threads, n_samples) 是线程级基准的核心执行单元：
源码路径：benchmarks/bench_hist_gradient_boosting_threading.py - one_run()（111-193行）
def one_run(n_threads, n_samples):
X_train = X_train_[:n_samples] # ① 数据切片
X_test = X_test_[:n_samples]
y_train = y_train_[:n_samples]
y_test = y_test_[:n_samples]
if sample_weight is not None:
sample_weight_train = sample_weight_train_[:n_samples]
else:
sample_weight_train = None
assert X_train.shape[0] == n_samples
assert X_test.shape[0] == n_samples
print("Fitting a sklearn model...")
tic = time()
est = sklearn.base.clone(sklearn_est) # ② 克隆基础估计器（隔离状态）
with threadpool_limits(n_threads, user_api="openmp"): # ③ 关键：在 OpenMP 级别限制线程数
est.fit(X_train, y_train, sample_weight=sample_weight_train) # ④ sklearn 训练
sklearn_fit_duration = time() - tic
tic = time()
sklearn_score = est.score(X_test, y_test) # ⑤ sklearn 评分
sklearn_score_duration = time() - tic
print("score: {:.4f}".format(sklearn_score))
print("fit duration: {:.3f}s,".format(sklearn_fit_duration))
print("score duration: {:.3f}s,".format(sklearn_score_duration))
lightgbm_score = None # ⑥ 初始化 LightGBM 指标占位
lightgbm_fit_duration = None
lightgbm_score_duration = None
if args.lightgbm: # ⑦ 可选：LightGBM 对比
print("Fitting a LightGBM model...")
lightgbm_est = get_equivalent_estimator(est, lib="lightgbm", n_classes=args.n_classes)
lightgbm_est.set_params(num_threads=n_threads) # ⑧ 设置 LightGBM 线程数
tic = time()
lightgbm_est.fit(X_train, y_train, sample_weight=sample_weight_train)
lightgbm_fit_duration = time() - tic
tic = time()
lightgbm_score = lightgbm_est.score(X_test, y_test)
lightgbm_score_duration = time() - tic
print("score: {:.4f}".format(lightgbm_score))
print("fit duration: {:.3f}s,".format(lightgbm_fit_duration))
print("score duration: {:.3f}s,".format(lightgbm_score_duration))
xgb_score = None # ⑨ 初始化 XGBoost 指标占位
xgb_fit_duration = None
xgb_score_duration = None
if args.xgboost: # ⑩ 可选：XGBoost 对比
print("Fitting an XGBoost model...")
xgb_est = get_equivalent_estimator(est, lib="xgboost", n_classes=args.n_classes)
xgb_est.set_params(nthread=n_threads) # ⑪ 设置 XGBoost 线程数
tic = time()
xgb_est.fit(X_train, y_train, sample_weight=sample_weight_train)
xgb_fit_duration = time() - tic
tic = time()
xgb_score = xgb_est.score(X_test, y_test)
xgb_score_duration = time() - tic
print("score: {:.4f}".format(xgb_score))
print("fit duration: {:.3f}s,".format(xgb_fit_duration))
print("score duration: {:.3f}s,".format(xgb_score_duration))
cat_score = None # ⑫ 初始化 CatBoost 指标占位
cat_fit_duration = None
cat_score_duration = None
if args.catboost: # ⑬ 可选：CatBoost 对比
print("Fitting a CatBoost model...")
cat_est = get_equivalent_estimator(est, lib="catboost", n_classes=args.n_classes)
cat_est.set_params(thread_count=n_threads) # ⑭ 设置 CatBoost 线程数
tic = time()
cat_est.fit(X_train, y_train, sample_weight=sample_weight_train)
cat_fit_duration = time() - tic
tic = time()
cat_score = cat_est.score(X_test, y_test)
cat_score_duration = time() - tic
print("score: {:.4f}".format(cat_score))
print("fit duration: {:.3f}s,".format(cat_fit_duration))
print("score duration: {:.3f}s,".format(cat_score_duration))
return ( # ⑮ 返回 12 项指标
sklearn_score,
sklearn_fit_duration,
sklearn_score_duration,
lightgbm_score,
lightgbm_fit_duration,
lightgbm_score_duration,
xgb_score,
xgb_fit_duration,
xgb_score_duration,
cat_score,
cat_fit_duration,
cat_score_duration,
)
这段代码是线程基准的"核心执行单元"。注意第 ③ 行的 threadpool_limits 上下文管理器，它来自 threadpoolctl，在 OpenMP 级别强制限制线程数，确保 sklearn（依赖 OpenMP 并行）在指定线程数下运行。对于其他库，脚本通过 set_params 显式设置各库的线程参数：LightGBM 用 num_threads、XGBoost 用 nthread、CatBoost 用 thread_count。这种"双管齐下"的线程控制策略，避免了环境变量干扰，实现了精准的并行扩展性对比。
多线程基准的整体控制流如下图所示：
源码路径：benchmarks/bench_hist_gradient_boosting_threading.py - __main__（195-280行）
max_threads = os.cpu_count() # ① 获取 CPU 核心数
n_threads_list = [2**i for i in range(8) if (2**i) < max_threads] # ② 生成 1,2,4,8... 线程列表
n_threads_list.append(max_threads) # ③ 添加最大线程数
sklearn_scores = [] # ④ 初始化各库指标结果列表
sklearn_fit_durations = []
sklearn_score_durations = []
lightgbm_scores = []
lightgbm_fit_durations = []
lightgbm_score_durations = []
xgb_scores = []
xgb_fit_durations = []
xgb_score_durations = []
cat_scores = []
cat_fit_durations = []
cat_score_durations = []
for n_threads in n_threads_list: # ⑤ 遍历所有线程配置
print(f"n_threads: {n_threads}")
(
sklearn_score,
sklearn_fit_duration,
sklearn_score_duration,
lightgbm_score,
lightgbm_fit_duration,
lightgbm_score_duration,
xgb_score,
xgb_fit_duration,
xgb_score_duration,
cat_score,
cat_fit_duration,
cat_score_duration,
) = one_run(n_threads, n_samples)
for scores, score in ( # ⑥ 将 one_run 返回值按库分别收集
(sklearn_scores, sklearn_score),
(sklearn_fit_durations, sklearn_fit_duration),
(sklearn_score_durations, sklearn_score_duration),
(lightgbm_scores, lightgbm_score),
(lightgbm_fit_durations, lightgbm_fit_duration),
(lightgbm_score_durations, lightgbm_score_duration),
(xgb_scores, xgb_score),
(xgb_fit_durations, xgb_fit_duration),
(xgb_score_durations, xgb_score_duration),
(cat_scores, cat_score),
(cat_fit_durations, cat_fit_duration),
(cat_score_durations, cat_score_duration),
):
scores.append(score)
if args.plot or args.plot_filename: # ⑦ 可视化扩展曲线
import matplotlib
import matplotlib.pyplot as plt
fig, axs = plt.subplots(2, figsize=(12, 12)) # ⑧ 创建 2 行子图
# 绘制各库的 fit/score 耗时随线程数变化的曲线
线程列表采用 2**i 等比序列，覆盖从单线程到全核心的范围，便于绘制对数刻度的扩展性曲线。可视化部分绘制两行子图：上行是 fit 耗时，下行是 score 耗时，横轴为线程数（对数刻度），多条曲线对应不同库。注意第 ④ 行的 12 个空列表与第 ⑥ 行的元组循环——这是一种简洁的"zip-like"解构模式，将 one_run 返回的 12 元结果按对应库分别追加到 12 个独立列表中，避免手动逐个 append。
81.7 分类特征与 Adult 数据集专项 —— 应对"非数值世界"的挑战
bench_hist_gradient_boosting_adult.py 是真实异构表格数据的"越野赛"，使用 Adult Census Income 数据集（数据 ID 179），对比 sklearn 的 HistGB 与 LightGBM 在处理原生分类特征时的性能差异。
主流程 __main__ 部分整合了 fetch_openml、ColumnTransformer 预处理、OrdinalEncoder 编码与原生分类特征支持：
源码路径：benchmarks/bench_hist_gradient_boosting_adult.py - __main__（37-89行）
data = fetch_openml(data_id=179, as_frame=True) # ① 从 OpenML 拉取 Adult 数据
X, y = data.data, data.target
cat_columns = make_column_selector(dtype_include="category")(X) # ② 识别类别型列
preprocessing = make_column_transformer(
(OrdinalEncoder(), cat_columns), # ③ 用 OrdinalEncoder 编码类别
remainder="passthrough", # ④ 数值列原样透传
verbose_feature_names_out=False, # ⑤ 关闭特征名前缀噪声
)
X = pd.DataFrame( # ⑥ 重新封装为 DataFrame
preprocessing.fit_transform(X),
columns=preprocessing.get_feature_names_out(),
)
n_classes = len(np.unique(y)) # ⑦ 统计类别数
n_features = X.shape[1] # ⑧ 统计特征数
n_categorical_features = len(cat_columns) # ⑨ 类别特征数
n_numerical_features = n_features - n_categorical_features # ⑩ 数值特征数
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
is_categorical = [True] * n_categorical_features + [False] * n_numerical_features # ⑪ 构建布尔掩码
est = HistGradientBoostingClassifier(
loss="log_loss",
learning_rate=lr,
max_iter=n_trees,
max_bins=max_bins,
max_leaf_nodes=n_leaf_nodes,
categorical_features=is_categorical, # ⑫ 告知 sklearn 哪些是分类特征
early_stopping=False,
random_state=0,
verbose=verbose,
)
fit(est, X_train, y_train, "sklearn") # ⑬ 训练 sklearn
predict(est, X_test, y_test) # ⑭ 评估
if args.lightgbm: # ⑮ 可选：LightGBM 对比
est = get_equivalent_estimator(est, lib="lightgbm", n_classes=n_classes)
est.set_params(max_cat_to_onehot=1) # dont use OHE # ⑯ 禁用 One-Hot
categorical_features = [
f_idx for (f_idx, is_cat) in enumerate(is_categorical) if is_cat
]
fit(est, X_train, y_train, "lightgbm", categorical_feature=categorical_features)
predict(est, X_test, y_test)
这段代码展示了真实异构表格数据的基准流程。fetch_openml(data_id=179) 拉取 Adult 数据集（约 48842 条记录，14 维特征，包含 8 个类别特征与 6 个数值特征）。关键技巧是 make_column_transformer + OrdinalEncoder 的组合：只对类别列做序数编码，数值列保持原样透传，避免了 One-Hot 编码带来的维度爆炸。is_categorical 布尔掩码告知 HistGB 哪些特征是分类的，使其能在分箱时采用最优的"分类特征分箱策略"（如对 8 种类别最多生成 8 个分箱而非 255 个）。对于 LightGBM，对应设置 categorical_feature 参数列表与 max_cat_to_onehot=1 禁用 One-Hot，保持两库"原生分类特征"处理方式的一致性。
Adult 基准的"异构数据预处理 → sklearn 原生分类 vs LightGBM 分类"对比流程如下图所示：
fit() 与 predict() 在该脚本中也是通用计时封装：
源码路径：benchmarks/bench_hist_gradient_boosting_adult.py - fit()（27-33行）
def fit(est, data_train, target_train, libname, **fit_params):
print(f"Fitting a {libname} model...") # ① 输出库标识
tic = time() # ② 启动计时
est.fit(data_train, target_train, **fit_params) # ③ 支持额外 fit 参数（如 categorical_feature）
toc = time()
print(f"fitted in {toc - tic:.3f}s")
源码路径：benchmarks/bench_hist_gradient_boosting_adult.py - predict()（35-44行）
def predict(est, data_test, target_test):
if args.no_predict:
return
tic = time()
predicted_test = est.predict(data_test)
predicted_proba_test = est.predict_proba(data_test)
toc = time()
roc_auc = roc_auc_score(target_test, predicted_proba_test[:, 1])
acc = accuracy_score(target_test, predicted_test)
print(f"predicted in {toc - tic:.3f}s, ROC AUC: {roc_auc:.4f}, ACC: {acc:.4f}")
fit() 增加了 **fit_params 形参，用于传递 LightGBM 的 categorical_feature 参数（它需要在 fit 时传入而非构造时设置）。predict() 保持与 HIGGS 脚本相同的 ROC AUC + Accuracy 双指标输出。
fit 与 predict 在 Adult 脚本中的协作流程如下图所示：
81.8 纯分类特征专项基准 —— 序数编码的"极限压力测试"
bench_hist_gradient_boosting_categorical_only.py 是"极限障碍赛"，将连续特征通过 KBinsDiscretizer 离散化为序数编码，所有特征标记为 categorical=True，测试全分类场景下 HistGB 与 LightGBM 的表现。
主流程 __main__ 部分展示了纯分类数据的构造与基准流程：
源码路径：benchmarks/bench_hist_gradient_boosting_categorical_only.py - __main__（37-87行）
X, y = make_classification(n_samples=n_samples, n_features=n_features, random_state=0) # ① 生成连续特征
X = KBinsDiscretizer(n_bins=n_categories, encode="ordinal").fit_transform(X) # ② 离散化为序数编码
print(f"Number of features: {n_features}") # ③ 输出特征数
print(f"Number of samples: {n_samples}") # ④ 输出样本数
is_categorical = [True] * n_features # ⑤ 所有特征标记为分类
est = HistGradientBoostingClassifier(
loss="log_loss",
learning_rate=lr,
max_iter=n_trees,
max_bins=max_bins,
max_leaf_nodes=n_leaf_nodes,
categorical_features=is_categorical, # ⑥ 全分类特征
early_stopping=False,
random_state=0,
verbose=verbose,
)
fit(est, X, y, "sklearn") # ⑦ sklearn 训练
predict(est, X) # ⑧ sklearn 预测（无 train_test_split）
if args.lightgbm: # ⑨ LightGBM 对比
est = get_equivalent_estimator(est, lib="lightgbm", n_classes=2)
est.set_params(max_cat_to_onehot=1) # dont use OHE
categorical_features = list(range(n_features)) # ⑩ 所有特征索引
fit(est, X, y, "lightgbm", categorical_feature=categorical_features)
predict(est, X)
这段代码构建了"全分类特征"场景：先 make_classification 生成连续数值数据，再用 KBinsDiscretizer(n_bins=20, encode="ordinal") 离散化为 20 档的序数编码，最后将所有特征标记为分类。这种设计测试 HistGB 与 LightGBM 在"无连续数值、全分箱"场景下的鲁棒性，考察库对纯离散特征的处理能力。
注意此脚本的 predict() 做了简化：仅计时，不计算 ROC/AUC。因为将连续特征离散化为序数后，原本的"序数关系"被破坏（"第 5 档"在物理上不一定大于"第 3 档"），评估预测准确性没有实际意义。
纯分类基准的"连续→离散→全分类"构造与对比流程如下图所示：
源码路径：benchmarks/bench_hist_gradient_boosting_categorical_only.py - fit()（27-33行）
def fit(est, data_train, target_train, libname, **fit_params):
print(f"Fitting a {libname} model...") # ① 输出库标识
tic = time() # ② 启动计时
est.fit(data_train, target_train, **fit_params) # ③ 支持额外 fit 参数
toc = time()
print(f"fitted in {toc - tic:.3f}s")
源码路径：benchmarks/bench_hist_gradient_boosting_categorical_only.py - predict()（35-41行）
def predict(est, data_test):
# We don't report accuracy or ROC because the dataset doesn't really make
# sense: we treat ordered features as un-ordered categories.
if args.no_predict: # ① 可选：跳过预测阶段
return
tic = time() # ② 启动预测计时
est.predict(data_test) # ③ 仅执行 predict
toc = time() # ④ 结束计时
print(f"predicted in {toc - tic:.3f}s") # ⑤ 仅输出耗时，不计算指标
fit() 保留了 **fit_params 形参以支持 LightGBM 的 categorical_feature 参数传递。predict() 简化至极——只做预测计时，不计算 ROC/AUC 或 Accuracy，因为纯分类场景下指标无物理意义。注释明确说明：ordered features 被当作 un-ordered categories 使用，破坏了原始数据的物理含义，因此放弃准确率评估以聚焦于"全分类场景下的训练/推理速度对比"。
fit 与 predict 在纯分类脚本中的极简协作流程如下图所示：
81.9 大规模真实数据集基准 —— HIGGS 玻色子的"重量级对决"
bench_hist_gradient_boosting_higgsboson.py 是"重量级对决"脚本，使用 UCI 机器学习仓库中的 HIGGS 数据集（2.6GB、1100万样本、28 维特征），考察各库在大规模数据下的内存效率与并行扩展能力。
load_data() 函数负责数据下载、缓存与解析：
源码路径：benchmarks/bench_hist_gradient_boosting_higgsboson.py - load_data()（44-58行）
@m.cache # ① joblib 内存缓存装饰器
def load_data():
filename = os.path.join(HERE, URL.rsplit("/", 1)[-1]) # ② 拼接本地路径
if not os.path.exists(filename): # ③ 若本地不存在
print(f"Downloading {URL} to {filename} (2.6 GB)...") # ④ 输出下载信息
urlretrieve(URL, filename) # ⑤ 从 UCI 下载
print("done.")
print(f"Parsing {filename}...") # ⑥ 输出解析信息
tic = time()
with GzipFile(filename) as f: # ⑦ 打开 gzip 压缩文件
df = pd.read_csv(f, header=None, dtype=np.float32) # ⑧ 用 float32 节省内存
toc = time()
print(f"Loaded {df.values.nbytes / 1e9:0.3f} GB in {toc - tic:0.3f}s")
return df
这段代码是数据加载的"核心入口"：@m.cache 装饰器（其中 m = Memory(location=args.cache_loc, mmap_mode="r")）缓存 DataFrame 解析结果，避免每次运行都重新下载/解析 2.6GB 的巨大文件。mmap_mode="r" 让缓存文件以内存映射方式读取，启动时几乎不占用内存。dtype=np.float32 将 64 位浮点降为 32 位，直接将内存占用减半。
load_data 的下载与解析流程如下图所示：
源码路径：benchmarks/bench_hist_gradient_boosting_higgsboson.py - __main__（60-124行）
df = load_data()
target = df.values[:, 0] # ① 第 0 列是标签
data = np.ascontiguousarray(df.values[:, 1:]) # ② 强制 C-contiguous 内存布局
data_train, data_test, target_train, target_test = train_test_split(
data, target, test_size=0.2, random_state=0
)
n_classes = len(np.unique(target))
if subsample is not None: # ③ 可选：子采样
data_train, target_train = data_train[:subsample], target_train[:subsample]
n_samples, n_features = data_train.shape
print(f"Training set with {n_samples} records with {n_features} features.")
if args.no_interactions: # ④ 可选：禁用特征交互
interaction_cst = [[i] for i in range(n_features)] # ⑤ 每特征独立成组
else:
interaction_cst = None
est = HistGradientBoostingClassifier(
loss="log_loss",
learning_rate=lr,
max_iter=n_trees,
max_bins=max_bins,
max_leaf_nodes=n_leaf_nodes,
early_stopping=False,
random_state=0,
verbose=1,
interaction_cst=interaction_cst, # ⑥ 设置交互约束
max_features=max_features, # ⑦ 特征采样比例
)
fit(est, data_train, target_train, "sklearn") # ⑧ sklearn 训练
predict(est, data_test, target_test) # ⑨ sklearn 评估
if args.lightgbm: # ⑩ 可选：LightGBM 对比
est = get_equivalent_estimator(est, lib="lightgbm", n_classes=n_classes)
fit(est, data_train, target_train, "lightgbm")
predict(est, data_test, target_test)
if args.xgboost: # ⑪ 可选：XGBoost 对比
est = get_equivalent_estimator(est, lib="xgboost", n_classes=n_classes)
fit(est, data_train, target_train, "xgboost")
predict(est, data_test, target_test)
if args.catboost: # ⑫ 可选：CatBoost 对比
est = get_equivalent_estimator(est, lib="catboost", n_classes=n_classes)
fit(est, data_train, target_train, "catboost")
predict(est, data_test, target_test)
HIGGS 基准的完整数据加载→切分→子采样→交互约束→多库对比流程如下图所示：
fit() 与 predict() 在 HIGGS 脚本中的协作流程如下图所示：
在大规模数据集场景下，每个工程配置项都不是"锦上添花"，而是决定基准测试能否跑起来的关键开关。以下表格列出了 HIGGS 基准中几个关键工程配置项的作用与加速原理——其中 mmap_mode="r" 的内存映射机制、Histogram 算法本身的分箱并行性、np.ascontiguousarray 强制 C-contiguous 的缓存友好性三者协同作用，才能让 2.6GB 的全量数据加载与训练成为可能。
| 配置项 | 作用 | 加速原理 |
|--------|------|----------|
| interaction_cst = [[i] for i in range(n_features)] | 禁用特征交互，每特征独立成组 | 分裂时无需考虑特征组合，减少候选分裂点；HistGB 默认会在分裂节点上尝试任意特征对（feature interaction），对 28 维特征而言，禁用后分裂评估从 O(n²) 降至 O(n)，且每个分箱直方图独立构建，并行度大幅提升 |
| subsample | 子采样前 N 条样本 | 减少数据量，快速验证流程；HIGGS 全量 880 万条训练样本完整训练一次需数十分钟，子采样到 10 万条可将单次运行压缩到秒级 |
| max_features | 每棵树仅采样部分特征 | 降低单树计算量，借鉴 Random Forest 思路；max_features=0.5 让每棵树的直方图构建只覆盖一半特征，并行任务量减半，且增强泛化能力 |
| np.ascontiguousarray | 强制 C-contiguous | 满足 Cython 内核的内存布局要求；DataFrame 转出的数组默认是 Fortran-contiguous，Cython 内核遍历会触发大量缓存未命中，强制转换后按行访问变为顺序访问 |
| mmap_mode="r" | 内存映射读取 | 减少内存占用，避免重复加载；缓存文件以 mmap 形式访问，多个进程可共享同一份物理内存，且不占用 RSS，直到实际访问页面才触发缺页中断 |
81.10 设计中的取舍
为什么用 threadpool_limits 而非环境变量控制线程数？ 环境变量（如 OMP_NUM_THREADS）需要在脚本启动前设置，且对所有 OpenMP 库生效，粒度粗、灵活性差。threadpool_limits 是 threadpoolctl 提供的上下文管理器，能在代码块级别精准控制线程数，且支持分别控制 OpenMP、MKL、BLIS 等不同线程池，实现"同一脚本测试不同线程配置"的精细需求。代价是需要为每个库单独设置线程参数（num_threads/nthread/thread_count），增加了一些样板代码。
为什么 bench_hist_gradient_boosting_higgsboson.py 禁用特征交互？ HIGGS 数据集有 28 个物理特征，interaction_cst=[[i] for i in range(n_features)] 让每特征独立成组，禁用了"任意两特征交互"的复杂性。这是一种性能对齐策略：HistGB 支持交互约束而 LightGBM/XGBoost/CatBoost 在某些配置下不支持任意交互，将 HistGB 配置简化到其他库的"最大公约数"水平，让对比更公平。代价是损失了部分建模能力，但在大规模数据基准中，速度收益远大于精度损失。
为什么 categorical_only 基准不计算 ROC/AUC？ 该脚本用 KBinsDiscretizer 将连续特征离散化后再视为分类特征，破坏了原始特征的物理序数关系（如"年龄 5 档"并不必然大于"年龄 3 档"）。在这种情况下计算 ROC/AUC 没有实际意义，模型只是机械地"记住"了分箱编号与标签的对应关系。简化 predict 为纯计时，让基准聚焦于"全分类场景下的训练/推理速度对比"。
为什么用 joblib.Memory 缓存而非手动管理临时文件？ HIGGS 数据集原始压缩包 2.6GB，解析为 DataFrame 后内存占用相近，若每次运行都重新下载与解析，将浪费数分钟甚至数十分钟的 IO 时间。joblib.Memory 装饰器以"hash 函数签名→磁盘缓存"机制自动管理缓存文件：首次调用执行实际下载/解析并缓存到磁盘（默认 /tmp），后续调用直接读取磁盘上的 pickle 文件并通过 mmap_mode="r" 以内存映射方式反序列化。相较于手动管理临时文件，该方案无需关心缓存清理、命名冲突、过期失效等问题，且 Memory(location=args.cache_loc) 支持自定义缓存目录，便于在集群环境下共享缓存。
为什么 HIGGS 脚本使用 float32 而非默认 float64？ 28 维 × 1100 万样本的 DataFrame 若用 float64 存储需约 2.5GB 内存，float32 直接减半至约 1.25GB，这在内存受限或需要多进程并行的环境下是决定性的优化。HistGB 的直方图分箱算法对数值精度不敏感（255 个分箱下，float32 的精度足以区分相邻分箱），因此精度损失可忽略。
为什么 Adult 基准用 ColumnTransformer + OrdinalEncoder 而非 OneHotEncoder？ Adult 数据集有 8 个类别特征，部分特征（如 occupation）有 14+ 种类别，One-Hot 编码会将特征维度从 14 膨胀到 100+，既浪费内存也增加树模型寻找有效分裂的难度。OrdinalEncoder 将每个类别映射为 0/1/2/... 的整数，配合 HistGB 的 categorical_features 布尔掩码，模型会采用专门的"分类特征分箱策略"（如 Fisher 1958 提出的最优划分），在不损失信息的前提下大幅降低计算复杂度。
81.11 动手练习
-
阅读并对比基准脚本的数据生成策略
- 对比
bench_hist_gradient_boosting.py、bench_hist_gradient_boosting_threading.py、bench_hist_gradient_boosting_categorical_only.py三个脚本的数据生成方式：
-
前两者使用
make_classification/make_regression生成连续数值特征 -
后者使用
KBinsDiscretizer将连续特征离散化为序数编码
回答问题：
-
为什么线程扩展性测试 (
bench_hist_gradient_boosting_threading.py) 将n_informative设为n_features // 2而标准测试设为n_classes？ -
纯分类特征基准为何不计算 ROC AUC 和 Accuracy？
-
make_classification的n_clusters_per_class=1和n_informative参数如何影响数据分布与模型训练难度？
- 对比
-
深入理解 get_equivalent_estimator 参数映射机制
- 阅读
sklearn/ensemble/_hist_gradient_boosting/utils.py中的get_equivalent_estimator源码（非提供源码，需自行查阅），重点关注：
-
learning_rate、max_iter、max_leaf_nodes、max_bins等核心参数的映射关系 -
categorical_features布尔数组如何转换为 LightGBM 的categorical_feature索引列表、XGBoost 的enable_categorical、CatBoost 的cat_features -
interaction_cst约束在三个目标库中是否支持及如何映射
回答问题：
-
为什么
bench_hist_gradient_boosting_threading.py中调用get_equivalent_estimator后还需手动set_params(num_threads=...)等线程参数？ -
loss='log_loss'映射到 LightGBM 的objective='binary'/multiclass，XGBoost 的objective='binary:logistic'/multi:softprob，CatBoost 的loss_function='Logloss'/MultiClass，这种映射是否总是精确等价？
- 阅读
-
设计自定义基准：对比 HistGB 与 LightGBM 在 GPU 上的训练性能
- 参考
bench_hist_gradient_boosting_threading.py的结构，设计一个新的基准脚本bench_hist_gpu.py，要求：
-
检测环境是否支持 GPU (LightGBM
device='gpu', XGBoosttree_method='gpu_hist', CatBoosttask_type='GPU')，HistGB 目前无原生 GPU 支持 -
复用
one_run结构，但移除threadpool_limits，改为对比 CPU (sklearn) vs GPU (LightGBM/XGBoost/CatBoost) -
增加
--gpu-id参数支持多 GPU 选择 -
记录显存占用峰值（可选，需
nvidia-smi或pynvml）
回答问题：
-
HistGB 在 CPU 上何时能击败 GPU 上的 LightGBM/XGBoost？（考虑数据传输开销、数据规模阈值）
-
如何修改
get_equivalent_estimator调用以启用目标库的 GPU 模式？
- 参考
-
分析 HIGGS 基准中的内存与缓存策略
- 阅读
bench_hist_gradient_boosting_higgsboson.py的load_data函数与主流程：
-
@m.cache装饰器缓存的是什么？mmap_mode='r'如何避免重复内存映射开销？ -
np.ascontiguousarray(df.values[:, 1:])为什么必要？（提示：HistGB/Cython 内核要求 C-contiguous 内存布局） -
subsample参数如何实现"快速调试模式"？
回答问题：
-
如果不使用
joblib.Memory缓存，每次运行基准测试需重复哪些耗时操作？ -
interaction_cst=[[i] for i in range(n_features)]禁用特征交互后，HistGB 的计算复杂度如何变化？为何能加速训练？ -
对比
max_features=1.0(默认) 与max_features=0.5(类似 Random Forest)，后者如何影响分箱直方图构建与并行度？
- 阅读
81.12 本章小结
这一章中我们深入了梯度提升库的基准测试世界，从标准赛道到多核挑战赛，从合成数据到 2.6GB 的 HIGGS 真实大数据，全面揭示了 HistGradientBoosting 与 LightGBM/XGBoost/CatBoost 四大选手的竞技场设计。首先我们通过 bench_hist_gradient_boosting.py 理解了 one_run 原子操作单元与"训练→计时→评分→计时"的标准化流程；接着剖析了 get_equivalent_estimator 这一跨库翻译官如何实现参数等价映射；然后通过 bench_hist_gradient_boosting_threading.py 学习了 threadpool_limits 控制 OpenMP 线程数与并行扩展性评测方法；接着通过 bench_hist_gradient_boosting_adult.py 与 bench_hist_gradient_boosting_categorical_only.py 理解了真实异构数据与纯分类特征场景的基准设计；最后通过 bench_hist_gradient_boosting_higgsboson.py 掌握了大规模数据集下的 joblib.Memory 缓存、内存映射、交互约束等高级工程技巧。
本章我们一起学习了以下概念：
| 概念 | 解释 |
|------|------|
| one_run() | 基准测试的原子操作单元：固定数据规模/线程数下的"训练→计时→评分→计时"标准化流程 |
| get_equivalent_estimator() | sklearn.ensemble._hist_gradient_boosting.utils 中的跨库参数翻译层，实现 HistGB 参数到 LightGBM/XGBoost/CatBoost 的等价映射 |
| threadpool_limits | 来自 threadpoolctl，在 OpenMP 级别强制限制线程数，消除环境变量干扰，实现精准的并行扩展性对比 |
| HIGGS 基准 | 大规模真实数据集（1100万样本）基准，包含 joblib 缓存、子采样、特征交互约束（interaction_cst）等工程化配置 |
| Adult 数据集基准 | 真实异构表格数据基准，使用 ColumnTransformer + OrdinalEncoder 预处理，对比 sklearn 原生 categorical_features 与 LightGBM categorical_feature 参数 |
| 纯分类特征基准 | KBinsDiscretizer 离散化合成数据，所有特征标记为 categorical，测试全分类场景下 HistGB 与 LightGBM（max_cat_to_onehot=1）的表现 |
| 交互约束 (interaction_cst) | HistGB 独有参数，限定特征交互方式；HIGGS 基准中用 [[i] for i in range(n_features)] 禁用交互，加速训练并对齐其他库行为 |
下一章中，我们将学习 独立基准脚本：PCA 与 SVD 专题——探秘"矩阵分解的效率光谱"，深入 PCA 求解器全景、KernelPCA 扩展性、随机 SVD 的精度与速度博弈以及增量 PCA 在流式数据下的降维表现，从另一个维度理解机器学习库的"性能艺术"。
第 82 章 —— 独立基准脚本：PCA 与 SVD 专题 —— 探秘"矩阵分解的效率光谱"
82.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解 PCA 多求解器在不同数据规模下的性能边界与 'auto' 启发式决策逻辑
-
掌握 KernelPCA 在样本量与成分数两个维度上的扩展性特征及近似求解器的适用场景
-
深入理解随机 SVD 中幂迭代次数与归一化策略对数值稳定性、精度与速度的三重权衡
-
掌握增量 PCA 与批量 PCA 在流式数据场景下的时空权衡与重构误差控制
-
学会设计严谨的基准实验：中位数计时抵抗抖动、正确性守门、分块计算大矩阵范数、缓存与复现性保障
82.2 生活类比
想象你是一位城市规划师，面对一座由数据构成的"数字大楼"。整栋大楼的结构信息都存储在一个巨大的矩阵里，而矩阵分解就是你要拆解这座大楼、提炼其骨架结构的过程。不同场景需要不同的拆解工具：当你需要在极短时间内估算出大楼的承重墙分布（主成分），你会选择抽芯取样的快速勘探法；当你需要为抗震鉴定出具权威报告，你会选择逐砖检查的完整拆解法；当你面对的是一座还在不断加盖新楼层的"活体大楼"（流式数据），你需要流水线式的分批处理。基准测试就像是建筑评测机构，用秒表、标准砝码、体检表和分装箱子给每种拆解方案出具权威报告。全量 SVD（full/dense）相当于"整体拆解法"，将整栋大楼一次性拆成砖块（奇异向量），精确但极其耗时、占场地（内存），适合小楼（小矩阵）。随机 SVD（randomized）则像"抽样勘探法"，只钻取几个芯样（随机投影）推断整体结构，速度快、省场地，配合"幂迭代"（多钻几次）和"归一化"（校准仪器）提高精度，适合摩天大楼（大矩阵）。ARPACK（arpack）则像"定向爆破法"，只炸开需要的楼层（前 k 个成分），中等楼层数高效，但想炸太多层反而慢。协方差特征分解（covariance_eigh）相当于"先压缩再拆"，先把大楼压扁成平面图（协方差矩阵），再拆解，极宽矩阵（特征>>样本）有奇效，但压扁过程本身可能撑爆仓库。增量 PCA（IncrementalPCA）则是"分批拆解流水线"，大楼太大放不下仓库，分批运进来拆、存零件、再组装，牺牲一点精度（近似误差）换取"无限高楼"（流式数据）的处理能力。KernelPCA 则像"非线性镜像重建"，先把大楼投影到高维镜像空间（核矩阵），再在镜像里拆解，镜像大小随样本量平方增长，必须用"抽样勘探"（randomized）或"定向爆破"（arpack）才能应对大规模镜像。
82.3 源码地图
benchmarks/bench_pca_solvers.py
├── __main__ (line 1-200) # PCA 求解器全景基准：数据生成、中位数计时、剪枝策略、CSV 导出、双维度可视化
├── measure_one() (line 45-65) # 单次基准测量：重复执行至累计 0.5s 取中位数，config_context 关闭 NaN 检查
benchmarks/bench_plot_svd.py
├── compute_bench() (line 20-60) # SciPy SVD vs randomized_svd 3D 曲面对比：低秩矩阵生成、网格扫描、surface 绘图
├── __main__ (line 62-100) # 入口：参数网格定义、结果收集、3D 绘图
benchmarks/bench_kernel_pca_solvers_time_vs_n_samples.py
├── __main__ (line 1-180) # KernelPCA 样本扩展性：make_circles 数据、dense/arpack/randomized 三求解器、assert 正确性守门、误差棒绘图
benchmarks/bench_kernel_pca_solvers_time_vs_n_components.py
├── __main__ (line 1-180) # KernelPCA 成分数扩展性：固定样本、指数网格 n_components、arpack 仅小成分测试、对数坐标可视化
benchmarks/bench_plot_randomized_svd.py
├── bench_a() (line 180-280) # 实验 A：固定数据扫描 n_iter 与归一化策略，时间 vs 范数误差
├── bench_b() (line 280-350) # 实验 B：可控秩矩阵，变化 n_components 与真实秩关系，研究最优 n_iter
├── bench_c() (line 350-420) # 实验 C：跨真实数据集对比 sklearn 默认配置与 fbpca
├── scalable_frobenius_norm_discrepancy() (line 150-180) # 分块计算 Frobenius 范数，避免 OOM
├── norm_diff() (line 130-150) # 谱范数(ARPACK svds)与 Frobenius 范数统一接口
├── svd_timing() (line 200-230) # 统一计时封装：randomized_svd 与 fbpca.pca 不同参数约定
├── get_data() (line 90-150) # 多数据源统一获取：fetch_*、手动下载(CIFAR/SVHN)、合成数据
├── plot_time_vs_s() (line 60-80) # 时间 vs 范数折线图绘制
├── scatter_time_vs_s() (line 80-110) # 时间 vs 范数散点图绘制
├── plot_power_iter_vs_s() (line 110-130) # n_iter vs 范数曲线绘制
├── handle_missing_dataset() (line 85-90) # 手动下载数据集缺失检查
├── unpickle() (line 80-85) # CIFAR 批次文件解包
├── __main__ (line 420-450) # 入口：三大实验串联、绘图展示
benchmarks/bench_plot_incremental_pca.py
├── benchmark() (line 30-50) # 单模型基准：fit/transform/inverse_transform 计时与 MAE 重构误差
├── fixed_batch_size_comparison() (line 90-120) # 固定批次、变 n_components：PCA vs IncrementalPCA
├── variable_batch_size_comparison() (line 120-160) # 固定成分、变批次：引入 randomized PCA 基线
├── plot_results() (line 20-25) # 通用曲线绘制
├── plot_feature_times() (line 50-60) # 成分数-时间对比图
├── plot_feature_errors() (line 60-70) # 成分数-误差对比图
├── plot_batch_times() (line 70-80) # 批次-时间对比图
├── plot_batch_errors() (line 80-90) # 批次-误差对比图
├── __main__ (line 160-190) # 入口：LFW 数据加载、标准化、两大实验调用、绘图展示
82.4 PCA 求解器全景对比 —— 从 full 到 randomized 的"效率光谱"
scikit-learn 的 PCA 提供 full、covariance_eigh、arpack、randomized、auto 五种 SVD 求解器，性能随数据规模与 n_components 变化剧烈。auto 启发式需在经验数据上验证其决策边界是否最优。我们通过中位数计时（measure_one 重复执行至累计 0.5 秒）抵抗抖动，输出 CSV 供离线分析。
核心实验设计围绕三个基准维度（100、1000、10000），生成 (n_samples, n_features) 组合，覆盖"高宽/宽高"矩阵。在双维度可视化上，我们分别固定 n_features 画 time vs n_samples，固定 n_samples 画 time vs n_features。同时测 fit 与 fit_transform，n_components=2 与 50 代表"极少/中等"成分数。
关键剪枝策略包括：covariance_eigh 在 n_features > 5000 时跳过（内存/时间不可接受）；arpack 与 full 在 log10(data.size) > 7 时跳过（极慢）。config_context(assume_finite=True) 关闭 NaN 检查以减少开销。
bench_plot_svd.py 的 3D 曲面对比直接展示 SciPy 标准 svd 与 randomized_svd (n_iter=0/3)，使用 make_low_rank_matrix 生成带有效秩与长尾的合成数据，模拟真实谱衰减。3D 曲面图（plot_surface）直观展示 n_samples × n_features 平面上的耗时地形。
82.4.1 measure_one 中位数计时实现
源码路径：benchmarks/bench_pca_solvers.py - measure_one()（45-65行）
def measure_one(data, n_components, solver, method_name="fit"):
print(
f"Benchmarking {solver=!r}, {n_components=}, {method_name=!r} on data with"
f" shape {data.shape}"
)
# ① 构造 PCA 估计器，固定随机种子保证可复现
pca = PCA(n_components=n_components, svd_solver=solver, random_state=0)
timings = [] # 存储每次执行的耗时
elapsed = 0 # 累计耗时
method = getattr(pca, method_name) # 获取 fit 或 fit_transform 方法
# ② config_context 临时关闭 NaN 检查（assume_finite=True），减少检查开销
with config_context(assume_finite=True):
# ③ 循环执行直至总耗时达到 0.5 秒
while elapsed < 0.5:
tic = perf_counter() # 记录开始时间
method(data) # 执行 fit/fit_transform
duration = perf_counter() - tic # 计算本次耗时
timings.append(duration) # 追加到耗时列表
elapsed += duration # 累加总耗时
# ④ 返回所有耗时样本的中位数，抗抖动
return np.median(timings)
这段代码定义了 measure_one 函数，实现单次基准测量。它通过 while 循环累计执行至 0.5 秒，然后取所有单次耗时的中位数，能够有效抵抗单次抖动。中位数比平均值更稳健，因为计时噪声常常出现长尾分布。
82.4.2 主循环：数据网格生成与剪枝策略
源码路径：benchmarks/bench_pca_solvers.py - __main__（1-200行）
# 第 82 章 —— 三个基准维度，围绕它们生成数据形状网格
REF_DIMS = [100, 1000, 10_000]
data_shapes = []
for ref_dim in REF_DIMS:
# 扩展 (ref_dim, n_features) 与 (n_features, ref_dim) 两种方向
data_shapes.extend([(ref_dim, 10**i) for i in range(1, 8 - int(log10(ref_dim)))])
data_shapes.extend([(ref_dim, 3 * 10**i) for i in range(1, 8 - int(log10(ref_dim)))])
data_shapes.extend([(10**i, ref_dim) for i in range(1, 8 - int(log10(ref_dim)))])
data_shapes.extend([(3 * 10**i, ref_dim) for i in range(1, 8 - int(log10(ref_dim)))])
# 第 82 章 —— 去重并排序
data_shapes = sorted(set(data_shapes))
print("Generating test datasets...")
rng = np.random.default_rng(0)
# 第 82 章 —— 生成对应形状的高斯随机矩阵
datasets = [rng.normal(size=shape) for shape in data_shapes]
# 第 82 章 —— 遍历所有数据集 × 成分数 × 方法的笛卡尔积
SOLVERS = ["full", "covariance_eigh", "arpack", "randomized", "auto"]
measurements = []
for data, n_components, method_name in itertools.product(
datasets, [2, 50], ["fit", "fit_transform"]
):
# 跳过 n_components 大于数据形状最小维度的情况
if n_components >= min(data.shape):
continue
for solver in SOLVERS:
# 剪枝1：covariance_eigh 在特征数过大时内存不可接受
if solver == "covariance_eigh" and data.shape[1] > 5000:
continue
# 剪枝2：arpack 与 full 在矩阵元素总数过大时过慢
if solver in ["arpack", "full"] and log10(data.size) > 7:
continue
time = measure_one(data, n_components, solver, method_name=method_name)
measurements.append(
{
"n_components": n_components,
"n_samples": data.shape[0],
"n_features": data.shape[1],
"time": time,
"solver": solver,
"method_name": method_name,
}
)
# 第 82 章 —— 保存为 CSV，便于后续分析
measurements = pd.DataFrame(measurements)
measurements.to_csv("bench_pca_solvers.csv", index=False)
这段代码实现了 PCA 求解器全景基准的主循环。三个核心设计要点是：（1）REF_DIMS 围绕 100/1000/10000 三个数量级扩展，覆盖"高瘦/矮胖"两类矩阵；（2）剪枝策略跳过明显不可行的组合，避免长时间空转；（3）结果以 DataFrame 存储并导出 CSV，方便后续 pandas 查询与可视化。
82.4.3 3D 曲面绘制：SciPy SVD vs randomized_svd
源码路径：benchmarks/bench_plot_svd.py - compute_bench()（1-100行）
def compute_bench(samples_range, features_range, n_iter=3, rank=50):
it = 0
results = defaultdict(lambda: []) # 每个算法名对应一个耗时列表
max_it = len(samples_range) * len(features_range)
# 网格扫描：每个 (n_samples, n_features) 组合
for n_samples in samples_range:
for n_features in features_range:
it += 1
print("====================")
print("Iteration %03d of %03d" % (it, max_it))
print("====================")
# 生成低秩矩阵：有效秩 50，尾部强度 0.2，模拟谱衰减
X = make_low_rank_matrix(
n_samples, n_features, effective_rank=rank, tail_strength=0.2
)
gc.collect() # 清理上一轮内存，避免累计
# 基准1：SciPy 精确 SVD（full_matrices=False 节省内存）
print("benchmarking scipy svd: ")
tstart = time()
svd(X, full_matrices=False)
results["scipy svd"].append(time() - tstart)
gc.collect()
# 基准2：randomized_svd 无幂迭代
print("benchmarking scikit-learn randomized_svd: n_iter=0")
tstart = time()
randomized_svd(X, rank, n_iter=0)
results["scikit-learn randomized_svd (n_iter=0)"].append(time() - tstart)
gc.collect()
# 基准3：randomized_svd 带幂迭代（默认 n_iter=4）
print("benchmarking scikit-learn randomized_svd: n_iter=%d " % n_iter)
tstart = time()
randomized_svd(X, rank, n_iter=n_iter)
results["scikit-learn randomized_svd (n_iter=%d)" % n_iter].append(
time() - tstart
)
return results
这段代码定义了 compute_bench 函数，对比 SciPy 精确 SVD 与 sklearn 的 randomized_svd。它使用 make_low_rank_matrix 生成具有指定有效秩和长尾强度的合成数据，模拟真实世界谱衰减场景。3D 曲面图（plot_surface）会直观展示 n_samples × n_features 平面上的耗时地形，让随机算法的加速优势一目了然。
82.5 KernelPCA 扩展性实验 —— 样本与维度的"双重考验"
KernelPCA 提供了 dense、arpack、randomized 三种 eigen_solver，在样本量与目标成分数两个维度上的扩展性差异巨大。我们通过两个正交实验揭示这些差异。
82.5.1 样本扩展性：n_samples 从 101 到 4000
bench_kernel_pca_solvers_time_vs_n_samples.py 固定 n_components=100、n_features=2，让 n_samples 从 101 增至 4000，观察三种求解器的时间复杂度分化。使用 make_circles 生成非线性可分数据，保证核矩阵有非平凡谱结构。每个配置重复 n_iter=3 次取均值±标准差，误差棒可视化。
关键工程实践是正确性守门：assert_array_almost_equal(np.abs(a_pred), np.abs(ref_pred)) 确保近似解在符号不定性下与 dense 解一致。奇异向量符号本身是不定的（如果 v 是奇异向量，-v 也是），直接比较会失败，取绝对值后再比较是标准做法。
源码路径：benchmarks/bench_kernel_pca_solvers_time_vs_n_samples.py - __main__（1-180行）
# 第 82 章 —— 1- 设计实验：n_samples 网格
min_n_samples, max_n_samples = 101, 4000
n_samples_grid_size = 4
n_samples_range = [
min_n_samples
+ np.floor((x / (n_samples_grid_size - 1)) * (max_n_samples - min_n_samples))
for x in range(0, n_samples_grid_size)
]
n_components = 100
n_iter = 3
include_arpack = False # arpack 在样本量大时过慢，默认关闭
# 第 82 章 —— 2- 生成非线性数据：make_circles 制造同心圆结构
n_features = 2
X, y = make_circles(n_samples=max_n_samples, factor=0.3, noise=0.05, random_state=0)
# 第 82 章 —— 3- 基准测试循环
ref_time = np.empty((len(n_samples_range), n_iter)) * np.nan
a_time = np.empty((len(n_samples_range), n_iter)) * np.nan
r_time = np.empty((len(n_samples_range), n_iter)) * np.nan
for j, n_samples in enumerate(n_samples_range):
n_samples = int(n_samples)
print("Performing kPCA with n_samples = %i" % n_samples)
X_train = X[:n_samples, :]
X_test = X_train
# A- 参考（dense）求解器，全特征分解
print(" - dense")
for i in range(n_iter):
start_time = time.perf_counter()
ref_pred = (
KernelPCA(n_components, eigen_solver="dense").fit(X_train).transform(X_test)
)
ref_time[j, i] = time.perf_counter() - start_time
# B- arpack（可选）
if include_arpack:
print(" - arpack")
for i in range(n_iter):
start_time = time.perf_counter()
a_pred = (
KernelPCA(n_components, eigen_solver="arpack").fit(X_train).transform(X_test)
)
a_time[j, i] = time.perf_counter() - start_time
# 正确性守门：符号不定性下取绝对值比较
assert_array_almost_equal(np.abs(a_pred), np.abs(ref_pred))
# C- randomized 近似求解器
print(" - randomized")
for i in range(n_iter):
start_time = time.perf_counter()
r_pred = (
KernelPCA(n_components, eigen_solver="randomized").fit(X_train).transform(X_test)
)
r_time[j, i] = time.perf_counter() - start_time
# 正确性守门
assert_array_almost_equal(np.abs(r_pred), np.abs(ref_pred))
这段代码实现了 KernelPCA 样本扩展性基准。核心观察是：dense 求解器随 n_samples 三次方增长（因为核矩阵是 N×N 矩阵，特征分解复杂度 O(N³)）；randomized 在大样本/中等成分数区间呈近似线性扩展；arpack 是中间态，仅在 n_components < 100 时高效。
82.5.2 成分数扩展性：n_components 从 1 到 1999
bench_kernel_pca_solvers_time_vs_n_components.py 固定 n_train=2000、n_features=2，让 n_components 从 1 指数增长至 1999，验证小 n_components 下 randomized 的优势。注意 n_components 网格使用 np.exp(np.log(max_n_compo)) 的对数指数采样，而非线性网格，这是为了在双对数坐标上获得均匀的采样点。
源码路径：benchmarks/bench_kernel_pca_solvers_time_vs_n_components.py - __main__（1-180行）
# 第 82 章 —— 1- 设计实验：指数网格
n_train, n_test = 2000, 1000
max_n_compo = 1999
n_compo_grid_size = 10
# 第 82 章 —— 对数均匀网格：在 [1, 1999] 上指数采样
n_compo_range = [
np.round(np.exp((x / (n_compo_grid_size - 1)) * np.log(max_n_compo)))
for x in range(0, n_compo_grid_size)
]
n_iter = 3
arpack_all = False # 仅在小成分时启用 arpack
# 第 82 章 —— 2- 生成数据
n_features = 2
X, y = make_circles(
n_samples=(n_train + n_test), factor=0.3, noise=0.05, random_state=0
)
X_train, X_test = X[:n_train, :], X[n_train:, :]
# 第 82 章 —— 3- 循环主体（与 n_samples 基准类似，但 arpack 仅在 n_components < 100 时测试）
for j, n_components in enumerate(n_compo_range):
n_components = int(n_components)
print("Performing kPCA with n_components = %i" % n_components)
# A- dense
print(" - dense solver")
for i in range(n_iter):
start_time = time.perf_counter()
ref_pred = (
KernelPCA(n_components, eigen_solver="dense").fit(X_train).transform(X_test)
)
ref_time[j, i] = time.perf_counter() - start_time
# B- arpack（仅小成分数）
if arpack_all or n_components < 100:
print(" - arpack solver")
for i in range(n_iter):
start_time = time.perf_counter()
a_pred = (
KernelPCA(n_components, eigen_solver="arpack").fit(X_train).transform(X_test)
)
a_time[j, i] = time.perf_counter() - start_time
assert_array_almost_equal(np.abs(a_pred), np.abs(ref_pred))
# C- randomized
print(" - randomized solver")
for i in range(n_iter):
start_time = time.perf_counter()
r_pred = (
KernelPCA(n_components, eigen_solver="randomized").fit(X_train).transform(X_test)
)
r_time[j, i] = time.perf_counter() - start_time
assert_array_almost_equal(np.abs(r_pred), np.abs(ref_pred))
# 第 82 章 —— 4- 绘图：对数 x 轴
ax.set_xscale("log")
ax.set_xlim(1, max(n_compo_range) * 1.1)
这段代码实现了成分数扩展性基准。关键设计是 arpack_all 开关与 n_components < 100 条件——因为 ARPACK 在大成分数时退化为全谱计算，耗时爆炸。对数 x 轴让指数增长的成分数与线性增长的耗时在同一张图上呈现。
82.6 随机 SVD 深入探秘 —— 近似算法的"精度与速度博弈"
随机 SVD 的精度受三大因素影响：幂迭代次数 n_iter、归一化策略 power_iteration_normalizer、过采样数 n_oversamples。本节通过三大实验量化这些参数的影响。在深入实验之前，我们先介绍脚本中复用度最高的三组辅助函数：norm_diff 统一谱范数与 Frobenius 范数的计算接口；svd_timing 屏蔽 sklearn 与 fbpca 的参数约定差异；get_data 统一管理多数据源的获取与缺失检查。
82.6.1 辅助函数：范数计算、计时封装与数据获取
源码路径：benchmarks/bench_plot_randomized_svd.py - norm_diff()（130-150行）、svd_timing()（200-230行）、get_data()（90-150行）
def norm_diff(A, norm=2, msg=True, random_state=None):
"""
Compute the norm diff with the original matrix, when randomized
SVD is called with *params.
norm: 2 => spectral; 'fro' => Frobenius
"""
if msg:
print("... computing %s norm ..." % norm)
if norm == 2:
# s = sp.linalg.norm(A, ord=2) # slow
# ARPACK 仅迭代求最大奇异值，比 sp.linalg.norm(A, ord=2) 快几个数量级
v0 = _init_arpack_v0(min(A.shape), random_state)
value = sp.sparse.linalg.svds(A, k=1, return_singular_vectors=False, v0=v0)
else:
if sp.sparse.issparse(A):
value = sp.sparse.linalg.norm(A, ord=norm)
else:
value = sp.linalg.norm(A, ord=norm)
return value
def svd_timing(
X, n_comps, n_iter, n_oversamples, power_iteration_normalizer="auto", method=None
):
"""
Measure time for decomposition
"""
print("... running SVD ...")
if method != "fbpca":
gc.collect()
t0 = time()
# sklearn 约定：n_oversamples 是单独的过采样数
U, mu, V = randomized_svd(
X,
n_comps,
n_oversamples=n_oversamples,
n_iter=n_iter,
power_iteration_normalizer=power_iteration_normalizer,
random_state=random_state,
transpose=False,
)
call_time = time() - t0
else:
gc.collect()
t0 = time()
# fbpca 约定：l 是过采样后的总列数（l = n_oversamples + n_comps）
U, mu, V = fbpca.pca(
X, n_comps, raw=True, n_iter=n_iter, l=n_oversamples + n_comps
)
call_time = time() - t0
return U, mu, V, call_time
def get_data(dataset_name):
print("Getting dataset: %s" % dataset_name)
# 统一封装 fetch_*、make_*、手动下载(CIFAR/SVHN)三大类数据源
if dataset_name == "lfw_people":
X = fetch_lfw_people().data
elif dataset_name == "20newsgroups":
X = fetch_20newsgroups_vectorized().data[:, :100000]
elif dataset_name == "olivetti_faces":
X = fetch_olivetti_faces().data
elif dataset_name == "rcv1":
X = fetch_rcv1().data
elif dataset_name == "CIFAR":
if handle_missing_dataset(CIFAR_FOLDER) == 0:
return # 目录不存在时优雅跳过
X1 = [unpickle("%sdata_batch_%d" % (CIFAR_FOLDER, i + 1)) for i in range(5)]
X = np.vstack(X1)
del X1
elif dataset_name == "SVHN":
if handle_missing_dataset(SVHN_FOLDER) == 0:
return
X1 = sp.io.loadmat("%strain_32x32.mat" % SVHN_FOLDER)["X"]
X2 = [X1[:, :, :, i].reshape(32 * 32 * 3) for i in range(X1.shape[3])]
X = np.vstack(X2)
del X1
del X2
elif dataset_name == "low rank matrix":
X = make_low_rank_matrix(
n_samples=500,
n_features=int(1e4),
effective_rank=100,
tail_strength=0.5,
random_state=random_state,
)
elif dataset_name == "uncorrelated matrix":
X, _ = make_sparse_uncorrelated(
n_samples=500, n_features=10000, random_state=random_state
)
else:
X = fetch_openml(dataset_name).data
return X
这三个辅助函数承担了脚本中绝大部分横向职责。norm_diff 的核心设计是：当 norm=2 时调用 svds(k=1) 替代 sp.linalg.norm(A, ord=2)，精确谱范数需要全 SVD 复杂度 O(N³)，而 ARPACK 只迭代找最大奇异值，复杂度 O(N² × 迭代次数)，大矩阵上快几个数量级。svd_timing 的关键设计在于屏蔽 sklearn 与 fbpca 的 API 差异：sklearn 的 n_oversamples 是单独的过采样数，而 fbpca 的 l 参数是过采样后的总列数（l = n_oversamples + n_comps），因此在 fbpca 分支要做参数映射。get_data 的核心设计是缺失检查的优雅处理——当 CIFAR 或 SVHN 的本地目录不存在时，handle_missing_dataset 返回 0，get_data 返回 None，外层循环 if X is None: continue 跳过该数据集而不是崩溃中断整个基准。
82.6.2 实验 A：固定数据扫描 n_iter 与归一化策略
bench_a(X, dataset_name, power_iter, n_oversamples, n_comps) 固定数据集，扫描 n_iter=0..5 与归一化策略 none/LU/QR/auto，绘制 时间 vs 范数误差 折线图，评估幂迭代收益与数值稳定性权衡。
源码路径：benchmarks/bench_plot_randomized_svd.py - bench_a()（180-280行）
def bench_a(X, dataset_name, power_iter, n_oversamples, n_comps):
all_time = defaultdict(list)
# 谱范数（默认关闭，太慢）
if enable_spectral_norm:
all_spectral = defaultdict(list)
X_spectral_norm = norm_diff(X, norm=2, msg=False, random_state=0)
# Frobenius 范数（默认使用，速度更快）
all_frobenius = defaultdict(list)
X_fro_norm = norm_diff(X, norm="fro", msg=False)
# 扫描 n_iter 与归一化策略
for pi in power_iter:
for pm in ["none", "LU", "QR"]:
print("n_iter = %d on sklearn - %s" % (pi, pm))
# svd_timing 返回 (U, s, V, time)，下方计算相对 Frobenius 误差
U, s, V, time = svd_timing(
X,
n_comps,
n_iter=pi,
power_iteration_normalizer=pm,
n_oversamples=n_oversamples,
)
label = "sklearn - %s" % pm
all_time[label].append(time)
# 分块计算 Frobenius 范数误差，并归一化为相对误差
f = scalable_frobenius_norm_discrepancy(X, U, s, V)
all_frobenius[label].append(f / X_fro_norm)
# 与 fbpca 对比（如果可用）
if fbpca_available:
print("n_iter = %d on fbca" % (pi))
U, s, V, time = svd_timing(
X,
n_comps,
n_iter=pi,
power_iteration_normalizer=pm,
n_oversamples=n_oversamples,
method="fbpca",
)
label = "fbpca"
all_time[label].append(time)
f = scalable_frobenius_norm_discrepancy(X, U, s, V)
all_frobenius[label].append(f / X_fro_norm)
# 绘图
if enable_spectral_norm:
title = "%s: spectral norm diff vs running time" % (dataset_name)
plot_time_vs_s(all_time, all_spectral, power_iter, title)
title = "%s: Frobenius norm diff vs running time" % (dataset_name)
plot_time_vs_s(all_time, all_frobenius, power_iter, title)
这段代码实现了实验 A，扫描幂迭代次数与归一化策略。关键观察：随着 n_iter 增加，误差下降但耗时上升；QR 归一化最稳健但计算量大；LU 比 QR 快且精度相当；auto 在 n_iter<=2 时使用 none，否则使用 LU，是经验最优解。
82.6.3 分块 Frobenius 范数：scalable_frobenius_norm_discrepancy
当重构矩阵过大（>4GB）无法装入内存时，必须分块计算范数。
源码路径：benchmarks/bench_plot_randomized_svd.py - scalable_frobenius_norm_discrepancy()（150-180行）
MAX_MEMORY = int(4e9) # 4GB 内存阈值
def scalable_frobenius_norm_discrepancy(X, U, s, V):
# 如果输入是稀疏但不大，或不是稀疏，可一次性重建
if not sp.sparse.issparse(X) or (
X.shape[0] * X.shape[1] * X.dtype.itemsize < MAX_MEMORY
):
A = X - U.dot(np.diag(s).dot(V)) # 一次性重构误差矩阵
return norm_diff(A, norm="fro")
# 否则按 batch_size 1000 分块累加范数平方
print("... computing fro norm by batches...")
batch_size = 1000
Vhat = np.diag(s).dot(V)
cum_norm = 0.0
for batch in gen_batches(X.shape[0], batch_size):
# 只取当前批次的行参与计算
M = X[batch, :] - U[batch, :].dot(Vhat)
cum_norm += norm_diff(M, norm="fro", msg=False)
# Frobenius 范数是各块范数平方和的平方根
return np.sqrt(cum_norm)
这段代码实现了分块 Frobenius 范数计算。它通过 gen_batches 生成行索引批次，每次只计算当前批次的误差矩阵范数平方并累加，最终开平方得到完整范数。这避免了 OOM 风险，使大规模稀疏矩阵的近似误差量化成为可能。
82.6.4 实验 B：可控秩矩阵与目标秩关系
源码路径：benchmarks/bench_plot_randomized_svd.py - bench_b()（280-350行）
def bench_b(power_list):
n_samples, n_features = 1000, 10000
data_params = {
"n_samples": n_samples,
"n_features": n_features,
"tail_strength": 0.7, # 较慢的谱衰减
"random_state": random_state,
}
dataset_name = "low rank matrix %d x %d" % (n_samples, n_features)
ranks = [10, 50, 100] # 三个真实秩
if enable_spectral_norm:
all_spectral = defaultdict(list)
all_frobenius = defaultdict(list)
for rank in ranks:
# 构造指定有效秩的矩阵
X = make_low_rank_matrix(effective_rank=rank, **data_params)
if enable_spectral_norm:
X_spectral_norm = norm_diff(X, norm=2, msg=False, random_state=0)
X_fro_norm = norm_diff(X, norm="fro", msg=False)
# 目标成分数取 rank/2, rank, 2*rank 三档
for n_comp in [int(rank / 2), rank, rank * 2]:
label = "rank=%d, n_comp=%d" % (rank, n_comp)
print(label)
for pi in power_list:
# LU 归一化在多次幂迭代下数值稳定
U, s, V, _ = svd_timing(
X,
n_comp,
n_iter=pi,
n_oversamples=2,
power_iteration_normalizer="LU",
)
f = scalable_frobenius_norm_discrepancy(X, U, s, V)
all_frobenius[label].append(f / X_fro_norm)
# 绘图：n_iter vs 范数
title = "%s: Frobenius norm diff vs n power iteration" % (dataset_name)
plot_power_iter_vs_s(power_iter, all_frobenius, title)
这段代码实现了实验 B，研究"目标成分数与矩阵真实秩的关系"如何影响最优 n_iter。当 n_comp 小于真实秩时，少量幂迭代即可收敛；当 n_comp 接近或大于真实秩时，可能需要更多幂迭代才能准确捕捉长尾结构。结论是 n_iter=2 几乎是普适的最优默认值。
82.6.5 实验 C：跨真实数据集对比 fbpca
bench_c(datasets, n_comps) 跨多个真实数据集（low rank matrix、lfw_people、olivetti_faces、20newsgroups、mnist_784、CIFAR、a3a、SVHN、uncorrelated matrix）对比 sklearn 默认配置与 fbpca。get_data() 统一封装 fetch_*、make_* 与手动下载（CIFAR/SVHN），缺失目录优雅跳过。
源码路径：benchmarks/bench_plot_randomized_svd.py - bench_c()（350-420行）
def bench_c(datasets, n_comps):
all_time = defaultdict(list)
all_frobenius = defaultdict(list)
for dataset_name in datasets:
X = get_data(dataset_name)
if X is None: # 数据集缺失时跳过
continue
X_fro_norm = norm_diff(X, norm="fro", msg=False)
n_comps = np.minimum(n_comps, np.min(X.shape))
# sklearn 默认配置：n_iter=2, n_oversamples=10
label = "sklearn"
print("%s %d x %d - %s" % (dataset_name, X.shape[0], X.shape[1], label))
U, s, V, time = svd_timing(X, n_comps, n_iter=2, n_oversamples=10, method=label)
all_time[label].append(time)
f = scalable_frobenius_norm_discrepancy(X, U, s, V)
all_frobenius[label].append(f / X_fro_norm)
# fbpca 默认配置：n_iter=2, l=n_oversamples+n_comps
if fbpca_available:
label = "fbpca"
print("%s %d x %d - %s" % (dataset_name, X.shape[0], X.shape[1], label))
U, s, V, time = svd_timing(
X, n_comps, n_iter=2, n_oversamples=2, method=label
)
all_time[label].append(time)
f = scalable_frobenius_norm_discrepancy(X, U, s, V)
all_frobenius[label].append(f / X_fro_norm)
if len(all_time) == 0:
raise ValueError("No tests ran. Aborting.")
title = "normalized Frobenius norm diff vs running time"
scatter_time_vs_s(all_time, all_frobenius, datasets, title)
这段代码实现了实验 C，跨真实数据集对比。注意 svd_timing 在 fbpca 分支使用 l=n_oversamples+n_comps，因为 fbpca 的 l 参数表示"过采样后的总列数"而非单独的过采样数，这是与 sklearn 约定的主要差异。跨数据集散点图直观显示 sklearn 默认配置在大多数场景下达到 fbpca 相当或更优的精度-速度平衡。
82.7 增量 PCA 与批次学习 —— 流式数据下的"降维战士"
IncrementalPCA 通过分批增量 SVD 实现"无限高楼"的流式数据降维能力。本节通过两个维度对比 PCA（全量批量）、IncrementalPCA（分批增量）与 PCA(svd_solver='randomized')（随机近似），揭示时空权衡与重构误差控制。
82.7.1 绘图辅助函数
源码路径：benchmarks/bench_plot_incremental_pca.py - plot_results()（20-25行）、plot_feature_times()（50-60行）、plot_feature_errors()（60-70行）、plot_batch_times()（70-80行）、plot_batch_errors()（80-90行）
def plot_results(X, y, label):
# 通用曲线绘制：仅 plt.plot 一条折线 + 散点 marker
plt.plot(X, y, label=label, marker="o")
def plot_feature_times(all_times, batch_size, all_components, data):
# 成分数-时间对比图：标题含数据规模，X 轴为成分数，Y 轴为时间
plt.figure()
plot_results(all_components, all_times["pca"], label="PCA")
plot_results(
all_components, all_times["ipca"], label="IncrementalPCA, bsize=%i" % batch_size
)
plt.legend(loc="upper left")
plt.suptitle(
"Algorithm runtime vs. n_components\n LFW, size %i x %i"
% data.shape
)
plt.xlabel("Number of components (out of max %i)" % data.shape[1])
plt.ylabel("Time (seconds)")
def plot_feature_errors(all_errors, batch_size, all_components, data):
# 成分数-误差对比图：Y 轴为 MAE 重构误差
plt.figure()
plot_results(all_components, all_errors["pca"], label="PCA")
plot_results(
all_components,
all_errors["ipca"],
label="IncrementalPCA, bsize=%i" % batch_size,
)
plt.legend(loc="lower left")
plt.suptitle("Algorithm error vs. n_components\nLFW, size %i x %i" % data.shape)
plt.xlabel("Number of components (out of max %i)" % data.shape[1])
plt.ylabel("Mean absolute error")
def plot_batch_times(all_times, n_features, all_batch_sizes, data):
# 批次-时间对比图：X 轴为批次大小，标题含 n_components
plt.figure()
plot_results(all_batch_sizes, all_times["pca"], label="PCA")
plot_results(all_batch_sizes, all_times["ipca"], label="IncrementalPCA")
plt.legend(loc="lower left")
plt.suptitle(
"Algorithm runtime vs. batch_size for n_components %i\n LFW,"
" size %i x %i" % (n_features, data.shape[0], data.shape[1])
)
plt.xlabel("Batch size")
plt.ylabel("Time (seconds)")
def plot_batch_errors(all_errors, n_features, all_batch_sizes, data):
# 批次-误差对比图：Y 轴为 MAE，标题含 n_components
plt.figure()
plot_results(all_batch_sizes, all_errors["pca"], label="PCA")
plot_results(all_batch_sizes, all_errors["ipca"], label="IncrementalPCA")
plt.legend(loc="lower left")
plt.suptitle(
"Algorithm error vs. batch_size for n_components %i\n LFW,"
" size %i x %i" % (n_features, data.shape[0], data.shape[1])
)
plt.xlabel("Batch size")
plt.ylabel("Mean absolute error")
这五个绘图函数遵循同一套命名约定：plot_results 是最底层的 plt.plot 封装；plot_feature_* 系列横轴为 n_components；plot_batch_* 系列横轴为 batch_size；*_times 系列纵轴为时间；*_errors 系列纵轴为 MAE 重构误差。统一的标题格式（含数据规模）、统一的图例位置（左上/左下）、统一的标签文案，让最终生成的 4×2 = 8 张子图具备视觉一致性，便于读者横向对比。
82.7.2 双指标基准：benchmark 函数
源码路径：benchmarks/bench_plot_incremental_pca.py - benchmark()（30-50行）
def benchmark(estimator, data):
gc.collect() # 计时前清理内存，避免上一轮分配干扰
print("Benching %s" % estimator)
t0 = time()
estimator.fit(data) # 计时 fit
training_time = time() - t0
data_t = estimator.transform(data) # transform
data_r = estimator.inverse_transform(data_t) # inverse_transform 回原空间
# 重构误差：MAE（平均绝对误差）而非 MSE
reconstruction_error = np.mean(np.abs(data - data_r))
return {"time": training_time, "error": reconstruction_error}
这段代码定义了 benchmark 函数，单模型基准测试入口。它使用 gc.collect() 在计时前清理垃圾，避免上一次迭代的内存分配残留干扰本次计时。重构误差用 MAE 而非 MSE，这是因为人脸图像的像素值差异用绝对值更直观，且 MAE 对极端噪声点更稳健。
82.7.3 固定批次大小、变化成分数
源码路径：benchmarks/bench_plot_incremental_pca.py - fixed_batch_size_comparison()（90-120行）
def fixed_batch_size_comparison(data):
# 5 个成分数点：从 n_features//10 到 n_features
all_features = [
i.astype(int) for i in np.linspace(data.shape[1] // 10, data.shape[1], num=5)
]
batch_size = 1000
all_times = defaultdict(list)
all_errors = defaultdict(list)
for n_components in all_features:
pca = PCA(n_components=n_components)
ipca = IncrementalPCA(n_components=n_components, batch_size=batch_size)
# 同时基准 PCA 和 IncrementalPCA
results_dict = {
k: benchmark(est, data) for k, est in [("pca", pca), ("ipca", ipca)]
}
for k in sorted(results_dict.keys()):
all_times[k].append(results_dict[k]["time"])
all_errors[k].append(results_dict[k]["error"])
plot_feature_times(all_times, batch_size, all_features, data)
plot_feature_errors(all_errors, batch_size, all_features, data)
这段代码实现了固定批次大小、变化成分数的对比。linspace 在 5 个等间距点上采样成分数，覆盖从少量到全量成分的完整区间。绘图函数 plot_feature_times 和 plot_feature_errors 分别绘制时间-成分数与误差-成分数曲线。
82.7.4 固定成分数、变化批次大小（引入 randomized 基线）
源码路径：benchmarks/bench_plot_incremental_pca.py - variable_batch_size_comparison()（120-160行）
def variable_batch_size_comparison(data):
# 10 个批次大小点：从 n_samples//10 到 n_samples
batch_sizes = [
i.astype(int) for i in np.linspace(data.shape[0] // 10, data.shape[0], num=10)
]
# 4 个成分数点：每个成分数都画一张图
for n_components in [
i.astype(int) for i in np.linspace(data.shape[1] // 10, data.shape[1], num=4)
]:
all_times = defaultdict(list)
all_errors = defaultdict(list)
pca = PCA(n_components=n_components)
rpca = PCA(
n_components=n_components, svd_solver="randomized", random_state=1999
)
# 先跑 PCA 和 randomized PCA 作为水平基线（与 batch_size 无关）
results_dict = {
k: benchmark(est, data) for k, est in [("pca", pca), ("rpca", rpca)]
}
# 重复填充基线值，让其与 batch_size 数量一致，呈现为水平线
all_times["pca"].extend([results_dict["pca"]["time"]] * len(batch_sizes))
all_errors["pca"].extend([results_dict["pca"]["error"]] * len(batch_sizes))
all_times["rpca"].extend([results_dict["rpca"]["time"]] * len(batch_sizes))
all_errors["rpca"].extend([results_dict["rpca"]["error"]] * len(batch_sizes))
# 遍历批次大小，每次仅跑 IncrementalPCA
for batch_size in batch_sizes:
ipca = IncrementalPCA(n_components=n_components, batch_size=batch_size)
results_dict = {k: benchmark(est, data) for k, est in [("ipca", ipca)]}
all_times["ipca"].append(results_dict["ipca"]["time"])
all_errors["ipca"].append(results_dict["ipca"]["error"])
plot_batch_times(all_times, n_components, batch_sizes, data)
plot_batch_errors(all_errors, n_components, batch_sizes, data)
这段代码实现了固定成分数、变化批次大小的对比，并使用"水平基线"设计：PCA 与 randomized PCA 的时间/误差与 batch_size 无关，因此它们的值被复制到与 batch_size 等长的列表中，从而在曲线图上呈现为水平直线。这让 IncrementalPCA 与两个不可变基线的对比一目了然。
82.7.5 LFW 数据加载与标准化
源码路径：benchmarks/bench_plot_incremental_pca.py - __main__（160-190行）
faces = fetch_lfw_people(resize=0.2, min_faces_per_person=5)
# 第 82 章 —— 限制数据集到 5000 人（不关心具体是谁！）
X = faces.data[:5000]
n_samples, h, w = faces.images.shape
n_features = X.shape[1]
# 第 82 章 —— 标准化关键步骤
X -= X.mean(axis=0)
X /= X.std(axis=0)
fixed_batch_size_comparison(X)
variable_batch_size_comparison(X)
plt.show()
这段代码加载 LFW 人脸数据集并执行显式标准化。X -= X.mean(axis=0); X /= X.std(axis=0) 是不可变基线的关键预处理步骤，因为 IncrementalPCA 对特征尺度敏感。虽然 PCA 内部也会做均值中心化，但显式标准化能消除零方差特征的除零警告。
82.8 设计中的取舍
基准测试的设计是一门在准确性、可重复性与工程效率之间不断权衡的艺术。本节集中讨论这些脚本在五个关键设计点上所做的取舍。
Q1：为什么 measure_one 使用 while elapsed < 0.5 循环并取 np.median 而不是单次 timeit？
基准测试的单次执行时间通常只有几毫秒甚至更短，但操作系统调度、CPU 频率变化、内存分配抖动会带来显著的单次计时噪声。通过累计执行至 0.5 秒以上，我们能收集到数十甚至上百次耗时样本，取中位数能有效抵抗长尾分布噪声，比平均值更稳健——平均值会被一两个异常高的耗时样本拉偏，而中位数天然忽略极端值。此外，perf_counter 提供单调递增的高精度计时，比 time.time 更适合短时测量；timeit 模块虽然内置循环逻辑，但它的固定迭代次数策略不能根据单次耗时动态调整，可能在小矩阵上跑出数十万次无意义的重复执行。
Q2：config_context(assume_finite=True) 在这里的作用是什么？会影响数值正确性吗？
它的作用是临时关闭 check_array 中的有限性检查（NaN/Inf 检测），跳过这部分开销能减少 5-15% 的耗时。它不会影响数值正确性，因为合成数据是确定无 NaN 的；对于真实数据，如果用户已知数据干净，开启这个上下文是安全的。但若数据可能包含 NaN/Inf，必须先做预处理清洗，否则后续数值计算会出现难以诊断的 NaN 传播。这种"明确已知前提 → 跳过防御性检查"的模式，是基准测试追求"测算法本身而非测检查逻辑"的常见工程实践。
Q3：为什么 bench_kernel_pca_solvers_time_vs_n_components.py 中 arpack 仅在 n_components < 100 时测试？
ARPACK 的核心优势是"只计算指定数量的极端特征值"，但当请求的特征值数量接近矩阵维度时，迭代次数爆炸，效率退化至与 dense 相当甚至更差。n_components < 100 是经验阈值，在这个范围内 ARPACK 仍能保持迭代收敛速度。当核矩阵维度 N=2000 时，请求 1999 个特征值意味着几乎要算全谱，此时 ARPACK 不仅没有加速，反而因为额外的迭代开销比 dense 更慢。这与 bench_kernel_pca_solvers_time_vs_n_samples.py 中 arpack 默认关闭的取舍同理——在大样本/大成分数时 ARPACK 都是"高不成低不就"的中间态。
Q4：为什么 bench_plot_randomized_svd.py 使用 ARPACK 的 svds(k=1) 计算谱范数，而不是直接用 sp.linalg.norm(A, ord=2)？
精确谱范数需要计算所有奇异值并取最大值，即全 SVD 复杂度 O(N³)，对大矩阵完全不可承受。ARPACK svds(k=1) 只迭代找最大奇异值，复杂度 O(N² × 迭代次数)，在大矩阵上快几个数量级。这是基准测试设计中"用近似测量替代精确测量以换取可行性"的典型权衡——我们测的是"相对误差"而非绝对谱范数值，只要比较的基线一致，近似谱范数足够敏感地捕捉到不同参数配置的精度差异。
Q5：为什么 bench_plot_randomized_svd.py 实现了分块 Frobenius 范数？为什么 gc.collect() 放在计时前？
重构误差矩阵 A - U·diag(s)·V 的大小与原矩阵相当。当原矩阵超过 4GB 时，重构矩阵无法装入 RAM，必须按行分块累加范数平方，最终开方得到完整范数。gen_batches 是 sklearn 提供的批次生成工具，自动处理边界对齐。MAX_MEMORY = int(4e9) 这个阈值是经验设定，对应典型服务器的物理内存与对象头部开销的余量。gc.collect() 放在计时前的理由是：让本次测量从"干净状态"开始，避免上一次迭代的内存分配残留干扰本次的内存分配耗时，使计时更能反映真实算法性能。计时后清理则没有这种效果，因为本次计时已经结束。这种"测试隔离"思想与单元测试中的 setup/teardown 一脉相承。
82.9 动手练习
-
阅读 PCA 求解器基准核心逻辑
-
阅读
benchmarks/bench_pca_solvers.py中的measure_one函数与主循环：-
为什么使用
while elapsed < 0.5循环并取np.median而不是单次timeit？ -
config_context(assume_finite=True)在这里的作用是什么？会影响数值正确性吗？ -
剪枝条件
log10(data.size) > 7对应多大规模的矩阵？为什么arpack和full要剪枝而randomized不用？ -
data_shapes生成逻辑如何保证覆盖"高宽"与"宽高"矩阵？REF_DIMS的三个基准维度有何含义？
-
-
-
剖析 KernelPCA 扩展性实验的正确性守门
-
对比
bench_kernel_pca_solvers_time_vs_n_samples.py与bench_kernel_pca_solvers_time_vs_n_components.py：-
为什么用
assert_array_almost_equal(np.abs(pred), np.abs(ref_pred))而不是直接比较pred与ref_pred？奇异向量符号不定性如何影响比较？ -
arpack求解器为何在n_components < 100时测试，大成分数时跳过？其背后的算法复杂度特征是什么？ -
两个脚本中
n_samples_range与n_compo_range的生成方式有何不同（线性 vs 指数）？这反映了什么实验设计意图？
-
-
-
实现随机 SVD 范数计算的分块策略
-
阅读
benchmarks/bench_plot_randomized_svd.py中的scalable_frobenius_norm_discrepancy与norm_diff：-
MAX_MEMORY = 4e9(4GB) 如何决定是否分块？分块时batch_size=1000的选择依据是什么？ -
谱范数计算为何调用
_init_arpack_v0生成v0再用svds(k=1)？相比sp.linalg.norm(A, ord=2)有何优势？ -
fbpca.pca与randomized_svd在n_iter与l(oversamples) 参数约定上的差异是什么？svd_timing如何统一封装？ -
实验 B 中构造
effective_rank=10/50/100的低秩矩阵，n_comp取rank/2, rank, 2*rank分别代表什么几何意义？
-
-
-
设计增量 PCA 基准实验的对照组
-
分析
benchmarks/bench_plot_incremental_pca.py的实验设计：-
variable_batch_size_comparison中为何引入PCA(svd_solver='randomized')作为额外基线？它与IncrementalPCA的本质区别是什么？ -
LFW 数据为何必须显式
X -= X.mean(); X /= X.std()？PCA内部是否也做了标准化？为什么IncrementalPCA更敏感？ -
benchmark函数中gc.collect()的位置为何放在计时前？如果放在计时后会有什么后果？ -
重构误差用 MAE (
np.mean(np.abs(data - data_r))) 而非 MSE，这在评估人脸重构质量时有何考量？
-
-
82.10 本章小结
这一章我们深入剖析了 scikit-learn 中 PCA 与 SVD 相关的独立基准脚本，从求解器全景对比到 KernelPCA 扩展性，从随机 SVD 参数权衡到增量 PCA 流式学习，构建了一套完整的矩阵分解效率评测体系。首先，我们学习了 PCA 多求解器（full、covariance_eigh、arpack、randomized、auto）在不同数据规模下的性能边界与 'auto' 启发式决策逻辑。其次，我们掌握了 KernelPCA 在样本量与成分数两个维度上的扩展性特征及近似求解器的适用场景。接着，我们深入理解了随机 SVD 中幂迭代次数与归一化策略对数值稳定性、精度与速度的三重权衡。然后，我们学习了增量 PCA 与批量 PCA 在流式数据场景下的时空权衡与重构误差控制。最后，我们学会了设计严谨的基准实验：中位数计时抵抗抖动、正确性守门、分块计算大矩阵范数、缓存与复现性保障。
为方便后续查阅与复习，下表将本章涉及的核心脚本、关键函数与对应的设计要点进行了系统汇总：
| 概念 | 解释 |
|------|------|
| PCA 求解器全景 (bench_pca_solvers.py) | five solvers (full, covariance_eigh, arpack, randomized, auto) 在 n_samples/n_features 网格上中位数计时，剪枝策略与双维度 log-log 可视化指导 auto 启发式 |
| SVD 3D 曲面对比 (bench_plot_svd.py) | SciPy 精确 SVD 与 randomized_svd (n_iter=0/3) 在低秩合成矩阵上的耗时地形图，直观展示随机算法加速优势 |
| KernelPCA 样本扩展性 (bench_kernel_pca_solvers_time_vs_n_samples.py) | 固定 n_components=100，n_samples 101→4000，dense O(N³) 爆炸 vs randomized 近似线性，arpack 中间态，assert 守门数值一致性 |
| KernelPCA 成分数扩展性 (bench_kernel_pca_solvers_time_vs_n_components.py) | 固定 n_train=2000，n_components 1→1999 指数网格，随机求解器小成分极快，arpack 仅 <100 测试，对数坐标揭示复杂度分水岭 |
| 随机 SVD 幂迭代深度探秘 (bench_plot_randomized_svd.py) | 三大实验：A-幂迭代/归一化策略权衡；B-矩阵真实秩与目标秩关系；C-实战数据集对比 fbpca。分块 Frobenius 范数、ARPACK 谱范数、LU/QR/auto 归一化策略 |
| 增量 PCA vs 批量 PCA (bench_plot_incremental_pca.py) | LFW 人脸数据，双维度对比：固定 batch_size 变 n_components、固定 n_components 变 batch_size。引入 randomized PCA 基线，MAE 重构误差量化近似损失，标准化预处理关键 |
下一章中，我们将学习聚类与最近邻专题的独立基准脚本，深入 KMeans 快速变体、层次聚类、最近邻搜索与 Ward 算法的性能特征，理解它们在数据规模与维度变化下的扩展性差异。
第 83 章 —— 独立基准脚本：聚类与最近邻专题 —— 丈量"空间划分算法"
83.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解基准测试脚本在 scikit-learn 性能评估体系中的定位与作用
-
掌握 KMeans 与 MiniBatchKMeans 在不同数据规模、维度与批次大小下的性能权衡分析方法
-
理解层次聚类四种链接策略的计算复杂度差异与可视化对比技巧
-
掌握最近邻搜索中 kd_tree、ball_tree 与 brute 三种算法在构建与查询阶段的性能分相分析
-
理解 Ward 算法在 scikit-learn 与 SciPy 两种实现下的跨库性能基准设计
-
学会设计多因子正交实验、对数网格扫描、分相计时等性能基准测试工程方法
-
掌握 3D 曲面图、分面折线图、堆积柱状图、比率热力图等性能数据可视化技术
83.2 生活类比
想象这四个基准脚本是算法性能的"奥林匹克测试场"，每一行代码都是精心设计的比赛项目，每个可视化图表都是裁判手中的秒表和计时工具。具体来看，KMeans 对决 MiniBatchKMeans 就像全程马拉松与间歇训练的较量：全批次 KMeans 像马拉松运动员一次性跑完全程，每一步都经过深思熟虑，它要求对全量数据反复计算质心与距离，收敛稳定但计算成本高昂；MiniBatchKMeans 像间歇跑者分段冲刺，用小批次数据快速迭代，每一次"冲刺"只利用一小部分样本更新质心位置，以近似路径逼近最优解。批次大小就是"冲刺距离"，太短会导致频繁的批次切换开销（累死人），太长则退化为全批次收敛慢（变成马拉松）。继续来看层次聚类的较量——同一座山（数据），四种攀岩路线策略耗时天差地别：Single（单链接）像"贪心独木桥"，只找最近的两个点，优点是速度快，缺点是容易形成链状效应；Complete（全链接）像"稳扎稳打"，确保合并的两簇中任意两点都最远，安全但计算量大；Average（平均链接）像"民主协商"，用所有点对的平均距离做决策；Ward 像"结构工程师"，严格最小化合并带来的方差增加，每一步都需评估簇内方差的变化。再看最近邻搜索的三角对决——三种找人策略各有适用场景：kd_tree 像"按地址分区查户口"，先把城市划分成网格区域，查询时只需在相关区域搜索；ball_tree 像"画圈圈找人"，用超球体包围数据点，快速排除圈外区域；brute 像"挨家挨户敲门"，最简单直接，但样本量大时效率低下。索引构建阶段是"画地图"，查询阶段是"按图索骥"，不同地形（维度/样本量）适合不同策略。最后是Ward 双实现对比——同款引擎两套调校：scikit-learn 版本像"多核涡轮增压"（Cython 并行加速），SciPy 版本像"经典自然吸气"（单线程 C 实现）。两者实现的是同一套 Ward 最小方差合并的数学原理，工程实现差异造就了性能地图上的"蓝红分界线"。小规模数据时 Python 调用开销占比高，全批次 C 实现反而更快；大规模数据时 Cython 并行化优势显现，sklearn 弯道超车；中等地形下两库势均力敌，存在清晰的"等高线"决策边界。
83.3 源码地图
benchmarks/bench_plot_fastkmeans.py
├── compute_bench() # 核心基准：KMeans vs MiniBatchKMeans 二维网格扫描 (n_samples, n_features)
├── compute_bench_2() # 专项基准：MiniBatchKMeans 批次大小敏感性分析
├── __main__ # 可视化入口：3D 曲面图 + 2D 折线图
benchmarks/bench_plot_hierarchical.py
├── compute_bench() # 核心基准：AgglomerativeClustering 四链接横评 二维网格扫描
├── __main__ # 可视化入口：分面折线图 (2x2 子图矩阵)
benchmarks/bench_plot_neighbors.py
├── get_data() # 数据工厂：生成 dense 随机数据或 digits 真实数据
├── barplot_neighbors() # 核心基准：三因子正交实验 (N, D, k) 分相计时 (fit/kneighbors)
├── __main__ # 可视化入口：堆积柱状图 (构建/查询) + 对数坐标 + 双数据集重复
benchmarks/bench_plot_ward.py
├── __main__ # 核心基准+可视化：对数网格扫描 + 比率热力图 + 等高线决策边界
83.4 KMeans 对决 MiniBatchKMeans —— 批处理带来的"速度红利"
聚类算法的性能基准测试是理解算法工程特性的窗口。KMeans 作为最经典的聚类算法，其全批次 Lloyd 实现与 MiniBatchKMeans 的小批次随机优化形成鲜明对比。这个基准脚本通过二维网格扫描，揭示数据规模与特征维度对两种算法的联合影响。在样本与特征的双重维度上，全批次算法与小批次算法的性能边界并非线性变化——当样本数从 50 跃升到 150、特征数从 150 跃升到 50000 时，KMeans 的每次迭代成本呈双线性增长，而 MiniBatchKMeans 的增量更新机制则可能在某些区域反超。理解这一曲面拓扑，是为新数据选择合适聚类策略的关键。
83.4.1 二维网格扫描：样本与特征的联合影响
compute_bench() 函数通过双重 for 循环遍历 n_samples × n_features 的笛卡尔积，对每个数据规模组合同时运行 KMeans 和 MiniBatchKMeans，记录两者的训练耗时与惯性指标（聚类质量），最终返回一个嵌套字典包含四个指标序列。
# 第 83 章 —— benchmarks/bench_plot_fastkmeans.py - compute_bench() (第13-56行)
def compute_bench(samples_range, features_range):
it = 0 # 迭代计数器
results = defaultdict(lambda: []) # 自动初始化列表的结果容器
chunk = 100 # MiniBatchKMeans 批次大小固定为 100
max_it = len(samples_range) * len(features_range) # 总迭代次数
# 嵌套循环遍历样本数与特征数的所有组合
for n_samples in samples_range:
for n_features in features_range:
it += 1
# 打印进度条，便于长时间运行时监控
print("==============================")
print("Iteration %03d of %03d" % (it, max_it))
print("==============================")
print()
# 生成整数特征矩阵，范围 [-50, 51)，模拟通用稠密数据
data = nr.randint(-50, 51, (n_samples, n_features))
# 基准测试 1：标准 KMeans
print("K-Means")
tstart = time()
kmeans = KMeans(init="k-means++", n_clusters=10).fit(data) # k-means++ 初始化，n_clusters=10 固定
delta = time() - tstart
print("Speed: %0.3fs" % delta)
print("Inertia: %0.5f" % kmeans.inertia_)
print()
# 记录 KMeans 的速度与质量指标
results["kmeans_speed"].append(delta)
results["kmeans_quality"].append(kmeans.inertia_)
# 基准测试 2：MiniBatchKMeans
print("Fast K-Means")
# let's prepare the data in small chunks
mbkmeans = MiniBatchKMeans(
init="k-means++", n_clusters=10, batch_size=chunk # batch_size=chunk 控制变量
)
tstart = time()
mbkmeans.fit(data)
delta = time() - tstart
print("Speed: %0.3fs" % delta)
print("Inertia: %f" % mbkmeans.inertia_)
print()
print()
# 记录 MiniBatchKMeans 的速度与质量指标
results["MiniBatchKMeans Speed"].append(delta)
results["MiniBatchKMeans Quality"].append(mbkmeans.inertia_)
return results
逐行解析：
-
第 14 行：
defaultdict(lambda: [])创建自动初始化列表的结果容器，避免手动检查键存在性的繁琐。 -
第 15 行：固定
chunk = 100作为 MiniBatchKMeans 的批次大小，这是实验的控制变量。 -
第 20-21 行：嵌套循环覆盖
n_samples × n_features所有组合，构建完整的二维性能拓扑图。 -
第 28 行：
nr.randint(-50, 51, (n_samples, n_features))生成整型随机数据，避免浮点精度差异干扰结果。 -
第 32 行：
KMeans(init="k-means++", n_clusters=10)使用 k-means++ 初始化，n_clusters=10 固定，消除初始化差异干扰。 -
第 37-38 行：记录两种算法的
speed（训练耗时）和quality（惯性/聚类质量）指标。
这段代码展示了控制变量法的精髓：固定算法参数（init、n_clusters）、固定批次大小（chunk=100），仅变化数据规模，观察性能曲面的形态。
为什么用整数矩阵？ 整型数据确保距离计算的可复现性，避免浮点精度在不同平台间的微小差异。同时，对称分布（-50 到 51）使 Ward 的方差最小化准则有意义的输入范围。
83.4.2 批次大小敏感性分析
compute_bench_2() 函数聚焦单一变量，固定极高维度（50000 特征）和 8 个高斯簇的合成数据，仅变化 MiniBatchKMeans 的 batch_size 参数，量化批次大小对收敛速度与聚类质量的非线性影响。
# 第 83 章 —— benchmarks/bench_plot_fastkmeans.py - compute_bench_2() (第58-94行)
def compute_bench_2(chunks):
results = defaultdict(lambda: []) # 自动初始化列表的结果容器
n_features = 50000 # 固定极高维度，隔离批次大小的影响
# 构造 8 个高斯簇的合成数据，模拟结构化分布
means = np.array(
[
[1, 1], # 第 1 个簇中心
[-1, -1], # 第 2 个簇中心
[1, -1], # 第 3 个簇中心
[-1, 1], # 第 4 个簇中心
[0.5, 0.5], # 第 5 个簇中心
[0.75, -0.5], # 第 6 个簇中心
[-1, 0.75], # 第 7 个簇中心
[1, 0], # 第 8 个簇中心
]
)
X = np.empty((0, 2)) # 初始化空数据矩阵，准备逐簇拼接
for i in range(8):
# 每个簇生成 n_features 个点，加 0.8*randn 引入高斯噪声
X = np.r_[X, means[i] + 0.8 * np.random.randn(n_features, 2)]
max_it = len(chunks) # 单维扫描，总迭代次数等于批次数
it = 0
# 单维扫描：仅变化批次大小
for chunk in chunks:
it += 1
print("==============================")
print("Iteration %03d of %03d" % (it, max_it))
print("==============================")
print()
print("Fast K-Means")
tstart = time()
# batch_size=chunk 是实验自变量，n_clusters=8 匹配真实簇数
mbkmeans = MiniBatchKMeans(init="k-means++", n_clusters=8, batch_size=chunk)
mbkmeans.fit(X) # 执行小批次 KMeans 拟合
delta = time() - tstart # 计算训练耗时
print("Speed: %0.3fs" % delta)
print("Inertia: %0.3fs" % mbkmeans.inertia_)
print()
results["MiniBatchKMeans Speed"].append(delta)
results["MiniBatchKMeans Quality"].append(mbkmeans.inertia_)
return results
逐行解析：
-
第 60 行：固定
n_features = 50000，创造一个极高维度的数据场景，确保批次切换开销成为瓶颈。 -
第 61-66 行：构造 8 个高斯簇的合成数据，每个簇有 n_features 个样本，总计约 40 万点。
-
第 67-68 行：数据矩阵沿行方向拼接（
np.r_），最终形状为(8*n_features, 2)。 -
第 79 行：
batch_size=chunk是实验的自变量，从 500 到 10000 变化。 -
第 81 行：
n_clusters=8匹配数据中的真实簇数，使聚类任务有意义。
这个实验揭示了批次大小与收敛速度的非线性关系：过小的批次增加迭代次数和切换开销，过大的批次退化为全批次失去小批次的优势。
83.4.3 3D 曲面图可视化
__main__ 入口运行两组基准实验，使用 3D 曲面图展示 Speed 与 Quality 随样本/特征规模的联合变化，下半部分用 2D 折线图展示批次大小对收敛的单调影响。
# 第 83 章 —— benchmarks/bench_plot_fastkmeans.py - __main__ (第96-146行)
if __name__ == "__main__":
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d # register the 3d projection # noqa: F401
# 定义参数扫描范围
samples_range = np.linspace(50, 150, 5).astype(int) # 5 个样本数点，50~150
features_range = np.linspace(150, 50000, 5).astype(int) # 5 个特征数点，150~50000
chunks = np.linspace(500, 10000, 15).astype(int) # 15 个批次大小点
# 运行两组基准实验
results = compute_bench(samples_range, features_range)
results_2 = compute_bench_2(chunks)
# 计算全局上限，确保不同算法曲面可比
max_time = max(
[max(i) for i in [t for (label, t) in results.items() if "speed" in label]]
)
max_inertia = max(
[max(i) for i in [t for (label, t) in results.items() if "speed" not in label]]
)
# 创建 2x2 子图布局
fig = plt.figure("scikit-learn K-Means benchmark results")
for c, (label, timings) in zip("brcy", sorted(results.items())):
if "speed" in label:
# 左上：Speed 的 3D 曲面
ax = fig.add_subplot(2, 2, 1, projection="3d")
ax.set_zlim3d(0.0, max_time * 1.1) # 统一 Z 轴上限
else:
# 右上：Quality 的 3D 曲面
ax = fig.add_subplot(2, 2, 2, projection="3d")
ax.set_zlim3d(0.0, max_inertia * 1.1)
# 构建网格用于 3D 曲面绑图
X, Y = np.meshgrid(samples_range, features_range)
Z = np.asarray(timings).reshape(
samples_range.shape[0], features_range.shape[0]
)
ax.plot_surface(X, Y, Z.T, cstride=1, rstride=1, color=c, alpha=0.5)
ax.set_xlabel("n_samples")
ax.set_ylabel("n_features")
# 下半部分：批次大小的 2D 折线图
i = 0
for c, (label, timings) in zip("br", sorted(results_2.items())):
i += 1
ax = fig.add_subplot(2, 2, i + 2)
y = np.asarray(timings)
ax.plot(chunks, y, color=c, alpha=0.8)
ax.set_xlabel("Chunks")
ax.set_ylabel(label)
plt.show()
逐行解析：
-
第 107-109 行：定义三组实验参数，samples_range 和 features_range 用作 2D 网格扫描，chunks 用于单因素分析。
-
第 115-119 行：计算所有速度/质量指标的最大值，作为 3D 曲面图的统一 Z 轴上限，确保不同算法可比。
-
第 121-135 行：遍历 results 字典，根据标签区分 speed 和 quality，创建对应的 3D 子图。
-
第 132-135 行：
reshape将一维计时数组转换为与网格匹配的 2D 数组，Z.T转置确保 X/Y 轴对应正确。 -
第 137-146 行：为 results_2（批次大小实验）创建 2D 折线图，补充批次敏感性的单调趋势展示。
3D 曲面图的核心设计：
这个基准的核心价值在于揭示两个维度的权衡：样本数增长主要影响迭代成本，特征数增长主要影响每次距离计算的成本。3D 曲面图将这种联合影响可视化，帮助你在面对新数据时快速判断哪种算法更合适。
83.5 层次聚类基准 —— 四种链接策略的"攀爬竞赛"
层次聚类的四种链接策略（single、average、complete、ward）代表了不同的合并哲学。这个基准脚本通过固定算法框架、仅切换链接参数的设计，精确隔离链接策略对计算复杂度的影响。在 scikit-learn 中，四种链接共用 AgglomerativeClustering 的同一套底层框架，仅在合并准则层面存在差异，这使得横向对比成为可能。
83.5.1 单算法多链接横评
compute_bench() 函数在固定 n_clusters=10 的前提下，对四种链接策略进行全组合计时，最终返回以链接名称为键、计时列表为值的字典。
# 第 83 章 —— benchmarks/bench_plot_hierarchical.py - compute_bench() (第10-44行)
def compute_bench(samples_range, features_range):
it = 0
results = defaultdict(lambda: [])
max_it = len(samples_range) * len(features_range)
for n_samples in samples_range:
for n_features in features_range:
it += 1
print("==============================")
print("Iteration %03d of %03d" % (it, max_it))
print("n_samples %05d; n_features %02d" % (n_samples, n_features))
print("==============================")
print()
# 生成整数特征矩阵，范围 [-50, 51)
data = nr.randint(-50, 51, (n_samples, n_features))
# 遍历四种链接策略
for linkage in ("single", "average", "complete", "ward"):
print(linkage.capitalize())
tstart = time()
# 固定 n_clusters=10，仅切换 linkage 参数隔离策略差异
AgglomerativeClustering(linkage=linkage, n_clusters=10).fit(data)
delta = time() - tstart
print("Speed: %0.3fs" % delta)
print()
results[linkage].append(delta)
return results
逐行解析：
-
第 29 行：
nr.randint(-50, 51, (n_samples, n_features))生成整型稠密矩阵。 -
第 31-39 行：核心设计——固定
AgglomerativeClustering(n_clusters=10)，仅在第 31 行的 for 循环中切换linkage参数。 -
第 36 行：
time.time()获取壁钟时间（wall-clock time），包含 Python/Cython 调用开销和算法计算时间。 -
第 38 行：
results[linkage].append(delta)按链接策略累积计时结果。
为什么用整数矩阵？ 整型数据确保距离计算的可复现性，避免浮点精度在不同平台间的微小差异。同时，对称分布（-50 到 51）使 Ward 的方差最小化准则有意义的输入范围。
83.5.2 分面折线图可视化
__main__ 入口创建 2x2 分面子图矩阵，每个子图对应一个 n_features 值，子图内绘制四种链接策略随样本数变化的计时曲线。
# 第 83 章 —— benchmarks/bench_plot_hierarchical.py - __main__ (第46-82行)
if __name__ == "__main__":
import matplotlib.pyplot as plt
# 定义扫描范围
samples_range = np.linspace(1000, 15000, 8).astype(int) # 1000~15000 样本
features_range = np.array([2, 10, 20, 50]) # 固定 4 个维度值作为分面
results = compute_bench(samples_range, features_range)
max_time = max([max(i) for i in [t for (label, t) in results.items()]])
# 从 tab10 调色板取前 4 色，确保打印和色盲友好
colors = plt.get_cmap("tab10")(np.linspace(0, 1, 10))[:4]
lines = {linkage: None for linkage in results.keys()}
# 2x2 分面子图矩阵
fig, axs = plt.subplots(2, 2, sharex=True, sharey=True)
fig.suptitle("Scikit-learn agglomerative clustering benchmark results", fontsize=16)
for c, (label, timings) in zip(colors, sorted(results.items())):
# 将一维计时重塑为 (n_samples, n_features) 的 2D 矩阵
timing_by_samples = np.asarray(timings).reshape(
samples_range.shape[0], features_range.shape[0]
)
# 在每个子图中绘制四条链接曲线
for n in range(timing_by_samples.shape[1]):
ax = axs.flatten()[n] # 展平为一维索引
(lines[label],) = ax.plot(
samples_range, timing_by_samples[:, n], color=c, label=label
)
ax.set_title("n_features = %d" % features_range[n])
if n >= 2:
ax.set_xlabel("n_samples")
if n % 2 == 0:
ax.set_ylabel("time (s)")
# 右侧统一图例，避免子图内图例遮挡数据
fig.subplots_adjust(right=0.8)
fig.legend(
[lines[link] for link in sorted(results.keys())],
sorted(results.keys()),
loc="center right",
fontsize=8,
)
plt.show()
逐行解析：
-
第 52 行：
samples_range = np.linspace(1000, 15000, 8)生成 8 个样本量点（1000~15000）。 -
第 53 行：
features_range = np.array([2, 10, 20, 50])固定 4 个维度值，作为分面的基础。 -
第 59 行：
plt.subplots(2, 2, sharex=True, sharey=True)创建 2x2 分面，共享坐标轴确保视觉对比准确。 -
第 65-67 行：重塑计时数组，每列对应一个
n_features值，每行对应一个n_samples值。 -
第 69-77 行：在每个子图中绘制该维度下四种链接的计时曲线。
-
第 80-81 行：右侧统一图例，
loc="center right"避免子图内图例与数据重叠。
分面图的核心设计：
设计意图解读：分面维度（n_features）展示了维度增长如何改变链接策略的相对优劣——例如 Ward 在高维下可能退化更快，因为需要计算更多维度的方差。
83.6 最近邻搜索算法 —— kd_tree、ball_tree 与 brute 的"三角对决"
最近邻搜索是机器学习的基石算法，其性能直接决定了 KNN 分类、异常检测等下游任务效率。这个基准通过三因子正交实验设计，精确隔离样本数 N、维度 D、邻居数 k 对三种搜索算法的独立影响。
83.6.1 数据工厂：合成与真实数据双源
get_data() 函数统一封装两类数据源：dense 模式生成均匀随机矩阵用于最坏/平均情况性能测试，digits 模式加载手写数字数据集用于流形结构数据下的加速效果验证。
# 第 83 章 —— benchmarks/bench_plot_neighbors.py - get_data() (第13-28行)
def get_data(N, D, dataset="dense"):
"""生成用于基准测试的数据
Parameters
----------
N : int
样本数量
D : int
特征维度
dataset : str
数据类型，"dense" 生成随机数据，"digits" 使用真实手写数字
Returns
-------
X : ndarray
形状为 (N, D) 的数据矩阵
"""
if dataset == "dense":
np.random.seed(0) # 固定随机种子，确保可复现
return np.random.random((N, D)) # 生成稠密随机数据
elif dataset == "digits":
X, _ = datasets.load_digits(return_X_y=True) # 加载手写数字数据集
i = np.argsort(X[0])[::-1] # 按第一行方差降序排列特征
X = X[:, i]
return X[:N, :D] # 截取前 N 个样本和前 D 个特征
else:
raise ValueError("invalid dataset: %s" % dataset) # 非法数据集类型
逐行解析：
-
第 19 行：
np.random.seed(0)固定随机种子，确保 dense 模式下结果完全可复现。 -
第 23 行：
datasets.load_digits()加载 scikit-learn 内置的手写数字数据集（1797 样本，64 特征）。 -
第 24 行：
np.argsort(X[0])[::-1]按第一行（第一个样本）的方差降序排列特征，确保前 D 列包含最重要特征。 -
第 25 行：截取
[:N, :D]，允许 D > 64 时优雅截断（load_digits只返回 64 特征）。
双数据源的设计意图：dense 模式测试最坏/平均情况性能，digits 模式测试流形结构数据下的加速效果。
83.6.2 三因子正交实验：分相计时的核心逻辑
barplot_neighbors() 函数通过三个独立的单因子扫描实验，分别变化 N、D、k，每次实验内对三种算法进行分相计时（fit 与 kneighbors 分开），最终返回六组结果字典（每个因子 × 构建/查询）。
# 第 83 章 —— benchmarks/bench_plot_neighbors.py - barplot_neighbors() 实验部分
def barplot_neighbors(
Nrange=2 ** np.arange(1, 11), # 2, 4, 8, ..., 1024
Drange=2 ** np.arange(7), # 1, 2, 4, 8, 16, 32, 64
krange=2 ** np.arange(10), # 2, 4, 8, ..., 512
N=1000, D=64, k=5, # 基准值（控制变量）
leaf_size=30, # 树算法的叶子大小参数
dataset="digits", # 数据集类型
):
algorithms = ("kd_tree", "brute", "ball_tree") # 三种搜索算法
fiducial_values = {"N": N, "D": D, "k": k} # 基准值字典
# ------------------------------------------------------------
# 实验 1：变化 N，固定 D=64, k=5
N_results_build = {alg: np.zeros(len(Nrange)) for alg in algorithms} # 构建时间容器
N_results_query = {alg: np.zeros(len(Nrange)) for alg in algorithms} # 查询时间容器
for i, NN in enumerate(Nrange):
print("N = %i (%i out of %i)" % (NN, i + 1, len(Nrange)))
X = get_data(NN, D, dataset) # 生成对应样本量的数据
for algorithm in algorithms:
# 创建搜索器：n_neighbors 取 NN 与 k 的最小值防越界
nbrs = neighbors.NearestNeighbors(
n_neighbors=min(NN, k), algorithm=algorithm, leaf_size=leaf_size
)
# 分相计时：fit 与 kneighbors 分开
t0 = time()
nbrs.fit(X) # 阶段 1：索引构建
t1 = time()
nbrs.kneighbors(X) # 阶段 2：查询执行
t2 = time()
N_results_build[algorithm][i] = t1 - t0 # 索引构建耗时
N_results_query[algorithm][i] = t2 - t1 # 查询耗时
# ------------------------------------------------------------
# 实验 2：变化 D，固定 N=1000, k=5
D_results_build = {alg: np.zeros(len(Drange)) for alg in algorithms}
D_results_query = {alg: np.zeros(len(Drange)) for alg in algorithms}
for i, DD in enumerate(Drange):
print("D = %i (%i out of %i)" % (DD, i + 1, len(Drange)))
X = get_data(N, DD, dataset) # 生成对应维度的数据
for algorithm in algorithms:
nbrs = neighbors.NearestNeighbors(
n_neighbors=k, algorithm=algorithm, leaf_size=leaf_size
)
t0 = time()
nbrs.fit(X)
t1 = time()
nbrs.kneighbors(X)
t2 = time()
D_results_build[algorithm][i] = t1 - t0
D_results_query[algorithm][i] = t2 - t1
# ------------------------------------------------------------
# 实验 3：变化 k，固定 N=1000, D=64
k_results_build = {alg: np.zeros(len(krange)) for alg in algorithms}
k_results_query = {alg: np.zeros(len(krange)) for alg in algorithms}
X = get_data(N, DD, dataset) # 复用实验 2 的最终数据
for i, kk in enumerate(krange):
print("k = %i (%i out of %i)" % (kk, i + 1, len(krange)))
for algorithm in algorithms:
nbrs = neighbors.NearestNeighbors(
n_neighbors=kk, algorithm=algorithm, leaf_size=leaf_size
)
t0 = time()
nbrs.fit(X)
t1 = time()
nbrs.kneighbors(X)
t2 = time()
k_results_build[algorithm][i] = t1 - t0
k_results_query[algorithm][i] = t2 - t1
逐行解析：
-
第 32-38 行：定义三个维度的变化范围，使用
2 ** np.arange()生成指数级序列（对数均匀采样）。 -
第 39-41 行：定义基准值 N=1000, D=64, k=5，这是控制变量实验中的固定参数。
-
第 47-48 行：
neighbors.NearestNeighbors创建搜索器对象，指定算法和叶子大小。 -
第 50-58 行：分相计时——
fit()构建索引的时间与kneighbors()执行查询的时间分开记录，这是理解算法架构的关键。 -
第 60 行：
n_neighbors=min(NN, k)防止请求的邻居数超过样本数。
三因子正交实验设计：
83.6.3 堆积柱状图可视化：构建与查询的分相展示
barplot_neighbors() 的可视化部分将六组计时结果绘制为三个 3 行 1 列的子图，每个子图对应一个变化因子，子图内以堆积柱状图展示三种算法在每个因子取值下的构建（红）与查询（蓝）耗时。
# 第 83 章 —— benchmarks/bench_plot_neighbors.py - barplot_neighbors() 可视化部分
plt.figure(figsize=(8, 11)) # 创建 8x11 英寸画布容纳三个子图
# 遍历三组实验数据：每个元组包含子图位置、变化因子、构建结果、查询结果
for sbplt, vals, quantity, build_time, query_time in [
(311, Nrange, "N", N_results_build, N_results_query),
(312, Drange, "D", D_results_build, D_results_query),
(313, krange, "k", k_results_build, k_results_query),
]:
ax = plt.subplot(sbplt, yscale="log") # 创建子图，Y 轴对数刻度
plt.grid(True) # 显示网格线辅助读数
tick_vals = [] # X 轴刻度值收集器
tick_labels = [] # X 轴刻度标签收集器
# 计算动态基线，避免柱从零开始导致 log(0) 问题
bottom = 10 ** np.min(
[min(np.floor(np.log10(build_time[alg]))) for alg in algorithms]
)
for i, alg in enumerate(algorithms):
# 为每个算法的柱子分配 X 偏移位置，三组柱子并列
xvals = 0.1 + i * (1 + len(vals)) + np.arange(len(vals))
width = 0.8 # 柱宽
# 红柱：构建时间（从 bottom 基线开始）
c_bar = plt.bar(xvals, build_time[alg] - bottom, width, bottom, color="r")
# 蓝柱：查询时间（从构建时间顶部开始堆积）
q_bar = plt.bar(xvals, query_time[alg], width, build_time[alg], color="b")
# 收集 X 轴刻度
tick_vals += list(xvals + 0.5 * width)
tick_labels += ["%i" % val for val in vals]
# 子图内左上角标注算法名称，半透明白色背景框
plt.text(
(i + 0.02) / len(algorithms),
0.98,
alg,
transform=ax.transAxes,
ha="left",
va="top",
bbox=dict(facecolor="w", edgecolor="w", alpha=0.5),
)
plt.ylabel("Time (s)") # Y 轴标签
# X 轴固定刻度为实验值
ax.xaxis.set_major_locator(ticker.FixedLocator(tick_vals))
ax.xaxis.set_major_formatter(ticker.FixedFormatter(tick_labels))
# 旋转 X 轴标签防止重叠
for label in ax.get_xticklabels():
label.set_rotation(-90)
label.set_fontsize(10)
# 构建右侧标题："Varying N/D/k"
title_string = "Varying %s" % quantity
# 构建右侧描述：固定参数列表，如 "D = 64, k = 5"
descr_string = ""
for s in "NDk":
if s == quantity:
pass # 跳过变化因子本身
else:
descr_string += "%s = %i, " % (s, fiducial_values[s])
descr_string = descr_string[:-2] # 去掉末尾逗号
# 右侧旋转 -90° 绘制标题
plt.text(
1.01, 0.5, title_string,
transform=ax.transAxes, rotation=-90,
ha="left", va="center", fontsize=20,
)
# 右侧旋转 -90° 绘制固定参数描述
plt.text(
0.99, 0.5, descr_string,
transform=ax.transAxes, rotation=-90,
ha="right", va="center",
)
# 全局总标题标识数据集类型
plt.gcf().suptitle("%s data set" % dataset.capitalize(), fontsize=16)
# 全局图例标识红蓝柱含义
plt.figlegend((c_bar, q_bar), ("construction", "N-point query"), "upper right")
if __name__ == "__main__":
barplot_neighbors(dataset="digits") # digits 数据集基准
barplot_neighbors(dataset="dense") # dense 数据集基准
plt.show()
逐行解析：
-
第 139 行：
yscale="log"使用对数坐标，容纳从毫秒到秒的跨数量级差异。 -
第 144-146 行：计算动态基线
bottom = 10 ** np.min(...)，确保红柱（构建时间）从非零位置开始，避免log(0)导致的问题。 -
第 148-152 行：堆积柱状图核心——红柱的
bottom参数设为bottom，蓝柱的bottom参数设为build_time[alg]，实现红下蓝上的堆积效果。 -
第 154-156 行：收集 X 轴刻度值与标签。
-
第 158-167 行：在子图内标注算法名称，使用半透明白色背景框避免遮挡数据。
-
第 171-173 行：设置 X 轴为固定刻度（对应实验值）。
-
第 176-178 行：旋转 90° 防止标签重叠。
-
第 181 行：构建"变化 X"的标题字符串。
-
第 184-191 行：构建"固定参数"的描述字符串（如 "D = 64, k = 5"）。
-
第 194-211 行：在右侧旋转 -90° 绘制标题与描述字符串。
-
第 214 行：设置全局总标题。
-
第 217 行：全局图例标识红蓝柱含义。
-
第 220-224 行：
__main__入口对 digits 与 dense 两个数据集分别运行基准并显示图表。
堆积柱状图的设计意图：红色代表索引构建开销，蓝色代表查询开销。总高度 = 端到端延迟。通过颜色对比，可以直观判断算法属于"建图贵查路快"（树算法）还是"建图免费查路贵"（暴力算法）。
83.7 Ward 算法双实现对比 —— sklearn 与 SciPy 的"同源较量"
当同一算法存在多种实现时，跨库基准是评估工程优化差异的最佳方式。Ward 最小方差合并准则在 scikit-learn 和 SciPy 中都有实现，两者数学原理完全相同，差异仅在工程实现层面。
83.7.1 对数网格扫描：指数级规模覆盖
bench_plot_ward.py 的 __main__ 部分在 7×9 的对数网格上对两种实现进行计时，计算速度比 ratio = scikits_time / scipy_time，构建完整的"性能地形图"。
# 第 83 章 —— benchmarks/bench_plot_ward.py - __main__ (第13-51行)
"""
Benchmark scikit-learn's Ward implement compared to SciPy's
"""
import time
import matplotlib.pyplot as plt
import numpy as np
from scipy.cluster import hierarchy
from sklearn.cluster import AgglomerativeClustering
ward = AgglomerativeClustering(n_clusters=3, linkage="ward") # sklearn Ward 实例
# 第 83 章 —— 对数均匀采样，覆盖 3~1000 样本和 10~3000 特征
n_samples = np.logspace(0.5, 3, 9) # [3, 10, 17, 31, 56, 100, 177, 316, 562]
n_features = np.logspace(1, 3.5, 7) # 7 个值，范围 10 ~ 3162
N_samples, N_features = np.meshgrid(n_samples, n_features) # 构造 7×9 网格
scikits_time = np.zeros(N_samples.shape) # sklearn 计时矩阵
scipy_time = np.zeros(N_samples.shape) # scipy 计时矩阵
# 第 83 章 —— 嵌套循环：对数网格的每个 (n, p) 组合
for i, n in enumerate(n_samples):
for j, p in enumerate(n_features):
# 生成 i.i.d. 标准正态数据
X = np.random.normal(size=(n, p))
# scikit-learn 计时
t0 = time.time()
ward.fit(X)
scikits_time[j, i] = time.time() - t0
# SciPy 计时
t0 = time.time()
hierarchy.ward(X)
scipy_time[j, i] = time.time() - t0
ratio = scikits_time / scipy_time # 计算速度比
浙公网安备 33010602011771号