符号回归问题:GP简介

符号回归是GP中最为人熟知的问题之一(参见 参考)。它通常被用作新算法的调优问题,但也广泛应用于其他回归方法可能不适用的现实生活分布中。从概念上讲,这是一个简单的问题,因此它为DEAP中的GP框架提供了一个很好的入门示例。

所有符号回归问题都使用任意数据分布,并尝试用最准确的符号公式来拟合数据。通常,像 RMSE(均方根误差)或 MSE(均方误差)这样的度量被用来衡量个体的适应度。

在这个例子中,我们使用了一个经典的分布,即四次多项式 \((x^4 + x^3 + x^2 + x)\),这是一个一维分布。在范围 [-1, 1] 内生成了 20 个等距点,并用于评估适应度。

创建基本元素集

GP 程序中最关键的方面之一是原语集的选择。它们应该是良好的构建块,以便进化能够成功。在这个问题中,我们使用了一组经典的原语,即基本的算术函数:


# Define new functions
def protectedDiv(left, right):
    try:
        return left / right
    except ZeroDivisionError:
        return 1

pset = gp.PrimitiveSet("MAIN", 1)
pset.addPrimitive(operator.add, 2)
pset.addPrimitive(operator.sub, 2)
pset.addPrimitive(operator.mul, 2)
pset.addPrimitive(protectedDiv, 2)
pset.addPrimitive(operator.neg, 1)
pset.addPrimitive(math.cos, 1)

对除法的重新定义是为了防止零除错误(这会导致程序崩溃)。其他函数只是从 Python operator 模块的映射。函数后面的数字是原语的 arity,即它所接受的条目数。

在最后一行,我们声明了一个 Ephemeral 常量。这是一种特殊终端类型,它没有固定值。当程序将一个临时常量终端附加到树时,它包含的函数被执行,其结果作为常量终端插入。在这种情况下,这些常量终端可以取值 -1、0 或 1。

第二个参数是 PrimitiveSet 的输入数量。这里,由于我们只有一个一维回归问题,所以只有一个输入,但它可以有你想要的任意多个。默认情况下,这些输入被命名为“ARGx”,其中“x”是一个数字,但你可以轻松地重命名它们:

pset.addPrimitive(math.sin, 1)

创建者

作为一个进化程序,符号回归需要(至少)两种对象类型:一个包含基因型和适应度的个体。我们可以很容易地用创建者来创建它们:

pset.renameArguments(ARG0='x')

第一行创建了适应度对象(这是一个最小化问题,所以权重是负的)。weights 参数必须是一个权重的可迭代对象,即使只有一个适应度度量。第二行创建了个体对象本身。非常直接,我们可以看到它将基于一棵树,并向其中添加一个适应度。如果出于任何原因,用户想要添加任何其他属性(例如,一个将保存个体的文件),只需将该属性的任何类型添加到这一行即可。在此声明之后,任何生成的个体都将包含这些想要的属性。

工具箱

现在,我们想要注册一些特定于进化过程的参数。在DEAP中,这是通过工具箱完成的:

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)

def evalSymbReg(individual, points):
    # Transform the tree expression in a callable function
    func = toolbox.compile(expr=individual)
    # Evaluate the mean squared error between the expression
    # and the real function : x**4 + x**3 + x**2 + x
    sqerrors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points)
    return math.fsum(sqerrors) / len(points),

toolbox.register("evaluate", evalSymbReg, points=[x/10. for x in range(-10,10)])
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)

首先,创建一个工具箱实例(在某些问题类型如协同进化中,您可能需要创建多个工具箱)。然后,我们可以注册任何参数。前几行注册了如何创建个体(通过调用之前定义的原始集的 gp.genHalfAndHalf),以及如何创建种群(通过重复个体初始化)。

我们现在可以介绍评估函数,该函数将接收一个个体作为输入,并返回相应的适应度。此函数使用先前定义的 compile 函数将个体转换为其可执行形式——即一个程序。之后,评估只是简单的数学运算,其中评估个体产生的值与实际值之间的差异被平方并求和以计算均方误差(MSE),该误差作为个体的适应度返回。

警告

即使适应度只包含一个度量,请记住 DEAP 将其存储为可迭代对象。了解这一点后,您可以理解为什么评估函数必须返回一个元组值(即使它是一个 1-元组):

def evalSymbReg(individual, points):
    # Transform the tree expression in a callable function
    func = toolbox.compile(expr=individual)
    # Evaluate the mean squared error between the expression
    # and the real function : x**4 + x**3 + x**2 + x
    sqerrors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points)
    return math.fsum(sqerrors) / len(points),

仅返回值会产生奇怪的行为和错误,因为选择和统计函数依赖于适应度始终是一个可迭代对象这一事实。

之后,我们注册评估函数。我们还选择了选择方法(大小为3的锦标赛)、交配方法(在所有节点上均匀概率的单点交叉)和变异方法(均匀概率变异,可能会将一个新的完整子树附加到一个节点上)。

然后,我们装饰 mate 和 mutate 方法以限制生成个体的高度。这是为了避免遗传编程的一个重要缺点:膨胀。Koza 在他的遗传编程书中建议使用最大深度为 17。

此时,任何能够访问工具箱实例的结构也将能够访问所有这些已注册的参数。当然,用户可以根据自己的需求注册其他参数。

统计

虽然可选,统计数据在进化编程中通常很有用。DEAP 提供了一个简单的类,可以处理大部分“枯燥的工作”。在这种情况下,我们希望计算个体适应度和大小的均值、标准差、最小值和最大值。为此,我们将使用一个 MultiStatistics 对象。

    hof = tools.HallOfFame(1)

    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
    stats_size = tools.Statistics(len)
    mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size)
    mstats.register("avg", numpy.mean)
    mstats.register("std", numpy.std)

请注意,一个简单的 Statistics 对象可以被使用,就像在前面的例子中,当只需要对单一键进行统计时。

启动进化

此时,DEAP 已经拥有开始进化过程所需的所有信息,但还没有进行初始化。我们可以通过创建种群,然后调用一个完整的算法来开始进化。在这种情况下,我们将使用 eaSimple()

    random.seed(318)

    mstats.register("max", numpy.max)

名人堂是一个特定的结构,包含*n*个最佳个体(这里,仅包含最佳的一个)。

完整的 examples/gp/symbreg

参考

John R. Koza, “遗传编程:通过自然选择对计算机进行编程”, MIT Press, 1992, 第162-169页.