尧图网页设计 
简介基于遗传编程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)、 (x0)、 (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, arity4) pset.addPrimitive(np.add, 2, nameadd) pset.addPrimitive(np.subtract, 2, namesub) pset.addPrimitive(np.multiply, 2, namemul) pset.addPrimitive(protected_div, 2) pset.addPrimitive(np.sin, 1) pset.addPrimitive(np.cos, 1) pset.addEphemeralConstant(rand_const, lambda: random.uniform(-1, 1))这里arity4表示有四个输入变量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_divprotected_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 很容易出现(x0x1)*((x0-x0)x0)这类冗余结构树长惩罚的强度由一个parsimony_coefficient控制通常设成 0.0010.01 之间。另一个思路是把“模型复杂度”作为第二目标用 NSGA-II 做多目标演化但实现和工作量会明显增加评分上的收益需要权衡。此外如果训练数据分成了 train/validate/testGP 内部只应使用训练集计算适应度验证集负责判断是否早停测试集留到最后报告泛化误差。因为符号回归很容易过拟合树一旦变深几乎可以记住所有训练样本包括噪声。我见过一个项目在训练集 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, arityn_features) pset.addPrimitive(np.add, 2, nameadd) pset.addPrimitive(np.subtract, 2, namesub) pset.addPrimitive(np.multiply, 2, namemul) 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(**{fARG{i}: fx{i} for i in range(n_features)}) creator.create(FitnessMin, base.Fitness, weights(-1.0,)) creator.create(Individual, gp.PrimitiveTree, fitnesscreator.FitnessMin)代码前几行很直白protected_div处理除零PSet的定义决定了可搜索的算子集合。weights(-1.0,)表示单目标最小化如果你的增强版本引入两个目标可以写成weights(-1.0, -0.1)第二维对应复杂度但那样选择算子也要切换成selNSGA2。评估函数里需要把个体编译成可调用函数。DEAP 的toolbox.compile(exprindividual)会把PrimitiveTree转成一个普通 Python 可调用对象输入是n_features个位置参数。然后对全部训练样本做列式计算这里不要用 Python 显式循环先转成 NumPy 数组再向量化否则种群规模一大每代评估会慢到让你怀疑人生。def eval_sr(individual, X, y, parsimony_coef0.005): func toolbox.compile(exprindividual) # 向量化预测每个个体一次算完所有样本 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的tournsize3表示每次从 3 个个体里选 1 个既保证选择压力又不会让最优个体垄断下一代。交叉用cxOnePoint它交换两棵树的随机子树比cxBlend更容易保持树结构变异用mutUniform将选中的节点替换成一个随机的子树。toolbox base.Toolbox() toolbox.register(expr_init, gp.genHalfAndHalf, psetpset, 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, XX, yy) toolbox.register(select, tools.selTournament, tournsize3) toolbox.register(mate, gp.cxOnePoint) toolbox.register(expr_mut, gp.genFull, min_0, max_2) toolbox.register(mutate, gp.mutUniform, exprtoolbox.expr_mut, psetpset) POP_SIZE 300 N_GEN 50 CX_PB 0.7 MUT_PB 0.2 pop toolbox.population(nPOP_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(fGen {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 个个体通常 210不经过交叉和变异直接复制到下一代。它保证了最优解是单调不增的同时让你在记录适应度曲线时能画出一条持续下降的线。在 DEAP 中完整做法是def elitism_replace(pop, offspring, k5): # 从当前种群选出最好的 k 个替换 offspring 中适应度最差的 k 个 best tools.selBest(pop, k) worst_index sorted( range(len(offspring)), keylambda i: offspring[i].fitness.values[0], reverseTrue )[:k] for idx, ind in zip(worst_index, best): offspring[idx] tools.clone(ind)k太小不起作用太大又会让种群过快收敛。我的经验是k max(2, int(0.02 * POP_SIZE))300 的种群取 56。要注意替换的对象是 offspring 中适应度最差的个体而不是随机个体否则会破坏已经形成的优良结构。精英主义虽然简单但增强报告里记得要放同一参数下加与不加精英主义的两条适应度曲线评分者一眼就能看出它带来的收敛差异。4.2 多样性保护与抗膨胀深度惩罚之外的软硬约束单独的精英主义会让种群多样性快速流失典型表现是第 20 代以后所有个体都共享同一个子树骨架只是在外围加一点无关节点。常见做法有三种用selTournament的基础上增加selDoubleTournament它同时把树长也作为锦标赛比较指标给每代的重复个体做去重只保留一个副本或者引入岛屿模型把种群分成多个子岛屿每隔几代交换几个个体。第一种在 DEAP 内置最容易实现第三种适合你有并行计算资源时用代码量不大但对参数敏感。膨胀bloat是符号回归最顽固的问题。即使适应度树长惩罚已经加上树还是会在成千上万的代中缓慢变长。我的处理方式是双保险初始化深度max_depth4变异时限制生成子树深度max_2并且在每一代结束后对最优个体做一次“剪枝”如果一个节点的某个子树不影响该个体在训练集上的输出就把它替换成常数或变量。这一条虽然不是通用算法但结合可视化可以给报告提供大量素材。参数边界见下表参数常见范围代表性风险种群大小 POP_SIZE1001000过小早熟过大每代评估成本线性上升交叉概率 CX_PB0.60.9过高破坏好结构过低搜索停滞变异概率 MUT_PB0.10.4过高退化为随机搜索过低容易陷在局部最优锦标赛规模25越大选择压力越大越小多样性越好parsimony_coef0.0010.01过大偏向短公式过小膨胀失控树深度限制初始 35变异 13太浅无法表达复合函数太深导致不可解释和过拟合这张表里最容易被忽略的是“锦标赛规模”。不少初学者调到 7 或 8期望更快收敛结果最优解在第 10 代就锁死在一个局部结构里。另一个容易踩的坑是交叉概率和变异概率一起加起来超过 1。它们确实是独立事件但如果你对同一个后代既交叉又变异概率上会叠加很多变化最终种群的破坏性操作多于建设性操作。更稳妥的配置是CX_PB0.7, MUT_PB0.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(funsupported 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 模型的预测面数据点和预测面放在同一张图上比单纯画三维散点更清晰。画这些曲线时训练误差和验证误差用左轴、平均树深用右轴一眼就能看出膨胀是从哪一代开始的。本文还有配套的精品资源点击获取 符号回归实践：遗传编程自动发现公式（Python实现）
简介基于遗传编程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)、 (x0)、 (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, arity4) pset.addPrimitive(np.add, 2, nameadd) pset.addPrimitive(np.subtract, 2, namesub) pset.addPrimitive(np.multiply, 2, namemul) pset.addPrimitive(protected_div, 2) pset.addPrimitive(np.sin, 1) pset.addPrimitive(np.cos, 1) pset.addEphemeralConstant(rand_const, lambda: random.uniform(-1, 1))这里arity4表示有四个输入变量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_divprotected_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 很容易出现(x0x1)*((x0-x0)x0)这类冗余结构树长惩罚的强度由一个parsimony_coefficient控制通常设成 0.0010.01 之间。另一个思路是把“模型复杂度”作为第二目标用 NSGA-II 做多目标演化但实现和工作量会明显增加评分上的收益需要权衡。此外如果训练数据分成了 train/validate/testGP 内部只应使用训练集计算适应度验证集负责判断是否早停测试集留到最后报告泛化误差。因为符号回归很容易过拟合树一旦变深几乎可以记住所有训练样本包括噪声。我见过一个项目在训练集 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, arityn_features) pset.addPrimitive(np.add, 2, nameadd) pset.addPrimitive(np.subtract, 2, namesub) pset.addPrimitive(np.multiply, 2, namemul) 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(**{fARG{i}: fx{i} for i in range(n_features)}) creator.create(FitnessMin, base.Fitness, weights(-1.0,)) creator.create(Individual, gp.PrimitiveTree, fitnesscreator.FitnessMin)代码前几行很直白protected_div处理除零PSet的定义决定了可搜索的算子集合。weights(-1.0,)表示单目标最小化如果你的增强版本引入两个目标可以写成weights(-1.0, -0.1)第二维对应复杂度但那样选择算子也要切换成selNSGA2。评估函数里需要把个体编译成可调用函数。DEAP 的toolbox.compile(exprindividual)会把PrimitiveTree转成一个普通 Python 可调用对象输入是n_features个位置参数。然后对全部训练样本做列式计算这里不要用 Python 显式循环先转成 NumPy 数组再向量化否则种群规模一大每代评估会慢到让你怀疑人生。def eval_sr(individual, X, y, parsimony_coef0.005): func toolbox.compile(exprindividual) # 向量化预测每个个体一次算完所有样本 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的tournsize3表示每次从 3 个个体里选 1 个既保证选择压力又不会让最优个体垄断下一代。交叉用cxOnePoint它交换两棵树的随机子树比cxBlend更容易保持树结构变异用mutUniform将选中的节点替换成一个随机的子树。toolbox base.Toolbox() toolbox.register(expr_init, gp.genHalfAndHalf, psetpset, 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, XX, yy) toolbox.register(select, tools.selTournament, tournsize3) toolbox.register(mate, gp.cxOnePoint) toolbox.register(expr_mut, gp.genFull, min_0, max_2) toolbox.register(mutate, gp.mutUniform, exprtoolbox.expr_mut, psetpset) POP_SIZE 300 N_GEN 50 CX_PB 0.7 MUT_PB 0.2 pop toolbox.population(nPOP_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(fGen {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 个个体通常 210不经过交叉和变异直接复制到下一代。它保证了最优解是单调不增的同时让你在记录适应度曲线时能画出一条持续下降的线。在 DEAP 中完整做法是def elitism_replace(pop, offspring, k5): # 从当前种群选出最好的 k 个替换 offspring 中适应度最差的 k 个 best tools.selBest(pop, k) worst_index sorted( range(len(offspring)), keylambda i: offspring[i].fitness.values[0], reverseTrue )[:k] for idx, ind in zip(worst_index, best): offspring[idx] tools.clone(ind)k太小不起作用太大又会让种群过快收敛。我的经验是k max(2, int(0.02 * POP_SIZE))300 的种群取 56。要注意替换的对象是 offspring 中适应度最差的个体而不是随机个体否则会破坏已经形成的优良结构。精英主义虽然简单但增强报告里记得要放同一参数下加与不加精英主义的两条适应度曲线评分者一眼就能看出它带来的收敛差异。4.2 多样性保护与抗膨胀深度惩罚之外的软硬约束单独的精英主义会让种群多样性快速流失典型表现是第 20 代以后所有个体都共享同一个子树骨架只是在外围加一点无关节点。常见做法有三种用selTournament的基础上增加selDoubleTournament它同时把树长也作为锦标赛比较指标给每代的重复个体做去重只保留一个副本或者引入岛屿模型把种群分成多个子岛屿每隔几代交换几个个体。第一种在 DEAP 内置最容易实现第三种适合你有并行计算资源时用代码量不大但对参数敏感。膨胀bloat是符号回归最顽固的问题。即使适应度树长惩罚已经加上树还是会在成千上万的代中缓慢变长。我的处理方式是双保险初始化深度max_depth4变异时限制生成子树深度max_2并且在每一代结束后对最优个体做一次“剪枝”如果一个节点的某个子树不影响该个体在训练集上的输出就把它替换成常数或变量。这一条虽然不是通用算法但结合可视化可以给报告提供大量素材。参数边界见下表参数常见范围代表性风险种群大小 POP_SIZE1001000过小早熟过大每代评估成本线性上升交叉概率 CX_PB0.60.9过高破坏好结构过低搜索停滞变异概率 MUT_PB0.10.4过高退化为随机搜索过低容易陷在局部最优锦标赛规模25越大选择压力越大越小多样性越好parsimony_coef0.0010.01过大偏向短公式过小膨胀失控树深度限制初始 35变异 13太浅无法表达复合函数太深导致不可解释和过拟合这张表里最容易被忽略的是“锦标赛规模”。不少初学者调到 7 或 8期望更快收敛结果最优解在第 10 代就锁死在一个局部结构里。另一个容易踩的坑是交叉概率和变异概率一起加起来超过 1。它们确实是独立事件但如果你对同一个后代既交叉又变异概率上会叠加很多变化最终种群的破坏性操作多于建设性操作。更稳妥的配置是CX_PB0.7, MUT_PB0.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(funsupported 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 模型的预测面数据点和预测面放在同一张图上比单纯画三维散点更清晰。画这些曲线时训练误差和验证误差用左轴、平均树深用右轴一眼就能看出膨胀是从哪一代开始的。本文还有配套的精品资源点击获取