遗传算法与深度学习实战（10）——遗传编程详解与实现

0. 前言

我们已经使用 DEAP 开发了各种问题的遗传算法 (Genetic Algorithms, GA) 解决方案。接下来，我们将使用 DEAP 来探索 GA 的子集，遗传编程 (Genetic Programming, GP)。GP 遵循与 GA 的相同原则，并采用许多相同的遗传算子。GA 和 GP 之间的关键区别在于基因或染色体的结构以及如何评估适应度，遗传编程和基因表达式编程 (Gene Expression Programming, GEP) 可以用于解决各种自动化和控制问题。在本节中，我们将介绍如何使用遗传编程来解决回归问题。GEP 也可以应用于从优化到搜索的其他问题。但回归和深度学习 ( Deep learning, DL) 更相关，我们可以使用类似 DL 的方式解决相同的问题。

1. 遗传编程原理

遗传编程 (Genetic Programming, GP) 是遗传算法的一种特殊形式。在这种特殊形式下，候选解(或个体)旨在找到最适合我们目的的解决方案，即实际的计算机程序。换句话说，当我们应用 GP 时，开发计算机程序的目的是找到一种在执行特定任务方面表现出色的程序。
GP 使我们能够更专注于如何解决问题，通过 GP，不仅能够搜索新颖的解决方案，并且开发可以用来推导这些解决方案的数学函数或程序代码，这些函数可以被重复使用或用于更好地理解特定问题。

1.1 基因表达式编程

在基因表达式编程 (Gene Expression Programming, GEP) 中，每个基因代表一个运算符、常量或输入，整个染色体或基因序列表示一个表达式树，其中运算符可以代表简单的加法、减法，也可以代表复杂的程序函数，常量和输入/参数代表单个标量值或更复杂的数组和张量。
下图展示了一个 GP 可能使用的基因序列。在图中，可以看到运算符和输入/常量的顺序形成了一个表达式树，可以将其作为方程式，然后使用数学规则进行计算。方程的输出可以与给定目标值进行比较，以确定个体的误差或适应度。

图中的基因序列显示了运算符、输入/参数和常量如何形成一个不均匀的叶节点表达式树。通过查看基因序列，可以看到第一个运算符是 $\ast$，它映射到树的根节点。接下来的两个基因扩展到第一级节点，表示为 $+$ 和 \sqrt {\ \ \ \ }      ​ 运算符。序列中接下来的两个基因映射到 $+$ 节点。随后，$-$ 运算符附加到 \sqrt {\ \ \ \ }      ​ 节点，最后几个基因映射到底层节点。
每个节点形成的子节点数取决于运算符的顺序。运算符的顺序可以是一元、二元、三元或 n 元的。为了简单起见，本节仅使用了二元运算符 $\ast$、$-$ 和 $+$ 以及一元运算符 \sqrt {\ \ \ \ }      ​。输入和常量没有顺序，始终表示表达式树中的叶节点。
通过表达式树，可以计算出方程结果，其结果代表输出，将输出与目标值进行比较，它们之间的差值就表示误差或适应度，此问题的目标是将适应度或误差减少到最小值。

1.2 遗传算子

通过分割每个父级的一部分并在父级之间交换分割的部分，从两个父级中创建了两个后代。类似地，用于树表示的交叉运算符可以从每个父级分离一个子树(一个分支或一组分支)，并在父级之间交换分离的分支以创建后代树，如下图所示：

同样，可以通过在候选解决方案中选择一个子树并将其替换为随机生成的子树来实现旨在将随机变化引入单个个体的突变算子。

2. 使用遗传编程解决回归问题

接下来，我们使用基因表达式编程 (Gene Expression Programming, GEP) 来推导方程，解决多变量回归问题，目标是使该方程在给定多个输入值的情况下成功地回归或预测输出值。

(1) 导入所需的库与数据集后，定义表示个体的一组基因。在代码中，可以看到三种不同类型的基因定义：运算符、常量和输入/参数：

