企业官网建设流程全解析

做数据建模这行久了，会遇到一类很具体的需求：对方不要一份预测结果，要一个能抄在纸上、能贴进工艺卡、能塞进单片机里跑的公式。几年前我接过一个能耗标定的活，输入是六个现场采集的物理量，输出是单位产出的能耗值。我先上了梯度提升树，交叉验证的误差压到了很低，交付时对方问了一句“这个模型具体是怎么算的”，我打开特征重要性图讲了十分钟，对方的表情说明了一切——他们要的不是一个黑箱有多准，而是一条能写在设备铭牌旁边的经验式。后来我把方向转到了符号回归，用的就是gplearn，一个把遗传编程包装成 scikit-learn 风格的 Python 库。它做的事情一句话能说清：让计算机自己进化出一段数学表达式，这段表达式既能拟合你的数据，又能被人类读懂。

这篇内容我想写给三类人看。第一类是手上有回归任务、但被要求“结果必须可解释到公式级别”的从业者；第二类是做特征工程做到瓶颈、想找一种自动化构造非线性组合特征的工程师；第三类是对遗传编程好奇、翻过文档但被那一堆参数劝退的朋友。gplearn 的官方 API 文档其实写得不算差，但它不会告诉你哪些参数是真正影响结果的、哪些坑是第一次跑必然会踩的、以及为什么你跑出来的公式长得像一团乱麻。这些只能靠自己一段一段试出来，我踩过的部分在这里摊开讲。

1. 当一个“能抄在纸上的公式”成为硬需求，我才认真用起 gplearn

1.1 树模型给不了的东西：符号表达式

先说清楚符号回归到底在解决什么问题。常规回归模型是在一个固定的函数形式里找参数，比如线性回归锁定的是 $y = w_0 + w_1x_1 + \dots + w_nx_n$，你调的是那组权重；多项式回归锁定的是各次幂的组合；神经网络锁定的是层数和激活函数。模型的形式是人定的，机器只负责填数。符号回归换个思路：函数形式本身也是待搜索的变量，搜索空间是从常数、特征变量和一组基本运算（加减乘除、开方、对数、三角）能拼出来的所有表达式。

这就带来一个很有意思的性质——搜索空间是离散且无限的。你不能用梯度下降去优化“该不该在这里加一个平方项”，因为这是一个结构性决策，不是连续参数。遗传编程正是为这类离散结构搜索设计的：维护一个表达式种群，用适应度筛选，用交叉和变异产生新个体，一代代迭代下去。gplearn 就是这套思路的工程实现。

我常拿它和一个比喻对照：调参像是在给定户型里摆家具，符号回归则连户型一起重新设计。好处是可能找到你完全没想到的结构，坏处是搜索空间太大，容易搜出一堆虽然拟合得好但物理上荒谬的式子。这个矛盾贯穿了使用 gplearn 的全过程。

1.2 gplearn 的定位：遗传编程 + scikit-learn 的接口习惯

gplearn 的接口设计明显在向 scikit-learn 靠拢，fit、predict、transform、get_params、set_params一应俱全，能直接塞进Pipeline，也能用GridSearchCV做搜索。这个设计选择对实际工作帮助很大，因为大部分人的数据预处理链路已经用 scikit-learn 搭好了，接进来不需要重写。

它主要提供三个估计器：

三个里面我用得最多的是前两个，SymbolicTransformer属于“想到了就用一下”的工具。分类器的加入时间比回归晚不少，实际用下来它在特征维度不高、决策边界接近某个简单表达式时表现还行，但不要指望它能替代梯度提升树在复杂分类任务上的地位。

内部实现上，gplearn 用 numpy 数组表示每个表达式，节点存函数编号和子节点索引，执行时按拓扑序遍历。这个实现方式决定了它的性能特征：单条表达式的求值很快，但种群规模乘以代数之后总计算量会迅速膨胀，后面讲参数时会详细算这笔账。

1.3 什么样的任务适合交给它，什么样的不适合

这一点必须先说，因为用错场景是最大的时间浪费。适合的典型情形有这么几类：

