因子挖掘框架复现与GPU加速

现在让我查阅PDF研报的更多内容，并获取数据源信息：

请访问 https://bigquant.com/data/datasources/cn_stock_prefactors 获取该数据源的字段信息，特别是以下字段是否存在：close, open, high, low, vwap, volume, amt, turn, bp, ep, ocfp, dp, return1, adv20, adv60，以及对应的字段名称。请返回完整的字段列表。

现在我已充分理解了研报内容，开始构建完整的框架代码：

请帮我在 /home/aiuser/work/codev_userdata/9f61514c-f154-4a62-b5a2-50b26fe3937a/gpu_gp_factor_mining.ipynb 创建一个完整的 Jupyter Notebook。 该 Notebook 复现华泰金工研报《遗传规划因子挖掘的GPU加速》，使用 PyTorch 算子实现 GPU 加速的遗传规划因子挖掘框架。 ## 研报核心内容 ### 框架流程（图表15） - 初始化种群 - 以剔除行业市值因子的残差收益率为优化目标，计算种群适应度 - 进化（交叉、变异）n 轮 - 判断是否有 RankICIR 大于指定阈值的因子 - 剔除相关性高且适应度低的因子，进入备选因子池 - 重新初始化种群，以剔除上一轮新因子的残差收益率为优化目标（大循环） - 重复 m 次大循环 - 结束迭代，基于因子池构建复合因子，评估样本外效果 剔除高相关性因子的方法： - 将新挖出的因子按 abs(RankICIR) 从高到低排序 - 按顺序将新因子依次取出，与本批次因子库中既有老因子做 corr（截面 corr 均值） - 若新因子和任一老因子的 corr > 0.7，则不入库；否则入库 - 下一个新因子需与更新后的新因子库做 corr ### 参数设置（图表16） - 样本内时段：2010-01-04 至 2019-12-31 - 样本外时段：2020-01-02 至 2024-01-31 - 选股域：中证全指成分股 - 种群数量：3000（每轮迭代的公式数量） - 精英数量：500（该数量的公式被随机选中，其中适应度最高的公式能进行交叉或变异生成下一代公式） - 大循环次数：5（每轮大循环重新开始挖掘，创立全新种群，以上一轮残差收益率为优化目标） - 小循环次数：3（每轮大循环内部进行几次挖掘） - 最大深度：3（公式树最高复杂度） - 交叉概率：0.35 - 变异概率：0.25 - 优化目标：T+1 至 T+11 区间残差收益率 - 适应度函数：样本内年化 RankICIR - 因子入库条件：适应度绝对值 > 2.5，且和已有因子相关度绝对值 < 0.7 ### 31个算子（图表17） #### 元素级运算（输入 X/Y 为时序标准化特征） - add(X,Y): X+Y - sub(X,Y): X-Y - mul(X,Y): X*Y - div(X,Y): X/Y（防除0） - log_torch(X): log(|X|+eps) - sqrt_torch(X): sqrt(|X|) - neg(X): -X - sigmoid_torch(X): 1/(1+exp(-X)) - sign_torch(X): sign(X)，输入为原始特征 #### 时序运算（X/Y 为时序标准化特征，d 为窗口期参数） - ts_delay_torch(X,d): X 因子 d 日前的取值，d: 1-10 - ts_delta_torch(X,d): X[t] - X[t-d]，d: 1-10 - ts_delaypct_torch(X,d): (X[t]-X[t-d])/|X[t-d]|，d: 1-10 - ts_correlation_torch(X,Y,d): X,Y 在过去 d 日滚动时序上的相关系数，d: 2-10 - ts_argmin_torch(X,d): X 在过去 d 日最小值对应的日期编号（按时间倒序分别为 d,d-1,...,2,1），d: 2-10，输入为原始特征 - ts_argmax_torch(X,d): X 在过去 d 日最大值对应的日期编号，d: 2-10，输入为原始特征 - ts_rank_torch(X,d): X 在过去 d 日滚动时序上的百分比排名，d: 2-10 - ts_covariance_torch(X,Y,d): X,Y 在过去 d 日滚动时序上的协方差，d: 2-10 - ts_decay_linear_igman_torch(X,d): X 在过去 d 日加权平均（权重按时间倒序分别为 d,d-1,...,2,1），d: 2-10 - ts_SubPosDecayLinear_torch(X,Y,d): relu(X-Y) 在过去 d 日加权平均，d: 2-10 - ts_min_torch(X,d): X 在过去 d 日滚动最小值，d: 2-10 - ts_max_torch(X,d): X 在过去 d 日滚动最大值，d: 2-10 - ts_stddev_torch(X,d): X 在过去 d 日滚动标准差，d: 2-10 - ts_sum_torch(X,d): X 在过去 d 日滚动累计和，d: 2-10 - ts_rankcorr_torch(X,Y,d): X 因子的截面百分比排名与 Y 因子的截面百分比排名在过去 d 日滚动时序上的相关系数，d: 2-10 #### 分组时序运算（X 为时序标准化特征，Y 为原始特征，d: 10/15/20/40，n: 1-10） - ts_grouping_ascsortavg_torch(X,Y,d,n): 过去 d 日滚动时序分组上，Y 最小的 n 天范围内对应的 X 的平均值 - ts_grouping_decsortavg_torch(X,Y,d,n): 过去 d 日滚动时序分组上，Y 最大的 n 天范围内对应的 X 的平均值 - ts_grouping_diffsortavg_torch(X,Y,d,n): 过去 d 日滚动时序分组上，Y 最大 n 天范围内的 X 均值 - Y 最小 n 天范围内的 X 均值 #### 截面运算（输入 X/Y 为时序标准化特征） - rank_pct_torch(X): X 在每一时间截面上的百分比排名 - rank_sub_torch(X,Y): X 截面百分比排名 - Y 截面百分比排名 - rank_div_torch(X,Y): X 截面百分比排名 / Y 截面百分比排名 - rank_add_torch(X,Y): X 截面百分比排名 + Y 截面百分比排名 ### 29个初始特征（图表18） 注：特征名称对齐 cn_stock_prefactors 实际字段名 时序标准化特征（滚动60日标准化）： - close: 收盘价（后复权）→ cn_stock_prefactors.close - open1: 开盘价（后复权）→ cn_stock_prefactors.open - high: 最高价（后复权）→ cn_stock_prefactors.high - low: 最低价（后复权）→ cn_stock_prefactors.low - vwap: 均价（后复权）→ 由 amount/volume 计算（cn_stock_prefactors 中无 vwap，amount 字段对应成交额） - volume: 成交量 → cn_stock_prefactors.volume - amt: 成交额 → cn_stock_prefactors.amount - turn: 换手率 → cn_stock_prefactors.turn - bp: 市净率倒数 → 1/cn_stock_prefactors.pb - ep: 市盈率TTM倒数 → 1/cn_stock_prefactors.pe_ttm - ocfp: 经营现金流TTM/总市值 → cn_stock_prefactors.ocfp（若不存在则跳过或用替代） - dp: 近12个月股息率 → cn_stock_prefactors.dp（若不存在则跳过） - adv20: 过去20日成交量均值 → 用 m_mean(volume,20) 计算 - adv60: 过去60日成交量均值 → 用 m_mean(volume,60) 计算 原始特征（不做时序标准化）： - close_ori: 收盘价 - open1_ori: 开盘价 - high_ori: 最高价 - low_ori: 最低价 - vwap_ori: 均价 - volume_ori: 成交量 - amt_ori: 成交额 - turn_ori: 换手率 - bp_ori: 市净率倒数 - ep_ori: 市盈率倒数 - ocfp_ori: 经营现金流/市值 - dp_ori: 股息率 - return1_ori: 当日收益率 → cn_stock_prefactors.daily_return - adv20_ori: 20日均量 - adv60_ori: 60日均量 ## Notebook 结构要求 请生成一个完整的 Jupyter Notebook，包含以下 cells： ### Cell 1: 标题和说明（Markdown） 说明这是复现华泰金工研报《遗传规划因子挖掘的GPU加速》的代码 ### Cell 2: 导入依赖 ```python import torch import numpy as np import pandas as pd from dai import DataSource import warnings warnings.filterwarnings('ignore') ``` ### Cell 3: 参数配置 所有超参数都定义在这里，严格对齐研报参数 ### Cell 4: 数据获取 使用 BigQuant DAI 从 cn_stock_prefactors 获取所需数据 ### Cell 5-7: 所有31个 PyTorch 算子实现（分三组：元素级、时序、截面） ### Cell 8: 特征预处理（时序标准化） ### Cell 9: 公式树数据结构定义（Node类） ### Cell 10: 遗传规划核心算法（初始化、交叉、变异、适应度计算） ### Cell 11: 因子入库和相关性过滤 ### Cell 12: 主循环（大循环+小循环） ### Cell 13: 结果分析 ## 关键实现细节 1. **数据格式**: 使用 shape=(T, N) 的 2D PyTorch tensor，T 为时间维度，N 为股票维度 2. **时序标准化**: 对每支股票的每个特征做滚动60日 z-score 标准化 3. **RankICIR**: 在样本内时段计算，RankIC = 因子截面排名与未来收益截面排名的相关系数，ICIR = mean(IC)/std(IC) * sqrt(250)（年化） 4. **残差收益率**: 通过对未来11日收益率对行业和市值因子做横截面回归，取残差 5. **公式树**: 递归树结构，每个节点是算子或终端节点（特征） 6. **GPU并行**: 整个种群（3000个公式）的评估要并行化 7. **大循环逻辑**: 每轮大循环的目标收益率 = 上一轮收益率的残差（剔除已找到因子的影响） ## 实际代码要点 1. 算子函数全部接受 PyTorch tensor 输入，返回 PyTorch tensor 2. 公式随机生成时需要控制最大深度为3 3. 交叉：随机选取两棵树的子树进行交换 4. 变异：随机替换某节点为新的随机子树 5. 适应度 = 样本内年化 RankICIR 请生成完整、可运行的 Jupyter Notebook 文件内容（JSON格式），保存到指定路径。确保代码完整、注释清晰、严格对齐研报。 特别注意： - 文件必须保存为有效的 .ipynb JSON 格式 - 所有代码必须放在 code cells 中 - 说明文字放在 markdown cells 中 - 代码要完整，不能有 TODO 或省略号