import operator
import math
import randomimport numpy as npfrom deap import algorithms
from deap import base
from deap import creator
from deap import tools
from deap import gpfrom sklearn.datasets import load_boston
x, y = load_boston(return_X_y=True)
X = np.swapaxes(x,0,1)
inputs = X.shape[0]import pandas as pdboston_dataset = load_boston()
boston = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)
boston['MEDV'] = boston_dataset.target
boston.head()pset = gp.PrimitiveSet("MAIN", inputs)
pset.addPrimitive(np.add, 2, name="vadd")
pset.addPrimitive(np.subtract, 2, name="vsub")
pset.addPrimitive(np.multiply, 2, name="vmul")
pset.addPrimitive(protectedDiv, 2)
pset.addPrimitive(np.negative, 1, name="vneg")
pset.addPrimitive(np.cos, 1, name="vcos")
pset.addPrimitive(np.sin, 1, name="vsin")
pset.addEphemeralConstant("rand101", lambda: random.randint(-1,1))

(2) 创建 toolbox 之后，设置表达式树评估器。在本节中，我们使用一个开箱即用的表达式树生成器 genHalfandHalf，它有 50% 的概率使用两种不同形式的树之一：

creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)toolbox = base.Toolbox()
toolbox.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=2)
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("compile", gp.compile, pset=pset)toolbox.register("evaluate", evalSymbReg)
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)

在表达式树生成中，我们可以假设树是使用两个基本规则生成的。一个规则假设树中的所有叶节点都位于同一层次上，即均匀的；另一条规则假设叶节点可以是不均匀的。上一小节示例是一个不均匀叶节点表达式树的示例，使用 genHalfAndHalf 函数能够生成这两种形式的树。

(3) 编写 plot_expression() 函数，可视化方程表达式。使用 plot_expression() 函数绘制表达式树，输出最佳个体。在本节的问题中，我们希望适应度最小化或接近零。使用 plot_expression() 函数绘制的表达式树使用 Fruchterman-Reingold 算法定位节点，Networkx 提供了多种布局算法：

import matplotlib.pyplot as plt
import networkx as nxdef plot_expression(individual):options = {"node_size": 500, "alpha": 0.8}nodes, edges, labels = gp.graph(individual)g = nx.Graph()g.add_nodes_from(nodes)g.add_edges_from(edges)pos = nx.spring_layout(g)  nx.draw_networkx_nodes(g, pos, **options)nx.draw_networkx_edges(g, pos, width=1.0, alpha=0.5)nx.draw_networkx_labels(g, pos, labels, font_size=9, font_color='k')  plt.show()

GP 演化的目标是重建方程，计算目标值。在本节中，我们使用的是波士顿房价数据集，也可以将相同的原理应用于由输入特征 x 和目标输出 y 定义的结构化数据来解决其他形式的回归问题。

(4) 定义函数 protected_div()，替代常用的除法运算符，这样做是为了避免在使用 NumPy 时遇到除零错误，该函数来定义表达式基元集的除法运算符：

def protectedDiv(left, right):with np.errstate(divide='ignore',invalid='ignore'):x = np.divide(left, right)if isinstance(x, np.ndarray):x[np.isinf(x)] = 1x[np.isnan(x)] = 1elif np.isinf(x) or np.isnan(x):x = 1return x

(5) 定义适应度函数 evalSymbReg()，通过传入每个输入来评估编译表达式与值之间的误差，NumPy 能够一次处理所有样本数据行( 10000 行)，以输出总误差：

def evalSymbReg(individual):  # Transform the tree expression in a callable functionfunc = toolbox.compile(expr=individual)    # Evaluate the sum of squared difference between the expression      diff = math.sqrt(np.sum((func(*X.tolist()) - y)**2)) # unfold list of inputs x using *return diff,

(6) 运行演化过程：