• 需要交付公式的场景。工艺参数标定、设备能耗模型、传感器补偿系数、材料配比经验式，这类任务最终的产出物就是一条式子，模型精度反而不是第一优先级。
• 样本量不大但特征物理含义明确的场景。几百到几千条样本、特征在 5 到 20 维之间，这是 gplearn 最舒服的区间。特征太多它会陷入维度灾难，样本太多则算力吃不消。
• 想验证某个领域经验公式是否还有改进空间。你可以把它已有的公式作为先验，观察进化出来的结构是否接近。
• 自动构造交互特征。用SymbolicTransformer生成一批非线性组合，再喂给线性模型或者树模型。

反过来，不适合的情形同样明确。图像、文本、高维稀疏特征直接上 gplearn 基本是自找麻烦，因为进化的搜索空间随特征数指数膨胀。目标函数噪声极大、信噪比极低的场景也不合适，遗传编程很容易把噪声拟合成复杂结构。还有一类常见误用是把它当成“自动调参工具”，指望它帮你找到一个比现有模型更准的黑箱——如果你的目标只是精度，梯度提升树在绝大多数表格数据上都更省事。

提示：判断标准可以简化成一句话——如果你的交付物里必须出现一条人类能读懂的数学表达式，gplearn 值得投入时间；如果不需要，先考虑常规模型。

2. 动笔写 fit 之前，先把数据和算子集这两件事定死

2.1 安装与版本兼容的几个细节

安装本身没难度，pip install gplearn就行。但版本兼容有两处我踩过：

第一处是 numpy 版本。gplearn 内部大量使用np.float64和数组切片操作，numpy 2.x 早期版本下部分老版本 gplearn 会报类型相关的警告甚至错误。稳妥做法是把 gplearn 升到当前最新版，或者把 numpy 锁在 1.26 这一线。

第二处是 scikit-learn 版本。因为 gplearn 继承了BaseEstimator和RegressorMixin，sklearn 大版本变更时偶尔会出现_validate_data相关的兼容问题。我现在的习惯是在项目里固定一组版本，写进requirements.txt：

pip install "numpy<2" "scikit-learn>=1.3,<1.6" "gplearn>=0.4.2" "joblib>=1.3"

装完做个三十秒的冒烟测试，确认能导入、能跑通一个极小的例子，比后面排查环境问题划算得多：

import numpy as np from gplearn.genetic import SymbolicRegressor rng = np.random.RandomState(42) X = rng.uniform(-3, 3, size=(200, 2)) y = X[:, 0] ** 2 + X[:, 1] * 2.0 est = SymbolicRegressor( population_size=500, generations=10, function_set=('add', 'sub', 'mul', 'div'), random_state=42, verbose=0, ) est.fit(X, y) print(est._program)

如果这段跑出来打印了一个形如add(mul(X0, X0), mul(2.0, X1))的表达式，环境就没问题了。注意verbose=1会打印每一代的最优适应度，调试时很有用，但如果放在批处理脚本里记得关掉，否则日志会刷屏。

2.2 输入矩阵的三条硬性前提

gplearn 对输入数据有几个隐含要求，文档里提得不显眼，但不满足会出各种诡异结果。

第一条，不能有缺失值。它不会帮你做任何插补，遇到 NaN 会直接在表达式求值时把整条程序污染成 NaN，然后这个个体被淘汰，你可能观察到种群适应度长期停在初始水平。跑之前务必X.isna().sum()检查一遍。

第二条，特征最好先标准化。这条是重点。为什么？因为 gplearn 的默认算子集里包含add和sub。假如你有两个特征，一个是价格（几千到几万），一个是比例（0 到 1），那么add(X0, X1)这个组合的方差几乎完全由 X0 主导，进化过程会发现把它加进来对拟合有帮助，于是大量个体都往这个方向跑，浪费了本该用于探索真正交互结构的预算。更麻烦的是量纲混在一起时，生成的公式在物理上讲不通，交付时会被人一眼看出问题。

我一般统一做 Z-score 或者 min-max，具体看特征分布。标准化带来的一个副作用是表达式里的系数需要做反变换，这一点后面在公式交付那节会讲怎么处理。

第三条，样本要足够覆盖特征空间。遗传编程是靠适应度筛选的，如果某个特征在训练集里只在很小的范围内变化，进化出的表达式在这个特征上的行为基本是外推，验证时容易翻车。我习惯在训练前对每个特征画个直方图，看看分布是不是严重偏斜，偏斜的做对数变换。

2.3 function_set 不是越多越好：算子选择的量纲逻辑

