简介:基于遗传编程(GP)实现符号回归的资料包,面向计算机科学研究者、机器学习爱好者以及初学进化计算的技术从业者,重点解决非线性函数发现、数学表达式建模等复杂问题,强调从基础实现到方法优化再到成果展示的完整链路。压缩包内为1个docx文档,体积仅13KB,但内容安排紧凑:不仅给出符号回归GP的基础任务与增强/修改方向,也明确了精英主义、选择算子、启发式策略等可加分项,并对类型化问题、可视化呈现、统计分析和LaTeX报告提出了评分标准。文档以Python作为实现语言,还提供了GitHub上的CSV数据样本与参考文献,可直接用于实测。目前已有118人学习。读者按步骤推进,既能独立搭建适用于不同问题的GP系统,也能掌握高维数据可视化、算法性能评估和高水平科技写作等实用技能,适合课程作业与自主实战。
1. 符号回归不是“调参”,而是让进化算法自己找公式
拿到一份 CSV,最后一列是因变量,前几列是特征,你并不知道数据背后的生成函数——这是符号回归的典型起手式。传统做法是先假设 y = βX + ε,再用最小二乘去估 β;但如果真实关系带 sin、乘积、分段逻辑,线性假设很快暴露出模型偏差。遗传编程(GP)把候选公式编码成表达式树,靠选择、交叉、变异在函数空间中搜索,让模型结构自行涌现。这个任务虽然看起来像是机器学习入门作业,实际跑起来后你会遇到除零、符号爆炸、早熟和过拟合这些教科书不会写给你的问题。下文先讲编码,再给出可运行的 Python 实现,最后聊增强方案和能拿去报告的收尾技巧。
2. 符号回归问题的编码:树形表达、原始集与适应度压力
在动手写遗传编程之前,先想清楚搜索空间长什么样。符号回归的目标不是确定一组系数,而是找到一个从输入向量到输出标量的函数 f(x0, x1, …, xn)。GP 将候选函数表示为一棵表达式树,树的内部节点是运算符,叶子节点是变量或常数;每个个体对应一棵树,也对应一个完整的候选模型。这种编码方式的好处是:它天然覆盖复合函数、嵌套结构和任意深度,不需要像神经网络那样预先指定层数或激活类型。代价是搜索空间巨大,且存在大量语义等价的冗余表达式,比如 (x*1)、 (x+0)、 (x-x)。因此,后面的适应度函数和遗传算子都需要对冗余与膨胀做约束。
2.1 树形编码:一棵表达式树就是一个候选模型
以二维输入 x0、x1 为例,(x0 + x1) * x0 可以表达为根节点mul、左子树add(x0,x1)、右子树x0。树越深,能表达的复杂度越高,但搜索难度和计算成本也随之上升。GP 初始化通常有两个经典算法:full 和 grow。full 生成所有叶子深度一致的树,grow 允许叶子出现在不同深度;实际使用中会用genHalfAndHalf折中,一半 full、一半 grow。这个策略在 DEAP 里用genHalfAndHalf一行解决,但理解它有助于你解释为什么初始种群里已经有一部分“长得像”合理公式的个体。
# 用 DEAP 定义原始集(PrimitiveSet) pset = gp.PrimitiveSet("MAIN", arity=4) pset.addPrimitive(np.add, 2, name="add") pset.addPrimitive(np.subtract, 2, name="sub") pset.addPrimitive(np.multiply, 2, name="mul") pset.addPrimitive(protected_div, 2) pset.addPrimitive(np.sin, 1) pset.addPrimitive(np.cos, 1) pset.addEphemeralConstant("rand_const", lambda: random.uniform(-1, 1))这里arity=4表示有四个输入变量,DEAP 默认把它们命名成ARG0到ARG3;随后你可以用renameArguments改成x0到x3,方便呈现。protected_div是自定义函数,不能直接用/,否则遇到除零会产生inf或nan,进而让适应度比较失效。至于addEphemeralConstant,它会在每次生成个体时采样一个随机常数,搜索过程中常数是固定的,但同一个原始集里可以演化出多个不同常数。这是符号回归里很重要的一步:如果不加入随机常数,最终公式就只能由变量和算子组成,系数耦合在结构里,很难搜索出例如2.7 * sin(x0)这样带缩放因子的大数。
2.2 原始集与终止集的选择决定了搜索空间
原始集和终止集直接划定“GP 可能发现什么”。你放入sin、cos,它就有机会发现周期性;放入exp和log,就能覆盖指数增长或对数衰减关系;不放任何保护算子,就会因为异常值把种群污染掉。下面这张表是我在处理回归数据时常用的一组配置,也是这份课程数据(来自cs4XX-EvolutionaryComputation仓库)里足够通用的起点。
| 类别 | 元素 | 说明 |
|---|---|---|
| 二元算子 | add, sub, mul, protected_div | protected_div 在 |
| 一元函数 | sin, cos, exp, log | 视数据特征选择;log 也需要保护 |
| 终止符 | x0, x1, x2, x3 | 输入特征,数量由 CSV 列数决定 |
| 随机常数 | uniform(-1, 1) | 常数采样范围影响系数搜索效率;范围过大会让搜索变慢 |
| 深度限制 | 初始化 max_=3,变异 max_=2~4 | 防止一开始就生成不可读的深树 |
这不是一个“加上去就更好”的清单。原始集越大,搜索空间指数级变大;原始集太小,可能漏掉真实函数。比如一个由x1*x2 + x3生成的数据,如果没有乘法算子,GP 再进化多少代也找不到正确答案。在做任务时,可以先对不同实例做特征相关分析,看数据是全为正、有无负值、是否周期性波动,再决定是否把sin/cos加回去。对多元实例,我一般先把add/sub/mul/protected_div放进去跑一轮,然后用sympy把最优个体打印出来,判断是否缺少某类算子,再决定是否扩充原始集。
2.3 适应度函数要处理异常、复杂度和多目标
适应度函数在 GP 里不只是误差度量,它还承担着引导搜索方向的作用。最直接的是均方根误差(RMSE):每个个体在训练集所有样本上做完前向计算,求预测值与真实值的 RMSE,越小越好。如果数据有噪声,RMSE 是合理选择,因为它对离群点比 MAE 更敏感,能把“偏离大片点”的个体压下去。除了误差,还要加一个树长惩罚项,否则 GP 很容易出现(x0+x1)*((x0-x0)+x0)这类冗余结构;树长惩罚的强度由一个parsimony_coefficient控制,通常设成 0.001~0.01 之间。另一个思路是把“模型复杂度”作为第二目标,用 NSGA-II 做多目标演化,但实现和工作量会明显增加,评分上的收益需要权衡。
此外,如果训练数据分成了 train/validate/test,GP 内部只应使用训练集计算适应度,验证集负责判断是否早停,测试集留到最后报告泛化误差。因为符号回归很容易过拟合:树一旦变深,几乎可以记住所有训练样本,包括噪声。我见过一个项目在训练集 RMSE 降到 0.01,但测试集误差高达 0.8,问题就出在没有对复杂度做正则化。实际报告里,至少要在适应度函数里同步记录 len(individual),并在每代输出最优个体的训练误差、验证误差、树长,这样才好展示“增长-过拟合”的拐点。
3. 用 Python 实现 GP 求解器:进化循环写清楚,剩下交给运气
第2章的表只是定义了搜索空间,真正让 GP 跑起来的是“初始种群 → 评估 → 选择 → 交叉变异 → 替换”这个循环。这里我直接给出一个基于 DEAP 的可运行版本,只保留核心部分,适合作为课程作业的起点。选择 DEAP 不是因为仓库里的实现不可用,而是它把树表示、交叉变异原语都封装好了,你可以把精力放在算子策略和后续增强上;如果你希望完全从零实现,本质也是拆开toolbox.mate和toolbox.mutate,自己在PrimitiveTree节点列表上做交换。
3.1 初始化与评估函数
先安装依赖:pip install numpy pandas deap matplotlib sympy scikit-learn。DEAP 的gp模块负责树的生成、编译和遗传操作,base和creator负责适应度与个体类型的运行时构造。个体类型要在主模块顶部定义,因为 DEAP 的creator会在运行时创建一个类;如果你把它放进函数里,多进程并行时会重新创建类导致冲突。
import random import numpy as np import pandas as pd from deap import base, creator, tools, gp # 保护除法:分母绝对值小于阈值时返回 1.0 def protected_div(x, y): if abs(y) < 1e-6: return 1.0 return x / y # 读取数据:假设最后一列是 y data = pd.read_csv("instance.csv") X = data.iloc[:, :-1].values y = data.iloc[:, -1].values n_features = X.shape[1] pset = gp.PrimitiveSet("MAIN", arity=n_features) pset.addPrimitive(np.add, 2, name="add") pset.addPrimitive(np.subtract, 2, name="sub") pset.addPrimitive(np.multiply, 2, name="mul") pset.addPrimitive(protected_div, 2) pset.addPrimitive(np.sin, 1) pset.addPrimitive(np.cos, 1) pset.addEphemeralConstant("rand_const", lambda: random.uniform(-1, 1)) pset.renameArguments(**{f"ARG{i}": f"x{i}" for i in range(n_features)}) creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)代码前几行很直白:protected_div处理除零;PSet的定义决定了可搜索的算子集合。weights=(-1.0,)表示单目标最小化;如果你的增强版本引入两个目标,可以写成weights=(-1.0, -0.1),第二维对应复杂度,但那样选择算子也要切换成selNSGA2。
评估函数里需要把个体编译成可调用函数。DEAP 的toolbox.compile(expr=individual)会把PrimitiveTree转成一个普通 Python 可调用对象,输入是n_features个位置参数。然后对全部训练样本做列式计算:这里不要用 Python 显式循环,先转成 NumPy 数组再向量化,否则种群规模一大,每代评估会慢到让你怀疑人生。
def eval_sr(individual, X, y, parsimony_coef=0.005): func = toolbox.compile(expr=individual) # 向量化预测:每个个体一次算完所有样本 y_pred = np.array([func(*row) for row in X]) # 处理 inf/nan:如果出现非法值,直接给一个很大的惩罚 if not np.all(np.isfinite(y_pred)): return (1e12,) rmse = np.sqrt(np.mean((y_pred - y) ** 2)) return (rmse + parsimony_coef * len(individual),)len(individual)在 DEAP 中返回树的节点总数,用它作为结构复杂度惩罚项。如果预测值出现inf或nan,直接返回1e12而不是报错,这样这些个体会在锦标赛选择中被快速淘汰,但不会让程序崩溃。parsimony_coef需要反复试,太大个体会倾向用x0这类短表达式,太小又是深度膨胀;后面第4章会给一个经验范围。
3.2 进化循环与关键参数
接入遗传算子后,主循环并不长。selTournament的tournsize=3表示每次从 3 个个体里选 1 个,既保证选择压力,又不会让最优个体垄断下一代。交叉用cxOnePoint,它交换两棵树的随机子树,比cxBlend更容易保持树结构;变异用mutUniform,将选中的节点替换成一个随机的子树。
toolbox = base.Toolbox() toolbox.register("expr_init", gp.genHalfAndHalf, pset=pset, min_=1, max_=3) toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr_init) toolbox.register("population", tools.initRepeat, list, toolbox.individual) toolbox.register("evaluate", eval_sr, X=X, y=y) toolbox.register("select", tools.selTournament, tournsize=3) toolbox.register("mate", gp.cxOnePoint) toolbox.register("expr_mut", gp.genFull, min_=0, max_=2) toolbox.register("mutate", gp.mutUniform, expr=toolbox.expr_mut, pset=pset) POP_SIZE = 300 N_GEN = 50 CX_PB = 0.7 MUT_PB = 0.2 pop = toolbox.population(n=POP_SIZE) hall_of_fame = tools.HallOfFame(1) stats = tools.Statistics(lambda ind: ind.fitness.values[0]) stats.register("min", np.min) stats.register("avg", np.mean) for gen in range(N_GEN): # 下一代先从选择压力选出候选 offspring = toolbox.select(pop, POP_SIZE) offspring = list(map(toolbox.clone, offspring)) # 两两交叉 for child1, child2 in zip(offspring[::2], offspring[1::2]): if random.random() < CX_PB: toolbox.mate(child1, child2) del child1.fitness.values, child2.fitness.values # 个体变异 for mutant in offspring: if random.random() < MUT_PB: toolbox.mutate(mutant) del mutant.fitness.values # 重新评估变化过的个体 invalid_ind = [ind for ind in offspring if not ind.fitness.valid] fitnesses = map(toolbox.evaluate, invalid_ind) for ind, fit in zip(invalid_ind, fitnesses): ind.fitness.values = fit # 精英保留:把上一代最优个体原样送回下一代 best_prev = tools.selBest(pop, 1)[0] offspring[-1] = tools.clone(best_prev) pop[:] = offspring hall_of_fame.update(pop) record = stats.compile(pop) print(f"Gen {gen}: min={record['min']:.6f}, avg={record['avg']:.6f}")逻辑说明:第一,cxpb和mutpb的概率是互相独立的,也就是说同一个个体可能先交叉再变异;DEAP 文档里也允许这样搭配,只要不把两个概率加总为 1。第二,交叉/变异后必须删除fitness.values,否则 DEAP 会认为该个体已经评估过,跳过重新评估,这是最常见的一个坑。第三,我在这里用了一种最简单的精英主义:把上一代最优个体放在后代的最后一个位置。它不参与交叉变异,但会覆盖掉后代中最后一个个体;缺点是减少了种群的多样性,后面第4章会讨论更稳定的前 k 个精英保留。
注意:交叉或变异后必须删除
fitness.values,否则新个体继承旧适应度,进化曲线会变成一条假的水平线。
3.3 输出并检查最优个体
跑完 50 代后,从hall_of_fame里取最优个体,直接打印树对象可读性很差。建议用sympy做符号化简,并至少在验证集上重新算一次误差:
from sympy import simplify from deap import gp as deap_gp best = hall_of_fame[0] expr_str = deap_gp.stringify(best) print(expr_str) # 原始树字符串 simplified = simplify(expr_str) # 有时能化简,但不要抱太大期望 print(simplified)stringify会把树还原成中缀表达式字符串,simplify只能处理代数恒等式,对含sin/cos的表达式作用有限。更多时候,你需要直接从数据里挑几个样本,把预测和真实值画成散点对比图;如果拟合优度 R² 高于 0.99,再考虑把它作为最终答案。这里需要强调的是,GP 找到的“公式”和真实生成函数很可能在形式上完全不一样,但只要在新样本上预测精度符合要求,就说明它已经学到了数据的生成规律;不要追求形式上完全相同。
4. 增强 GP 时的三个有效方向:精英主义、多样性保护与早停防暴长
基础实现拿 4 分后,剩下的分数大致靠增强、可视化和报告。增强部分最容易出效果的是对搜索过程的干预,而不是换一套新算法。这里给出三个投入产出比最高的方向,以及实际调参时遇到的边界。
4.1 精英主义:让每一代的最优个体不再丢失
定义:每一代选择最优秀的 k 个个体(通常 2~10),不经过交叉和变异直接复制到下一代。它保证了最优解是单调不增的,同时让你在记录适应度曲线时能画出一条持续下降的线。在 DEAP 中,完整做法是:
def elitism_replace(pop, offspring, k=5): # 从当前种群选出最好的 k 个,替换 offspring 中适应度最差的 k 个 best = tools.selBest(pop, k) worst_index = sorted( range(len(offspring)), key=lambda i: offspring[i].fitness.values[0], reverse=True )[:k] for idx, ind in zip(worst_index, best): offspring[idx] = tools.clone(ind)k太小不起作用,太大又会让种群过快收敛。我的经验是k = max(2, int(0.02 * POP_SIZE)),300 的种群取 5~6。要注意替换的对象是 offspring 中适应度最差的个体,而不是随机个体;否则会破坏已经形成的优良结构。精英主义虽然简单,但增强报告里记得要放同一参数下加与不加精英主义的两条适应度曲线,评分者一眼就能看出它带来的收敛差异。
4.2 多样性保护与抗膨胀:深度惩罚之外的软硬约束
单独的精英主义会让种群多样性快速流失,典型表现是第 20 代以后所有个体都共享同一个子树骨架,只是在外围加一点无关节点。常见做法有三种:用selTournament的基础上增加selDoubleTournament,它同时把树长也作为锦标赛比较指标;给每代的重复个体做去重,只保留一个副本;或者引入岛屿模型,把种群分成多个子岛屿,每隔几代交换几个个体。第一种在 DEAP 内置,最容易实现;第三种适合你有并行计算资源时用,代码量不大但对参数敏感。
膨胀(bloat)是符号回归最顽固的问题。即使适应度树长惩罚已经加上,树还是会在成千上万的代中缓慢变长。我的处理方式是双保险:初始化深度max_depth=4,变异时限制生成子树深度max_=2,并且在每一代结束后,对最优个体做一次“剪枝”:如果一个节点的某个子树不影响该个体在训练集上的输出,就把它替换成常数或变量。这一条虽然不是通用算法,但结合可视化可以给报告提供大量素材。参数边界见下表:
| 参数 | 常见范围 | 代表性风险 |
|---|---|---|
| 种群大小 POP_SIZE | 100~1000 | 过小早熟,过大每代评估成本线性上升 |
| 交叉概率 CX_PB | 0.6~0.9 | 过高破坏好结构,过低搜索停滞 |
| 变异概率 MUT_PB | 0.1~0.4 | 过高退化为随机搜索,过低容易陷在局部最优 |
| 锦标赛规模 | 2~5 | 越大选择压力越大,越小多样性越好 |
| parsimony_coef | 0.001~0.01 | 过大偏向短公式,过小膨胀失控 |
| 树深度限制 | 初始 3~5,变异 1~3 | 太浅无法表达复合函数,太深导致不可解释和过拟合 |
这张表里最容易被忽略的是“锦标赛规模”。不少初学者调到 7 或 8,期望更快收敛,结果最优解在第 10 代就锁死在一个局部结构里。另一个容易踩的坑是交叉概率和变异概率一起加起来超过 1。它们确实是独立事件,但如果你对同一个后代既交叉又变异,概率上会叠加很多变化,最终种群的破坏性操作多于建设性操作。更稳妥的配置是CX_PB=0.7, MUT_PB=0.2,剩下 10% 的个体由精英保留和直接复制进入下一代。
4.3 早停策略与统计检验
评估完每代后,连续 N 代最优适应度没有改善,就可以提前终止。DEAP 的HallOfFame只记录历史最优,你需要自己维护一个best_history列表,判断best_history[-5] == best_history[-1]。早停阈值过大会白跑算力,过小会把一个刚进入新搜索区域的种群误判为收敛。我建议在报告里把“完整跑 50 代”和“早停 5 代不动就停”放在一起,展示两者的训练/测试误差对比,这就是一种可接受的统计分析。更高阶的做法是用独立 t 检验比较两个随机种子的结果,但注意 GP 本身随机性很强,通常至少跑 10 个种子再算均值和标准差,不要只汇报一次运气。
5. 把 GP 结果变成能汇报的产出:符号表达式、可视化与类型化扩展
最后的产出部分,重点不是“再多跑几代”,而是把实验结果呈现成课程评分者不需要在代码里找的东西。按作业评分框架,LaTeX 报告、图表、统计分析和类型化问题可以叠加不少分数,但这些分数都建立在一个可工作的 GP 实现上。我从技巧层面讲两个我常用的生成方式。
5.1 用 sympy 生成 LaTeX 公式
报告中贴一棵树不是好主意,贴一张sympy渲染的公式才是。DEAP 树节点转sympy表达式可以用stringify后交给sympy.sympify,但含protected_div时符号简写会失效。更可控的办法是自己遍历树:把节点类型映射到sympy函数,叶子变量映射为Symbol("x0"),然后调用print_latex(expr)。例如:
from sympy import Symbol, latex, sin, cos, exp, log node_map = { "add": lambda x, y: x + y, "sub": lambda x, y: x - y, "mul": lambda x, y: x * y, "protected_div": lambda x, y: x / y, "sin": sin, "cos": cos, "exp": exp, "log": log, } def tree_to_sympy(expr): # expr 是 deap.gp.PrimitiveTree,这里做递归转换 stack = [] for node in expr: if node.arity == 0: if isinstance(node.value, str) and node.value.startswith("x"): stack.append(Symbol(node.value)) else: stack.append(node.value) else: args = [stack.pop() for _ in range(node.arity)][::-1] if node.name in node_map: stack.append(node_map[node.name](*args)) else: raise ValueError(f"unsupported node: {node.name}") return stack[-1]得到一个 sympy 表达式后,latex(expr)可以直接生成公式排版代码,放进 LaTeX 的 equation 环境。报告中再配合一幅“真实值 vs 预测值”散点图,就能让评分者快速相信你的 GP 确实工作正常。如果还想更直观,可用matplotlib画出树结构,但树容易过大;通常只画最优个体和几棵有代表性个体就够。
5.2 类型化 GP 的实战切入
类型化 GP 就是把原始集里的节点加上返回类型约束。比如树既要生成float特征,又需要布尔判断if x0 > 0.3时切换到另一条表达式,就需要让if_then_else接受一个布尔类型分支和两个 float 类型分支。这类问题在 UCI 数据集中很多,比如根据传感器数值判断设备状态,输出是分类标签。DEAP 的pset.addPrimitive本身不限制类型,但你可以通过定义不同类型的 PrimitiveSet 或在树的节点上做类型检查来模拟。更省力的方式是直接选一个回归类数据集,把目标值做离散化转换成布尔输出,然后给 GP 添加AND/OR/NOT算子,这样也能展示类型化能力。关键在于不要选用教程里的垃圾邮件检测,因为评分会排除它。
5.3 可视化进化过程的三条线
最后一种高性价比产出是画“训练误差、验证误差、树平均深度”三条曲线。横轴是代数,左纵轴是 RMSE,右纵轴是平均节点数;你会看到误差下降在某个代数后变缓,平均深度仍在上升,这就是膨胀的证据。配合四维数据可视化时,可以用scikit-learn的 PCA 把高维输入降到二维,再用等高线图画出 GP 模型的预测面,数据点和预测面放在同一张图上,比单纯画三维散点更清晰。画这些曲线时,训练误差和验证误差用左轴、平均树深用右轴,一眼就能看出膨胀是从哪一代开始的。
本文还有配套的精品资源,点击获取