Sklearn-源码解析-书-v1-0-三十六-
Sklearn 源码解析（书）v1.0（三十六）
逐行解析：
-
第 26-27 行：
np.logspace(0.5, 3, 9)生成 9 个在对数空间均匀分布的样本数点（从 √10 ≈ 3.16 到 10³=1000）。 -
第 28 行：
np.logspace(1, 3.5, 7)生成 7 个特征数点，在对数空间从 10 均匀扩展到约 3162。 -
第 29 行：
meshgrid创建 7×9 的参数网格，共 63 个实验点。 -
第 36 行：
np.random.normal(size=(n, p))生成 i.i.d. 标准正态数据，球对称分布使 Ward 合并顺序具有理论确定性。 -
第 39-42 行：分别对 sklearn 和 scipy 的 Ward 实现计时，记录到对应数组。
-
第 44 行：
ratio = scikits_time / scipy_time计算速度比，<1 表示 sklearn 更快。
对数网格的设计意图：对数空间均匀采样确保在各个数量级上都有实验点，避免线性采样时小规模被大规模"淹没"。63 个点的网格构建出完整的"性能地形图"。
83.7.2 比率热力图 + 等高线：跨库决策边界可视化
# 第 83 章 —— benchmarks/bench_plot_ward.py - 可视化 (第53-67行)
plt.figure("scikit-learn Ward's method benchmark results")
# 第 83 章 —— 对数色映射：对称展示加速/减速
plt.imshow(np.log(ratio), aspect="auto", origin="lower") # 对比率取对数后用 imshow 显示
plt.colorbar()
# 第 83 章 —— 绘制 ratio=1 的等值线（两库性能相等线）
plt.contour(
ratio,
levels=[1],
colors="k",
)
plt.yticks(range(len(n_features)), n_features.astype(int))
plt.ylabel("N features")
plt.xticks(range(len(n_samples)), n_samples.astype(int))
plt.xlabel("N samples")
plt.title("Scikit's time, in units of scipy time (log)")
plt.show()
逐行解析：
-
第 55 行：
plt.imshow(np.log(ratio), ...)对比率取对数后用 imshow 显示，对数变换使加速和减速在视觉上对称。 -
第 57-61 行：
plt.contour(ratio, levels=[1], colors="k")绘制ratio=1的等值线（黑色），精确划分 sklearn 胜（浅色）vs scipy 胜（深色）区域。 -
第 62-66 行：坐标轴标注原始样本数/特征数值（而非索引），保持工程直觉。
热力图 + 等高线的设计意图：
为了直观展示比率值与性能优势之间的关系，下表列出了不同区域比率所代表的含义：
| 区域 | ratio | log(ratio) | 颜色 | 含义 |
|------|-------|------------|------|------|
| 右上角 | < 1 | < 0 | 浅色 | sklearn 更快（并行化优势） |
| 左下角 | > 1 | > 0 | 深色 | scipy 更快（调用开销占比小） |
| 等值线 | = 1 | = 0 | 黑色 | 两库性能相等 |
结果解读规律：大规模高维数据下 sklearn 的 Cython 并行化优势明显，小规模低维数据下 scipy 的单线程 C 实现更高效（无 Python/Cython 调用开销）。
83.8 设计中的取舍
为什么使用壁钟时间而非 CPU 时间？ scikit-learn 大量使用 Cython 和底层 C/C++ 库，CPU 时间可能无法准确反映多线程并行化的实际收益。壁钟时间是最终用户体验的核心指标，直接回答"需要等待多久"。当 Cython 代码释放 GIL 后并行执行时，多核 CPU 时间可能远大于单线程壁钟时间，此时壁钟时间更能反映实际加速比。在多进程或分布式环境中，CPU 时间甚至可能不可用，因此壁钟时间成为跨平台、跨架构比较的统一标尺。基准测试的目标是回答"这个算法在真实部署场景中需要等待多久"，壁钟时间直接对应用户的等待体感，是评估工程实现差异的最终标尺。
为什么层次聚类基准固定 n_clusters=10？ n_clusters 的变化会影响合并操作的终止时机，但在基准设计中需要控制这个变量。固定为 10 使得链接策略的比较更加纯粹，避免簇数差异引入的混淆。如果允许 n_clusters 变化，结果将混合"算法本身的复杂度"与"合并终止时机"两个因素，无法单独评估链接策略对计算复杂度的影响。10 这个值既保证合并过程有足够步数展示链接差异，又不会让计算量过小导致计时精度不足。从统计学角度看，固定终止条件可以视为对四种链接策略施加相同的"问题规模"，使其计算量比较具有可比性。
为什么最近邻基准分开测量 fit 和 kneighbors？ 树算法的索引构建是"一次性投资"，查询是"多次回报"。这种分相计时揭示了算法的架构权衡：brute 的 fit 是 O(1)（无需预处理），但 query 是 O(N)；树算法的 fit 是 O(N log N) 或 O(N²)，但 query 可以达到 O(log N) 或 O(N/D)。在实际部署中，如果查询次数远大于构建次数（如在线服务场景），树算法的总成本可能更低；如果是离线批量查询（一次性处理），brute 可能更优。分开计时让用户能够根据应用场景选择最优算法。这种分相设计是软件性能分析中的经典方法，类似 profiling 中的 split-phase timing，能够精确识别"投资"与"回报"阶段的成本分配，帮助架构师做出符合访问模式的算法选型决策。
83.9 动手练习
-
扩展 KMeans 基准：引入 Elkan 算法与稀疏数据
-
修改
bench_plot_fastkmeans.py的compute_bench函数，增加KMeans(algorithm="elkan")作为第三个对比对象 -
增加稀疏矩阵数据生成分支 (使用
scipy.sparse.random)，对比稠密/稀疏数据下三种算法的相对性能 -
将 3D 曲面图升级为分组柱状图，更清晰展示三算法在稠密/稀疏两种数据下的 Speed/Quality 对比
-
思考问题：Elkan 算法利用三角不等式避免距离计算，在什么数据规模/维度下优势最大？稀疏数据下 MiniBatchKMeans 的批次处理机制是否仍有优势？
-
-
深度剖析最近邻算法的'维度灾难'临界点
-
基于
bench_plot_neighbors.py设计实验：固定 N=5000, k=10，将 D 扩展到 2^0 ~ 2^10 (1~1024 维) -
记录 kd_tree、ball_tree、brute 的查询耗时，绘制查询耗时随维度增长的对数坐标曲线
-
计算并标注 kd_tree/ball_tree 相对于 brute 的加速比 (speedup = brute_time / tree_time)
-
寻找加速比跌破 1.0 (即树算法不如暴力) 的'维度灾难'临界维度
-
思考问题：为什么高维下树算法会退化？结合'超球体体积集中在表面'和'空间分割失效'解释，leaf_size 参数如何影响这个临界点？
-
-
Ward 算法跨库基准的统计严谨性增强
-
改进
bench_plot_ward.py使其达到发布级基准标准：每个 (n, p) 组合重复运行 5 次，取中位数而非单次计时 (参考bench_pca_solvers.py的measure_one) -
增加预热运行 (warm-up) 消除 JIT/缓存冷启动影响
-
计算 95% 置信区间，在热力图上用半透明误差条或数值标注展示
-
增加内存峰值监控 (使用
tracemalloc或memory_profiler)，绘制内存比率热力图 -
思考问题：单次计时为何不可靠？OS 调度、GC、CPU 频率调节等噪声源有哪些？为什么 scikit-learn 的 Ward 在小规模数据上可能慢于 SciPy？
-
-
设计统一的基准测试可视化组件库
-
观察四个脚本的可视化代码，提取共性设计模式，实现一个可复用的
BenchmarkVisualizer类 -
统一的图表风格：配色方案、字体大小、网格、图例位置
-
通用的 3D 曲面图方法：
plot_surface_3d(X, Y, Z, xlabel, ylabel, zlabel, title) -
通用的分面折线图方法：
plot_facet_lines(results_dict, x_vals, facet_vals, xlabel, ylabel) -
通用的堆积柱状图方法：
plot_stacked_bars(categories, build_times, query_times, xlabel) -
通用的比率热力图方法：
plot_ratio_heatmap(ratio_matrix, x_vals, y_vals, xlabel, ylabel) -
将四个脚本的绘图代码重构为调用该组件库
-
思考问题：如何设计 API 兼容不同基准的数据结构 (defaultdict, nested dict, ndarray)？组件库应放在
benchmarks/utils/plotting.py还是asv_benchmarks/benchmarks/common.py？
-
83.10 本章小结
这一章中我们深入探索了 scikit-learn 基准测试体系中的聚类与最近邻专题，揭开了算法性能评估的工程方法论。通过四个独立基准脚本，我们学习了如何通过控制变量法、对数网格扫描、三因子正交实验等设计原则，构建可靠、可复现的性能基准实验。首先我们剖析了 KMeans 与 MiniBatchKMeans 的二维网格扫描，理解批次大小对收敛速度与质量的非线性权衡；其次我们探索了层次聚类四种链接策略的分面折线图可视化，揭示链接方式对计算复杂度的决定性影响；接着我们深入最近邻搜索的堆积柱状图设计，精确隔离样本量、维度、邻居数三个因子对 kd_tree、ball_tree、brute 算法的独立作用；最后我们学习了 Ward 算法跨库基准的比率热力图，直观展示 sklearn 与 SciPy 在不同数据规模下的胜负分界线。
本章我们一起学习了以下概念：
| 概念 | 解释 |
|------|------|
| compute_bench (fastkmeans) | 二维网格扫描揭示 KMeans/MiniBatchKMeans 随样本/特征规模的性能拓扑 |
| compute_bench_2 | 单维扫描批次大小，量化 MiniBatchKMeans 收敛速度与质量的权衡曲线 |
| 3D 曲面图可视化 | Speed/Quality 双指标 3D 曲面对比，统一 Z 轴上限保证跨算法可比性 |
| 层次聚类链接横评 | 固定算法切换 linkage 参数，分面折线图展示四策略随样本/维度的计时差异 |
| 三因子正交实验设计 | N/D/k 单因子变化法，隔离最近邻搜索中样本量、维度、邻居数的独立影响 |
| 分相计时 (fit vs query) | 索引构建与查询分开计时，揭示树算法'建图贵查路快' vs 暴力'建图免费查路贵'的架构权衡 |
| 堆积柱状图 + 对数坐标 | 红底蓝顶堆积柱+对数Y轴+动态基线，跨数量级展示构建/查询耗时占比 |
| 对数网格扫描 | 指数级步长覆盖 n_samples/n_features 3 个数量级，构建完整性能地形图 |
| 比率热力图 + 等高线 | sklearn_time/scipy_time 比率取对数着色，ratio=1 黑线精确划分两库胜负区域 |
| 工程基准测试方法论 | 控制变量法、随机种子固定、多轮重复、真实/合成数据双源验证、可视化即分析 |
下一章中，我们将继续深入独立基准脚本的世界，聚焦 GLM、Lasso、SGD 回归、稀疏化预测与 OMP/LARS 等线性模型专题，剖析正则化路径计算、Gram 矩阵预计算与坐标下降对回归性能的影响，探索稀疏表示如何在不损失精度的前提下压缩预测耗时。
第 84 章 —— 独立基准脚本：回归与稀疏专题 —— 剖析"线性模型的瘦身之道"
84.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解广义线性模型（Ridge、OLS、LassoLars）在不同维度下的性能对比方法
-
掌握坐标下降与最小角回归在Lasso路径计算上的效率差异
-
掌握SGD回归变体（普通、平均、ElasticNet、Ridge）的超参数调优与基准测试策略
-
理解模型稀疏化对推理加速的量化评估方法
-
掌握OMP与LARS在稀疏编码任务上的性能对比与Gram矩阵预计算影响分析
-
了解跨语言基准测试（sklearn vs glmnet）的工程实现细节
84.2 生活类比
把这一章的七个基准脚本想象成一场算法奥林匹克运动会的不同赛场，七场赛事从不同侧面刻画了 scikit-learn 中线性模型的工程表现。bench_glm.py 是一场短跑预赛——三种线性模型在随机跑道（数据）上同距离冲刺，看谁随维度增长最稳；bench_lasso.py 则是中长跑分组赛——Lasso（坐标下降）与 LassoLars（LARS）两种战术在固定跑道宽度（特征）与不同圈数（样本）的组合上反复比拼，验证两种路径算法在不同维度组合下的效率分水岭；bench_plot_lasso_path.py 是一次全程计时摄影——不只看终点，要拍全程轨迹（正则化路径），四种计时器（LARS/CD ±Gram）在样本-特征二维赛场上采集 25 个测点的耗时，绘制成 3D 耗时地形图；bench_glmnet.py 是跨国对抗赛——Python 阵营的 sklearn 与 R 阵营的 glmnet 在同一器材（数据/超参）上比速度、比精度（RMSE）、比动作标准度（系数差异），衡量两个生态系统在 Lasso 实现上的工程差异；bench_sgd_regression.py 是接力赛与变阵实验——四支队伍（ElasticNet/SGD/A-SGD/Ridge）在不同赛道长度（样本）和宽度（特征）上反复试错，目标是调配速度（学习率）与耐力（迭代）找到最优组合；bench_sparsify.py 是轻装上阵挑战赛——训练时背重包（稠密权重），比赛时扔掉零负重（sparsify），看能不能在稀疏跑道上跑得更快且成绩（R²）不降级；bench_plot_omp_lars.py 则是双人同跑 PK 赛——OMP 与 LARS 在稀疏编码赛道上，带/不带向导图（Gram 矩阵）各跑一遍，最终输出谁领先多少的热力地图。
84.3 源码地图
benchmarks/bench_glm.py
├── __main__ # GLM三模型基准测试入口
│ ├── 数据生成: 随机方阵 X, Y
│ ├── 计时器: datetime.now() 墙钟时间
│ ├── 模型训练: Ridge / LinearRegression / LassoLars
│ └── 可视化: Matplotlib 折线图对比
benchmarks/bench_lasso.py
├── compute_bench() # 核心基准计算函数
│ ├── 数据生成: make_regression (10% informative)
│ ├── 数据归一化: X /= sqrt(sum(X**2, axis=0))
│ ├── GC控制: gc.collect() 减少抖动
│ ├── 模型训练: Lasso / LassoLars (alpha=0.01)
│ └── 计时: time() 秒级精度
├── __main__ # 双场景实验入口
│ ├── 场景1: 固定特征(10) 变样本(100~1e6) + Gram预计算
│ ├── 场景2: 固定样本(2000) 变特征(500~3000) 无预计算
│ └── 可视化: 2x1 子图对比
benchmarks/bench_plot_lasso_path.py
├── compute_bench() # 四路径基准计算
│ ├── 数据生成: 低秩+粗尾分布 (effective_rank)
│ ├── Gram预计算: G=X.T@X, Xy=X.T@y
│ ├── 四路径对比: lars_path(±Gram) / lasso_path(±precompute)
│ ├── 结果聚合: defaultdict(list) -> 2D矩阵
│ └── GC/刷新: gc.collect() + sys.stdout.flush()
├── __main__ # 3D曲面可视化入口
│ ├── 网格: samples/features 10~2000 取5点
│ ├── 绘图: plot_surface 4子图
│ └── 归一化: Z轴统一至 max_time*1.1
benchmarks/bench_glmnet.py
├── rmse() # 均方根误差计算工具
├── bench() # 统一基准工厂函数
│ ├── 输入: factory, X, Y, X_test, Y_test, ref_coef
│ ├── 计时/预测/RMSE/系数差异
│ └── 返回: 训练耗时
├── __main__ # 跨语言对决入口
├── 实验1: 固定特征(1000) 变样本(500~10000)
├── 实验2: 固定样本(500) 变特征(100~2000)
├── 数据划分: 训练集前i*step, 测试集固定后1000
└── 双图输出: sklearn(蓝) vs glmnet(红)
benchmarks/bench_sgd_regression.py
├── __main__ # SGD回归全家桶基准
│ ├── 参数网格: samples(100~10k) x features(10/100/1000)
│ ├── 数据标准化: 训练集统计量 -> 测试集复用
│ ├── 四模型配置:
│ │ ├── ElasticNet (CD求解)
│ │ ├── SGDRegressor (invscaling, eta0=0.01, power_t=0.25)
│ │ ├── A-SGD (eta0=0.002, power_t=0.05, average=True)
│ │ └── Ridge (解析/迭代求解)
│ ├── 结果存储: (n_samples, n_features, 2) RMSE+Time
│ └── 可视化: m行2列子图矩阵
benchmarks/bench_sparsify.py
├── sparsity_ratio() # 稀疏度计算工具函数
├── benchmark_dense_predict() # 稠密预测微基准函数
├── benchmark_sparse_predict() # 稀疏预测微基准函数
├── score() # R²评分打印工具
├── __main__ # 稀疏化推理加速验证
│ ├── 数据构造: 5000x300 人工稀疏化输入/系数
│ ├── 模型: SGDRegressor(penalty='l1', alpha=0.2)
│ ├── 稀疏化: clf.sparsify() -> coef_ 变 csr_matrix
│ ├── 微基准函数: benchmark_dense/sparse_predict
│ │ ├── 300次循环预测
│ │ └── 适配 kernprof 逐行剖析
│ ├── 评估: R²无损 + 稀疏输入推理加速~30%
│ └── 稀疏度度量: nnz / (n_samples * n_features)
benchmarks/bench_plot_omp_lars.py
├── compute_bench() # OMP vs LARS 效率对比
│ ├── 数据生成: make_sparse_coded_signal
│ │ ├── n_samples=1, n_components=n_features
│ │ ├── n_features=n_samples (转置后对称)
│ │ └── n_nonzero_coefs = n_features//10
│ ├── 四组合: LARS/OMP x ±Gram预计算
│ ├── Gram预计算: G=X.T@X, Xy=X.T@y
│ ├── 结果: 时间比率矩阵 (LARS/OMP)
│ └── 对称扫描: samples/features 1000~5000 同步5点
├── __main__ # 热力图可视化入口
├── 双子图: 带Gram / 不带Gram
├── 配色: 中心1对称, >1 OMP快, <1 LARS快
└── 布局: colorbar水平置底
84.4 广义线性模型对比 —— Ridge、OLS 与 LassoLars 的"三角赛跑"
84.4.1 生活类比
把这场"三角赛跑"想象成一场体测三项全能——三个运动员（算法）要依次完成立定跳远（Ridge）、跳高（OLS）、撑杆跳（LassoLars），而维度就如同跳高架的逐级升高的横杆，从 500 厘米一路升到 20000 厘米。三项都是"跳"，但风格各异：Ridge 像穿着铅衣跳——闭式解保证步态稳定但负重大；OLS 像赤脚跳——伪逆解直接但关节磨损快；LassoLars 像走钢丝式的逐级逼近——每一步都在计算方向与幅度的微妙平衡。三条折线就像三个运动员的成绩单，谁在低维度轻盈谁在高维度吃力，一目了然。
核心任务在于使用随机方阵数据，对比 Ridge（岭回归）、OLS（普通最小二乘）与 LassoLars（最小角回归 Lasso）三种线性模型的训练耗时随维度增长的曲线。具体而言，维度从 500 递增至 20000，步长 500，共 40 个数据点。
关键实现细节包括以下几个方面：数据生成使用 X = np.random.randn(n_samples, n_features) 与 Y = np.random.randn(n_samples)，确保每次迭代数据独立同分布；计时器使用 datetime.now() 进行墙钟时间测量，单位为秒；模型配置上 Ridge(alpha=1.0)、LinearRegression()、LassoLars() 均使用默认求解器，fit_intercept=True（默认）；结果可视化用 Matplotlib 绘制折线图，红/绿/蓝分别对应 Ridge/OLS/LassoLars，横轴为特征维度，纵轴为训练时间。
源码路径：benchmarks/bench_glm.py - __main__（18-57行）
if __name__ == "__main__":
import matplotlib.pyplot as plt
n_iter = 40 # 总迭代次数，对应40个维度点
time_ridge = np.empty(n_iter) # 存储Ridge模型耗时
time_ols = np.empty(n_iter) # 存储OLS模型耗时
time_lasso = np.empty(n_iter) # 存储LassoLars模型耗时
dimensions = 500 * np.arange(1, n_iter + 1) # 维度从500到20000
for i in range(n_iter):
print("Iteration %s of %s" % (i, n_iter)) # 进度日志
n_samples, n_features = 10 * i + 3, 10 * i + 3 # 方阵维度同步增长
X = np.random.randn(n_samples, n_features) # 随机设计矩阵
Y = np.random.randn(n_samples) # 随机目标向量
start = datetime.now() # 墙钟计时起点
ridge = linear_model.Ridge(alpha=1.0)
ridge.fit(X, Y) # 训练Ridge（闭式解）
time_ridge[i] = (datetime.now() - start).total_seconds() # 记录耗时
start = datetime.now() # 重新计时OLS
ols = linear_model.LinearRegression()
ols.fit(X, Y) # 训练OLS（伪逆解）
time_ols[i] = (datetime.now() - start).total_seconds()
start = datetime.now() # 重新计时LassoLars
lasso = linear_model.LassoLars() # 实例化LassoLars模型
lasso.fit(X, Y) # 训练LassoLars（最小角回归）
time_lasso[i] = (datetime.now() - start).total_seconds()
plt.figure("scikit-learn GLM benchmark results") # 创建图形
plt.xlabel("Dimensions")
plt.ylabel("Time (s)")
plt.plot(dimensions, time_ridge, color="r") # 红色Ridge曲线
plt.plot(dimensions, time_ols, color="g") # 绿色OLS曲线
plt.plot(dimensions, time_lasso, color="b") # 蓝色LassoLars曲线
plt.legend(["Ridge", "OLS", "LassoLars"], loc="upper left")
plt.axis("tight")
plt.show()
这段代码定义了 GLM 三模型基准测试的主流程。它通过 datetime.now() 计时，依次训练 Ridge、OLS、LassoLars 三种线性模型，维度从 500 递增到 20000。最终用 Matplotlib 绘制红绿蓝三条曲线，横轴为维度，纵轴为训练耗时，直观展示三种算法在不同数据规模下的扩展性差异。注意此脚本使用 datetime.now() 而非 time.time()，前者精度为微秒级，更适合长时间运行的基准测试。
下面用 Mermaid 图展示三种 GLM 算法的计算复杂度差异：
84.5 Lasso 对决 LassoLars —— 坐标下降 vs 最小角回归的"路径之争"
84.5.1 生活类比
把这场"Lasso 对决 LassoLars"想象成两种登山策略的 PK——坐标下降（CD）像是一个"沿单条山脊逐级攀升"的策略，每一步只调整一个变量方向，适合平坦宽阔的山地（样本多、特征少）；最小角回归（LARS）则像是"沿等高线角度切入"的策略，每一步都在寻找当前残差与所有特征夹角的等分线，适合山势险峻（特征多）的小队伍。固定特征（10）变样本的实验像是在固定海拔的山地反复清点人数——CD 在人山人海中凭借"单兵作战"的优势能快速推进，LARS 在庞大队伍中重新计算等分角的开销逐渐成为瓶颈；固定样本（2000）变特征的实验则像是在固定人数的登山队中不断切换攀登工具——CD 在工具增多时依然能逐项处理，但 LARS 在特征维度膨胀时重新计算等分角的代价急剧上升。两次实验的交叉点正是两种算法的"效率分水岭"。
核心任务包含两大实验场景：固定特征数（10）变样本数（1001e6），固定样本数（2000）变特征数（5003000）。两个场景均仅 10% 特征为信息量特征（n_informative = nf // 10），模拟稀疏真实场景。数据预处理阶段按列归一化 X /= np.sqrt(np.sum(X**2, axis=0))，消除尺度影响。
关键实现细节覆盖以下几个层面：核心函数 compute_bench 通过双重循环遍历样本/特征组合，每轮强制 gc.collect() 减少 GC 抖动；两模型统一超参数为 alpha=0.01、fit_intercept=False，precompute 由外部控制（首实验 True，次实验 False）；计时精度使用 time() 秒级计时，记录 fit 端到端耗时；可视化采用两行子图，上行固定特征变样本，下行固定样本变特征，蓝/红线分别代表 Lasso/LassoLars。
源码路径：benchmarks/bench_lasso.py - compute_bench()（17-53行）
def compute_bench(alpha, n_samples, n_features, precompute):
lasso_results = [] # 收集Lasso(CD)耗时
lars_lasso_results = [] # 收集LassoLars(LARS)耗时
it = 0
for ns in n_samples: # 外层遍历样本数
for nf in n_features: # 内层遍历特征数
it += 1
print("==================")
print("Iteration %s of %s" % (it, max(len(n_samples), len(n_features))))
print("==================")
n_informative = nf // 10 # 仅10%特征为信息量
X, Y, coef_ = make_regression( # 合成稀疏回归数据
n_samples=ns,
n_features=nf,
n_informative=n_informative,
noise=0.1,
coef=True,
)
X /= np.sqrt(np.sum(X**2, axis=0)) # 按列归一化，消除尺度影响
gc.collect() # 强制GC，减少内存抖动
print("- benchmarking Lasso")
clf = Lasso(alpha=alpha, fit_intercept=False, precompute=precompute) # CD求解器
tstart = time() # 秒级计时起点
clf.fit(X, Y)
lasso_results.append(time() - tstart) # 记录Lasso耗时
gc.collect() # 再次GC，确保两次实验条件一致
print("- benchmarking LassoLars")
clf = LassoLars(alpha=alpha, fit_intercept=False, precompute=precompute) # LARS求解器
tstart = time()
clf.fit(X, Y)
lars_lasso_results.append(time() - tstart) # 记录LassoLars耗时
return lasso_results, lars_lasso_results # 返回两组耗时列表
这段代码是 Lasso 双模型基准的核心计算函数。它通过双重循环遍历样本与特征组合，使用 make_regression 生成仅 10% 特征有信息量的稀疏数据。每次训练前强制 gc.collect() 减少垃圾回收抖动，并按列归一化消除数据尺度对算法收敛速度的影响。Lasso 使用坐标下降（CD）求解器，LassoLars 使用最小角回归（LARS）求解器，两者共享 alpha=0.01 与 fit_intercept=False 设置。
源码路径：benchmarks/bench_lasso.py - __main__（55-97行）
if __name__ == "__main__":
import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso, LassoLars
alpha = 0.01 # 统一正则化参数
# ========== 场景1：固定特征变样本 ==========
n_features = 10 # 低维特征
list_n_samples = np.linspace(100, 1000000, 5).astype(int) # 样本从100到1e6
lasso_results, lars_lasso_results = compute_bench( # 启用Gram预计算
alpha, list_n_samples, [n_features], precompute=True
)
plt.figure("scikit-learn LASSO benchmark results")
plt.subplot(211) # 上子图：固定特征变样本
plt.plot(list_n_samples, lasso_results, "b-", label="Lasso") # 蓝线Lasso
plt.plot(list_n_samples, lars_lasso_results, "r-", label="LassoLars") # 红线LassoLars
plt.title("precomputed Gram matrix, %d features, alpha=%s" % (n_features, alpha))
plt.legend(loc="upper left")
plt.xlabel("number of samples")
plt.ylabel("Time (s)")
plt.axis("tight")
# ========== 场景2：固定样本变特征 ==========
n_samples = 2000 # 中等样本量
list_n_features = np.linspace(500, 3000, 5).astype(int) # 特征从500到3000
lasso_results, lars_lasso_results = compute_bench( # 不使用Gram预计算
alpha, [n_samples], list_n_features, precompute=False
)
plt.subplot(212) # 下子图：固定样本变特征
plt.plot(list_n_features, lasso_results, "b-", label="Lasso")
plt.plot(list_n_features, lars_lasso_results, "r-", label="LassoLars")
plt.title("%d samples, alpha=%s" % (n_samples, alpha))
plt.legend(loc="upper left")
plt.xlabel("number of features")
plt.ylabel("Time (s)")
plt.axis("tight")
plt.show()
这段代码是 Lasso 基准测试的实验入口与可视化部分。它设计了两个互补的实验场景——固定 10 维特征变化样本数（100 到 100 万）启用 Gram 预计算，以及固定 2000 样本变化特征数（500 到 3000）不预计算。两个子图共享 Lasso/LassoLars 对比框架，蓝线代表坐标下降，红线代表最小角回归，让读者直观看到两种算法在不同维度组合下的效率分水岭。
下面用 Mermaid 图展示双场景对比设计：
84.6 Lasso 正则化路径基准 —— LARS 与坐标下降的"全景扫描"
84.6.1 生活类比
把这场"全景扫描"想象成地形测绘任务——四组测绘队员（lars_path ±Gram、lasso_path ±precompute）被派到 25 块不同大小的土地上同时作业，每块地都有"样本数×特征数"的二维坐标。带 Gram 矩阵的队员如同自带地形数据库的测绘员，无需现场勘测就能定位每个特征的相关性；不带 Gram 的队员则需在每块地上现场重新计算内积，成本随土地面积急剧攀升。LARS 算法的队员每次移动都要更新完整的方向矩阵，因此自带数据库对其加速显著；坐标下降（CD）队员每次只调整单个坐标方向，查询需求小，预计算收益有限。最终绘制的 3D 曲面就是这片"样本-特征高原"的耗时地形图——蓝/青/红/黄四色曲面在同一海拔基准线（max_time * 1.1）上叠加，凹陷处即是算法的"性能洼地"。
核心任务是对比四种路径计算方式：lars_path（带/不带 Gram 矩阵）、lasso_path（带/不带 Gram 矩阵）。数据特征为低秩结构（effective_rank = min(n_samples, n_features) / 10）加粗尾分布，模拟真实高维数据。样本数与特征数均从 10 到 2000 取 5 个点，共 25 组组合，形成二维网格。
关键实现细节包括以下几个方面：Gram 矩阵预计算采用 G = np.dot(X.T, X) 与 Xy = np.dot(X.T, y)，传入 lars_path_gram 与 lasso_path(precompute=True)；结果聚合使用 defaultdict(list) 收集四组耗时，最终重塑为 (n_samples, n_features) 矩阵；3D 可视化使用 Matplotlib plot_surface 绘制四个曲面子图，Z 轴统一归一化至 max_time * 1.1，便于横向比较。
源码路径：benchmarks/bench_plot_lasso_path.py - compute_bench()（17-71行）
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("====================")
dataset_kwargs = {
"n_samples": n_samples,
"n_features": n_features,
"n_informative": n_features // 10, # 10%信息量特征
"effective_rank": min(n_samples, n_features) / 10, # 低秩结构
"bias": 0.0,
}
print("n_samples: %d" % n_samples)
print("n_features: %d" % n_features)
X, y = make_regression(**dataset_kwargs) # 合成低秩数据
gc.collect()
print("benchmarking lars_path (with Gram):", end="")
sys.stdout.flush() # 实时刷新输出
tstart = time()
G = np.dot(X.T, X) # 预计算Gram矩阵
Xy = np.dot(X.T, y)
lars_path_gram(Xy=Xy, Gram=G, n_samples=y.size, method="lasso") # Gram-LARS
delta = time() - tstart
print("%0.3fs" % delta)
results["lars_path (with Gram)"].append(delta)
gc.collect()
print("benchmarking lars_path (without Gram):", end="")
sys.stdout.flush()
tstart = time()
lars_path(X, y, method="lasso") # 原始LARS路径
delta = time() - tstart
print("%0.3fs" % delta)
results["lars_path (without Gram)"].append(delta)
gc.collect()
print("benchmarking lasso_path (with Gram):", end="")
sys.stdout.flush()
tstart = time()
lasso_path(X, y, precompute=True) # CD路径带Gram
delta = time() - tstart
print("%0.3fs" % delta)
results["lasso_path (with Gram)"].append(delta)
gc.collect()
print("benchmarking lasso_path (without Gram):", end="")
sys.stdout.flush()
tstart = time()
lasso_path(X, y, precompute=False) # CD路径不带Gram
delta = time() - tstart
print("%0.3fs" % delta)
results["lasso_path (without Gram)"].append(delta)
return results # 返回四组耗时字典
这段代码定义了 Lasso 正则化路径的四象限基准计算函数。它在样本-特征二维网格上扫描四种路径算法组合——lars_path 带/不带 Gram 与 lasso_path 带/不带 Gram。每次实验前 gc.collect() 减少内存抖动，实验中 sys.stdout.flush() 实时刷新进度。关键设计：通过预计算 Gram 矩阵 G = X.T @ X 与 Xy = X.T @ y，将 lars_path_gram 的输入从设计矩阵简化为二次型，显著加速高维场景。
源码路径：benchmarks/bench_plot_lasso_path.py - __main__（73-108行）
if __name__ == "__main__":
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d # 注册3D投影
samples_range = np.linspace(10, 2000, 5).astype(int) # 5个样本点
features_range = np.linspace(10, 2000, 5).astype(int) # 5个特征点
results = compute_bench(samples_range, features_range)
max_time = max(max(t) for t in results.values()) # 找到全局最大耗时
fig = plt.figure("scikit-learn Lasso path benchmark results")
i = 1
for c, (label, timings) in zip("bcry", sorted(results.items())): # 四色对应四算法
ax = fig.add_subplot(2, 2, i, projection="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.8) # 绘制3D曲面
ax.set_xlabel("n_samples")
ax.set_ylabel("n_features")
ax.set_zlabel("Time (s)")
ax.set_zlim3d(0.0, max_time * 1.1) # Z轴统一归一化
ax.set_title(label)
i += 1
plt.show()
这段代码是 Lasso 路径基准的 3D 可视化入口。它将 compute_bench 返回的耗时字典重塑为二维矩阵，绘制四个 3D 曲面。核心设计：Z 轴统一归一化至 max_time * 1.1，保证四个子图的 Z 轴范围一致，便于横向比较算法耗时地形。bcry 颜色序列对应蓝/青/红/黄四种路径算法的曲面，直观呈现它们在样本-特征空间的耗时分布。sorted(results.items()) 按字符串字母序排序四个 key，依次为 lars_path (with Gram) / lars_path (without Gram) / lasso_path (with Gram) / lasso_path (without Gram)，分别绑定蓝/青/红/黄四色。
下面用 Mermaid 图展示四路径基准的数据流与可视化链路：
84.7 glmnet 大对决 —— sklearn 与 R 生态的"跨语言擂台"
84.7.1 生活类比
把这场"跨语言擂台"想象成两位厨师的世界级厨艺比拼——一位来自意大利的 Python 厨师（sklearn），另一位来自法国的 R 厨师（glmnet），他们要在两个"厨房"（不同数据集）里比拼做同一道菜（Lasso）。实验 1 让两人同时做一道"千人宴"（样本数 500→10000，固定 1000 食材），看谁能在食客增多时仍保持稳定出餐速度；实验 2 改为"百人品鉴会"（固定 500 食客，食材从 100 增到 2000），看谁在菜品种类爆炸时仍能精准调味。三重评判标准——出餐速度（训练耗时）、菜的味道（测试 RMSE）、摆盘精度（系数与真实值偏差）——综合评定两位厨师谁更胜一筹。这场比赛的趣味在于：跨语言的实现差异（解释型 vs 编译型、API 设计风格、底层 BLAS 调用）都会成为胜负手。
核心任务是横向对比 sklearn.linear_model.Lasso 与 glmnet.elastic_net.Lasso 在相同数据、相同超参数（alpha=0.1）下的训练速度与预测精度。扫描双维度：固定 1000 特征变样本数（500~10000，步长 500），固定 500 样本变特征数（100~2000，步长 100）。评价指标涵盖训练耗时、测试集 RMSE、系数平均绝对差（相对真实系数 coef_）。
关键实现细节包括以下几点：统一基准函数 bench(factory, X, Y, X_test, Y_test, ref_coef) 通过注入模型工厂，消除代码重复；数据划分采用训练集前 i*step 样本，测试集固定后 1000 样本，保证测试集分布一致；强制 gc.collect() 前置，减少内存管理噪声；双图输出第一图随样本数增长，第二图随特征数增长，蓝/红线对应 sklearn/glmnet；辅助工具函数 rmse() 计算均方根误差，用于量化预测精度差异。
源码路径：benchmarks/bench_glmnet.py - rmse()（19-21行）
def rmse(a, b):
return np.sqrt(np.mean((a - b) ** 2)) # 均方根误差计算
这段代码定义了 RMSE 计算工具函数。它使用 NumPy 向量化操作计算预测值与真实值之间的均方根误差，作为跨语言基准的精度指标之一。
源码路径：benchmarks/bench_glmnet.py - bench()（23-33行）
def bench(factory, X, Y, X_test, Y_test, ref_coef):
gc.collect() # 强制GC，保证内存状态一致
tstart = time()
clf = factory(alpha=alpha).fit(X, Y) # 工厂注入模型实例化与训练
delta = time() - tstart
# stop time
print("duration: %0.3fs" % delta) # 输出训练耗时
print("rmse: %f" % rmse(Y_test, clf.predict(X_test))) # 输出测试RMSE
print("mean coef abs diff: %f" % abs(ref_coef - clf.coef_.ravel()).mean()) # 系数偏差
return delta # 返回耗时供上层聚合
这段代码定义了跨语言基准的统一工厂函数。它通过 factory(alpha=alpha).fit(X, Y) 模式支持不同框架的模型注入（sklearn 或 glmnet），消除了重复的计时与评估代码。关键设计：同时输出训练耗时、测试 RMSE 与系数平均绝对差三个指标，从速度、精度、参数一致性三个维度全面对比。
源码路径：benchmarks/bench_glmnet.py - __main__（35-107行）
if __name__ == "__main__":
# Delayed import of matplotlib.pyplot
import matplotlib.pyplot as plt
from glmnet.elastic_net import Lasso as GlmnetLasso # 延迟导入glmnet
from sklearn.linear_model import Lasso as ScikitLasso
scikit_results = [] # sklearn耗时列表
glmnet_results = [] # glmnet耗时列表
n = 20
step = 500
n_features = 1000 # 固定1000特征
n_informative = n_features // 10 # 100个信息量特征（整数除法）
n_test_samples = 1000 # 固定1000测试样本
for i in range(1, n + 1):
print("==================")
print("Iteration %s of %s" % (i, n))
print("==================")
X, Y, coef_ = make_regression( # 一次性生成足够多数据
n_samples=(i * step) + n_test_samples,
n_features=n_features,
noise=0.1,
n_informative=n_informative,
coef=True,
)
X_test = X[-n_test_samples:] # 测试集固定尾部1000样本
Y_test = Y[-n_test_samples:]
X = X[: (i * step)] # 训练集递增取前i*step样本
Y = Y[: (i * step)]
print("benchmarking scikit-learn: ")
scikit_results.append(bench(ScikitLasso, X, Y, X_test, Y_test, coef_)) # sklearn
print("benchmarking glmnet: ")
glmnet_results.append(bench(GlmnetLasso, X, Y, X_test, Y_test, coef_)) # glmnet
plt.clf()
xx = range(0, n * step, step)
plt.title("Lasso regression on sample dataset (%d features)" % n_features)
plt.plot(xx, scikit_results, "b-", label="scikit-learn") # 蓝线sklearn
plt.plot(xx, glmnet_results, "r-", label="glmnet") # 红线glmnet
plt.legend()
plt.xlabel("number of samples to classify")
plt.ylabel("Time (s)")
plt.show()
# ========== 实验2：固定样本变特征 ==========
scikit_results = []
glmnet_results = []
n = 20
step = 100
n_samples = 500 # 固定500样本
for i in range(1, n + 1):
print("==================")
print("Iteration %02d of %02d" % (i, n))
print("==================")
n_features = i * step # 特征从100递增到2000
n_informative = n_features // 10 # 整数除法避免浮点
X, Y, coef_ = make_regression(
n_samples=(i * step) + n_test_samples,
n_features=n_features,
noise=0.1,
n_informative=n_informative,
coef=True,
)
X_test = X[-n_test_samples:]
Y_test = Y[-n_test_samples:]
X = X[:n_samples] # 训练集固定前500样本
Y = Y[:n_samples]
print("benchmarking scikit-learn: ")
scikit_results.append(bench(ScikitLasso, X, Y, X_test, Y_test, coef_))
print("benchmarking glmnet: ")
glmnet_results.append(bench(GlmnetLasso, X, Y, X_test, Y_test, coef_))
xx = np.arange(100, 100 + n * step, step)
plt.figure("scikit-learn vs. glmnet benchmark results")
plt.title("Regression in high dimensional spaces (%d samples)" % n_samples)
plt.plot(xx, scikit_results, "b-", label="scikit-learn")
plt.plot(xx, glmnet_results, "r-", label="glmnet")
plt.legend()
plt.xlabel("number of features")
plt.ylabel("Time (s)")
plt.axis("tight")
plt.show()
这段代码是跨语言对决的双实验入口。它设计了实验 1：固定 1000 特征变样本数（500~10000）与实验 2：固定 500 样本变特征数（100~2000）两个互补的扫描维度。每次迭代生成足够大的数据矩阵，训练集递增切片，测试集固定尾部 1000 样本，保证测试分布一致。蓝线代表 sklearn，红线代表 glmnet，从速度维度对比两个生态系统的 Lasso 实现。
下面用 Mermaid 图展示跨语言对决的双实验数据流：
84.8 SGD 回归基准 —— 随机梯度下降的"稳定与速度"
84.8.1 生活类比
把这场"SGD 回归基准"想象成四个赛艇队的水上拉力赛——ElasticNet 是四人同步划桨队（坐标下降，全局协调但队员多就慢）；普通 SGD 是单人双桨队（效率高但遇风浪易偏航）；A-SGD 是单桨加稳定翼队（学习率更小、桨频更慢但加装了平均稳定器，长距离更稳）；Ridge 是带引擎的快艇队（闭式解瞬时启动，但引擎在复杂水域不如划桨船灵活）。每组船队要在三个"训练水域"（10/100/1000 维特征）反复横渡不同长度的河道（100~10000 样本），同时记录航行精度（RMSE）和完赛时间（训练时间）。最终 3×2 子图矩阵就像三块水域的"战绩板"，揭示每个队伍在不同水域的速度-精度权衡。
核心任务围绕四模型混战展开：ElasticNet（坐标下降）、SGDRegressor（普通 SGD）、SGDRegressor(average=...)（平均 SGD）、Ridge（解析/迭代求解）。实验构建三维参数网格：样本数（100~10000，5 点）× 特征数（10/100/1000）× 2 指标（RMSE、训练时间）。数据标准化阶段训练集按列/目标零均值单位方差，测试集复用训练集统计量。
关键实现细节包括以下几点：SGD 超参精心调校为 alpha=alpha/n_train 对齐正则化强度，learning_rate='invscaling'、eta0=0.01、power_t=0.25，A-SGD 进一步降低 eta0=0.002、power_t=0.05 并开启 average；结果存储采用 (n_samples, n_features, 2) 三维数组，最后一维 0=RMSE 1=Time；可视化布局为 m x 2 子图矩阵，每行对应一种特征数，左列 RMSE 曲线，右列耗时曲线，四色图例贯穿始终。
源码路径：benchmarks/bench_sgd_regression.py - __main__（22-115行）
if __name__ == "__main__":
list_n_samples = np.linspace(100, 10000, 5).astype(int) # 5个样本数
list_n_features = [10, 100, 1000] # 3个特征数
n_test = 1000
max_iter = 1000
noise = 0.1
alpha = 0.01
sgd_results = np.zeros((len(list_n_samples), len(list_n_features), 2)) # 三维结果数组
elnet_results = np.zeros((len(list_n_samples), len(list_n_features), 2))
ridge_results = np.zeros((len(list_n_samples), len(list_n_features), 2))
asgd_results = np.zeros((len(list_n_samples), len(list_n_features), 2))
for i, n_train in enumerate(list_n_samples):
for j, n_features in enumerate(list_n_features):
X, y, coef = make_regression( # 合成回归数据
n_samples=n_train + n_test,
n_features=n_features,
noise=noise,
coef=True,
)
X_train = X[:n_train] # 切分训练测试
y_train = y[:n_train]
X_test = X[n_train:]
y_test = y[n_train:]
print("=======================")
print("Round %d %d" % (i, j))
print("n_features:", n_features)
print("n_samples:", n_train)
# Shuffle data
idx = np.arange(n_train)
np.random.seed(13)
np.random.shuffle(idx) # 固定种子洗牌
X_train = X_train[idx]
y_train = y_train[idx]
std = X_train.std(axis=0) # 训练集统计量
mean = X_train.mean(axis=0)
X_train = (X_train - mean) / std # 标准化特征
X_test = (X_test - mean) / std # 测试集复用统计量
std = y_train.std(axis=0)
mean = y_train.mean(axis=0)
y_train = (y_train - mean) / std # 标准化目标
y_test = (y_test - mean) / std
# ========== ElasticNet 基准 ==========
gc.collect()
print("- benchmarking ElasticNet")
clf = ElasticNet(alpha=alpha, l1_ratio=0.5, fit_intercept=False) # CD求解
tstart = time()
clf.fit(X_train, y_train)
elnet_results[i, j, 0] = mean_squared_error(clf.predict(X_test), y_test) # 记录RMSE
elnet_results[i, j, 1] = time() - tstart # 记录耗时
# ========== SGD 基准 ==========
gc.collect()
print("- benchmarking SGD")
clf = SGDRegressor( # 普通SGD
alpha=alpha / n_train, # 正则化按样本数缩放
fit_intercept=False,
max_iter=max_iter,
learning_rate="invscaling", # 学习率衰减策略
eta0=0.01,
power_t=0.25,
tol=1e-3,
)
tstart = time()
clf.fit(X_train, y_train)
sgd_results[i, j, 0] = mean_squared_error(clf.predict(X_test), y_test)
sgd_results[i, j, 1] = time() - tstart
# ========== A-SGD 基准 ==========
gc.collect()
print("max_iter", max_iter)
print("- benchmarking A-SGD")
clf = SGDRegressor( # 平均SGD
alpha=alpha / n_train,
fit_intercept=False,
max_iter=max_iter,
learning_rate="invscaling",
eta0=0.002, # 更小的初始学习率
power_t=0.05, # 更慢的衰减
tol=1e-3,
average=(max_iter * n_train // 2), # 开启平均
)
tstart = time()
clf.fit(X_train, y_train)
asgd_results[i, j, 0] = mean_squared_error(clf.predict(X_test), y_test)
asgd_results[i, j, 1] = time() - tstart
# ========== Ridge 基准 ==========
gc.collect()
print("- benchmarking RidgeRegression")
clf = Ridge(alpha=alpha, fit_intercept=False) # 闭式/迭代解
tstart = time()
clf.fit(X_train, y_train)
ridge_results[i, j, 0] = mean_squared_error(clf.predict(X_test), y_test)
ridge_results[i, j, 1] = time() - tstart
# ========== 可视化部分 ==========
i = 0
m = len(list_n_features)
plt.figure("scikit-learn SGD regression benchmark results", figsize=(5 * 2, 4 * m))
for j in range(m):
plt.subplot(m, 2, i + 1)
plt.plot(list_n_samples, np.sqrt(elnet_results[:, j, 0]), label="ElasticNet")
plt.plot(list_n_samples, np.sqrt(sgd_results[:, j, 0]), label="SGDRegressor")
plt.plot(list_n_samples, np.sqrt(asgd_results[:, j, 0]), label="A-SGDRegressor")
plt.plot(list_n_samples, np.sqrt(ridge_results[:, j, 0]), label="Ridge")
plt.legend(prop={"size": 10})
plt.xlabel("n_train")
plt.ylabel("RMSE")
plt.title("Test error - %d features" % list_n_features[j])
i += 1
plt.subplot(m, 2, i + 1) # 时间子图
plt.plot(list_n_samples, np.sqrt(elnet_results[:, j, 1]), label="ElasticNet")
plt.plot(list_n_samples, np.sqrt(sgd_results[:, j, 1]), label="SGDRegressor")
plt.plot(list_n_samples, np.sqrt(asgd_results[:, j, 1]), label="A-SGDRegressor")
plt.plot(list_n_samples, np.sqrt(ridge_results[:, j, 1]), label="Ridge")
plt.legend(prop={"size": 10})
plt.xlabel("n_train")
plt.ylabel("Time [sec]")
plt.title("Training time - %d features" % list_n_features[j])
i += 1
plt.subplots_adjust(hspace=0.30)
plt.show()
这段代码是 SGD 回归全家桶基准的主流程。它构建了样本数 × 特征数 × 模型的三维评估空间，记录每个组合的 RMSE 与耗时。关键设计：SGD 使用 alpha=alpha/n_train 缩放正则化强度以对齐 SGD 目标函数 1/n * loss + alpha * penalty 的数学定义；A-SGD 进一步降低 eta0=0.002 与 power_t=0.05 并开启 average=True，通过更温和的学习率衰减与参数平均换取更稳定的收敛。最终绘制 m×2 子图矩阵，每行对应一种特征数，左列 RMSE、右列耗时，四条曲线对比 ElasticNet/SGD/A-SGD/Ridge 在精度与速度上的权衡。
下面用 Mermaid 图展示 SGD 超参调优的对比：
84.9 稀疏化预测探秘 —— 系数瘦身带来的"推理加速器"
84.9.1 生活类比
把这场"稀疏化加速"想象成行李减重对跑步成绩的优化——一位马拉松选手在训练时背着满载重物（稠密系数 ndarray），跑完后把背包里的空水瓶、空包装统统扔掉（sparsify()），只留下真正用得上的能量棒（csr_matrix 中存储的非零权重）。这场比赛的趣味在于：选手本身没换（模型未重训），只是行囊变轻了。训练时背重包是为了容纳所有可能的"装备"（潜在的非零权重），比赛时扔掉零负重则是利用数据本身的稀疏特性——后 17% 的输入样本全为零，根本用不上对应装备。微基准函数 benchmark_dense_predict 与 benchmark_sparse_predict 就像在同一条赛道上用精密计时器分别测试"背包跑"和"轻装跑"——300 次循环消除预热噪声，1774 μs/次 vs 1271 μs/次的对比直观揭示了"数据结构匹配输入特性"带来的工程红利。
核心任务是验证 SGDRegressor(penalty='l1') 训练后调用 sparsify() 将稠密系数转为稀疏 CSR 矩阵，对推理速度的提升。实验合成 5000×300 稠密矩阵并人工稀疏化（前 2500 样本保留，后 2500 置零），真实系数 300 维仅前 150 非零。对比指标涵盖模型稀疏度、测试集 R²、300 次 predict 循环的总耗时（密集 vs 稀疏输入）。
关键实现细节覆盖以下方面：稀疏度度量函数 sparsity_ratio() 通过 np.count_nonzero(X) / (n_samples * n_features) 计算；稀疏化流程 clf.sparsify() 将 coef_ 从 ndarray 转为 csr_matrix，随后 predict 自动走稀疏内核；微基准函数 benchmark_dense_predict / benchmark_sparse_predict 分离热路径，便于 kernprof 逐行剖析；评分工具函数 score() 封装 R² 计算与打印。关键结论是稀疏模型在稀疏输入上推理快约 30%（示例中 1774 → 1271 μs/次），且 R² 无损。
源码路径：benchmarks/bench_sparsify.py - sparsity_ratio()（31-32行）
def sparsity_ratio(X):
return np.count_nonzero(X) / float(n_samples * n_features) # 稀疏度 = 非零元素占比
这段代码定义了稀疏度度量工具函数。它通过非零元素数量除以矩阵总元素数，得到 0~1 之间的稀疏度比率，用于量化输入数据、真实系数与模型权重的稀疏程度。
源码路径：benchmarks/bench_sparsify.py - benchmark_dense_predict()（52-54行）
@profile # kernprof逐行剖析装饰器
def benchmark_dense_predict():
for _ in range(300): # 300次循环消除JIT/缓存预热噪声
clf.predict(X_test) # 稠密矩阵预测
这段代码是稠密预测微基准函数。它通过 @profile 装饰器适配 kernprof 逐行剖析工具，循环 300 次 predict 以获取稳定的耗时统计。
源码路径：benchmarks/bench_sparsify.py - benchmark_sparse_predict()（57-60行）
@profile
def benchmark_sparse_predict():
X_test_sparse = csr_matrix(X_test) # 转换为CSR稀疏矩阵
for _ in range(300):
clf.predict(X_test_sparse) # 稀疏矩阵预测走稀疏内核
这段代码是稀疏预测微基准函数。它在循环外预先将测试矩阵转为 CSR 稀疏格式，循环内调用 predict 触发稀疏推理路径。关键设计：转换操作放在循环外，避免每次迭代重复转换开销，准确测量推理热路径耗时。
源码路径：benchmarks/bench_sparsify.py - score()（62-64行）
def score(y_test, y_pred, case):
r2 = r2_score(y_test, y_pred) # R²评分
print("r^2 on test data (%s) : %f" % (case, r2)) # 打印场景化结果
这段代码是 R² 评分打印工具函数。它将 R² 计算与场景标签（如"dense model"/"sparse model"）一起输出，便于对比稀疏化前后的精度差异。
源码路径：benchmarks/bench_sparsify.py - __main__（1-67行）
np.random.seed(42) # 固定随机种子保证可复现
def sparsity_ratio(X):
return np.count_nonzero(X) / float(n_samples * n_features)
n_samples, n_features = 5000, 300 # 数据规模
X = np.random.randn(n_samples, n_features)
inds = np.arange(n_samples)
np.random.shuffle(inds)
X[inds[int(n_features / 1.2) :]] = 0 # 人工稀疏化输入（int(300/1.2)=250，前250样本保留，后2750置零）
print("input data sparsity: %f" % sparsity_ratio(X))
coef = 3 * np.random.randn(n_features)
inds = np.arange(n_features)
np.random.shuffle(inds)
coef[inds[n_features // 2 :]] = 0 # 人工稀疏化系数（后一半置零）
print("true coef sparsity: %f" % sparsity_ratio(coef))
y = np.dot(X, coef)
y += 0.01 * np.random.normal((n_samples,)) # 添加噪声
n_samples = X.shape[0]
X_train, y_train = X[: n_samples // 2], y[: n_samples // 2]
X_test, y_test = X[n_samples // 2 :], y[n_samples // 2 :] # 50/50切分
print("test data sparsity: %f" % sparsity_ratio(X_test))
###############################################################################
clf = SGDRegressor(penalty="l1", alpha=0.2, max_iter=2000, tol=None) # L1正则化SGD
clf.fit(X_train, y_train)
print("model sparsity: %f" % sparsity_ratio(clf.coef_))
def benchmark_dense_predict():
for _ in range(300):
clf.predict(X_test)
def benchmark_sparse_predict():
X_test_sparse = csr_matrix(X_test)
for _ in range(300):
clf.predict(X_test_sparse)
def score(y_test, y_pred, case):
r2 = r2_score(y_test, y_pred)
print("r^2 on test data (%s) : %f" % (case, r2))
score(y_test, clf.predict(X_test), "dense model") # 稠密模型评分
benchmark_dense_predict() # 稠密推理基准
clf.sparsify() # 核心：系数稀疏化
score(y_test, clf.predict(X_test), "sparse model") # 稀疏模型评分
benchmark_sparse_predict() # 稀疏推理基准
这段代码是稀疏化推理加速验证的主流程。它首先构造 5000×300 的人工稀疏化数据（输入通过 int(300/1.2)=250 的切片下标保留前 250 样本、置零后 2750 行，真实系数后一半置零），训练 L1 正则化的 SGD 回归模型。然后依次执行稠密推理、稀疏化、稀疏推理三个阶段。核心结论：clf.sparsify() 将 coef_ 从 ndarray 转为 csr_matrix 后，推理热路径走稀疏内核，在稀疏输入上耗时从 1774 μs/次降至 1271 μs/次（加速约 30%），且 R² 完全一致。
下面用 Mermaid 图展示稀疏化推理加速的流程：
84.10 OMP 与 LARS 对垒 —— 稀疏编码的"精度与速度赛跑"
84.10.1 生活类比
把这场"OMP 与 LARS 对垒"想象成两位侦探破案方式的 PK——OMP（正交匹配追踪）像是一位"按图索骥"的侦探，每一步只找与当前线索最相关的单个证据（选原子），然后用最小二乘法做投影修正，工具简单但偶尔会走弯路；LARS（最小角回归）则像是一位"同步多线并进"的侦探，每一步都要计算所有证据与案件夹角的等分线，并实时更新"破案矩阵"（Cholesky 分解），信息更全但工具负担重。带 Gram 矩阵的侦探自带"证据关联数据库"，无需现场重新核对证据间的关系，因此对 LARS 这种"全局协调型"侦探加速明显；而 OMP 这种"局部决策型"侦探，本身只需要查每个证据的自相关（Gram 对角），完整数据库收益有限。比率热力图就像两位侦探在不同案件规模下的"破案效率对照表"——蓝色表示 LARS 领先、红色表示 OMP 领先、对角线附近两者势均力敌，一图尽览算法相对优势的二维分布。
核心任务是对比 orthogonal_mp (OMP) 与 lars_path (LARS) 在稀疏编码任务上的计算效率。实验采用四种组合：LARS/OMP × 带/不带 Gram 矩阵，输出耗时比率热力图（time(LARS)/time(OMP)）。数据生成使用 make_sparse_coded_signal 生成字典 X (n_samples, n_components) 与稀疏码 y，n_nonzero_coefs = n_features // 10。
关键实现细节涵盖以下几个层面：维度扫描时样本数/特征数同步从 1000 到 5000 取 5 点，共 25 组，矩阵对称；Gram 预计算采用 G = X.T @ X、Xy = X.T @ y，传入 lars_path_gram 与 orthogonal_mp(precompute=True)；结果归一化直接计算比率矩阵，matshow 以 1 为中心对称配色，>1 表示 OMP 更快，<1 表示 LARS 更快；双子图并列布局，左图带 Gram，右图不带 Gram，配色条水平置底，直观展示 Gram 缓存对两算法的加速差异。
源码路径：benchmarks/bench_plot_omp_lars.py - compute_bench()（17-78行）
def compute_bench(samples_range, features_range):
it = 0
results = dict()
lars = np.empty((len(features_range), len(samples_range))) # LARS无Gram矩阵
lars_gram = lars.copy() # LARS带Gram矩阵
omp = lars.copy() # OMP无Gram矩阵
omp_gram = lars.copy() # OMP带Gram矩阵
max_it = len(samples_range) * len(samples_range)
for i_s, n_samples in enumerate(samples_range):
for i_f, n_features in enumerate(features_range):
it += 1
n_informative = n_features // 10
print("====================")
print("Iteration %03d of %03d" % (it, max_it))
print("====================")
dataset_kwargs = {
"n_samples": 1, # 信号长度固定为1
"n_components": n_features, # 字典原子数
"n_features": n_samples, # 转置后样本数=特征数
"n_nonzero_coefs": n_informative, # 10%非零系数
"random_state": 0,
}
print("n_samples: %d" % n_samples)
print("n_features: %d" % n_features)
y, X, _ = make_sparse_coded_signal(**dataset_kwargs)
X = np.asfortranarray(X.T) # 转置并Fortran连续
gc.collect()
print("benchmarking lars_path (with Gram):", end="")
sys.stdout.flush()
tstart = time()
G = np.dot(X.T, X) # 预计算Gram矩阵
Xy = np.dot(X.T, y)
lars_path_gram(Xy=Xy, Gram=G, n_samples=y.size, max_iter=n_informative)
delta = time() - tstart
print("%0.3fs" % delta)
lars_gram[i_f, i_s] = delta # 记录到矩阵
gc.collect()
print("benchmarking lars_path (without Gram):", end="")
sys.stdout.flush()
tstart = time()
lars_path(X, y, Gram=None, max_iter=n_informative) # 原始LARS
delta = time() - tstart
print("%0.3fs" % delta)
lars[i_f, i_s] = delta
gc.collect()
print("benchmarking orthogonal_mp (with Gram):", end="")
sys.stdout.flush()
tstart = time()
orthogonal_mp(X, y, precompute=True, n_nonzero_coefs=n_informative) # OMP带Gram
delta = time() - tstart
print("%0.3fs" % delta)
omp_gram[i_f, i_s] = delta
gc.collect()
print("benchmarking orthogonal_mp (without Gram):", end="")
sys.stdout.flush()
tstart = time()
orthogonal_mp(X, y, precompute=False, n_nonzero_coefs=n_informative) # OMP不带Gram
delta = time() - tstart
print("%0.3fs" % delta)
omp[i_f, i_s] = delta
results["time(LARS) / time(OMP)\n (w/ Gram)"] = lars_gram / omp_gram # 计算比率矩阵
results["time(LARS) / time(OMP)\n (w/o Gram)"] = lars / omp
return results
这段代码是 OMP vs LARS 效率对比的核心计算函数。它在样本数/特征数同步扫描的方形网格上，对四种算法组合计时后直接计算耗时比率。关键设计：make_sparse_coded_signal 中设置 n_samples=1, n_components=n_features, n_features=n_samples，通过转置使最终矩阵对称，方便比率热力图呈现对角线几何特征。
源码路径：benchmarks/bench_plot_omp_lars.py - __main__（80-113行）
if __name__ == "__main__":
samples_range = np.linspace(1000, 5000, 5).astype(int) # 5个样本点
features_range = np.linspace(1000, 5000, 5).astype(int) # 5个特征点
results = compute_bench(samples_range, features_range)
max_time = max(np.max(t) for t in results.values())
import matplotlib.pyplot as plt
fig = plt.figure("scikit-learn OMP vs. LARS benchmark results")
for i, (label, timings) in enumerate(sorted(results.items())):
ax = fig.add_subplot(1, 2, i + 1) # 双子图并列
vmax = max(1 - timings.min(), -1 + timings.max()) # 以1为中心对称
plt.matshow(timings, fignum=False, vmin=1 - vmax, vmax=1 + vmax) # 比率热力图
ax.set_xticklabels([""] + [str(each) for each in samples_range])
ax.set_yticklabels([""] + [str(each) for each in features_range])
plt.xlabel("n_samples")
plt.ylabel("n_features")
plt.title(label)
plt.subplots_adjust(0.1, 0.08, 0.96, 0.98, 0.4, 0.63)
ax = plt.axes([0.1, 0.08, 0.8, 0.06]) # 水平colorbar
plt.colorbar(cax=ax, orientation="horizontal")
plt.show()
这段代码是 OMP vs LARS 比率热力图的可视化入口。它将 compute_bench 返回的两个比率矩阵（带/不带 Gram）以 matshow 绘制为双子图热力图。关键设计：vmin=1-vmax, vmax=1+vmax 使配色以 1 为中心对称——比率 > 1 表示 OMP 更快（红色），比率 < 1 表示 LARS 更快（蓝色），直观展示两个算法在不同维度组合下的相对优势区域。
下面用 Mermaid 图展示 OMP vs LARS 的完整算法决策流：
84.11 设计中的取舍
为什么 bench_glm.py 与 bench_lasso.py 选用不同的计时函数？ 这源自精度与可读性的权衡。datetime.now() 返回带微秒精度的 datetime 对象，适合长时间运行的累计计时；time.time() 返回浮点秒数，适合微基准中的瞬时计时。本章中 GLM 单次实验可能耗时数十秒，datetime.now() 的结果更易读；而 Lasso 单次实验在毫秒到秒级，time.time() 精度更高且开销更低。这一选择并非随意，而是与实验时长量级匹配。
为什么 A-SGD 在 bench_sgd_regression.py 中使用更小的 eta0=0.002 与 power_t=0.05？ 这一组合背后是平均策略对学习率衰减敏感度的理论依据。平均 SGD 通过对历史参数取平均降低方差，但平均操作对学习率衰减更敏感——若衰减过快，平均窗口内的有效学习率过小，会导致收敛停滞。因此 A-SGD 使用更小的初始学习率与更慢的衰减指数，让平均窗口内的学习率保持温和下降，从而平衡收敛速度与稳定性。这是 SGD 平均策略的工程经验值，反映了偏差-方差权衡的实践智慧。
为什么 bench_sparsify.py 中稀疏化能加速推理却不损失精度？ 本质原因在于 L1 正则化训练后的 SGD 模型系数本就是稀疏的——大量权重为零。clf.sparsify() 仅仅是把这种潜在的稀疏结构显式化（ndarray → csr_matrix），并未改变模型本身。预测时 sklearn 自动检测 coef_ 类型并选择对应的矩阵-向量乘积内核——稀疏矩阵-稠密向量乘法跳过零元素的乘加运算，因此在输入也稀疏时实现 ~30% 加速，且数值结果完全一致。这一设计揭示了"数据结构匹配算法特性"在工程优化中的核心价值。
为什么 bench_plot_omp_lars.py 采用比率矩阵而非绝对耗时进行可视化？ 源于一种"消除系统偏差、放大相对差异"的设计哲学。比率矩阵消除了硬件与数据生成的差异，让 OMP vs LARS 的相对优势在不同维度组合下都直观可见。以 1 为中心对称的配色（vmin=1-vmax, vmax=1+vmax）让"谁更快"的语义直接映射到颜色——读者一眼就能定位到 OMP 优势区（红色）与 LARS 优势区（蓝色），无需逐一对照数值。这种可视化设计将定量基准转化为定性洞察。
为什么 OMP 的 precompute=True 仅缓存 Gram 对角而非完整 Gram？ 根源在于两算法的查询模式差异。OMP 算法每步只需计算残差与各原子的相关性，Gram 对角元素（每个原子的自相关）足以避免重复计算 O(n_samples) 的内积；而完整 Gram 矩阵在 OMP 中仅用于正交化投影的少量查询，缓存它收益有限。这与 LARS 算法形成鲜明对比——LARS 每步需要 O(n_features) 次完整 Gram 查询来更新 Cholesky 分解，预计算完整 Gram 能省去 O(n_samples*n_features^2) 的冗余内积。这就是热力图中 Gram 矩阵对 LARS 加速明显、对 OMP 加速有限的几何原因，也体现了"按算法特性定制缓存策略"的设计原则。
84.12 动手练习
-
分析基准脚本的实验设计缺陷与改进
-
阅读
benchmarks/bench_glm.py与benchmarks/bench_lasso.py的数据生成与计时逻辑：-
为什么
bench_glm.py使用datetime.now()而bench_lasso.py使用time.time()？精度差异对微秒级操作有何影响？ -
bench_lasso.py中gc.collect()放在每次fit前，而bench_glm.py无 GC 控制，这会如何影响结果可比性？ -
两脚本均未设置
np.random.seed，如何量化随机数据波动对实验结论的置信区间？ -
设计改进方案：统一计时器、固定随机种子、增加预热轮次、多次重复取中位数
-
-
-
复现并扩展 Lasso 路径基准的 Gram 矩阵效应分析
-
运行
benchmarks/bench_plot_lasso_path.py并观察 3D 曲面：-
解释为什么
lars_path带 Gram 矩阵在大特征/小样本时加速明显，而在大样本/小特征时可能减速？ -
对比
lasso_path(precompute=True)与lars_path_gram的数值稳定性差异（提示：Cholesky 分解条件数） -
修改脚本：增加
n_informative从 10% 到 50% 的扫描，观察稀疏度对四算法相对性能的影响 -
将
make_regression替换为真实数据集（如fetch_california_housing），验证合成数据结论的泛化性
-
-
-
设计 SGD 回归超参数敏感性基准实验
-
基于
benchmarks/bench_sgd_regression.py的 SGD 超参配置：-
为什么 A-SGD 使用更小的
eta0=0.002和power_t=0.05？推导平均策略对学习率衰减敏感度的理论依据 -
alpha=alpha/n_train的缩放是否符合 SGD 目标函数1/n * loss + alpha * penalty的数学定义？验证不缩放时的收敛行为 -
设计网格搜索实验：
learning_rate ∈ {constant, optimal, invscaling, adaptive}×eta0 ∈ {0.001, 0.01, 0.1}×power_t ∈ {0.1, 0.25, 0.5} -
增加
early_stopping=True与validation_fraction=0.1，对比固定max_iter=1000的泛化性能差异
-
-
-
量化模型稀疏化对端到端推理管线的加速比
-
深入
benchmarks/bench_sparsify.py的微基准设计：-
benchmark_dense_predict与benchmark_sparse_predict的循环次数(300)是否足够消除 JIT/缓存预热噪声？如何用timeit.repeat统计分布？ -
稀疏化前后
clf.coef_的内存占用对比（sys.getsizeof+coef_.data.nbytes），计算压缩率 -
扩展实验：将
X_test稀疏度从 2.7% 扫描至 0.1%~50%，绘制稀疏度-加速比曲线，寻找盈亏平衡点 -
对比
SGDRegressor(penalty='l1')与Lasso(alpha=...)稀疏化后的推理速度，解释坐标下降与 SGD 稀疏模式差异
-
-
-
解读 OMP 与 LARS 热力图背后的算法复杂度理论
-
分析
benchmarks/bench_plot_omp_lars.py的比率热力图：-
推导 OMP 与 LARS 的时间复杂度：OMP 每步 O(n_features) 选原子 + O(k^2) 更新投影；LARS 每步 O(n_features) 计算相关 + O(k^3) 更新方向
-
解释热力图中对角线附近(样本≈特征)比率接近 1，而非对称区域显著偏离的几何原因
-
Gram 矩阵预计算将 LARS 从 O(n_samples*n_features^2) 降至 O(n_features^3)，为何对 OMP 加速有限？（提示：OMP 的
precompute=True仅缓存 Gram 对角） -
修改
n_nonzero_coefs从n_features//10扫描至n_features//2，观察稀疏度对比率矩阵拓扑结构的相变
-
-
84.13 本章小结
这一章中我们深入剖析了 scikit-learn 中针对线性模型与稀疏编码的独立基准脚本，理解它们如何从速度、精度与内存多维度对比不同算法的工程表现。其次，我们解构了 bench_glm.py 通过随机方阵在 500-20000 维度上的耗时增长曲线；接着，我们分析了 bench_lasso.py 双场景下坐标下降与最小角回归的效率分水岭；然后，我们探索了 bench_plot_lasso_path.py 在样本-特征二维网格上的 3D 耗时地形图；之后，我们对比了 bench_glmnet.py 跨语言基准测试中的速度、RMSE 与系数一致性；进一步，我们审视了 bench_sgd_regression.py 在样本/特征/学习率三维空间下的 SGD 家族权衡；再后，我们验证了 bench_sparsify.py 中稀疏化推理 ~30% 加速且 R² 无损的工程结论；最后，我们解读了 bench_plot_omp_lars.py 的比率热力图如何揭示 Gram 矩阵对 LARS 与 OMP 的非对称加速效应。
下表对本章涉及的七个基准脚本的核心定位做一个系统梳理：
| 概念 | 解释 |
|------|------|
| bench_glm.py | Ridge/OLS/LassoLars 随维度线性增长的训练时间对比基线 |
| bench_lasso.py | Lasso(CD) 与 LassoLars(LARS) 在样本/特征维度扫描下的计算效率分水岭 |
| bench_plot_lasso_path.py | 四种正则化路径算法在样本-特征二维网格上的 3D 耗时地形图 |
| bench_glmnet.py | sklearn Lasso 与 glmnet Lasso 跨生态系统的速度/精度/系数一致性三维对决 |
| bench_sgd_regression.py | 四种线性求解器在样本量/特征数/学习率策略三维空间的 RMSE 与时间权衡全景 |
| bench_sparsify.py | L1 稀疏化后模型在稀疏输入上推理加速 ~30% 且 R² 无损的工程验证 |
| bench_plot_omp_lars.py | OMP 与 LARS 带/不带 Gram 矩阵的耗时比率热力图，揭示预计算对两算法的非对称加速效应 |
下一章中，我们将学习独立基准脚本：分类器与数据专题 —— 跨算法"大阅兵"。该章将涵盖 MNIST、Covertype、20 Newsgroups 等真实数据集上多算法的训练/推理耗时与错误率对比，以及决策树、在线 OCSVM 等专题基准，揭示 scikit-learn 在大规模真实数据上的性能画像。
第 85 章 —— 独立基准脚本：分类器与数据专题 —— 跨算法"大阅兵"
85.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解 scikit-learn 基准测试脚本的标准设计模式：数据加载、模型配置、计时评估、结果输出
-
掌握 MNIST、Covertype、20 Newsgroups 等经典数据集在基准测试中的预处理差异（归一化、稀疏格式保持、分类特征处理）
-
学会设计可复现的基准实验：统一随机种子注入、并行度控制、内存映射缓存复用
-
理解核近似（Nystroem/RBFSampler）如何将非线性 SVM 转化为线性问题以适配在线学习器（SGDOneClassSVM）
-
掌握决策树扩展性分析的双维度实验设计（样本量 vs 维度）及原始计时绘图方法
-
了解异常检测基准中 ROC/AUC 多种子平均、插值对齐与对数坐标可视化的工程实践
85.2 生活类比
想象基准测试是一场"算法奥林匹克大赛"。在这场盛会中，每一段代码都扮演着特定的角色。数据集 (MNIST / Covertype / 20 Newsgroups) 就是标准赛道——MNIST 像短跑跑道，特征稠密且均匀；Covertype 像山地越野赛道，有 58 万样本的漫长距离；20 Newsgroups 像障碍跑赛道，稀疏矩阵密布。三条赛道地形差异巨大，选手必须配备不同的装备。预处理就像运动员的热身与装备，其中 joblib.Memory 是更衣室储物柜——一次下载、多次复用，绝不重复换装；mmap_mode='r' 是零拷贝启动器，起跑时不必把整个数据集塞进内存。ESTIMATORS 字典就是参赛名单：朴素贝叶斯选手朴素而快速、树选手稳扎稳打、集成接力队人多势众、核近似+线性 SVM 是穿着"降维打击"装备的黑马。随机种子与 n_jobs 注入则是统一比赛规则——同一起跑线（seed）、同一风速辅助（并行度），保证公平竞技。计时器（time / datetime）是电子计时牌，分别记录起跑反应时间（fit）和冲刺速度（predict）。错误率 / 准确率 / AUC 是裁判评分标准，不同赛道用不同计分卡。bench_online_ocsvm 中的核近似最让人眼前一亮：Nystroem 把复杂的非线性赛道（RBF 核）铺平成线性高速公路，让 SGD 在线选手能跑赢传统 LibSVM 选手，这正是效率革命的精髓。最后，可视化柱状图与折线图就是赛后数据复盘大屏，对数坐标是放大镜，能看清跨数量级的时间差异。整场大赛的每一行代码，都是为了让选手在公平、透明、可复现的条件下同场竞技。
85.3 源码地图
benchmarks/bench_mnist.py
├── load_data() # 数据加载、归一化、分割、缓存 (joblib.Memory + mmap)
├── ESTIMATORS # 10 种分类器字典（含 Pipeline 组合）
├── __main__ # 参数解析、随机种子/n_jobs 注入、计时循环、结果排序输出
benchmarks/bench_covertype.py
├── load_data() # fetch_covtype、二值化目标、前 10 列标准化、固定分割、缓存
├── ESTIMATORS # 8 种模型：线性、树、集成、朴素贝叶斯
├── __main__ # 同 MNIST 流程，输出类别分布统计
benchmarks/bench_20newsgroups.py
├── ESTIMATORS # 文本分类典型模型（NB、LR、树集成、Adaboost、基线）
├── __main__ # 稀疏矩阵加载(csc/csr)、密度统计、准确率评估、无缓存/无 n_jobs 注入
benchmarks/bench_tree.py
├── bench_scikit_tree_classifier() # 分类树拟合+预测计时 (datetime + gc.collect)
├── bench_scikit_tree_regressor() # 回归树拟合+预测计时
├── __main__ # 双循环：样本量扩展(10k-100k) + 维度扩展(500-5k)，matplotlib 双子图绘制
benchmarks/bench_online_ocsvm.py
├── print_outlier_ratio() # 异常比例统计
├── autolabel_auc() # AUC 柱状图数值标注
├── autolabel_time() # 时间柱状图数值标注 (对数坐标)
├── __main__ # 5 数据集加载/预处理(LabelBinarizer/StandardScaler)、LibSVM vs Online(SGD+Nystroem)对比
│ ├── 多种子 ROC 曲线插值平均 (interp1d -> 统一 FPR 轴)
│ ├── AUC/训练/预测时间记录
│ └── 三张对数/线性柱状图可视化 (autolabel 标注)
85.4 MNIST 分类器基准 —— 手写数字上的"算法嘉年华"
85.4.1 数据加载与缓存机制
MNIST 基准以 joblib.Memory 配合 mmap_mode='r' 起手，实现了"一次下载、磁盘缓存、内存映射复用"的完整数据生命周期。让我们先看缓存与加载函数。
源码路径：benchmarks/bench_mnist.py - load_data()（48-73行）
memory = Memory(
os.path.join(get_data_home(), "mnist_benchmark_data"),
mmap_mode="r" # 内存映射只读模式，多个进程可零拷贝共享
)
@memory.cache # 函数级缓存装饰器：相同输入直接读磁盘，不重跑函数
def load_data(dtype=np.float32, order="F"):
"""Load the data, then cache and memmap the train/test split"""
# ① 打印加载提示，便于观察进度
print("Loading dataset...")
# ② fetch_openml 获取 MNIST 数据集
data = fetch_openml("mnist_784", as_frame=True)
# ③ check_array 统一数组类型与内存布局
X = check_array(data["data"], dtype=dtype, order=order)
y = data["target"]
# ④ 像素归一化：原始 0-255 整数除以 255 映射到 [0,1] 浮点
X = X / 255
# ⑤ 创建训练-测试分割（遵循 [Joachims, 2006] 文献设定）
print("Creating train-test split...")
n_train = 60000
X_train = X[:n_train]
y_train = y[:n_train]
X_test = X[n_train:]
y_test = y[n_train:]
# ⑥ 返回四个 NumPy 数组，joblib 将其以 mmap 形式持久化
return X_train, X_test, y_train, y_test
这段代码定义/展示了MNIST 数据的"零成本复用"加载管线：通过 Memory 缓存装饰器与 mmap_mode='r' 只读内存映射，首次运行下载并保存到磁盘缓存，后续运行直接内存映射，避免重复加载的 I/O 与内存开销。
mmap_mode='r' 是一种零拷贝机制：数据保留在磁盘上，通过虚拟内存按需分页加载，多个进程可以同时读取同一份数据而互不干扰。这对基准测试至关重要——我们可能需要跑多次实验，每次都从头读取 70,000 张图片会浪费大量时间。
85.4.2 ESTIMATORS 字典设计
接下来看参赛名单。MNIST 的 ESTIMATORS 字典集成了 10 种风格迥异的分类器。
源码路径：benchmarks/bench_mnist.py - ESTIMATORS（76-106行）
ESTIMATORS = {
# ① 基线模型：始终预测最频繁类
"dummy": DummyClassifier(),
# ② 单棵决策树
"CART": DecisionTreeClassifier(),
# ③ 极端随机树集成
"ExtraTrees": ExtraTreesClassifier(),
# ④ 经典随机森林
"RandomForest": RandomForestClassifier(),
# ⑤ Nystroem 核近似 + 线性 SVM 组合
# 顺序：核近似必须在 LinearSVC 之前——LinearSVC 只能处理线性可分问题，
# Nystroem 先把原始 784 维像素映射到 1000 维近似核特征空间，
# 这样 LinearSVC 才能在近似空间里"近似地"做非线性分类
"Nystroem-SVM": make_pipeline(
Nystroem(gamma=0.015, n_components=1000), # 1000 维近似特征
LinearSVC(C=100)
),
# ⑥ RBFSampler 随机傅里叶特征 + 线性 SVM
# 与 Nystroem 同理：RBFSampler 先做显式核映射，LinearSVC 跟在后面做线性分类
"SampledRBF-SVM": make_pipeline(
RBFSampler(gamma=0.015, n_components=1000), LinearSVC(C=100)
),
# ⑦ SAG 求解器逻辑回归（适合大样本）
"LogisticRegression-SAG": LogisticRegression(
solver="sag", tol=1e-1, C=1e4),
# ⑧ SAGA 求解器逻辑回归（SAG 的方差缩减改进版）
"LogisticRegression-SAGA": LogisticRegression(
solver="saga", tol=1e-1, C=1e4),
# ⑨ SGD 优化器的多层感知机
"MultilayerPerceptron": MLPClassifier(
hidden_layer_sizes=(100, 100), max_iter=400, alpha=1e-4,
solver="sgd", learning_rate_init=0.2, momentum=0.9,
verbose=1, tol=1e-4, random_state=1,
),
# ⑩ Adam 优化器的多层感知机
"MLP-adam": MLPClassifier(
hidden_layer_sizes=(100, 100), max_iter=400, alpha=1e-4,
solver="adam", learning_rate_init=0.001,
verbose=1, tol=1e-4, random_state=1,
),
}
这段代码定义/展示了MNIST 基准的 10 人参赛阵容：覆盖了从基线（dummy）、单棵树（CART）、集成（ExtraTrees/RandomForest）、核近似+线性 SVM（Nystroem/RBF-Sampler）、大规模线性模型（SAG/SAGA）到深度学习（MLP-SGD/Adam）的完整算法谱系。其中 Nystroem-SVM 与 SampledRBF-SVM 用 make_pipeline 将非线性核转化为线性可解问题，是 MNIST 高维数据上的重要对比点。
注意 Nystroem(gamma=0.015, n_components=1000) 与 RBFSampler(gamma=0.015, n_components=1000) 两者都把 784 维原始像素映射到 1000 维近似特征空间，但原理不同：Nystroem 基于 Nyström 低秩近似，RBFSampler 基于随机傅里叶特征。这种"非线性 → 线性"的桥梁让 LinearSVC 能处理原本需要核 SVM 的问题。核近似必须放在 LinearSVC 之前：因为 LinearSVC 只能画线性决策边界，必须先由 Nystroem/RBFSampler 在特征空间中"展开"非线性结构，LinearSVC 才能在展开后的近似空间完成分类；如果调换顺序，LinearSVC 会在原始像素空间直接画超平面，丢失非线性能力。
85.4.3 超参数统一注入与计时协议
最后看主流程，理解随机种子注入、n_jobs 控制与计时协议。
源码路径：benchmarks/bench_mnist.py - __main__（108-185行）
if __name__ == "__main__":
parser = argparse.ArgumentParser()
# ① 选择参评分类器子集
parser.add_argument(
"--classifiers", nargs="+", choices=ESTIMATORS, type=str,
default=["ExtraTrees", "Nystroem-SVM"],
)
# ② 并行度控制
parser.add_argument("--n-jobs", nargs="?", default=1, type=int)
# ③ 数据内存布局选择（Fortran / C）
parser.add_argument("--order", nargs="?", default="C", type=str,
choices=["F", "C"])
# ④ 随机种子统一注入
parser.add_argument("--random-seed", nargs="?", default=0, type=int)
args = vars(parser.parse_args())
# ⑤ 加载数据（首次会缓存，后续直接 mmap 复用）
X_train, X_test, y_train, y_test = load_data(order=args["order"])
# ⑥ 遍历每个分类器：注入参数 -> 计时训练 -> 计时预测 -> 计算错误率
error, train_time, test_time = {}, {}, {}
for name in sorted(args["classifiers"]):
print("Training %s ... " % name, end="")
estimator = ESTIMATORS[name]
estimator_params = estimator.get_params()
# ⑦ 自动注入 random_state（遍历所有以 random_state 结尾的参数名）
estimator.set_params(**{
p: args["random_seed"]
for p in estimator_params if p.endswith("random_state")
})
# ⑧ 仅当模型支持 n_jobs 时才注入并行度
if "n_jobs" in estimator_params:
estimator.set_params(n_jobs=args["n_jobs"])
# ⑨ 分别记录 fit 与 predict 耗时
time_start = time()
estimator.fit(X_train, y_train)
train_time[name] = time() - time_start
time_start = time()
y_pred = estimator.predict(X_test)
test_time[name] = time() - time_start
# ⑩ 使用 zero_one_loss 计算错误率作为评分标准
error[name] = zero_one_loss(y_test, y_pred)
print("done")
# ⑪ 按错误率升序输出格式化结果表
print("Classification performance:")
print("===========================")
print("{0: <24} {1: >10} {2: >11} {3: >12}".format(
"Classifier ", "train-time", "test-time", "error-rate"))
print("-" * 60)
for name in sorted(args["classifiers"], key=error.get):
print("{0: <23} {1: >10.2f}s {2: >10.2f}s {3: >12.4f}".format(
name, train_time[name], test_time[name], error[name]))
这段代码定义/展示了MNIST 基准的"裁判台"：参数解析 → 数据加载 → 遍历每个分类器注入随机种子与并行度 → 分段计时训练/预测 → 用 zero_one_loss 评估错误率 → 按错误率升序输出排名表。
随机种子注入的精妙之处在于自动遍历 get_params() 返回的所有参数名，凡是后缀为 random_state 的参数都被统一赋值。这种通用模式适用于任何 sklearn 估计器，包括 Pipeline 中的嵌套组件（虽然此处只检查顶层参数名，不会探测 Pipeline 内部 step__param 形式的嵌套参数）。
时间复杂度与基准结果的对比：典型输出中 MLP-adam 训练 53s、错误率 0.0224；Nystroem-SVM 训练 113s、错误率 0.0228；dummy 训练 0s、错误率 0.8973。这直观展示了"精度与速度的权衡"——没有任何单一模型在两个维度上同时最优。
下面用 Mermaid 图描绘整个 MNIST 基准的执行流程：
85.5 Covertype 森林覆盖类型基准 —— 多模型"山地越野赛"
85.5.1 大规模数据集加载与领域知识驱动预处理
Covertype 数据集有 58 万样本、54 维特征，挑战在于：前 10 列是数值型海拔/距离等连续特征，后 44 列是"植被类型"的二值分类特征。我们必须只对前 10 列标准化，否则会破坏分类特征的语义。
源码路径：benchmarks/bench_covertype.py - load_data()（56-85行）
@memory.cache
def load_data(dtype=np.float32, order="C", random_state=13):
"""Load the data, then cache and memmap the train/test split"""
print("Loading dataset...")
# ① fetch_covtype 获取数据集，开启 shuffle 与随机种子保证可复现
data = fetch_covtype(
download_if_missing=True, shuffle=True, random_state=random_state
)
X = check_array(data["data"], dtype=dtype, order=order)
# ② 目标二值化：预测"是否为类别 1 (Spruce/Fir)"，其余归为反类
y = (data["target"] != 1).astype(int)
print("Creating train-test split...")
# ③ 遵循文献 [Joachims, 2006] 的固定分割：522911 训练，其余测试
n_train = 522911
X_train = X[:n_train]
y_train = y[:n_train]
X_test = X[n_train:]
y_test = y[n_train:]
# ④ 关键领域知识：仅对前 10 个数值特征标准化
# Covertype 前 10 维是连续数值（海拔、距离等），后 44 维是 0/1 植被类型二值特征
# 标准化会让数值特征等权重贡献，后 44 维若一起标准化会破坏 0/1 语义
mean = X_train.mean(axis=0) # 全 54 维的均值
std = X_train.std(axis=0) # 全 54 维的标准差
# ⑤ 后 44 个分类特征强制 mean=0, std=1，避免被错误标准化
# 这种"先计算再覆盖"的写法非常巧妙：先全量计算，再对分类特征位重置为 no-op，
# 后续 (X - mean) / std 对后 44 列就退化为 (X - 0) / 1 = X，数学上保持原样
mean[10:] = 0.0
std[10:] = 1.0
# ⑥ 应用标准化：使用训练集的 mean/std 同步处理训练与测试
X_train = (X_train - mean) / std
X_test = (X_test - mean) / std
return X_train, X_test, y_train, y_test
这段代码定义/展示了Covertype 数据的"地形适配"预处理：二值化目标 → 固定分割 → 仅标准化前 10 列数值特征而保持后 44 列分类特征原样，体现了"领域知识驱动预处理"的设计哲学。领域知识驱动预处理：标准化是为了让连续型特征处于同一量级，避免某些特征因数值大而主导距离/梯度计算；二值特征本身取值就是 0/1，标准化会把它变成 (-1, 1) 区间，反而破坏其物理含义。
mean[10:] = 0.0; std[10:] = 1.0 是非常巧妙的写法：先计算所有特征的统计量，再把后 44 列的 mean/std 强制设为 0/1。由于 (X - 0) / 1 = X，这些特征在数学上保持原样，但通过标准化器接口保持了数据流的一致性。
85.5.2 ESTIMATORS 字典与计时流程
Covertype 的 ESTIMATORS 字典选出了 8 种适合大规模表格数据的模型：既有线性组，也有树组和概率组，覆盖面广但又突出"线性模型可扩展性"这一重点。
源码路径：benchmarks/bench_covertype.py - ESTIMATORS（88-97行）
ESTIMATORS = {
# ① 梯度提升树（250 棵树）
"GBRT": GradientBoostingClassifier(n_estimators=250),
# ② 极端随机树（20 棵）
"ExtraTrees": ExtraTreesClassifier(n_estimators=20),
# ③ 随机森林（20 棵）
"RandomForest": RandomForestClassifier(n_estimators=20),
# ④ 单棵决策树（min_samples_split=5 略放宽）
"CART": DecisionTreeClassifier(min_samples_split=5),
# ⑤ 随机梯度下降线性分类
"SGD": SGDClassifier(alpha=0.001),
# ⑥ 高斯朴素贝叶斯
"GaussianNB": GaussianNB(),
# ⑦ Liblinear 坐标下降线性 SVM
"liblinear": LinearSVC(
loss="l2", penalty="l2", C=1000, dual=False, tol=1e-3),
# ⑧ SAG 求解器逻辑回归（仅 2 轮迭代展示收敛速度）
"SAG": LogisticRegression(solver="sag", max_iter=2, C=1000),
}
这段代码定义/展示了Covertype 基准的 8 人参赛阵容：线性组（SGD / liblinear / SAG）、树组（CART / ExtraTrees / RandomForest / GBRT）、概率组（GaussianNB）。重点是线性模型训练速度差异巨大：SGD 仅需约 1 秒即可训练完 58 万样本，Liblinear 需要约 16 秒，SAG 仅 2 轮迭代就展示了不同求解器的收敛特性。
__main__ 流程与 MNIST 基本一致，但增加了类别分布统计：
源码路径：benchmarks/bench_covertype.py - __main__（106-185行）
if __name__ == "__main__":
parser = argparse.ArgumentParser()
# ① 命令行参数：参评分类器子集、并行度、内存布局、随机种子
parser.add_argument(
"--classifiers", nargs="+", choices=ESTIMATORS, type=str,
default=["liblinear", "GaussianNB", "SGD", "CART"],
)
parser.add_argument("--n-jobs", nargs="?", default=1, type=int)
parser.add_argument("--order", nargs="?", default="C", type=str,
choices=["F", "C"])
parser.add_argument("--random-seed", nargs="?", default=13, type=int)
args = vars(parser.parse_args())
# ② 加载数据（首次缓存，后续 mmap 复用）
X_train, X_test, y_train, y_test = load_data(
order=args["order"], random_state=args["random_seed"]
)
# ③ 输出类别分布统计：包含特征数、类别数、dtype、pos/neg 数与内存 MB
# pos/neg 让基准报告自带不平衡数据画像——Covertype 类别 1（Spruce/Fir）
# 占比仅约 24%，这种不平衡分布会显著影响模型的训练策略与评估指标
print("%s %d" % ("number of features:".ljust(25), X_train.shape[1]))
print("%s %d" % ("number of classes:".ljust(25), np.unique(y_train).size))
print("%s %s" % ("data type:".ljust(25), X_train.dtype))
# ④ 关键统计：训练样本数 + 正样本数 + 负样本数 + 内存 MB
# np.sum(y_train == 1) 统计正类（Spruce/Fir）的样本数
# int(X_train.nbytes / 1e6) 把字节数转换为 MB 单位
print(
"%s %d (pos=%d, neg=%d, size=%dMB)"
% (
"number of train samples:".ljust(25),
X_train.shape[0],
np.sum(y_train == 1),
np.sum(y_train == 0),
int(X_train.nbytes / 1e6),
)
)
print(
"%s %d (pos=%d, neg=%d, size=%dMB)"
% (
"number of test samples:".ljust(25),
X_test.shape[0],
np.sum(y_test == 1),
np.sum(y_test == 0),
int(X_test.nbytes / 1e6),
)
)
# ⑤ 遍历分类器：注入 random_state + n_jobs，分段计时，错误率评估
error, train_time, test_time = {}, {}, {}
for name in sorted(args["classifiers"]):
print("Training %s ... " % name, end="")
estimator = ESTIMATORS[name]
estimator_params = estimator.get_params()
estimator.set_params(**{
p: args["random_seed"]
for p in estimator_params if p.endswith("random_state")
})
if "n_jobs" in estimator_params:
estimator.set_params(n_jobs=args["n_jobs"])
# ⑥ 分段计时：fit 与 predict 各自独立记录
time_start = time()
estimator.fit(X_train, y_train)
train_time[name] = time() - time_start
time_start = time()
y_pred = estimator.predict(X_test)
test_time[name] = time() - time_start
error[name] = zero_one_loss(y_test, y_pred)
print("done")
# ⑦ 按错误率升序输出格式化结果表
print("Classification performance:")
print("===========================")
print("%s %s %s %s" % ("Classifier ", "train-time", "test-time", "error-rate"))
print("-" * 44)
for name in sorted(args["classifiers"], key=error.get):
print("%s %s %s %s" % (
name.ljust(12),
("%.4fs" % train_time[name]).center(10),
("%.4fs" % test_time[name]).center(10),
("%.4f" % error[name]).center(10),
))
这段代码定义/展示了Covertype 基准的"越野赛裁判台"：与 MNIST 主流程几乎一致，差异点在于默认分类器子集（liblinear/GaussianNB/SGD/CART 四种典型大规模表格学习模型）、默认随机种子为 13、以及在训练/测试样本统计中额外输出 pos（正样本数）与 neg（负样本数），让基准报告自带不平衡数据画像。Covertype 类别 1（Spruce/Fir）占比仅约 24%，这种不平衡分布会显著影响模型的训练策略与评估指标。
下面用表格对比 MNIST 与 Covertype 主流程的差异，揭示两个基准在数据规模、模型选择与报告信息上的不同侧重：
| 维度 | MNIST | Covertype |
|------|-------|-----------|
| 默认分类器 | ExtraTrees / Nystroem-SVM | liblinear / GaussianNB / SGD / CART |
| 默认随机种子 | 0 | 13 |
| 输出统计字段 | n_samples / dtype / 内存 MB | 增加 pos / neg 计数 |
| 数据规模 | 70k | 580k |
下面用 Mermaid 图描绘 Covertype 基准的整体执行流程：
85.6 20 Newsgroups 文本分类基准 —— 自然语言的"性能试验场"
85.6.1 ESTIMATORS 字典设计
20 Newsgroups 基准的 ESTIMATORS 字典专门为稀疏文本数据挑选模型：朴素贝叶斯（天然适配词频离散分布）、逻辑回归（线性模型在稀疏高维空间收敛快）、决策树集成（验证非线性能力上限）、AdaBoost（提升弱学习器的早期基线）、Dummy（基线对照）。
源码路径：benchmarks/bench_20newsgroups.py - ESTIMATORS（18-27行）
ESTIMATORS = {
# ① 基线：始终预测最频繁类，衡量其他模型的相对优势
"dummy": DummyClassifier(),
# ② 随机森林：限制 max_features='sqrt' 防高维过拟合
"random_forest": RandomForestClassifier(
max_features="sqrt", min_samples_split=10),
# ③ 极端随机树：同 RF 限制，对比随机性增强的影响
"extra_trees": ExtraTreesClassifier(
max_features="sqrt", min_samples_split=10),
# ④ 逻辑回归：线性模型在 100 万+ 词表维度的标准基线
"logistic_regression": LogisticRegression(),
# ⑤ 多项式朴素贝叶斯：稀疏离散特征的天然搭档
"naive_bayes": MultinomialNB(),
# ⑥ AdaBoost：10 棵树弱提升，早期经典集成基线
"adaboost": AdaBoostClassifier(n_estimators=10),
}
这段代码定义/展示了20 Newsgroups 基准的 6 人参赛阵容：dummy 提供下界基线，MultinomialNB 是稀疏词频的经典搭档，LogisticRegression 是线性模型标杆，RandomForest/ExtraTrees 验证非线性能力上限（但因高维耗时较长），AdaBoost 用 10 棵弱树提供早期集成基线。整体配置遵循"由简到繁、由稀疏专用到通用"的文本分类模型选择逻辑。
85.6.2 稀疏文本数据的高效加载
20 Newsgroups 基准与众不同：它没有 joblib 缓存，却精心区分了训练集与测试集的稀疏格式。
源码路径：benchmarks/bench_20newsgroups.py - __main__（35-112行）
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument(
"-e", "--estimators", nargs="+", required=True, choices=ESTIMATORS
)
args = vars(parser.parse_args())
data_train = fetch_20newsgroups_vectorized(subset="train")
data_test = fetch_20newsgroups_vectorized(subset="test")
# ① 训练集用 CSC 格式：列切片高效，适合按特征操作的拟合阶段
# CSC（Compressed Sparse Column）按列存储，访问单列或切片列时只需读取对应索引，
# 在 fit 阶段如果需要按词表维度做统计或特征选择，CSC 高效得多
X_train = check_array(
data_train.data, dtype=np.float32, accept_sparse="csc")
# ② 测试集用 CSR 格式：行切片高效，适合按样本操作的预测阶段
# CSR（Compressed Sparse Row）按行存储，访问单行或切片行时只需读取对应索引，
# 在 predict 阶段需要逐文档（逐行）查询，CSR 高效得多
X_test = check_array(
data_test.data, dtype=np.float32, accept_sparse="csr")
y_train = data_train.target
y_test = data_test.target
print("20 newsgroups")
print("=============")
# ③ 打印稀疏统计：形状、格式、密度（nnz / 总元素数）、dtype
# X_train.nnz / np.prod(X_train.shape) 即密度 = 非零元素数 / 总元素数
# 20 Newsgroups TF-IDF 矩阵的密度通常只有 0.1% 左右，意味着 99.9% 的元素都是 0
# 如果不小心转成稠密矩阵，内存占用会爆炸式增长
print(f"X_train.shape = {X_train.shape}")
print(f"X_train.format = {X_train.format}")
print(f"X_train.dtype = {X_train.dtype}")
print(f"X_train density = {X_train.nnz / np.prod(X_train.shape)}")
print(f"y_train {y_train.shape}")
print(f"X_test {X_test.shape}")
print(f"X_test.format = {X_test.format}")
print(f"X_test.dtype = {X_test.dtype}")
print(f"y_test {y_test.shape}")
print()
print("Classifier Training")
print("===================")
accuracy, train_time, test_time = {}, {}, {}
for name in sorted(args["estimators"]):
clf = ESTIMATORS[name]
# ④ 简化版种子注入：仅设 random_state（不遍历 get_params）
# 用 try-except 容忍部分模型无 random_state 参数
try:
clf.set_params(random_state=0)
except (TypeError, ValueError):
pass
print("Training %s ... " % name, end="")
t0 = time()
clf.fit(X_train, y_train)
train_time[name] = time() - t0
t0 = time()
y_pred = clf.predict(X_test)
test_time[name] = time() - t0
# ⑤ 使用 accuracy_score 而非 zero_one_loss
accuracy[name] = accuracy_score(y_test, y_pred)
print("done")
这段代码定义/展示了稀疏文本数据的"地形优化"加载策略：训练集 CSC（按列切高效，利于词表构建）、测试集 CSR（按行切高效，利于逐文档预测），accept_sparse 参数保持稀疏格式不转稠密。
accept_sparse 参数告诉 check_array 接受并保持稀疏矩阵格式——若不指定，默认会尝试转稠密，导致 1.8 万文档 × 100 万词表的 TF-IDF 矩阵在内存中直接爆掉；指定后，矩阵在后续拟合/预测全流程都保持稀疏表示，内存占用与计算量都按 nnz（非零元素数）而非总元素数计算。
X_train.nnz / np.prod(X_train.shape) 是密度（density）指标：nnz 是非零元素个数，除以总元素数得到非零比例。20 Newsgroups TF-IDF 矩阵的密度通常只有 0.1% 左右，意味着 99.9% 的元素都是 0。如果不小心转成稠密矩阵，内存占用会爆炸式增长。
为何不需要缓存？20newsgroups 数据集本身约 60MB，向量化过程只需几秒，每次重新生成完全可接受。相比 MNIST 的 70,000 张原始图像需要先 fetch_openml 再归一化再分割的复杂流程，这里的开销小得多。
85.6.3 简化的评估流程
源码路径：benchmarks/bench_20newsgroups.py - __main__（训练循环）
为何放弃统一注入 n_jobs？原因在于 20 Newsgroups 的模型组合中，AdaBoostClassifier 等模型的 base_estimator 是决策树而非直接接受 n_jobs 的并行接口；同时 MultinomialNB 不支持并行。统一注入 n_jobs 会破坏部分模型的实例化，因此采用"按模型特性分别处理"的简化策略。
accuracy_score 与 zero_one_loss 的关系很简单：accuracy = 1 - zero_one_loss，因为 zero_one_loss 计算的是错误率。两者本质相同，只是表述角度不同——准确率越高越好，错误率越低越好。
下面用 Mermaid 图描绘 20 Newsgroups 基准的执行流程：
85.7 决策树基准 —— 样本与维度的"成长烦恼"
85.7.1 决策树回归器的计时函数
bench_tree.py 中 bench_scikit_tree_regressor() 与 bench_scikit_tree_classifier() 几乎完全对称，唯一区别是把分类器换成 DecisionTreeRegressor。
源码路径：benchmarks/bench_tree.py - bench_scikit_tree_regressor()（39-53行）
def bench_scikit_tree_regressor(X, Y):
"""Benchmark with scikit-learn decision tree regressor"""
# ① 局部导入减少启动开销，回归场景才需要 sklearn.tree
from sklearn.tree import DecisionTreeRegressor
# ② 强制垃圾回收，减少 GC 抖动对微秒级计时的干扰
# Python 默认启用自动 GC，但 GC 暂停会引入不可预测的微秒级延迟
# 计时前显式 gc.collect() 把内存中待回收对象清空，让计时窗口内尽量不发生 GC
gc.collect()
# ③ 使用 datetime 微秒级精度开始计时
tstart = datetime.now()
# ④ 回归树实例化 + 同时拟合与预测（in-sample 评估耗时）
# 与分类器对称：分类 Gini 准则 vs 回归 MSE 准则，底层 _splitter.pyx 不同
clf = DecisionTreeRegressor()
clf.fit(X, Y).predict(X)
# ⑤ 计算 timedelta 转换为浮点秒
delta = datetime.now() - tstart
# ⑥ 追加到全局列表，与分类结果并列绘图
scikit_regressor_results.append(
delta.seconds + delta.microseconds / mu_second)
这段代码定义/展示了决策树回归器的"镜像计时"实现：gc.collect() 减少 GC 抖动、datetime.now() 微秒级精度、同时记录 fit 与 predict 总耗时，结果追加到全局列表。其设计哲学与分类计时函数完全一致，唯一不同点是分类树用 Gini/Entropy 分裂准则、回归树用 MSE（方差缩减）。
回归 vs 分类在底层 _splitter.pyx 的差异主要在分裂准则计算：分类 Gini 是 2 * p * (1-p)（p 为正类比例），回归 MSE 是 (y_left - y_mean_left)² + (y_right - y_mean_right)²。两者计算量相当，理论上回归略快（无需概率归一化），但差异并不显著。
85.7.2 双维度扩展性实验设计
bench_tree.py 是一份非常独特的脚本——它没有命令行参数、没有磁盘缓存、没有 ESTIMATORS 字典，完全专注于算法复杂度的二维分析。
源码路径：benchmarks/bench_tree.py - bench_scikit_tree_classifier()（22-36行）
def bench_scikit_tree_classifier(X, Y):
"""Benchmark with scikit-learn decision tree classifier"""
from sklearn.tree import DecisionTreeClassifier
# ① 强制垃圾回收，减少 GC 抖动对计时的干扰
# 同 regressor：计时窗口内不触发 GC，才能得到稳定微秒级结果
gc.collect()
# ② 使用 datetime 微秒级计时
# datetime.now() 返回带微秒精度的 datetime 对象，差值运算更自然
tstart = datetime.now()
# ③ 同时记录拟合与预测耗时
# 与 regressor 对称：fit 用 Gini/Entropy 准则，底层 _splitter.pyx 路径不同
clf = DecisionTreeClassifier()
clf.fit(X, Y).predict(X)
delta = datetime.now() - tstart
# ④ 将 timedelta 转换为浮点秒
scikit_classifier_results.append(
delta.seconds + delta.microseconds / mu_second)
这段代码定义/展示了决策树分类器的"原始计时"实现：gc.collect() 减少 GC 抖动、datetime.now() 微秒级精度、同时记录 fit 与 predict 总耗时，结果追加到全局列表。
datetime.now() 与 time.time() 的核心差异在于：前者返回带微秒精度的 datetime 对象，后者返回浮点秒。前者更易于做差值运算（delta = datetime.now() - tstart），后者更通用且跨平台一致。在微基准测试中，两者精度相当。
下面看主流程的双循环设计：
源码路径：benchmarks/bench_tree.py - __main__（56-115行）
if __name__ == "__main__":
print("============================================")
print("Warning: this is going to take a looong time")
print("============================================")
n = 10
step = 10000
n_samples = 10000
dim = 10
n_classes = 10
# ① 实验 1：固定维度，变样本量（耐力测试）
for i in range(n):
n_samples += step # 每次增加 10000
X = np.random.randn(n_samples, dim) # 10 维固定
Y = np.random.randint(0, n_classes, (n_samples,))
bench_scikit_tree_classifier(X, Y)
Y = np.random.randn(n_samples)
bench_scikit_tree_regressor(X, Y)
xx = range(0, n * step, step)
plt.figure("scikit-learn tree benchmark results")
plt.subplot(211) # ② 上子图：2x1 网格的第一格
plt.title("Learning with varying number of samples")
plt.plot(xx, scikit_classifier_results, "g-", label="classification")
plt.plot(xx, scikit_regressor_results, "r-", label="regression")
plt.legend(loc="upper left")
plt.xlabel("number of samples")
plt.ylabel("Time (s)")
# ③ 重置全局列表，准备第二轮实验
scikit_classifier_results = []
scikit_regressor_results = []
n = 10
step = 500
start_dim = 500
n_classes = 10
# ④ 实验 2：固定样本量，变维度（爆发力测试）
dim = start_dim
for i in range(0, n):
dim += step
X = np.random.randn(100, dim) # 样本数固定 100
Y = np.random.randint(0, n_classes, (100,))
bench_scikit_tree_classifier(X, Y)
Y = np.random.randn(100)
bench_scikit_tree_regressor(X, Y)
xx = np.arange(start_dim, start_dim + n * step, step)
plt.subplot(212) # ⑤ 下子图：2x1 网格的第二格
plt.title("Learning in high dimensional spaces")
plt.plot(xx, scikit_classifier_results, "g-", label="classification")
plt.plot(xx, scikit_regressor_results, "r-", label="regression")
plt.legend(loc="upper left")
plt.xlabel("number of dimensions")
plt.ylabel("Time (s)")
plt.axis("tight")
plt.show() # ⑥ 一次性弹出双子图
这段代码定义/展示了决策树的双维度扩展性实验：实验 1 固定 10 维、样本量从 10k 步进到 100k（耐力测试）；实验 2 固定 100 样本、维度从 500 步进到 5000（爆发力测试）。两个实验分别隔离样本量与维度对决策树构建时间的影响。
plt.subplot(211) 与 plt.subplot(212) 划分 2 行 1 列网格：第一行展示"样本量 vs 时间"耐力曲线，第二行展示"维度 vs 时间"爆发力曲线，绿线分类、红线回归。
下面用 Mermaid 图描绘决策树基准的双循环结构：
下面用表格总结五个脚本的设计哲学差异，从数据规模、预处理重点、模型组合到输出形式多维度对比：
| 脚本 | 数据规模 | 预处理重点 | 模型组合 | 输出形式 |
|------|---------|-----------|---------|---------|
| bench_mnist.py | 70k 样本 | 像素归一化 + mmap | 10 种（含 Pipeline） | 文本表格 |
| bench_covertype.py | 58 万样本 | 仅前 10 列标准化 | 8 种（线性/树） | 文本表格 |
| bench_20newsgroups.py | 1.8 万文档 | csc/csr 分离 | 6 种（文本专用） | 文本表格 |
| bench_tree.py | 合成数据 | 无（随机生成） | 仅决策树 2 种 | matplotlib 图 |
| bench_online_ocsvm.py | 多源混合 | LabelBinarizer + 标准化 | 2 种 OCSVM | 柱状图 |
85.8 在线 OCSVM 对决 LibSVM —— 异常检测的"效率革命"
85.8.1 异常比例统计辅助函数
bench_online_ocsvm.py 在每个数据集加载完毕后，都会调用 print_outlier_ratio() 输出该数据集的类别分布与异常占比。这是异常检测基准的关键预处理步骤：异常率直接决定后续 OneClassSVM 的 nu 参数选择。
源码路径：benchmarks/bench_online_ocsvm.py - print_outlier_ratio()（35-44行）
def print_outlier_ratio(y):
"""
Helper function to show the distinct value count of element in the target.
Useful indicator for the datasets used in bench_isolation_forest.py.
"""
# ① np.unique 返回排序后的唯一值及其出现次数
uniq, cnt = np.unique(y, return_counts=True)
print("----- Target count values: ")
# ② 逐行打印每个类别及其计数（处理 bytes/str/数字任意类型）
for u, c in zip(uniq, cnt):
print("------ %s -> %d occurrences" % (str(u), c))
# ③ 计算异常比例：最小计数 / 总样本数
# 对于 OCSVM 而言，"少的那一类"即为异常类
print("----- Outlier ratio: %.5f" % (np.min(cnt) / len(y)))
这段代码定义/展示了异常比例的可读化输出：np.unique(..., return_counts=True) 同时获取唯一值与频次，遍历打印每个类别的样本数，最后以最小计数除以总数得到异常率。对于 KDDCup99 这类极不平衡数据集，正常样本占 99% 以上、异常样本仅占 0.5% 左右，nu=0.05 的设定允许 OCSVM 把 5% 的训练数据标记为"远离决策边界"。
注意 str(u) 的设计意图：KDD 的 y 是 bytes 类型（如 b"normal."），而 forestcover 的 y 是整数（0/1），str() 统一了打印格式。
85.8.2 数据集加载与预处理
bench_online_ocsvm.py 涉及 5 个异构数据集，需要根据数据集名选择不同预处理路径。SA、SF 因含字符串型分类特征需 LabelBinarizer 独热编码，http/smtp/forestcover 则直接转 float。
源码路径：benchmarks/bench_online_ocsvm.py - __main__（47-115行）
datasets = ["http", "smtp", "SA", "SF", "forestcover"]
novelty_detection = False # 若为 True，训练集会被过滤为只含正常样本（半监督设置）
random_states = [42]
nu = 0.05
results_libsvm = np.empty((len(datasets), n_axis + 5))
results_online = np.empty((len(datasets), n_axis + 5))
for dat, dataset_name in enumerate(datasets):
print(dataset_name)
# ① KDDCup99 子集加载（http/smtp/SA/SF）
if dataset_name in ["http", "smtp", "SA", "SF"]:
dataset = fetch_kddcup99(
subset=dataset_name, shuffle=False, percent10=False, random_state=88
)
X = dataset.data
y = dataset.target
# ② forestcover 单独处理：只取类别 2 和 4 的样本
if dataset_name == "forestcover":
dataset = fetch_covtype(shuffle=False)
X = dataset.data
y = dataset.target
# 正常数据 = 属性 2，异常数据 = 属性 4
s = (y == 2) + (y == 4)
X = X[s, :]
y = y[s]
y = (y != 2).astype(int)
# ③ SF 子集特殊处理：第 1 列是字符串型分类特征
if dataset_name == "SF":
# 需要先把 X[:, 1] 转字符串再 LabelBinarizer，因为 KDD 的分类特征是 bytes
lb = LabelBinarizer()
x1 = lb.fit_transform(X[:, 1].astype(str))
# 把独热编码列拼回原矩阵，替代原字符串列
X = np.c_[X[:, :1], x1, X[:, 2:]]
y = (y != b"normal.").astype(int)
# ④ SA 子集特殊处理：第 1/2/3 列都是字符串型分类特征
if dataset_name == "SA":
lb = LabelBinarizer()
x1 = lb.fit_transform(X[:, 1].astype(str))
x2 = lb.fit_transform(X[:, 2].astype(str))
x3 = lb.fit_transform(X[:, 3].astype(str))
X = np.c_[X[:, :1], x1, x2, x3, X[:, 4:]]
y = (y != b"normal.").astype(int)
# ⑤ http/smtp 没有分类特征，只需把 y 的 bytes 转 int
if dataset_name in ["http", "smtp"]:
y = (y != b"normal.").astype(int)
print_outlier_ratio(y)
n_samples, n_features = np.shape(X)
if dataset_name == "SA": # LibSVM 在 SA 全量上太慢，只取 1/20 训练
n_samples_train = n_samples // 20
else:
n_samples_train = n_samples // 2
n_samples_test = n_samples - n_samples_train
print("n_train: ", n_samples_train)
print("n_features: ", n_features)
这段代码定义/展示了5 个异构数据集的"按需适配"预处理管线：KDDCup99 子集统一走 fetch_kddcup99，forestcover 单独走 fetch_covtype 并按类别 2/4 切片；SA/SF 子集因含字符串型分类特征必须用 LabelBinarizer 独热编码；http/smtp 只做 y 的 bytes→int 转换。SA 子集因样本量过大只取 1/20 作为训练集以避免 LibSVM 训练过慢。
LabelBinarizer 是 sklearn 中专门处理单列多类别特征的编码器：fit_transform 接收一列字符串/字节，输出对应的独热编码矩阵。SA 子集有 3 列字符串特征、SF 子集有 1 列，所以分别调用 3 次和 1 次 fit_transform，最后用 np.c_ 沿列轴拼接回去。
fetch_covtype(shuffle=False) 后通过 (y == 2) + (y == 4) 创建布尔掩码——这是 numpy 布尔索引的常见用法：先用比较运算得到 bool 数组，再用它过滤 X 和 y。y = (y != 2).astype(int) 把二分类标签 0/1 化：原类别 2（正常）变为 0、类别 4（异常）变为 1。
85.8.3 核近似赋能在线学习
bench_online_ocsvm.py 是本章最复杂的脚本，它不仅要做基准对比，还要展示核近似如何让在线学习成为可能。
源码路径：benchmarks/bench_online_ocsvm.py - __main__（140-170行）
gamma = 1 / n_features # ① OCSVM 的默认 gamma 参数（核宽度倒数）
for random_state in random_states:
X, y = shuffle(X, y, random_state=random_state)
X_train = X[:n_samples_train]
X_test = X[n_samples_train:]
y_train = y[:n_samples_train]
y_test = y[n_samples_train:]
std = StandardScaler()
# ② LibSVM 基线：原生 RBF 核 OCSVM（复杂度 O(n²~n³)）
print("----------- LibSVM OCSVM ------------")
ocsvm = OneClassSVM(kernel="rbf", gamma=gamma, nu=nu)
pipe_libsvm = make_pipeline(std, ocsvm)
tstart = time()
pipe_libsvm.fit(X_train)
fit_time_libsvm += time() - tstart
# ③ 在线方案：Nystroem 核近似 + SGDOneClassSVM（线性复杂度）
print("----------- Online OCSVM ------------")
nystroem = Nystroem(gamma=gamma, random_state=random_state)
online_ocsvm = SGDOneClassSVM(nu=nu, random_state=random_state)
pipe_online = make_pipeline(std, nystroem, online_ocsvm)
tstart = time()
pipe_online.fit(X_train)
fit_time_online += time() - tstart
这段代码定义/展示了核近似 vs 原生核 SVM 的"效率对决"：LibSVM 基线用原生 RBF 核 OCSVM（O(n²~n³) 复杂度），在线方案用 Nystroem 显式核近似把数据映射到 ~1000 维特征空间，再丢给 SGDOneClassSVM 做线性训练（O(n) 复杂度）。两者共用 StandardScaler 预处理，make_pipeline 串联整个流程。
Nystroem 的关键参数 gamma=1/n_features 是 sklearn OneClassSVM 的默认值，它定义了 RBF 核 exp(-gamma * ||x-y||²) 的宽度倒数。在 100 维特征空间中，gamma=0.01 意味着 10 个单位的欧氏距离才会让核值衰减到 exp(-1)≈0.37，这与原 OCSVM 的默认行为保持一致。
85.8.4 ROC 曲线插值平均
多种子实验需要把不同随机种子下的 ROC 曲线合并。直接对 AUC 取平均会丢失曲线形状信息，这里采用插值对齐的更精细做法。
源码路径：benchmarks/bench_online_ocsvm.py - __main__（170-200行）
# 第 85 章 —— ① 统一 FPR 轴：1000 个等距点
n_axis = 1000
x_axis = np.linspace(0, 1, n_axis)
# 第 85 章 —— ② 计算当前种子的 ROC 曲线
tstart = time()
scoring = -pipe_libsvm.decision_function(X_test) # 分数越低越正常
predict_time_libsvm += time() - tstart
fpr_libsvm_, tpr_libsvm_, _ = roc_curve(y_test, scoring)
# 第 85 章 —— ③ 用 interp1d 把当前曲线插值到统一 x_axis
f_libsvm = interp1d(fpr_libsvm_, tpr_libsvm_)
# 第 85 章 —— ④ 累加 TPR（后面再除以种子数取平均）
tpr_libsvm += f_libsvm(x_axis)
# 第 85 章 —— ... online 模型同理 ...
# 第 85 章 —— ⑤ 对所有种子取平均
tpr_libsvm /= len(random_states)
tpr_libsvm[0] = 0.0 # 起点必须为 0（数学约束）
auc_libsvm = auc(x_axis, tpr_libsvm) # 用统一轴计算平均 AUC
这段代码定义/展示了多种子 ROC 曲线的"插值对齐"平均法：先构造统一 FPR 轴（1000 个等距点），用 interp1d 把每条 ROC 曲线插值到该轴上累加 TPR，最后除以种子数取平均。tpr[0] = 0.0 是数学约束（FPR=0 时 TPR 必须为 0），保证曲线在原点开始。
为何不直接平均 AUC 数值？因为 AUC 是一个标量，多种子平均后我们只能得到一个均值，无法得到置信区间或标准差。而插值后的平均 ROC 曲线保留了完整的曲线形状信息——我们可以从中观察到模型在不同误报率区间的行为，比如"在 FPR=0.1 时平均 TPR=0.8"。
85.8.5 AUC 柱状图数值标注
可视化部分由三组函数组成：autolabel_auc 处理 AUC（线性坐标），autolabel_time 处理时间（对数坐标）。两者仅在格式串上略有差异——AUC 用 "%.3f"（保留 3 位小数，因 AUC 在 0-1 之间）、时间用 "%.1f"（保留 1 位小数，因时间跨数量级）。
源码路径：benchmarks/bench_online_ocsvm.py - autolabel_auc()（149-163行）
def autolabel_auc(rects, ax):
"""Attach a text label above each bar displaying its height."""
for rect in rects:
height = rect.get_height()
# ① 在柱顶 1.05 倍高度处添加文本标注（柱顶上方 5%）
ax.text(
rect.get_x() + rect.get_width() / 2.0,
1.05 * height,
# ② 浮点 3 位小数（AUC 范围 0-1，3 位足够区分）
"%.3f" % height,
ha="center", va="bottom",
)
这段代码定义/展示了AUC 柱状图的数值标注：autolabel_auc 与 autolabel_time 设计完全相同，唯一差异是格式串——前者用 "%.3f"（AUC 精度敏感），后者用 "%.1f"（时间跨数量级，整数部分已具区分力）。两者都用 1.05 * height 把文字放在柱顶 5% 上方，避免与柱体重叠。
85.8.6 可视化与数值标注
三张柱状图分别展示 AUC、训练时间、预测时间。其中时间柱状图使用对数坐标，因为 LibSVM 与 Online SVM 跨数量级差异巨大。
源码路径：benchmarks/bench_online_ocsvm.py - autolabel_time()（165-179行）
def autolabel_time(rects, ax):
"""Attach a text label above each bar displaying its height."""
for rect in rects:
height = rect.get_height()
# ① 在柱顶 1.05 倍高度处添加文本标注
ax.text(
rect.get_x() + rect.get_width() / 2.0,
1.05 * height,
"%.1f" % height, # ② 浮点 1 位小数
ha="center", va="bottom",
)
源码路径：benchmarks/bench_online_ocsvm.py - __main__（260-300行）
fig, ax = plt.subplots(figsize=(15, 8))
ax.set_ylabel("Training time (sec) - Log scale")
ax.set_yscale("log") # ① 对数 y 坐标：跨数量级时间差异
rect_libsvm = ax.bar(ind, fit_time_libsvm_all, color="r", width=width)
rect_online = ax.bar(ind + width, fit_time_online_all, color="y", width=width)
ax.legend((rect_libsvm[0], rect_online[0]), ("LibSVM", "Online SVM"))
# 第 85 章 —— ② x 轴标签包含数据集名 + 样本量 + 维度
ax.set_xticks(ind + width / 2)
ax.set_xticklabels(x_tickslabels) # 每个标签形如 "http\n$n=488,398$\n$d=3$"
autolabel_time(rect_libsvm, ax)
autolabel_time(rect_online, ax)
plt.show()
这段代码定义/展示了对数坐标柱状图与数值标注：训练/预测时间跨数量级时，set_yscale("log") 是必要选择；autolabel_time 在每根柱顶部 5% 高度处添加数值文本，让数据可读性大幅提升；x 轴标签组合数据集名、训练样本数、特征维度，一图三信息。
下面用 Mermaid 图描绘异常检测基准的整体流程：
85.9 设计中的取舍
为什么 MNIST/Covertype 使用 joblib.Memory 缓存，而 20 Newsgroups 与 bench_tree 不需要？
joblib.Memory 的核心价值是"避免重复下载与解析的开销"。MNIST 数据下载约几十 MB、解析像素耗时数秒，Covertype 数据更是接近 11 MB 的 CSV；如果每次跑实验都重新下载，会显著拖慢迭代节奏。相反，20 Newsgroups 的 TF-IDF 数据已经预向量化、约 60 MB，加载只需几秒，每次重新生成完全可接受；bench_tree 用 np.random.randn 现场合成数据，根本没有"复用"的需求。缓存不是越多越好，要看重复生成的开销与缓存本身的复杂度是否匹配。
为何 bench_online_ocsvm 选择 Nystroem 而不是直接用 SGDOneClassSVM？
SGDOneClassSVM 本身是线性模型，只能在原始特征空间画线性决策边界，无法处理非线性可分的异常检测问题。现实中 KDDCup99 等异常数据往往需要非线性边界。Nystroem 提供了 RBF 核的显式近似——把数据先映射到 ~1000 维的近似核空间，再丢给线性 SGDOneClassSVM 处理。这是"线性模型的效率 + 非线性模型的表达力"的完美结合，代价是 Nystroem 本身的映射计算与超参数 n_components 的选择。
为何决策树基准使用 datetime.now() 而其他脚本都用 time.time()？
本质上两者精度相当（都是微秒级），区别在于 API 风格：datetime.now() 返回带日期的对象，做差值运算（delta = datetime.now() - tstart）非常自然；time.time() 返回浮点秒，做差值需要手动管理参考点。bench_tree 写成约 2010 年代初的风格，那个时期 datetime 微秒级计时是较常见的写法。对当代基准而言，time.perf_counter() 是更好的选择——它返回单调时钟，不受系统时间调整影响。
85.10 动手练习
-
对比五大基准脚本的数据加载与缓存策略差异
-
哪些脚本使用了
joblib.Memory缓存、哪些没有？背后的依据是什么？（MNIST 与 Covertype 因为下载/解析开销大而启用缓存；20 Newsgroups 数据已预向量化故无需缓存；bench_tree 的合成数据根本没有缓存价值） -
bench_20newsgroups.py为何不需要缓存装饰器？其稀疏矩阵格式（csc/csr）分离策略针对什么场景？（稀疏矩阵格式分离是为了在拟合阶段高效列操作、在预测阶段高效行操作，并避免稀疏转稠密的内存爆炸） -
bench_tree.py的数据来源是什么？为何不适合也不需要缓存？（np.random.randn现场合成数据，不适合也不需要缓存） -
bench_online_ocsvm.py处理哪 5 类数据集？哪些需要LabelBinarizer独热编码预处理，为什么？（http/smtp/SA/SF/forestcover 共 5 个；KDDCup99 的 SA/SF 子集因含字符串型分类特征而需要LabelBinarizer独热编码）
-
-
分析 ESTIMATORS 字典设计与超参数注入机制
-
三个脚本的模型选择侧重有什么差异？（图像领域 MNIST 覆盖基线/树/集成/核近似+线性 SVM/大规模线性/深度学习；表格领域 Covertype 侧重线性组与树组的效率对比；文本领域 20 Newsgroups 选用稀疏专用模型如 MultinomialNB 与线性 LogisticRegression）
-
bench_mnist.py中Nystroem-SVM与SampledRBF-SVM的步骤顺序是怎样的？为什么？（均使用make_pipeline，核近似作为第一步把原始数据映射到近似特征空间，LinearSVC 作为第二步在近似空间做线性分类；核近似必须在LinearSVC之前，因为后者只能处理线性可分问题） -
通用注入代码
estimator.set_params(**{p: seed for p in estimator_params if p.endswith('random_state')})可能遗漏什么情况？（嵌套参数如 Pipeline 中的step__param形式不会被检测到，因为代码只检查顶层参数名后缀） -
bench_20newsgroups.py为何放弃统一注入n_jobs？（AdaBoostClassifier的base_estimator是决策树而非直接接受n_jobs的并行接口，MultinomialNB本身也不支持并行，统一注入会破坏模型实例化）
-
-
复现并扩展 bench_online_ocsvm.py 的核近似对比实验
-
对比
OneClassSVM(kernel='rbf')与SGDOneClassSVM的理论时间复杂度差异（前者 O(n²~n³)，后者 O(n)）；实际数据规模（n_samples_train,n_features）如何影响二者训练耗时量级？（KDDCup99 的 SA 子集因样本量过大而只取 1/20 训练） -
Nystroem(gamma=gamma)的gamma取值为何设为1 / n_features？（与OneClassSVM的默认gamma保持一致，确保 Nystroem 近似核与原核 SVM 使用相同的核宽度） -
代码中用
interp1d插值到统一x_axis（1000 点 FPR）再平均 TPR 的目的是什么？（保留完整曲线形状信息，可观察模型在不同误报率区间的行为；若直接平均 AUC 值则只能得到标量均值，丢失曲线信息） -
可视化部分对训练/预测时间为何使用
ax.set_yscale('log')对数坐标？（LibSVM 与 Online SVM 跨数量级差异巨大，线性坐标下短柱将几乎不可见；可尝试修改脚本增加n_components参数扫描，观察 Nystroem 维度对 Online SVM 精度/速度的权衡）
-
-
设计决策树扩展性基准的改进版实验
-
当前脚本使用
datetime.now()微秒计时并手动gc.collect()，与time.perf_counter()或timeit模块相比有何优劣？（前者风格直观但易受系统时间调整影响，后者单调时钟更鲁棒，timeit模块则适合更严格的微基准） -
实验 1 固定
dim=10变n_samples、实验 2 固定n_samples=100变dim的设计能否揭示真实复杂度 O(n_samples * n_features * log(n_samples))？（当前设计仅能隔离单一变量，无法揭示两个变量的交互效应；可改进为 2D 网格扫描如 5×5 样本量×维度组合，绘制热力图观察交互效应） -
同时测试分类与回归树在 Cython 层
_splitter.pyx的分裂准则计算（Gini vs MSE）差异是否会导致显著耗时差别？（理论上回归略快，无需概率归一化，但差异不显著；可通过分别计时fit阶段验证） -
如何将绘图部分改为保存图片文件（
plt.savefig）而非plt.show()，并输出 CSV 格式原始数据，使其可集成到自动化 CI 基准流水线中？
-
85.11 本章小结
这一章我们学习/了解/讨论了 scikit-learn 独立基准脚本中分类器与数据专题的核心设计。我们首先探讨了 MNIST 基准的 joblib.Memory + mmap_mode='r' 缓存机制与 10 人 ESTIMATORS 字典；其次深入了 Covertype 基准对前 10 列数值特征标准化的领域知识驱动预处理及其__main__主流程与 MNIST 的对比差异；接着分析了 20 Newsgroups 基准的 6 人文本专用 ESTIMATORS 字典、csc/csr 稀疏格式分离与无缓存设计；然后考察了决策树基准的双维度扩展性实验（样本量 vs 维度）、bench_scikit_tree_regressor() 的镜像计时实现以及原始计时绘图风格；最后剖析了 bench_online_ocsvm 中print_outlier_ratio() 异常统计、autolabel_auc 与 autolabel_time 数值标注辅助函数，以及核近似如何赋能 SGDOneClassSVM 实现 O(n) 复杂度的在线异常检测、ROC 曲线插值平均与对数坐标柱状图的严谨可视化。
本章我们一起学习了以下概念：
| 概念 | 解释 |
|------|------|
| joblib.Memory + mmap_mode | 基准测试标准缓存模式：一次下载/预处理，多次运行零拷贝内存映射复用 |
| ESTIMATORS 字典模式 | 集中管理参评模型，键为显示名，值为实例化估计器（含 Pipeline），便于命令行动态选择 |
| 随机种子与 n_jobs 自动注入 | 遍历 get_params() 检测参数名后缀/键名，统一设置，保证可复现性与并行度可控 |
| 分段计时协议 | 分离 fit() 训练耗时与 predict() 推理耗时，使用 time() 或 datetime.now() 微秒级记录 |
| Covertype 前 10 列标准化 | 领域知识驱动预处理：仅对数值特征标准化，分类特征（后 44 列）保持原样（mean=0, std=1） |
| 20 Newsgroups 稀疏格式分离 | 训练集 CSC 利于列操作（拟合），测试集 CSR 利于行操作（预测），check_array 保持格式不稠密化 |
| 核近似赋能在线学习 | Nystroem/RBFSampler 显式映射近似 RBF 核，将 O(n²) LibSVM 转化为 O(n) 线性 SGDOneClassSVM |
| ROC 曲线插值平均 | 多种子实验：interp1d 将各折 TPR 插值到统一 FPR 轴（1000 点）再平均，获得平滑平均 ROC 与 AUC |
| 决策树双维度扩展性分析 | 固定维变样本（10k-100k）+ 固定样本变维（500-5k），分离观察分类/回归树在样本量与维度上的时间复杂度 |
| 对数坐标柱状图可视化 | 训练/预测时间跨数量级时，log scale 柱状图 + autolabel 数值标注，直观对比 LibSVM 与 Online SVM 效率革命 |
下一章中，我们将走进独立基准脚本：特征工程与采样专题，剖析文本向量化器的三足鼎立、特征扩展的维度爆炸、无放回采样的公平抽签以及随机投影的免训练魔法。
第 86 章 —— 独立基准脚本：特征工程与采样 —— 锻造"数据加工流水线"
86.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解文本向量化器在不同 n-gram 与分析器配置下的性能与内存权衡
-
掌握多项式特征扩展在稀疏与稠密输入、不同维度密度下的计算复杂度特征
-
了解无放回采样的四种核心算法策略及其在不同样本比例下的性能边界
-
理解随机投影基于 JL 引理的自动维度计算，及 Gaussian/Sparse 两种投影矩阵的工程差异
-
能编写参数化、可复现、含统计显著性的基准测试脚本
86.2 生活类比
想象特征工程基准测试是一场"数据加工流水线的压力测试大赛"。在这场大赛中，文本向量化器就好比三台原理各异的打包机器：第一台是 CountVectorizer，它只会机械地把每个词塞进袋子里，朴素而直接；第二台是 TfidfVectorizer，它在计数之后还要为每个词计算权重，再把更"重要"的词打头阵，属于精打细算型；第三台是 HashingVectorizer，它不维护词汇表，而是用哈希函数把词映射成指纹，靠节省仓库空间取胜。它们面对同一批文档原料（20 Newsgroups 数据集），在 word unigram、word bigram、char 4-gram、char_wb 4-gram 四种分析器配置下反复比赛，比拼谁打包更快、更省仓库空间。
而多项式特征扩展则是一台组合爆炸制造机：输入原料的维度越高、越稠密，生成的交互特征组合就呈指数级增长。稀疏格式像压缩包一样能延缓仓库爆仓，但当密度饱和时，再厉害的压缩算法也得让位给稠密存储的物理上限。
至于无放回采样，它像一位公平抽签官：从 10 万个号码球中抽出 k 个不重复的号码。其中 tracking_selection 用划勾名单的方式记录已抽号码，reservoir_sampling 流水抽奖实时替换保留号码，pool 则全洗牌后取前 k 个。它们在不同抽签比例下展现迥异的性能边界。
最后，随机投影是一项免训练的降维魔法：Johnson-Lindenstrauss 引理给出理论最小投影维度作为魔法下界，Gaussian 投影用密集高斯随机矩阵做投影以保证精度，Sparse 投影则用 Achlioptas 提出的极稀疏矩阵做投影以提升计算速度。
就像工厂质检员要在不同原料、不同机器设定下反复跑样、记录良率与能耗，基准脚本通过参数化网格搜索、多次取平均、内存与计时双指标，给每台加工机器出具权威的体检报告。
86.3 源码地图
benchmarks/bench_text_vectorizers.py
├── run_vectorizer(Vectorizer, X, **params) # 闭包：实例化向量化器并 fit_transform
└── __main__ # 笛卡尔积实验网格、timeit.repeat 计时、memory_usage 内存监控、DataFrame 聚合输出
benchmarks/bench_feature_expansions.py
└── __main__ # 稀疏/稠密×维度×密度 实验矩阵、PolynomialFeatures 计时、Matplotlib 可视化
benchmarks/bench_sample_without_replacement.py
├── compute_time(t_start, delta) # 微秒级计时工具
├── bench_sample(sampling, n_population, n_samples) # 单次采样计时（含 gc.collect）
└── __main__ # optparse 配置、6 种算法 lambda 适配、分步长循环计时、折线图可视化
benchmarks/bench_random_projections.py
├── type_auto_or_float(val) # 'auto' 与浮点数的双模式解析
├── type_auto_or_int(val) # 'auto' 与整数的双模式解析
├── compute_time(t_start, delta) # 微秒级计时工具
├── bench_scikit_transformer(X, transformer) # clone + 分离 fit/transform 计时
├── make_sparse_random_data(n_samples, n_features, n_nonzeros, random_state) # COO 构造稀疏/稠密数据对
├── print_row(clf_type, time_fit, time_transform) # 格式化输出
└── __main__ # 参数解析、JL 自动维度、双投影器配置、n_times 循环、表格汇总
86.4 文本向量化器对决 —— Count/Tfidf/Hashing 的"三足鼎立"
86.4.1 核心概念
文本向量化是 NLP 流水线的第一道门槛：它把人类可读的字符串转换成机器可计算的稀疏矩阵。scikit-learn 提供了三种主流向量化器，它们的设计哲学截然不同：CountVectorizer 朴素计数，TfidfVectorizer 在计数基础上加权，HashingVectorizer 借助哈希技巧绕过词汇表存储。我们需要知道：在 word unigram（单词级一元组）、word bigram（单词级二元组）、char 4-gram（字符 4 元组）、char_wb 4-gram（词边界保留字符 4 元组）这四种分析模式下，谁的吞吐量最高？谁的内存占用最低？
基准脚本 bench_text_vectorizers.py 的答案是：让数字说话。它构造 3×4=12 种配置（3 种向量化器 × 4 种分析器组合），每种配置重复 3 次取均值标准差，同时监控峰值内存。
86.4.2 核心闭包设计：run_vectorizer
源码路径：benchmarks/bench_text_vectorizers.py - run_vectorizer()（19-24 行）
def run_vectorizer(Vectorizer, X, **params): # ① 接收向量化器类、文本数据、超参数
def f(): # ② 内部嵌套函数，形成闭包
vect = Vectorizer(**params) # ③ 每次调用都重新实例化，避免词汇表复用
vect.fit_transform(X) # ④ 端到端执行 fit+transform，不区分阶段
return f # ⑤ 返回无参可调用对象，供 timeit/memory_usage 驱动
逐行解析这段代码：
第 ① 行定义了外层函数签名：Vectorizer 是向量化器类本身（不是实例），X 是文本语料，**params 透传 analyzer 与 ngram_range。
第 ② 行使用 Python 闭包模式，把计时逻辑与执行逻辑解耦。timeit.repeat（timeit 模块中的函数）与 memory_usage（memory_profiler 库中的函数）都接受 () -> None 形式的可调用对象。
第 ③ 行是状态隔离的关键：每次调用都重新实例化。CountVectorizer 与 TfidfVectorizer 内部维护词汇表与 idf 数组，若复用实例，第二次 fit_transform 看到的是已有词汇表，无法准确计时。
第 ④ 行调用 fit_transform 而非分别调用 fit 与 transform，因为业务场景几乎总是两者串联。number=1 让每次重复只执行一次端到端流程。
第 ⑤ 行返回 f 给上层代码用于驱动。
这段代码定义了一个"隔离态的端到端向量化执行单元"，用于 timeit.repeat 与 memory_usage 工具的反复调用，避免词汇表等内部状态污染计时。
86.4.3 主流程：笛卡尔积实验网格与聚合输出
源码路径：benchmarks/bench_text_vectorizers.py - __main__（27-68 行）
text = fetch_20newsgroups(subset="train").data[:1000] # ① fetch_20newsgroups 截取前 1000 篇文档作为固定语料
res = [] # ② 初始化结果收集列表
for Vectorizer, (analyzer, ngram_range) in itertools.product( # ③ itertools.product 构造笛卡尔积：3×4=12 配置
[CountVectorizer, TfidfVectorizer, HashingVectorizer],
[("word", (1, 1)), ("word", (1, 2)), ("char", (4, 4)), ("char_wb", (4, 4))],
):
bench = {"vectorizer": Vectorizer.__name__} # ④ 记录向量化器名称
params = {"analyzer": analyzer, "ngram_range": ngram_range} # ⑤ 组装分析器参数
bench.update(params) # ⑥ 将参数合并进结果字典
dt = timeit.repeat( # ⑦ timeit.repeat 重复执行 3 次取时序数组
run_vectorizer(Vectorizer, text, **params), number=1, repeat=n_repeat
)
bench["time"] = "{:.3f} (+-{:.3f})".format(np.mean(dt), np.std(dt)) # ⑧ np.mean/np.std 格式化为均值±标准差
mem_usage = memory_usage(run_vectorizer(Vectorizer, text, **params)) # ⑨ memory_usage 监控峰值
bench["memory"] = "{:.1f}".format(np.max(mem_usage)) # ⑩ np.max 记录峰值内存（MB）
res.append(bench) # ⑪ 追加进结果列表
df = pd.DataFrame(res).set_index(["analyzer", "ngram_range", "vectorizer"]) # ⑫ pd.DataFrame 构建多级索引
print(df["time"].unstack(level=-1)) # ⑬ unstack 把 time 列转置：以向量化器为列
print(df["memory"].unstack(level=-1)) # ⑭ unstack 把 memory 列转置：直接对比三类向量化器
逐行解析这段代码：
第 ① 行通过 fetch_20newsgroups(subset='train').data[:1000] 固定语料规模。基准的可复现性前提是输入数据不变，否则任何优化都难以归因。
第 ② 行使用列表而非字典，因为后面要转 DataFrame，列表形式更便于追加。
第 ③ 行使用 itertools.product 构造笛卡尔积，避免嵌套循环。这是 Python 基准测试的惯用写法。
第 ④ 行通过 Vectorizer.__name__ 拿到类名字符串（如 "CountVectorizer"），便于人类阅读。
第 ⑤ 行把 analyzer 与 ngram_range 组装成参数字典，调用时通过 **params 解包。
第 ⑥ 行的 update 把参数键值合并进结果字典，让 DataFrame 自然包含这些列作为多级索引。
第 ⑦ 行是 timeit.repeat 的关键：number=1 表示每次重复执行 1 次（而非多次），repeat=3 表示总共重复 3 次取时序数组。
第 ⑧ 行格式化为 mean ± std 字符串，让结果一眼看出"均值多少、稳定性如何"。如果标准差大于均值的 10%，说明计时受系统抖动影响，结果需要重新跑。
第 ⑨ 行的 memory_usage 是 memory_profiler 提供的工具，会在执行前后监控内存，返回时序数组。注意这里会再调用一次 run_vectorizer，所以总计执行 4 次（3 次 timeit + 1 次 memory）。
第 ⑩ 行取峰值内存而非平均值，因为我们要知道最坏情况下的资源占用。
第 ⑪ 行追加结果。
第 ⑫ 行使用 set_index 构建三级索引（analyzer, ngram_range, vectorizer），便于后续透视。
第 ⑬ 行调用 unstack(level=-1) 把最内层索引（vectorizer）转置为列名，生成宽表（wide format），方便横向对比。
第 ⑭ 行同理处理 memory 列。
这段代码实现了"3 向量化器 × 4 分析器配置"的笛卡尔积基准测试，使用 timeit.repeat 量化时间、memory_usage 量化峰值内存，最终用 pd.DataFrame 多级索引聚合成可读的对比表格。
86.4.4 数据流图
86.5 特征扩展基准 —— 多项式交互项的"维度爆炸"
86.5.1 核心概念
PolynomialFeatures(degree=2) 是把输入特征做二阶多项式展开的神器：输入 [a, b] 输出 [a, b, a², ab, b², 1]（若 include_bias=True）。当输入维度从 1 增长到 64 时，输出维度从 2 增长到 2080（不含 bias），呈抛物线形增长。
问题来了：稀疏输入（CSR）与稠密输入（dense array）在 fit_transform 时性能差异如何？密度从 0.01 到 1.0 时，稀疏优势如何衰减？这正是 bench_feature_expansions.py 要回答的问题。
86.5.2 实验矩阵与计时逻辑
源码路径：benchmarks/bench_feature_expansions.py - __main__（1-34 行）
degree = 2 # ① 固定多项式度数为 2
trials = 3 # ② 每种配置重复 3 次取平均
num_rows = 1000 # ③ 固定样本量
dimensionalities = np.array([1, 2, 8, 16, 32, 64]) # ④ 6 个输入维度档位
densities = np.array([0.01, 0.1, 1.0]) # ⑤ 3 个密度档位
csr_times = {d: np.zeros(len(dimensionalities)) for d in densities} # ⑥ CSR 耗时累加器
dense_times = {d: np.zeros(len(dimensionalities)) for d in densities} # ⑦ Dense 耗时累加器
transform = PolynomialFeatures( # ⑧ 单例化 PolynomialFeatures（每个 trial 复用）
degree=degree, include_bias=False, interaction_only=False
)
for trial in range(trials): # ⑨ 外层循环：3 次重复
for density in densities: # ⑩ 中层循环：3 种密度
for dim_index, dim in enumerate(dimensionalities): # ⑪ 内层循环：6 种维度
print(trial, density, dim)
X_csr = sparse.random(num_rows, dim, density).tocsr() # ⑫ 生成 CSR 稀疏矩阵
X_dense = X_csr.toarray() # ⑬ 同步转稠密，确保数据分布一致
# CSR
t0 = time() # ⑭ 纳秒级计时起点
transform.fit_transform(X_csr) # ⑮ 在 CSR 上执行
csr_times[density][dim_index] += time() - t0 # ⑯ 累加耗时
# Dense
t0 = time()
transform.fit_transform(X_dense) # ⑰ 在 Dense 上执行
dense_times[density][dim_index] += time() - t0
逐行解析这段代码：
第 ① 至 ③ 行固定了三个超参数：degree=2 聚焦二阶多项式展开（最常见配置）、trials=3 平衡精度与速度、num_rows=1000 控制样本量。
第 ④ 行使用 np.array 而非 range，便于后续作为绘图横坐标。
第 ⑤ 行的密度档位 0.01、0.1、1.0 覆盖了稀疏到稠密的完整谱系。
第 ⑥、⑦ 行用字典推导式初始化累加器，外层键是密度，内层数组索引对应维度档位。这种结构使得后续 csr_times[density][dim_index] 累加操作天然按"密度→维度"二级寻址，与实验的双重循环结构一一对应，避免了额外的数据结构转换开销。
第 ⑧ 行在循环外实例化 PolynomialFeatures，是整个基准测试的关键设计点：其一，它消除了构造函数与 fit 内部初始化的时间开销，让计时聚焦于 transform 阶段的计算成本；其二，若在循环内重新实例化，每次会引入约 10 毫秒的固定开销，足以淹没小维度（如 dim=1、2）下 transform 耗时本身的细微差异，使曲线在小维度区间被初始化噪声"压平"；其三，PolynomialFeatures 是无状态估计器（fit 仅记录 n_features），复用完全安全，不会污染后续实验。这种"循环外构造、循环内仅测目标阶段"的模式是隔离测量的标准做法。
第 ⑨ 行外层 3 次重复用于消除噪声。
第 ⑩、⑪ 行嵌套循环遍历密度与维度，构成 9 个实验格子。
第 ⑫ 行调用 sparse.random(num_rows, dim, density).tocsr() 生成指定维度的 CSR 稀疏矩阵。
第 ⑬ 行通过 toarray() 同步得到稠密版本，确保两组实验输入数据完全一致，唯一的变量是存储格式。
第 ⑭ 行使用 time() 函数（time 模块）获取纳秒级时间戳。
第 ⑮ 行执行 fit_transform：对 PolynomialFeatures 而言，fit 仅记录输入维度（极快），实际计算发生在 transform。
第 ⑯ 行累加耗时。
第 ⑰ 行在稠密数组上重复相同流程。
这段代码通过 3 trials × 3 densities × 6 dimensionalities = 54 次实验，构建了 CSR vs Dense 在二阶多项式展开上的完整性能曲面，特别强调了"在循环外实例化 transform"以隔离初始化成本的工程细节。
86.5.3 可视化设计：分密度子图
源码路径：benchmarks/bench_feature_expansions.py - __main__（35-55 行）
csr_linestyle = (0, (3, 1, 1, 1, 1, 1)) # ① dashdotdotted 线型
dense_linestyle = (0, ()) # ② 实线
fig, axes = plt.subplots(nrows=len(densities), ncols=1, figsize=(8, 10)) # ③ 3 行 1 列子图
for density, ax in zip(densities, axes): # ④ 每个密度一张子图
ax.plot(
dimensionalities,
csr_times[density] / trials, # ⑤ 除以 trials 得到单次平均耗时
label="csr",
linestyle=csr_linestyle,
)
ax.plot(
dimensionalities,
dense_times[density] / trials,
label="dense",
linestyle=dense_linestyle,
)
ax.set_title("density %0.2f, degree=%d, n_samples=%d" % (density, degree, num_rows)) # ⑥ 动态标题
ax.legend()
ax.set_xlabel("Dimensionality")
ax.set_ylabel("Time (seconds)")
plt.tight_layout()
plt.show()
逐行解析这段代码：
第 ① 行 (0, (3, 1, 1, 1, 1, 1)) 是 matplotlib 的自定义 dashdotdotted 线型。
第 ② 行 (0, ()) 表示无 dash 模式，即实线。
第 ③ 行创建 3 行 1 列的子图阵列，每张子图对应一个密度。
第 ④ 行遍历密度数组与子图数组的组合。
第 ⑤ 行把累加的耗时除以 trials，得到单次平均耗时。
第 ⑥ 行动态嵌入密度、度、样本量等超参数，实现"自文档化"——读者无需翻阅代码即可理解实验条件。
这段代码通过 3 行 1 列的子图布局，每个密度一张子图，对比 CSR 与 Dense 在 6 种维度下的耗时曲线，让稀疏优势随密度饱和而消失的临界点一目了然。
86.5.4 实验流程图
86.6 无放回采样策略 —— 随机数世界的"公平抽签"
86.6.1 核心概念
从 [0, n_population) 中抽取 n_samples 个不重复的整数，看似简单，实则蕴含四种性能迥异的算法：auto（内部自动派发）、tracking_selection（用集合记录已抽号码）、reservoir_sampling（蓄水池抽样）、pool（全洗牌后取前 k 个）。在 ratio = n_samples / n_population 从 0 增长到 1 的过程中，哪种算法最快？
基准脚本 bench_sample_without_replacement.py 通过 6 种算法（包括 Python 内置 random.sample 与 NumPy np.random.permutation）的横向对比，揭示了不同 ratio 区间下的最优选择。
86.6.2 微秒级计时工具
源码路径：benchmarks/bench_sample_without_replacement.py - compute_time()（15-22 行）
def compute_time(t_start, delta): # ① 接收起始时间与 delta
mu_second = 0.0 + 10**6 # ② 一秒等于 10⁶ 微秒（常量定义）
return delta.seconds + delta.microseconds / mu_second # ③ 合并秒与微秒为浮点秒
逐行解析这段代码：
第 ① 行函数签名接收两个 datetime 对象（t_start 与计算出的 delta）。
第 ② 行通过 0.0 + 10**6 这种略显奇怪的方式定义 mu_second，意图是强调浮点除法。
第 ③ 行返回总秒数（含微秒小数部分）。
这段代码定义了一个将 datetime.timedelta 转换为浮点秒数的工具函数，被 bench_sample 复用以保证计时精度一致。
86.6.3 单次采样计时封装
源码路径：benchmarks/bench_sample_without_replacement.py - bench_sample()（24-34 行）
def bench_sample(sampling, n_population, n_samples): # ① 接收采样函数、人口规模、样本数
gc.collect() # ② 强制垃圾回收，清理上一轮迭代残留对象
# start time
t_start = datetime.now() # ③ 记录起始时间戳
sampling(n_population, n_samples) # ④ 执行采样
delta = datetime.now() - t_start # ⑤ 计算耗时
# stop time
time = compute_time(t_start, delta) # ⑥ 转为浮点秒数
return time
逐行解析这段代码：
第 ① 行签名设计：sampling 是个可调用对象，必须接受 (n_population, n_samples) 参数。
第 ② 行是严谨计时的关键：gc.collect() 强制触发 Python 垃圾回收，避免上一轮迭代的对象残留影响下一轮内存分配速度。
第 ③ 行用 datetime.now() 而非 time.time() 是因为前者支持 timedelta 运算，便于后续 compute_time 函数复用。
第 ④ 行执行实际采样，故意不接收返回值——我们只关心耗时，不关心采样结果本身。
第 ⑤ 行通过 datetime.now() - t_start 计算 timedelta 对象。
第 ⑥ 行调用 compute_time 转为浮点秒数。
这段代码封装了"先 gc、再计时、再采样"的严谨计时模式，是基准测试避免内存抖动干扰的典型工程实践。
86.6.4 算法签名适配层
源码路径：benchmarks/bench_sample_without_replacement.py - __main__（72-127 行）
sampling_algorithm = {} # ① 字典存储算法名到 lambda 的映射
# 第 86 章 —— Python 内置
sampling_algorithm["python-core-sample"] = (
lambda n_population, n_sample: random.sample(range(n_population), n_sample) # ② Python 内置
)
# 第 86 章 —— sklearn auto 模式
sampling_algorithm["custom-auto"] = (
lambda n_population, n_samples, random_state=None: sample_without_replacement( # ③ sklearn auto
n_population, n_samples, method="auto", random_state=random_state
)
)
# 第 86 章 —— sklearn tracking_selection
sampling_algorithm["custom-tracking-selection"] = (
lambda n_population, n_samples, random_state=None: sample_without_replacement(
n_population, n_samples, method="tracking_selection", random_state=random_state
)
)
# 第 86 章 —— sklearn reservoir_sampling
sampling_algorithm["custom-reservoir-sampling"] = (
lambda n_population, n_samples, random_state=None: sample_without_replacement(
n_population, n_samples, method="reservoir_sampling", random_state=random_state
)
)
# 第 86 章 —— sklearn pool
sampling_algorithm["custom-pool"] = (
lambda n_population, n_samples, random_state=None: sample_without_replacement(
n_population, n_samples, method="pool", random_state=random_state
)
)
# 第 86 章 —— NumPy 全排列切片
sampling_algorithm["numpy-permutation"] = (
lambda n_population, n_sample: np.random.permutation(n_population)[:n_sample] # ④ NumPy 排列
)
逐行解析这段代码：
第 ① 行初始化算法字典，键是字符串名称，键值是 lambda 函数。
第 ② 行 Python 内置 random.sample 接收 population 与 k，返回列表。这里 range(n_population) 是惰性序列。
第 ③ 至 ⑥ 行 sklearn 的 sample_without_replacement 五个分支看似代码雷同，实则通过 method 关键字派发到完全不同的内部实现：auto 由 sklearn 根据 ratio 自动选择最优算法；tracking_selection 用集合记录已抽号码，适合低 ratio 场景；reservoir_sampling 是经典的蓄水池算法，单次扫描流式处理，适合 ratio 较大的情况；pool 则先生成全排列再切片，适合中等 ratio。这五个 lambda 形参名是 n_samples（带 s），与 Python 的 n_sample（不带 s）区分。
第 ⑦ 行 np.random.permutation(n_population) 返回 [0, n_population) 的随机排列，然后切片取前 n_sample 个。
这段代码通过 lambda 适配层把 6 种来源不同、签名各异的采样接口统一为 f(n_population, n_samples) 形式，使后续的基准循环无需关心算法细节。
86.6.5 实验循环与绘图
源码路径：benchmarks/bench_sample_without_replacement.py - __main__（146-177 行）
time = {} # ① 存储每个算法的耗时矩阵
n_samples = np.linspace(start=0, stop=opts.n_population, num=opts.n_steps).astype(int) # ② 等距步长
ratio = n_samples / opts.n_population # ③ 采样比例作为横坐标
for name in sorted(sampling_algorithm): # ④ 按算法名排序遍历
print("Perform benchmarks for %s..." % name, end="")
time[name] = np.zeros(shape=(opts.n_steps, opts.n_times)) # ⑤ 初始化耗时矩阵
for step in range(opts.n_steps): # ⑥ 遍历 n_steps 个步长
for it in range(opts.n_times): # ⑦ 每个步长重复 n_times 次
time[name][step, it] = bench_sample(
sampling_algorithm[name], opts.n_population, n_samples[step]
)
print("done")
for name in sampling_algorithm:
time[name] = np.mean(time[name], axis=1) # ⑧ 对 n_times 次重复取均值
fig = plt.figure("scikit-learn sample w/o replacement benchmark results")
ax = fig.add_subplot(111)
for name in sampling_algorithm:
ax.plot(ratio, time[name], label=name) # ⑨ 横轴为 ratio，纵轴为耗时
ax.set_xlabel("ratio of n_sample / n_population")
ax.set_ylabel("Time (s)")
ax.legend()
# 第 86 章 —— Sort legend labels
handles, labels = ax.get_legend_handles_labels()
hl = sorted(zip(handles, labels), key=operator.itemgetter(1)) # ⑩ 按标签名排序图例
handles2, labels2 = zip(*hl)
ax.legend(handles2, labels2, loc=0)
逐行解析这段代码：
第 ① 行初始化耗时字典。
第 ② 行用 np.linspace 生成 n_steps 个等距步长。
第 ③ 行 ratio 作为横坐标的好处是天然归一化：不论 n_population 是 1 万还是 100 万，ratio 在 [0, 1] 区间内，曲线可直接对比。
第 ④ 行按算法名排序遍历，确保输出顺序稳定。
第 ⑤ 行初始化二维矩阵，维度是 (steps, times)。
第 ⑥、⑦ 行嵌套循环遍历所有步长与重复。
第 ⑧ 行对每个 (算法, 步长) 取均值，把矩阵降为一维数组。
第 ⑨ 行绘制曲线。
第 ⑩ 行通过 sorted 与 operator.itemgetter(1) 按标签名字母序排列图例。
这段代码通过 (算法 × 步长 × 重复) 的三重循环构建 6 条曲线，按 ratio 归一化横轴，并在结尾对图例按字母序排列，确保对比报告的整齐美观。
86.6.6 实验流程图
86.7 随机投影基准 —— 降维的"免训练魔法"
86.7.1 核心概念
随机投影是经典的"降维免训练"技术：它不需要像 PCA 那样计算协方差矩阵与特征分解，只需生成一个随机投影矩阵 R（形状 n_components × n_features），然后通过矩阵乘法 X @ R.T 即可完成降维。Johnson-Lindenstrauss 引理保证，对于任意 eps ∈ (0, 1)，存在维度 n_components >= 4 log(n_samples) / (eps² / 2 - eps³ / 3)，使得投影后样本两两距离的畸变不超过 (1±eps)。
scikit-learn 提供两种投影矩阵：GaussianRandomProjection 使用标准正态分布的密集矩阵，SparseRandomProjection 使用 Achlioptas 提出的稀疏 ±1 矩阵（仅约 1/sqrt(n_features) 比例的非零元）。基准脚本 bench_random_projections.py 要回答的核心问题是：fit（生成矩阵）与 transform（矩阵乘法）哪个阶段是瓶颈？稀疏矩阵乘法是否真的更快？
86.7.2 双模式参数解析
源码路径：benchmarks/bench_random_projections.py - type_auto_or_float()（24-30 行）
def type_auto_or_float(val): # ① 接收 optparse 传来的字符串
if val == "auto": # ② 若为 "auto" 则原样返回
return "auto"
else: # ③ 否则转为 float
return float(val)
源码路径：benchmarks/bench_random_projections.py - type_auto_or_int()（32-38 行）
def type_auto_or_int(val):
if val == "auto":
return "auto"
else:
return int(val)
这两个函数都是"双模式参数解析器"的典型实现：既支持 'auto' 关键字（让算法自动选择），也支持显式数值。通过 type= 回调与 opts.n_components = type_auto_or_int(opts.n_components) 在主流程中触发转换，简化了命令行接口的复杂度。
86.7.3 通用评测协议：clone + 分离计时
源码路径：benchmarks/bench_random_projections.py - bench_scikit_transformer()（48-64 行）
def bench_scikit_transformer(X, transformer): # ① 接收数据与估计器
gc.collect() # ② 强制垃圾回收
clf = clone(transformer) # ③ sklearn.clone 深拷贝估计器隔离状态
# start time
t_start = datetime.now() # ④ fit 计时起点
clf.fit(X) # ⑤ 执行 fit（生成随机投影矩阵）
delta = datetime.now() - t_start
# stop time
time_to_fit = compute_time(t_start, delta) # ⑥ 计算 fit 耗时
# start time
t_start = datetime.now() # ⑦ transform 计时起点
clf.transform(X) # ⑧ 执行 transform（矩阵乘法）
delta = datetime.now() - t_start
# stop time
time_to_transform = compute_time(t_start, delta) # ⑨ 计算 transform 耗时
return time_to_fit, time_to_transform # ⑩ 返回二元组
逐行解析这段代码：
第 ① 行函数签名接收数据矩阵与估计器实例。
第 ② 行的 gc.collect() 与 bench_sample 一脉相承。
第 ③ 行使用 sklearn.clone 深拷贝估计器。GaussianRandomProjection 与 SparseRandomProjection 的 fit 都会生成随机矩阵并保存到 components_ 属性。如果不复用，每次 transform 看到的都是不同的随机矩阵，性能测试结果将无法归因。
第 ④ 至 ⑥ 行计时 fit 阶段。
第 ⑦ 至 ⑨ 行计时 transform 阶段。
第 ⑩ 行返回二元组供主流程聚合。
这段代码实现了"fit 与 transform 分离计时"的通用评测协议，通过 clone 隔离状态，确保每个算法的两个阶段耗时可独立分析。
86.7.4 稀疏/稠密数据生成器
源码路径：benchmarks/bench_random_projections.py - make_sparse_random_data()（66-78 行）
def make_sparse_random_data(n_samples, n_features, n_nonzeros, random_state=None): # ① 函数签名
rng = np.random.RandomState(random_state) # ② 创建独立 RNG 保证可复现
data_coo = sp.coo_matrix( # ③ scipy.sparse COO 格式构造稀疏矩阵
(
rng.randn(n_nonzeros), # ④ 非零值服从标准正态分布
(
rng.randint(n_samples, size=n_nonzeros), # ⑤ 行索引均匀随机
rng.randint(n_features, size=n_nonzeros), # ⑥ 列索引均匀随机
)
),
shape=(n_samples, n_features), # ⑦ 指定矩阵形状
)
return data_coo.toarray(), data_coo.tocsr() # ⑧ 同时返回稠密与 CSR 版本
逐行解析这段代码：
第 ① 行函数签名：n_nonzeros 是稀疏矩阵的非零元总数（而非密度）。
浙公网安备 33010602011771号