function_set是新手最容易“塞满”的参数。看着有十几个内置函数，手一抖全填进去，结果就是搜索空间炸开，收敛变慢，最后跑出来一个sin(log(abs(div(X3, X5))))这种完全没法解释的东西。我现在的原则是从最小的算子集开始，按需增加。

内置算子分两类。二元算子有add、sub、mul、div、max、min，一元算子有sqrt、log、abs、neg、inv、sin、cos、tan、sig。其中几个是带保护机制的：

保护机制让程序不会因为一次除零就整个崩掉，这是好事；但它同时意味着错误被静默吞掉了。这是后面要单独讲的坑之一。

选择算子的实用建议：

• 只在明确知道目标关系包含周期性时加sin、cos。加进去之后基本一定会被进化过程滥用，因为三角函数在某些区间内几乎可以拟合任何平滑曲线。
• 物理量通常非负，可以放心用sqrt；如果是已经中心化过的特征，sqrt会大量触发保护返回 0，等于白占一个搜索方向。
• div几乎必开，比值关系在工程里太常见了。
• 如果已知目标单调，考虑把abs去掉，它能构造出对称结构但会引入不可导点。
• max、min用得少，但在有分段特性的场景有价值，比如“取两个限制因素中更紧的那个”。

我还试过一个技巧：如果对目标关系有先验，比如知道一定是若干个变量的乘积形式，那就把add、sub从算子集里去掉，只留mul、div，强制进化搜索乘性结构。这种“缩小搜索空间换收敛速度”的做法，在特征少的时候效果很明显。

另外 gplearn 支持自定义一元和二元函数，通过_Function包装：

import numpy as np from gplearn.functions import _Function from gplearn.genetic import SymbolicRegressor def _square(x): return np.square(x) square = _Function(function=_square, name='square', arity=1) est = SymbolicRegressor( population_size=1000, generations=20, function_set=('add', 'sub', 'mul', 'div', square), random_state=0, )

自定义函数对收敛的帮助有时很大——把“平方”这种高频出现的结构直接作为一个原子操作，等于告诉搜索器“这条路是通的”，省下用mul(X0, X0)两步拼出来的成本。注意要求是向量化函数，接收和返回都是 numpy 数组。

3. 三个核心参数决定成败：种群、代数与简约系数

3.1 population_size 与 generations 的算力换算法

这两个参数一起决定了总计算量，而 gplearn 的计算量是可以粗算的。基本公式是：

总求值次数 ≈ population_size × generations

每一次求值要在max_samples指定比例的样本上执行一遍表达式树。默认population_size=1000、generations=20，就是两万次求值。听起来不多，但单次求值的成本取决于表达式深度和样本量。假设你有 2000 条样本、平均树深 10 层，那单次求值大约是几十次 numpy 数组运算，两万次乘起来就是几十万次数组操作——单机跑十几分钟到半小时很正常。

population_size控制的是每代内的多样性。太小（比如 200）会导致种群过早收敛到局部最优，十几代之后所有个体长得差不多；太大（比如 20000）则每代进化缓慢，同样算力下代数不够，一样得不到好结果。我的经验区间是1000 到 5000，特征维度高、搜索空间大的时候往 5000 走，特征只有四五个的时候 1000 就够了。

generations控制迭代深度。太少看不到收敛，太多则后期基本停滞，纯浪费算力。判断方法是打开verbose=1，观察每一代最优适应度的变化：

• 前 5 到 10 代通常快速下降；
• 之后进入平台期，偶尔有突破；
• 如果连续 10 代以上完全没有改善，可以停了。

我一般先设generations=20跑一轮摸底，看收敛曲线在第几代开始走平，再决定要不要加到 50 甚至 100。盲目设generations=500是新手最常见的浪费——你可能要等好几个小时，最后拿到的公式和跑到 30 代时一模一样。

还有一个stopping_criteria参数，当某代最优个体的适应度低于这个阈值时提前停止。这个参数要小心用：它比较的是metric计算出的适应度值，如果你用默认的'pearson'（注意 gplearn 会把它转成 1 减去相关系数作为适应度，越小越好），那阈值设 0.01 意味着要求相关系数达到 0.99，很多真实数据根本达不到，等于这个参数永远不生效。反而是设得太宽松（比如 0.9）会让它跑两代就停。我通常先关掉它（默认值 0.0 相当于永不触发），等摸清数据能达到的水平后再设。

