0. 前言
地球上的生命存在共生关系，数百万个物种相互依存以求生存，我们可以用协同进化描述这种关系。当我们尝试解决复杂问题时，可以利用进化方法来模拟协同进化。在本节中，我们将使用波士顿房地产市场的样本结构化数据集，利用协同进化解决回归问题。
1. 协同进化原理
1.1 协同进化理论基础
协同进化 (Coevolving) 是指不同生物种群之间相互影响和适应的过程，这种相互作用可以导致它们彼此之间的进化变化。在生态系统中，不同种群之间的相互作用可以促进它们的进化，这种相互作用包括共生、竞争、捕食等。例如，共生关系中的两个物种可以通过相互依赖来进化，竞争中的物种可以通过适应性变化来提高竞争力，捕食者和猎物之间的进化则可能导致双方的适应性变化。
此外，协同进化还可以涉及到物种与其环境之间的相互作用，例如某些植物对特定的授粉者进化以实现更有效的传粉。这种进化方式强调了生物之间复杂而密切的关系，这些关系不仅塑造了物种的结构和功能，也在演化的过程中起到了重要作用。
协同进化是生物进化过程中重要的一部分，强调了物种之间相互作用和依赖如何推动生物进化的变化和多样性。
1.2 问题描述
波士顿房地产数据集包含 13 个特征列，用于预测波士顿市场上房屋的价格。我们可以单独使用遗传编程来尝试推导出一个方程，但结果会过于复杂。在本节中，我们将结合遗传编程 (Genetic Programming, GP) 与遗传算法 (Genetic Algorithms, GA)，使用协同进化方法，以期得到一个更简单的输出。
2. 协同进化实现
接下来，我们将使用协同进化(同时使用 GP 和 GA )解决解决波士顿房价预测问题。
2.1 代码实现
(1) 导入的波士顿房地产数据集，波士顿房地产数据集可以从 sklearn.datasets 中加载，返回特征 x 以及目标 y。然后交换轴以适应行、特征格式和提取的输入数量，输入数量定义了推导方程中参数的数量：
import operator
import math
import random
import time
import numpy as np
from deap import algorithms
from deap import base
from deap import creator
from deap import tools
from deap import gp
from IPython.display import clear_output
#GA
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("IndGA", list, fitness=creator.FitnessMax)
#GP
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)
from 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 pd
boston_dataset = load_boston()
boston = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)
boston.head()
在使用遗传编程解决回归问题一节中，将方程演化到最小适应度分数低于 135 的结果如下图所示，得到的方程表明，并非所有特征都被用于推导方程，只有 ARG5、ARG7、ARG10、ARG11 和 ARG12 是相关的。这也意味着 GP 解决方案通过忽略不相关的特征自动执行特征选择。
(2) 接下来，同时使用遗传编程 (Genetic Programming, GP) 与遗传算法 (Genetic Algorithms, GA) 方法微调推导方程。首先，构建两个 toolbox——一个用于 GA 种群，另一个用于 GP：
toolbox_ga = base.Toolbox()
toolbox_ga.register("float", random.uniform, -1, 1)
toolbox_ga.register("individual", tools.initRepeat, creator.IndGA, toolbox_ga.float, inputs)
toolbox_ga.register("population", tools.initRepeat, list, toolbox_ga.individual)
toolbox_ga.register("select", tools.selTournament, tournsize=3)
toolbox_ga.register("mate", tools.cxTwoPoint)
toolbox_ga.register("mutate", tools.mutGaussian, mu=0, sigma=0.01, indpb=0.05)
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)] = 1
x[np.isnan(x)] = 1
elif np.isinf(x) or np.isnan(x):
x = 1
return x
#@title Create Set of Operators
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.randrange(-1000,1000))
toolbox_gp = base.Toolbox()
toolbox_gp.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=2)
toolbox_gp.register("individual", tools.initIterate, creator.Individual, toolbox_gp.expr)
toolbox_gp.register("population", tools.initRepeat, list, toolbox_gp.individual)
toolbox_gp.register("compile", gp.compile, pset=pset)
(3) 在协同进化中，GP 求解器用于构建一个推导方程。GA 求解器与其并行地协同进化一个大小等于输入/特征数的标量基因序列。对于波士顿房地产数据集，输入数量等于 13，每个 GA 标量基因值用于缩放输入到方程中的特征，应用于评估函数 evalSymbReg() 中。使用评估函数 evalSymbReg() 评估适应度时，传递两个个体(包括一个 GP 个体和一个 GA 个体)，每次使用这个函数评估时，都是使用 GA 和 GP 种群的两个个体进行评估的：
def evalSymbReg(individual, points):
# Transform the tree expression in a callable function
func = toolbox_gp.compile(expr=individual)
# Evaluate the mean squared error between the expression
p = np.expand_dims(points, axis=1)
x = X * np.asarray(p)
diff = math.sqrt(np.sum((func(*x.tolist()) - y)**2))
return diff,
toolbox_gp.register("evaluate", evalSymbReg, points=[x/10. for x in range(-10,10)])
toolbox_gp.register("select", tools.selTournament, tournsize=3)
toolbox_gp.register("mate", gp.cxOnePoint)
toolbox_gp.register("expr_mut", gp.genFull, min_=0, max_=2)
toolbox_gp.register("mutate", gp.mutUniform, expr=toolbox_gp.expr_mut, pset=pset)
toolbox_gp.decorate("mate", gp.staticLimit(key=operator.attrgetter("height"), max_value=5))
toolbox_gp.decorate("mutate", gp.staticLimit(key=operator.attrgetter("height"), max_value=5))
import matplotlib.pyplot as plt
import networkx as nx
def 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()
(4) 通常，在协同进化的情况下，我们不希望两种方法以相同的速率演化。例如，在本节示例中，我们设定 GA 种群(标量基因)的进化速度比 GP 种群更快，这使得 GA 方法能够更快的微调 GP 方程中使用的参数的比例或权重，本质上，这使得 GA 方法可以微调方程以获得更好的拟合。在本节中，我们设置在 GA 和 GP 之间的进化速率为 10:1：
GA_GEN, GP_GEN, BASE_POP = 1, 10, 10000
pop_ga = toolbox_ga.population(n=BASE_POP*GA_GEN)
pop_gp = toolbox_gp.population(n=BASE_POP*GP_GEN)
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)
logbook = tools.Logbook()
logbook.header = "gen", "type", "evals", "std", "min", "avg", "max"
best_ga = tools.selRandom(pop_ga, 1)[0]
best_gp = tools.selRandom(pop_gp, 1)[0]
for ind in pop_gp:
ind.fitness.values = toolbox_gp.evaluate(ind, points=best_ga)
for ind in pop_ga:
ind.fitness.values = toolbox_gp.evaluate(best_gp, points=ind)
record = stats.compile(pop_ga)
logbook.record(gen=0, type='ga', evals=len(pop_ga), **record)
record = stats.compile(pop_gp)
logbook.record(gen=0, type='gp', evals=len(pop_gp), **record)
(5) 在演化过程中，控制 GA 和 GP 种群的进化频率(使用 GA_GEN 控制)：
print(logbook.stream)
CXPB, MUTPB, NGEN = 0.5, 0.2, 100
done = False
# Begin the evolution
for g in range(1, NGEN):
if done: break
if (g+1) % GA_GEN == 0:
# Select and clone the offspring
off_ga = toolbox_ga.select(pop_ga, len(pop_ga))
off_ga = [toolbox_ga.clone(ind) for ind in off_ga]
# Apply crossover and mutation
for ind1, ind2 in zip(off_ga[::2], off_ga[1::2]):
if random.random() < CXPB:
toolbox_ga.mate(ind1, ind2)
del ind1.fitness.values
del ind2.fitness.values
for ind in off_ga:
if random.random() < MUTPB:
toolbox_ga.mutate(ind)
del ind.fitness.values
pop_ga = off_ga
if (g+1) % GP_GEN == 0:
off_gp = toolbox_gp.select(pop_gp, len(pop_gp))
off_gp = [toolbox_gp.clone(ind) for ind in off_gp]
for ind1, ind2 in zip(off_gp[::2], off_gp[1::2]):
if random.random() < CXPB:
toolbox_gp.mate(ind1, ind2)
del ind1.fitness.values
del ind2.fitness.values
for ind in off_gp:
if random.random() < MUTPB:
toolbox_gp.mutate(ind)
del ind.fitness.values
# Replace the old population by the offspring
pop_gp = off_gp
# Evaluate the individuals with an invalid fitness
for ind in pop_gp:
ind.fitness.values = toolbox_gp.evaluate(ind, points=best_ga)
for ind in pop_ga:
ind.fitness.values = toolbox_gp.evaluate(best_gp, points=ind)
record = stats.compile(pop_ga)
best_ga_min = record["min"]
logbook.record(gen=g, type='ga', evals=len(pop_ga), **record)
record = stats.compile(pop_gp)
best_gp_min = record["min"]
logbook.record(gen=g, type='gp', evals=len(pop_gp), **record)
#print(logbook.stream)
print(record)
best_ga = tools.selBest(pop_ga, 1)[0]
best_gp = tools.selBest(pop_gp, 1)[0]
if best_gp_min < 135:
done = True
print("Sovled")
if (g+1) % 1 == 0:
clear_output()
print(logbook.stream)
plot_expression(best_gp)
output = ["%.2f" % e for e in best_ga]
print(f"Best GA {output}")
time.sleep(1)
print("Best individual GA is %s, %s" % (best_ga, best_ga.fitness.values))
print("Best individual GP is %s, %s" % (best_gp, best_gp.fitness.values))
plot_expression(best_gp)
评估适应度时，函数需要来自两个种群( GA 和 GP )的个体。这意味着评估个体的代码与每个种群相关联。当我们评估 GA 或 GP 种群的适应度时，我们使用另一个种群的最佳表现者，这种简化近似了每个种群的最高适应度。另一种方法是循环遍历两个种群并测试每个组合，但这种方法时间复杂度很高，因此在本节中我们使用简化方法。
下图显示了最终得到的解决方案，假设适应度低于 135。如图所示，推导方程已经被大大简化。最佳的 GA 个体分配了用于修改方程输入的标量，并且可以看到最终方程如何应用输入权重。
最终的推导方程显示了预测波士顿房地产市场价格的可能解决方案。如果我们查看结果方程及波士顿房地产特征，可以看到 ARG5 (NOX, 氧化氮浓度)、ARG12 (LSTAT, 低收入人口比例)被识别为主要使用的特征。
使用遗传编程求解波士顿房价问题得到的解决方案中，可以看到 ARG5、ARG12 与 ARG7、ARG10 和 ARG11 被认为是重要特征。协同进化的解决方案能够通过传递到方程中的输入进一步减少特征的权重，这可能导致方程过于简化，但通过这种方法，我们可以确定波士顿房地产数据集中 NOX 和 LSTAT 特征之间的关键相关性。
2.2 结果分析
当 GP 演化时，经常创建过于复杂的方程。长时间运行的演化甚至可能因表达树的大小而中断。将 GA 与 GP 结合使用可以对方程进行微调。然而，这可能导致问题的过度简化。数据集的大小是另一个问题，因为波士顿房地产数据集只有 500 多个数据样本，在大多数情况下，使用更多数据可以获得更好的结果。
虽然，本节得到的协同进化解决方案并不完美。它确实能够解决 GP 中经常出现的复杂问题。最终答案可能过于简单，因此忽略了其他关键特征而未能推断出足够准确的结果。然而，它足够简单，可以使用简单计算快速得到结果。我们在后续学习中将使用其他协同进化解决方案来平衡应用于深度学习 ( Deep learning, DL) 的多种形式的进化计算 (Evolutionary Computation, EC)，协同进化可以将多种形式的 EC 结合在一起，以解决复杂问题。
小结
协同进化是指确定两个或多个个体种群来解决特定问题的独特任务的算法，协同进化可以通过最小化和缩放推导方程中的特征来找到复杂的解决方案。在本节中，我们使用了协同进化(结合遗传编程和遗传算法)在波士顿房地产数据集上解决回归问题。
系列链接
遗传算法与深度学习实战（1）——进化深度学习
遗传算法与深度学习实战（2）——生命模拟及其应用
遗传算法与深度学习实战（3）——生命模拟与进化论
遗传算法与深度学习实战（4）——遗传算法详解与实现
遗传算法与深度学习实战（5）——遗传算法中常用遗传算子
遗传算法与深度学习实战（6）——遗传算法框架DEAP
遗传算法与深度学习实战（7）——DEAP框架初体验
遗传算法与深度学习实战（8）——使用遗传算法解决N皇后问题
遗传算法与深度学习实战（9）——使用遗传算法解决旅行商问题
遗传算法与深度学习实战（10）——使用遗传算法重建图像
遗传算法与深度学习实战（11）——遗传编程详解与实现
遗传算法与深度学习实战（12）——粒子群优化详解与实现