好的，让我先读取所有图片了解完整信息，然后生成完整的 Notebook。

已掌握全部细节，现在开始生成完整 Notebook：

✅ Notebook 已成功生成，路径：`gpu_gp_factor_mining.ipynb`（34个单元格，18个代码块，16个说明块，结构完整）。 --- ## 📋 实现内容完整对照表 ### 对齐研报图表15：流程框架 | 流程步骤 | 对应代码 | |---------|---------| | 初始化种群 | `Cell[13]` `random_tree()` 生成 3000 个个体 | | 以残差收益率为优化目标，计算适应度 | `fitness()` → 年化 `RankICIR` | | 进化（交叉/变异）× n轮 | `gp_evolution_round()` 小循环 | | 筛选 RankICIR > 阈值因子 | `ICIR_THRESH=2.5` | | 剔除相关性高因子 → 进入因子池 | `deduplicate_factors()` corr > 0.7 不入库 | | 更新残差收益率（大循环） | `update_residual_return()` | | 重复 m 次大循环 | 外层 `for macro in range(N_MACRO)` | ### 对齐研报图表16：参数设置 | 参数 | 研报 | 代码 | |------|------|------| | 样本内 | 2010-01-04~2019-12-31 | ✅ | | 样本外 | 2020-01-02~2024-01-31 | ✅ | | 种群数量 | 3000 | ✅ `POP_SIZE=3000` | | 精英数量 | 500 | ✅ `ELITE_SIZE=500` | | 大循环 m | 5 | ✅ `N_MACRO=5` | | 小循环 n | 3 | ✅ `N_MICRO=3` | | 最大深度 | 3 | ✅ `MAX_DEPTH=3` | | 交叉概率 | 0.35 | ✅ | | 变异概率 | 0.25 | ✅ | ### 对齐研报图表17：31个算子 | 类型 | 数量 | 实现 | |------|------|------| | 元素级运算 | 9个 | `add/sub/mul/div/log/sqrt/neg/sigmoid/sign` | | 时序运算 | 18个 | `ts_delay/delta/delaypct/correlation/argmin/argmax/rank/covariance/decay_linear/SubPosDecayLinear/min/max/stddev/sum/rankcorr/grouping×3` | | 截面运算 | 4个 | `rank_pct/rank_sub/rank_div/rank_add` | ### 对齐研报图表18：29个初始特征 - `close/open1/high/low/vwap/volume/amt/turn/bp/ep/ocfp/dp/return1/adv20/adv60` × 原始 + 时序标准化（滚动60日z-score） > ⚠️ **风险提示**：本框架为研报复现，历史回测存在过拟合、幸存者偏差等风险，建议配合走样本外、滚动回测等方法验证因子有效性，不构成任何投资建议。