3.2 parsimony_coefficient：给公式膨胀的刹车片

这个参数我认为是 gplearn 里最值得花时间调的。它的作用是给适应度加一个复杂度惩罚项：

实际适应度 = 原始适应度 + parsimony_coefficient × 表达式长度

表达式越长，惩罚越大。没有这个惩罚，遗传编程会不可抑制地膨胀——因为更长的表达式总能拟合得更好一点点，而交叉操作又天然倾向于产生更长的后代。跑几十代之后你会得到一棵深度十几层、包含几十个节点的巨型表达式，在训练集上拟合完美，在验证集上一塌糊涂。

parsimony_coefficient的默认值是 0.001，这个值对于不同的问题差异很大。取值逻辑是相对metric的量级：如果你用mse，误差可能上百甚至上千，那 0.001 的系数相对于长度惩罚几乎可以忽略；如果你用pearson，适应度在 0 到 2 之间，那 0.001 就相对明显了。所以我一般先看原始适应度的量级，再定系数。

'auto'这个选项值得单独说。它会根据每代种群中个体长度与适应度的协方差动态计算惩罚系数，本质上是让惩罚强度自适应种群状态。我实测下来它在多数任务上给出的结果与手调的合理值接近，适合作为探索阶段的默认值，但如果最终要交付，还是建议手动调几档对比一下。

调这个参数的最佳观察方式是画出精度-复杂度权衡曲线：固定其他参数，扫几个系数值，记录每个值下的验证集误差和公式节点数。理想的结果应该落在拐点上——再增大系数，误差明显上升，而节点数下降不多。

3.3 交叉变异概率与锦标赛选择：进化到底有没有真的在发生

gplearn 提供了三种遗传操作，通过概率参数控制：

est = SymbolicRegressor( population_size=2000, generations=30, p_crossover=0.7, p_subtree_mutation=0.1, p_hoist_mutation=0.05, p_point_mutation=0.1, p_point_replace=0.05, tournament_size=20, max_samples=0.9, parsimony_coefficient=0.001, random_state=42, n_jobs=4, verbose=1, )

p_crossover是交叉概率，从两个父代各取一段子树互换。这是探索的主力，负责组合已有结构。p_subtree_mutation是子树变异，随机生成一棵新子树替换掉原树的一个分支，负责引入全新的结构。p_hoist_mutation是把某个子树提上来替代整棵树，作用是缩短表达式，我觉得它和 parsimony 系数是互补的两道刹车。p_point_mutation是点变异，随机替换树上的单个节点（换函数或换变量），负责局部微调。

这五个概率之和必须为 1，否则会报错。默认值（0.9 / 0.01 / 0.01 / 0.01）偏向交叉，探索性强但收敛慢。我常用的一组是 0.7 / 0.1 / 0.05 / 0.1，把变异的比重提上来一些。理由是：当初始种群质量不高时，纯靠交叉很难跳出局部区域，需要足够的变异来维持多样性。

tournament_size是锦标赛选择的规模，默认 20。每代选出新父代时随机抽这么多候选，挑适应度最好的那个。这个值越大，选择压力越强，收敛越快但也越容易早熟。种群规模是 1000 的时候，20 差不多是 2% 的采样比例；如果种群规模只有 200，20 就占了 10%，压力偏大，我一般按种群规模的 1% 到 3% 来设。

p_point_replace是点变异时每个节点被替换的概率，通常保持默认 0.05 就好，调它的收益不大。

有一件事值得观察：如果跑了十几代之后你会发现种群中适应度最好的若干个体结构高度相似，说明多样性已经耗尽。gplearn 不会直接暴露种群结构，但你可以通过est._programs拿到每一代的种群列表，统计一下不同表达式字符串的数量：

import numpy as np last_gen = est._programs[-1] unique_programs = {str(p) for p in last_gen if p is not None} print(f"末代种群规模 {len(last_gen)}，去重后 {len(unique_programs)}")

如果去重后的数量不到种群规模的二十分之一，基本可以判定早熟了，该考虑加大变异概率或者缩小tournament_size。

4. 一次从零到出公式的完整实操复盘

4.1 构造一个可验证的模拟数据集

真实验证过 gplearn 能力的做法是先造一个你知道答案的数据集，看它能不能把答案找回来。我用一个含四个变量、带噪声、含一个交互项和一个比值项的目标函数：