def eaSimple(population, toolbox, cxpb, mutpb, ngen, stats=None, halloffame=None):  logbook = tools.Logbook()logbook.header = ['gen', 'nevals'] + (stats.fields if stats else [])# Evaluate the individuals with an invalid fitnessinvalid_ind = [ind for ind in population if not ind.fitness.valid]fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)for ind, fit in zip(invalid_ind, fitnesses):ind.fitness.values = fitif halloffame is not None:halloffame.update(population)record = stats.compile(population) if stats else {}logbook.record(gen=0, nevals=len(invalid_ind), **record)    print(logbook.stream)done = False# Begin the generational processfor gen in range(1, ngen + 1):if done: return# Select the next generation individualsoffspring = toolbox.select(population, len(population))offspring = [toolbox.clone(ind) for ind in offspring]# Apply crossover and mutation on the offspringfor i in range(1, len(offspring), 2):if random.random() < cxpb:offspring[i - 1], offspring[i] = toolbox.mate(offspring[i - 1],offspring[i])del offspring[i - 1].fitness.values, offspring[i].fitness.valuesfor i in range(len(offspring)):if random.random() < mutpb:offspring[i], = toolbox.mutate(offspring[i])del offspring[i].fitness.values# Evaluate the individuals with an invalid fitnessinvalid_ind = [ind for ind in offspring if not ind.fitness.valid]fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)for ind, fit in zip(invalid_ind, fitnesses):ind.fitness.values = fit  if fit[0] <= 135:print("Solved")done = True# Update the hall of fame with the generated individualsif halloffame is not None:halloffame.update(offspring)            plot_expression(halloffame[0])  print(halloffame[0]) func = toolbox.compile(expr=halloffame[0]) rows = x.shape[0]random_indices = np.random.choice(rows, size=1, replace=False)rd = x[random_indices, :] test = func(*rd[0].tolist()) print(test)   # Replace the current population by the offspringpopulation[:] = offspring# Append the current generation statistics to the logbookrecord = stats.compile(population) if stats else {}logbook.record(gen=gen, nevals=len(invalid_ind), **record)print(logbook.stream)     pop = toolbox.population(n=10000)
hof = tools.HallOfFame(1)
stats = tools.Statistics(lambda ind: ind.fitness.values)
stats.register("avg", np.mean)
stats.register("std", np.std)
stats.register("min", np.min)
stats.register("max", np.max)eaSimple(pop, toolbox, 0.5, 0.1, 400, stats, halloffame=hof)

下图显示了一个演化的输出表达式树，其误差(适应度)低于 1。这个表达式树图是使用 plot_expression() 函数中定义的 network 创建的。当评估这个表达式树时，可以看到产生的方程和结果代码与用于计算 y 的原始函数相匹配。由于表达式树引入了另一个受保护的除法运算符，这导致最后一项与实际方程式有所不同。从数学上讲，求解后的输出表达式树与用于生成 y 的原始方程相匹配。

GP 提供了使用表达式树生成方程或实际编程代码的能力。这是因为，从根本上讲，所有的编程代码都可以表示为一个表达式树，其中像 if 语句这样的语句是一个布尔运算符，它接受二进制输入或复杂函数，表示为一个具有单个返回值的 n 元运算符。
GEP 的优势在于可以生成实际的数学函数或编程代码，以便进行优化和重用。然而，GEP 也可能会生成过于复杂的函数或代码，这可能导致解决方案无法使用。如果输入大于四个，可以看到在演化过程中生成的表达式树的复杂性极大的增加了。
可以通过完成以下问题进一步理解如何使用遗传编程解决回归问题：

• 更改目标函数，然后重新运行代码
• 删除一些运算符，然后重新运行
• 修改遗传算法的运算符、交叉、选择或突变参数，然后重新运行

小结

在本节中，我们使用基因表达式编程 (Gene Expression Programming, GEP) 来推导方程，解决多变量回归问题，目标是使该方程在给定多个输入值的情况下成功地回归或预测输出值。本节仅使用预先输入到目标方程式的随机输入来验证结果。然而，该方法同样可以用来执行回归，类似于在深度学习 (Deep learning, DL) 中使用的方式。

系列链接

遗传算法与深度学习实战（1）——进化深度学习
遗传算法与深度学习实战（2）——生命模拟及其应用
遗传算法与深度学习实战（3）——生命模拟与进化论
遗传算法与深度学习实战（4）——遗传算法详解与实现
遗传算法与深度学习实战（5）——遗传算法框架DEAP
遗传算法与深度学习实战（6）——DEAP框架初体验
遗传算法与深度学习实战（7）——使用遗传算法解决N皇后问题
遗传算法与深度学习实战（8）——使用遗传算法解决旅行商问题
遗传算法与深度学习实战（9）——使用遗传算法重建图像