让我先读取现有代码，定位性能瓶颈所在：

已找到根本原因，现在直接定位慢的三处代码并完整重写： **🔴 性能瓶颈根因（3处）：** 1. **所有滚动算子**（`ts_min/max/std/sum/correlation` 等）用 `for t in range(T)` Python循环，T≈2400次，每次触发一次GPU同步 → **最慢** 2. **`compute_rank_ic`** 用 Python `for t in range(T)` 逐截面计算 → **第二慢** 3. **`rank_pct_torch`** 用 Python `for t in range(T)` 逐截面排名 → **第三慢**

现在读取完整 notebook 的代码单元格内容，以便精确替换：

现在同步更新 `eval_tree` 和主循环里的 `compute_rank_ic` 引用，并加入**并行批量评估**逻辑：

验证通过，三个核心Cell已全部完成向量化。现在来解释一下修改内容： --- ## 🔴 根本原因分析 原代码卡在「小循环1/3 计算适应度」，30分钟无结果的原因是**三个叠加的 Python 循环**： ``` 3000棵树 × 2400天 × Python-for-t循环 = 720万次 Python→GPU 调度 每次都触发 CUDA 同步 → GPU 利用率 <5%，CPU 成为瓶颈 ``` --- ## ✅ 三处核心修复（Cell[12] / Cell[16] / Cell[20]） ### 修复1：Cell[12] — 滚动算子，`for t in range(T)` → `torch.unfold()` | 对比 | 旧代码 | 新代码 | |------|--------|--------| | 实现方式 | Python `for t in range(2400)` | `tensor.unfold(0, d, 1)` 一次性展开 | | GPU调度次数 | 2400次 | **1次** | | 涉及算子 | `ts_min/max/std/sum/corr/cov/rank/argmin/argmax/decay/grouping` | 全部 | ```python # 旧：Python 循环（极慢） for t in range(d-1, T): window = X[t-d+1:t+1] # 逐步切片 out[t] = window.min(0).values # 新：unfold 一次展开（快100-2400倍） def _unfold(X, d): # [T,N] → [T,N,d] padded = torch.cat([nan_pad, X], dim=0) return padded.unfold(0, d, 1) # 零拷贝视图，无Python循环 ``` ### 修复2：Cell[16] — `compute_rank_ic`，逐截面 → 全时间维度向量化 ```python # 旧：T=2400 次 Python 循环 for t in range(T): fv = factor_vals[t]; yv = Y[t] rf = fv.argsort().argsort() # [N] # 新：dim=1 整体向量化，一次完成所有截面 rf = f_fill.argsort(dim=1).argsort(dim=1).float() # [T, N] 全量 ic = (rfm * rym).sum(dim=1) / den # [T] 一步出结果 ``` ### 修复3：`rank_pct_torch`，逐截面 → 全向量化 ```python # 旧：for t in range(T) 循环 # 新： order = X_fill.argsort(dim=1).argsort(dim=1).float() # [T,N] 直接完成 ``` --- ## ⚡ 预期性能提升 | 阶段 | 旧耗时估算 | 新耗时估算 | 提速 | |------|-----------|-----------|------| | 单棵树评估 | ~3-5秒 | **10-50ms** | **60-300x** | | 每小循环3000棵 | ~2.5小时 | **5-25分钟** | 同上 | | 全程(15轮小循环) | **37小时+** | **1-6小时** | — | > ⚠️ **注意**：分组算子（`ts_grouping_*`，d=40时展开为 [2400,3000,40]=288MB）在 **RTX 4090（24GB显存）上没有问题**，如遇OOM可将 `D_CHOICES_GRP = [10, 15, 20]` 去掉40。