$$y = 2.5 \cdot x_1 + x_2 \cdot x_3 - \frac{x_4}{x_1 + 4} + \varepsilon$$

这里故意设计了几种不同结构：线性项、乘积交互、比值项，并且分母加了常数避免除零。

import numpy as np from gplearn.genetic import SymbolicRegressor from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler rng = np.random.RandomState(42) n = 1200 X = rng.uniform(0.5, 5.0, size=(n, 4)) y = 2.5 * X[:, 0] + X[:, 1] * X[:, 2] - X[:, 3] / (X[:, 0] + 4.0) y += rng.normal(0, 0.05, size=n) X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.25, random_state=42 ) scaler = StandardScaler() X_train_s = scaler.fit_transform(X_train) X_test_s = scaler.transform(X_test) est = SymbolicRegressor( population_size=3000, generations=30, function_set=('add', 'sub', 'mul', 'div'), metric='mse', parsimony_coefficient=0.0005, p_crossover=0.7, p_subtree_mutation=0.1, p_hoist_mutation=0.05, p_point_mutation=0.1, max_samples=0.9, feature_names=['x1', 'x2', 'x3', 'x4'], random_state=42, n_jobs=4, verbose=1, ) est.fit(X_train_s, y_train) print("最优表达式：") print(est._program) print(f"训练集 R2: {est.score(X_train_s, y_train):.4f}") print(f"测试集 R2: {est.score(X_test_s, y_test):.4f}") print(f"节点数: {est._program.length_} 树深: {est._program.depth_}")

跑这个例子的感受是：如果算子集只给四则运算、parsimony_coefficient设得合理，通常 20 代左右就能找到一个结构接近真实的表达式，测试集 R² 能到 0.98 以上。这算是个比较理想的情况，因为目标关系确实由这四个算子构成。

4.2 训练过程中的三件观察

第一件，verbose=1输出的适应度曲线不是单调下降的。你会看到某几代的最优值反而变差了。这不是 bug，是因为每代的“最优”是在新一代种群中选出来的，而新种群由上一代的采样重建，理论上被精英保留的个体应该保留下来——gplearn 确实做了精英保留，但parsimony_coefficient='auto'模式下惩罚系数逐代变化，会导致同一条程序的适应度在不同代之间不可比。所以看曲线时关注长期趋势就好，别纠结单代的抖动。

第二件，max_samples小于 1.0 时会出现oob_fitness_属性。这是袋外适应度，思路借自随机森林：每代对每个个体在未参与训练的样本子集上评估一次。这个值比训练适应度更能反映泛化能力，我在选最终表达式时会同时看两个——如果某个个体的训练适应度极好但袋外适应度很差，基本可以判定过拟合，弃之。

第三件，迭代后期会出现“结构等价但写法不同”的个体竞争。比如div(X0, X1)和mul(X0, inv(X1))在数学上等价，适应度几乎一样，但节点数差一个。这种情况下parsimony_coefficient会帮你选短的，不过前提是系数设得足够大。

4.3 把 program 对象翻译成能调用的函数

est.predict()能直接用，但交付场景需要的是独立的公式。est._program打印出来是类似add(mul(x2, x3), sub(mul(2.5, x1), div(x4, add(x1, 4.0))))这样的字符串。要把它变成可调用函数，最直接的方式是用eval构造：

import numpy as np from gplearn.functions import _function_map def compile_program(program, feature_names): """把 gplearn 的 program 对象编译成可调用函数""" expr = str(program) env = dict(_function_map) # 内置函数的 numpy 实现 env.pop('sqrt', None) # 按需清理，避免与 numpy 冲突 for idx, name in enumerate(feature_names): env[name] = f"X[:, {idx}]" # 用 eval 展开变量名，再构造 lambda expanded = expr for idx, name in enumerate(feature_names): expanded = expanded.replace(name, f"X{idx}") src = f"lambda X: {expanded}" return eval(src, {"__builtins__": {}}, dict(_function_map))

更稳妥的做法是直接使用 gplearn 内部提供的执行机制，因为_function_map里的函数签名和表达式里的调用形式完全一致：

from gplearn.functions import _function_map def build_callable(program, feature_names): expr = str(program) ns = {k: v for k, v in _function_map.items()} # 把特征名映射成从矩阵取列的 lambda for idx, name in enumerate(feature_names): ns[name] = (lambda i: (lambda X: X[:, i]))(idx) return eval(expr, {"__builtins__": {}}, ns) f = build_callable(est._program, ['x1', 'x2', 'x3', 'x4']) y_hat = f(X_test_s)

这里有个坑要提前说：_function_map是 gplearn 的内部对象，以下划线开头的模块属性在版本升级时可能变动。如果你要长期维护，最好自己写一份对应的实现，把div、log、sqrt、inv的保护行为原样复刻，否则线上跑出来的结果和训练时不一致。我见过一次线上预测和离线预测差的离谱，查到最后一比对，保护除零返回的值一个是 1.0 一个是 0.0。

另外，标准化之后的公式需要还原到原始量纲。如果用了 Z-score，还原方式是把每个x_i替换成 $(X_i - \mu_i)/\sigma_i$，然后展开多项式。这个展开在公式节点数超过 30 的时候会非常难看，实际交付时我宁愿用原始数据重新训练一遍，也不要交一个人工展开过的巨型公式。

5. 那些让公式看起来漂亮却毫无用处的坑

5.1 训练集 R² 0.99 而验证集崩盘：过拟合的三种形态

过拟合在符号回归里的表现和树模型不太一样，它有几种典型形态，识别方法各不相同。

第一种是结构过拟合。表达式里出现了大量细小常数，比如add(mul(1.037, X1), mul(-0.021, X2))这种。这些系数看起来精确，其实是进化过程在训练集噪声上做微调的结果。识别特征是常数位数的异常增长。解决办法是限制const_range，比如从(-1.0, 1.0)收窄到(-0.1, 0.1)，或者干脆设成None关掉常数生成——为什么关掉反而可能更好？因为关掉之后表达式只能由特征和算子组成，即使精度略降，结构往往更接近真实的物理关系。

第二种是深度过拟合。树深超过 15 层，节点数超过 50。gplearn 默认的init_depth=(2, 6)只控制初始深度，随着交叉和变异，树会一代代长起来，没有硬性上限。这类表达式在训练集上表现完美，但换个样本就崩。识别方式是看_program.depth_和测试集表现的相关性。

第三种是特征滥用。表达式里某个特征只出现一次，且系数极小，对整体输出几乎没影响，但让适应度提高了那么零点几个千分点。这种“装饰性节点”在节点数多的表达式里很常见，人工审查能发现，但这意味着你必须人工审查。

对抗过拟合的组合拳我总结成一张表：

其中max_samples的作用容易被低估。设成 0.8 意味着每个个体在 80% 随机子样本上评估，这天然引入了正则化，同时提供了oob_fitness_这个独立的泛化指标。我现在的默认配置基本都会开到 0.85 左右。

5.2 数值溢出与保护算子的静默失败

保护算子是个双刃剑，前面提过，这里展开讲。div在分母为 0 时返回 1.0，这看起来安全，但后果是这个错误完全不可见。我遇到过一个典型场景：目标值有明确的物理下界，进化出来的表达式里有个div(X2, sub(X1, X3))，而 X1 和 X3 在某些样本上恰好相等。结果那些样本的贡献被强行拉成 1.0，模型在这个区域完全失真，但因为样本量少，整体指标看不出来。

排查方法很土但有效：在测试集上检查残差的空间分布，看有没有局部异常峰。更主动的做法是手工算一遍表达式在每个样本上的中间值，找出接近保护触发条件的输入组合。gplearn 不提供单节点的中间输出，所以要自己用编译出的可调用函数分段计算。

另一个更隐蔽的问题是浮点溢出。inv(inv(inv(X1)))这种嵌套在 X1 接近 0 的时候会迅速变成天文数字，虽然是inv有保护，但三次嵌套中间层可能已经溢出成inf。我的应对是在编译出的可调用函数外面包一层：

def safe_callable(f, X, clip=1e6): with np.errstate(over='ignore', invalid='ignore', divide='ignore'): out = f(X) out = np.nan_to_num(out, nan=0.0, posinf=clip, neginf=-clip) return np.clip(out, -clip, clip)

这种处理当然是“打补丁”，但对交付场景来说，一个有界的错误输出比一个 NaN 更可控。

还有一个我一直没完全解决的问题：sqrt和log在输入为负时返回 0.0，这会让表达式的梯度信息在那一部分彻底消失。如果负数输入占比高，实际上相当于模型在那些区域是瞎猜的。处理方式是训练前确认特征的符号分布，全非负的特征才放心让sqrt、log进去。

5.3 特征共线性让进化原地打转

这个坑我吃了不止一次亏。当两个特征高度相关（相关系数超过 0.95）时，2.3 * X1和2.4 * X2的适应度几乎一样，进化过程会在两者之间反复横跳，消耗掉大量的搜索预算。表现出来就是适应度曲线在前几代快速下降后长期停滞，你换随机种子重跑，得到的是结构相似但换个变量的公式。

判断方法很简单，跑之前做个相关性矩阵，或者算一下 VIF。发现高相关就对，要么删一个，要么用 PCA 合成。我倾向于删，因为符号回归的价值就在于每个变量有物理含义，合成之后就丢了。

顺带说一个相关的问题：常数列和方差极小的特征要清理掉。如果某个特征在训练集里几乎是常数，进化过程会把它当成一个免费的常数项来用，add(mul(X1, 2.0), X5)里的 X5 其实只贡献了一个偏移量。这种表达式在训练集里看起来很合理，但换个时间段重新采集数据，X5 的方差变了，整个公式就失效。

6. 把它塞进真实流程：特征生成、公式交付与模型互补

6.1 用 SymbolicTransformer 造非线性特征

SymbolicTransformer的用途是生成一批表达式，把它们当作新特征列输出给下游模型。这个用法我认为被严重低估了——很多时候你不追求一个完美公式，只是想让线性模型能捕捉到一些非线性交互。

from gplearn.genetic import SymbolicTransformer function_set = ('add', 'sub', 'mul', 'div', 'sqrt', 'log', 'inv') gp = SymbolicTransformer( generations=20, population_size=2000, hall_of_fame=100, n_components=8, function_set=function_set, parsimony_coefficient=0.001, max_samples=0.85, metric='pearson', n_jobs=4, random_state=42, verbose=1, ) X_new = gp.fit_transform(X_train_s, y_train) print(f"新特征维度: {X_new.shape}")

关键参数是hall_of_fame和n_components。前者是进化结束后保留的最优个体数量，后者是从中挑出的、彼此相关性最低的若干表达式（保证新特征之间不冗余）。n_components必须小于等于hall_of_fame，这个约束不满足会报错。

实战心得是：n_components不要设太大，5 到 10 个足够。我试过设 20，结果下游模型提升有限，反而因为特征变多增加了过拟合风险。另外hall_of_fame设大一点（100 到 200）有助于n_components挑出真正互不相关的组合。

这个特征生成器的正确用法是把它放在 Pipeline 里，用交叉验证评估整体效果：

from sklearn.pipeline import Pipeline from sklearn.linear_model import RidgeCV from sklearn.preprocessing import StandardScaler pipe = Pipeline([ ('scaler', StandardScaler()), ('symbolic', gp), ('scaler2', StandardScaler()), ('ridge', RidgeCV(alphas=np.logspace(-3, 3, 13))), ]) pipe.fit(X_train, y_train) print(pipe.score(X_test, y_test))

注意gp是有监督的，fit_transform需要传y，所以它在 Pipeline 里必须放在能接触到标签的位置，这没问题，因为 Pipeline 的fit会把y传下去。

6.2 公式交付时的工程化处理

交付一条公式和交付一个模型文件是完全不同的两件事。我的做法是产出一份包含三部分的文档：表达式字符串、变量定义与单位、误差范围说明。表达式用 gplearn 的原始字符串形式给出，同时附上编译好的 Python 函数和一段验证脚本。

有一个细节很关键：公式里的常数要四舍五入到合理的有效位数。gplearn 生成的常数可能是1.0372914，这个精度在物理上是虚的，而且会让阅读者误以为这个系数很重要。我一般统一截断到三到四位有效数字，然后重新验证一遍误差变化——如果截断后误差上升超过 5%，说明这个系数确实在拟合噪声，应该回到训练阶段调大parsimony_coefficient。

还有一点，公式变量顺序要和业务口径对齐。gplearn 用X0、X1这样的占位名，如果训练时传了feature_names，表达式里会用真实名字，这大大减轻了交付时的对号入座工作。我强烈建议每次都传feature_names，尤其是特征超过五个的时候。

6.3 和主流模型放在一起比较的正确姿势

我常被问“gplearn 和梯度提升树哪个准”。这个问题本身就没问到点上，因为它们的优化目标不同。正确的比较方式是分三个维度：

外推这一行要重点说。树模型的预测值永远落在训练集目标值的区间内，出了这个区间只能给边界值。符号回归则可能给出合理的外推，也可能给出荒谬的外推。如果表达式里有div且分母可能趋近 0，外推会爆炸；如果表达式是线性或者低次多项式，外推就相对可靠。所以用 gplearn 做外推预测前，必须检查表达式结构，把含保护算子的分支单独标注出来。

我实际工作中的做法是两者并行：用梯度提升树给出精度基线，用 gplearn 给出可解释结构，然后比较两者的残差相关性。如果残差高度相关，说明符号回归抓住了主要结构，可以放心交付；如果残差几乎不相关，说明还有大量信息没被符号表达式捕捉，这时候要么加宽算子集，要么老老实实承认这个任务不适合符号回归。

7. 让它在有限算力里跑完：并行、缓存与调试手段

7.1 n_jobs 与 low_memory 的取舍

gplearn 的并行粒度是“每个个体一个任务”，用 joblib 分发。n_jobs设多少合适取决于两件事：一是单个任务的计算量，二是任务调度开销占比。当样本量大、表达式深时，单任务耗时长，调度开销可以忽略，n_jobs可以设到物理核心数。反过来，样本只有几百条时，单任务可能只花几毫秒，调度开销反而占大头，这时候多进程可能比单进程还慢。

我实测的经验是：样本量超过 1000 条、表达式平均深度超过 8 层时，n_jobs设成核心数的 0.75 倍收益最明显；轻量任务建议n_jobs=1老实跑。

low_memory是另一个值得开的开关。默认False时，gplearn 会保存历史每一代的完整种群对象，种群规模 5000、跑 50 代就是 25 万个对象常驻内存，很容易把内存吃满。设成True之后只保留当前代和上一代（因为交叉需要父代），内存占用能降一个数量级。代价是跑完之后你没法回溯历史种群做分析。我的习惯是调参阶段开low_memory=False方便分析，最终训练时开True。

7.2 warm_start 续跑与随机种子的复现性

warm_start=True允许从上一次fit的最终种群继续进化。这个功能在做长跑实验时非常有用——你可以先跑 20 代看看趋势，觉得有希望再加 30 代，而不用从头再来。用法很简单：

est = SymbolicRegressor(..., warm_start=True, random_state=42) est.fit(X_train_s, y_train) print(est._program, est.score(X_test_s, y_test)) est.set_params(generations=est.generations + 30) est.fit(X_train_s, y_train) # 接着上一代继续 print(est._program, est.score(X_test_s, y_test))

要注意set_params里generations是累加还是覆盖，取决于版本实现，稳妥做法是显式设成“当前已跑代数 + 新增代数”。另外 warm_start 之后_programs列表会继续追加，不会清空，分析时注意索引。

复现性方面，random_state只要固定了，单进程下的结果完全可复现。但切到多进程后要小心：gplearn 内部用check_random_state生成种子序列并按个体索引分配，理论上仍然可复现，但如果你的 numpy 版本和 joblib 后端有差异，随机数流的消耗顺序可能变化。我的做法是在实验记录里同时记下random_state、numpy 版本、gplearn 版本三样东西，缺一样都可能查不出差异来源。

7.3 快速定位“为什么它不收敛”的排查清单

跑了几十次之后我整理出一份排查顺序，基本上能覆盖九成的情况。按这个顺序检查，比盲目调参数高效得多：

还有一条不在这张表里但经常遇到：目标值的分布严重偏斜。如果 y 是长尾分布，用mse作为 metric 会让少数大值样本主导适应度，模型全力拟合那些极端点，普通样本反而拟合不好。处理方式是先对 y 做变换（对数、Box-Cox），训练完再反变换。这个技巧在能耗、金额、时长这类目标上几乎每次都要用。

我个人在实际操作中的体会是，gplearn 这类工具的价值不在于它能给出多高的精度，而在于它迫使我重新审视“数据里到底存在什么结构”。当你看着它一代代进化出各种各样的表达式时，会不断冒出“这个组合居然还有一定解释力”的念头，这些念头往往比最终那条公式更有价值。至于那些参数，调多了会发现真正需要反复斟酌的其实只有三个——算子集、简约系数和采样比例，剩下的调来调去差别都在噪声范围内。