关键词

遗传编程语法树Koza函数集终止符集
子树交叉子树变异符号回归自动函数发现程序演进
表达树树深度限制饱和度初始化方法适应度评估
闭合性充分性GP变体行为遗传多目标GP

摘要

遗传编程(Genetic Programming, GP)是由John Koza于1989年在其标志性著作《Genetic Programming: On the Programming of Computers by Means of Natural Selection》中系统提出的一种自动程序生成技术。与传统遗传算法固定长度编码不同,GP采用可变大小和形状的树形结构表示计算机程序,实现了程序空间的自动搜索与优化。这一突破性范式使计算机能够”自动编程”,已在符号回归、电路设计、游戏策略、自动化控制等领域展现出卓越的问题求解能力。


零、遗传编程入门:让程序自己进化

0.1 什么是遗传编程?

想象一下这个场景:你想让计算机自动写出一个能识别猫狗的程序,但你不知道该用什么特征、什么算法。传统做法是你得自己设计特征、搭神经网络、调参数,累死累活还不一定效果好。

遗传编程告诉你:别费劲了,让程序自己”进化”出来。

具体怎么玩?你先准备好一堆”零件”——比如加减乘除运算、变量x和y、一些随机常量。然后告诉GP:“这些零件随便你怎么组合,只要最后做出来的东西能解决问题就行。”

GP会:

  1. 随机用零件拼出一些”草稿程序”
  2. 测试每个草稿的效果(这就是”适应度”)
  3. 效果好的程序有更高几率”繁殖”,效果差的被淘汰
  4. “繁殖”过程中还有点”突变”,加点新花样
  5. 不断重复这几步,直到找到足够好的程序

整个过程就像养蛊——啊不,养一群程序,让它们互相竞争、进化,最后活下来的就是精英。

0.2 一个简单的例子

假设你想找到这样一个公式:输入x,输出接近sin(x) + x²。

用GP来做:

  1. 零件准备:函数集F = {+, -, *, /, sin, cos},终止符集T = {x, 随机常量}

  2. 生成草稿:随机组合出一些表达式树,比如:

    • 树A:sin(x) + x * x
    • 树B:cos(x) - x + x * x
    • 树C:sin(x) * x
  3. 评估适应度:用一堆x值测试,看每个公式算出来的结果和真实值差多少。误差越小,适应度越高。

  4. 进化:适应度高的留下繁殖下一代,过程中随机组合、变异。

  5. 迭代:重复评估和进化几十上百代。

几轮下来,你可能就得到了一个和sin(x) + x²非常接近的表达式。

这比你自己猜公式要靠谱多了——GP帮你穷举、筛选、进化,比人脑高效。

0.3 遗传编程 vs 遗传算法:为什么GP的”染色体”是树状的?

学过的同学可能会有疑问:遗传算法(GA)不是用二进制串或者实数向量当”染色体”吗?遗传编程怎么用树来表示?

好问题。

遗传算法处理的是参数优化问题——你知道问题是什么,也知道大概的解的结构,只是不知道具体的参数值。比如调度问题,你只是不知道每个任务的先后顺序怎么排最优。

遗传编程处理的是程序合成问题——你甚至不知道解应该是什么形式。可能是x+y,可能是sin(x)*cos(y),也可能是if x > 0 then x else -x

用固定长度的串来编码这种结构,编码长度本身就成了瓶颈。你得预先定好串有多长、每段代表什么意思,但程序的结构本身就是可变的。

所以GP用了树结构:

  • 树的形状可以自由变化,有的简单(只有几个节点),有的复杂(几十上百个节点)
  • 树的每个节点代表一个操作(函数调用)或者一个值(变量、常量)
  • 树的遍历可以解读成完整的程序

这种表示方法天然支持变长、变结构的解,非常适合程序合成这类问题。


一、遗传编程的理论基础

1.1 自动编程的愿景

传统编程要求人类程序员明确指定解决问题的每一个步骤。Koza提出的遗传编程旨在颠覆这一范式:给定问题的输入-输出示例或性能评估函数,让计算机通过进化过程自动生成解决问题的程序。

这一愿景的理论支撑来自两个核心原则:

闭合性原则(Closure):所有函数和终止符必须能够相互组合而不出错。例如,若函数集中包含返回布尔值的函数,则其他函数必须能够接受布尔输入。

充分性原则(Sufficiency):函数集和终止符集的组合必须能够表示任意复杂度的解空间。

1.2 与遗传算法的区别

特征遗传算法(GA)遗传编程(GP)
个体表示固定长度串可变大小树
基因含义编码决定自身即是程序
搜索空间离散、线性层次化、组合爆炸
应用领域参数优化程序合成

二、树形表示与语法结构

2.1 表达树的基本概念

GP中的个体是一棵有根有序树。每个内部节点代表一个函数(运算符),每个叶子节点代表一个终止符(变量或常量)。

数学表达:程序可形式化为从输入空间到输出空间的映射:

其中为输入变量空间,为输出空间。

2.2 表达树的构建示例

以数学表达式为例,其对应的表达树结构如下:

           ×
          / \
         +   -
        / \ / \
       x  3 y  2

中序遍历(标准数学表达式):

前序遍历(波兰表达式):

后序遍历(逆波兰表达式):

2.3 LISP的S表达式与树结构的渊源

你可能会好奇:为什么GP偏偏选了树结构?这要从Koza做GP时的编程环境说起。

1980年代末,Koza在斯坦福做研究时广泛使用LISP语言。LISP的核心数据结构就是S表达式(Symbolic Expression),本质上就是嵌套的列表和树。

比如这个LISP代码:

(* (+ x 3) (- y 2))

翻译成树就是:

       *
      / \
     +   -
    / \ / \
   x  3 y  2

所以GP选择树结构,有几分”历史偶然”的味道——Koza熟悉LISP,就顺手用了S表达式来表示程序。但这个选择意外地合适:树结构既能表达程序的层次结构,又便于定义交叉、变异等遗传操作。

后来有些研究者尝试用其他表示方法(比如线性编码、图结构),但树结构仍然是GP的主流选择。

2.4 树深度与复杂度控制

GP树的深度直接影响生成的程序复杂度:

  • 最大深度限制:防止”不可控生长”(bloat)现象
  • 初始深度控制:使用特定初始化方法控制初始种群

深度限制策略

class TreeNode:
    def __init__(self, value, depth=0):
        self.value = value
        self.depth = depth
        self.children = []
    
    def is_valid(self, max_depth):
        if self.depth > max_depth:
            return False
        return all(child.is_valid(max_depth) for child in self.children)
    
    def count_nodes(self):
        return 1 + sum(child.count_nodes() for child in self.children)
    
    def get_leaves(self):
        if not self.children:
            return [self.value]
        return [leaf for child in self.children for leaf in child.get_leaves()]
    
    def get_subtree(self, target_node):
        """获取以目标节点为根的子树"""
        if self == target_node:
            return self
        for child in self.children:
            result = child.get_subtree(target_node)
            if result:
                return result
        return None
    
    def replace_with(self, target_node, replacement):
        """将目标节点替换为新节点"""
        for i, child in enumerate(self.children):
            if child == target_node:
                self.children[i] = replacement
                return True
            if child.replace_with(target_node, replacement):
                return True
        return False

三、函数集与终止符集

3.1 函数集的设计

函数集应包含解决问题所需的所有运算符。典型的函数类别包括:

3.1.1 算术函数

function_set = {
    '+': lambda x, y: x + y,
    '-': lambda x, y: x - y,
    '*': lambda x, y: x * y,
    '/': lambda x, y: x / y if abs(y) > 1e-8 else 1,  # 保护性除法
    '%': lambda x, y: x % y if y != 0 else 0,
    '**': lambda x, y: x ** y if x > 0 else 0,
}

3.1.2 布尔函数

对于布尔问题,函数集通常包括:

boolean_functions = {
    'AND': lambda x, y: x and y,
    'OR': lambda x, y: x or y,
    'NOT': lambda x: not x,
    'IF': lambda cond, x, y: x if cond else y,
}

3.1.3 条件与循环函数

control_functions = {
    'IF-LESS': lambda x, y, z, w: z if x < y else w,
    'IF-EQUAL': lambda x, y, z, w: z if x == y else w,
    'PROGN2': lambda x, y: y,  # 顺序执行
    'PROGN3': lambda x, y, z: z,
}

3.2 终止符集的设计

终止符集提供程序的输入和常量:

terminal_set = {
    'variables': ['x', 'y', 'z', 't'],
    'constants': list(range(-10, 11)),
    'ephemeral_random': lambda: random.uniform(-1, 1),  # 临时随机常量
    'funcalls': ['sin', 'cos', 'exp', 'log', 'sqrt'],
}

Note

临时随机常量(Ephemeral Random Constants):在初始化时生成随机常量值,交叉和变异操作可能改变这些常量但不会产生新值。这是GP中处理数值常量的经典技术。

3.3 闭合性检查

def check_closure(function_set, terminal_set):
    """
    验证函数集的闭合性
    每个函数的所有参数位置必须能够接受任何函数的返回值
    """
    for name, func in function_set.items():
        arity = func.__code__.co_argcount
        for i in range(arity):
            # 检查第i个参数位置
            compatible = any(
                is_compatible(output_type, param_types[i])
                for func2 in function_set.values()
                for output_type in get_output_types(func2)
            )
            # 也必须接受终止符
            compatible |= any(
                is_compatible(t, param_types[i])
                for t in terminal_set
            )
            if not compatible:
                print(f"Closure violation: {name} at position {i}")
    return True

四、遗传操作算子

4.1 初始化方法

GP的初始化决定了种群的”基因池”丰富程度。如果初始种群太单一,进化就容易陷入局部最优。

4.1.1 完整方法(Full Method)

所有叶子节点达到指定深度,确保树的完全性:

def full_initialization(depth, function_set, terminal_set):
    if depth == 0:
        return random.choice(terminal_set)
    return TreeNode(
        random.choice(function_set),
        children=[full_initialization(depth-1, function_set, terminal_set) 
                  for _ in range(arity(random.choice(function_set)))]
    )

完整方法生成的树每一层都塞满了函数节点,看起来”整齐”但缺乏多样性。

4.1.2 生长方法(Grow Method)

允许叶子节点在任何深度出现,产生更自然的形状:

def grow_initialization(max_depth, function_set, terminal_set):
    if max_depth == 0 or random.random() < ratio_terminals:
        return random.choice(terminal_set)
    func = random.choice(function_set)
    return TreeNode(
        func,
        children=[grow_initialization(max_depth-1, function_set, terminal_set)
                  for _ in range(arity(func))]
    )

生长方法能生成各种形状的树:有的深、有的浅、有的”歪瓜裂枣”。这样种群更多样。

4.1.3 混合方法(Ramped Half-and-Half)

在实际应用中,我们通常混合使用完整法和生长法:

def ramped_initialization(pop_size, min_depth=2, max_depth=6, 
                          function_set, terminal_set):
    population = []
    depth_range = range(min_depth, max_depth + 1)
    
    for i in range(pop_size):
        depth = list(depth_range)[i % len(depth_range)]
        if i % 2 == 0:
            tree = full_initialization(depth, function_set, terminal_set)
        else:
            tree = grow_initialization(depth, function_set, terminal_set)
        population.append(tree)
    
    return population

Tip

Ramped Half-and-Half能产生多样化的树结构,是Koza推荐的标准初始化方法。

4.2 子树交叉(Crossover)

交叉是GP中最重要的遗传操作,相当于有性繁殖——两个”父母”各拿出一部分基因,组合成”孩子”。

4.2.1 标准子树交叉

随机选择两个父代树中的交叉点,交换对应的子树:

import copy
 
def subtree_crossover(parent1, parent2, max_depth=None):
    # 随机选择交叉点
    crossover_point1 = random.choice(get_all_nodes(parent1))
    crossover_point2 = random.choice(get_all_nodes(parent2))
    
    # 深度检查
    if max_depth:
        subtree_depth = get_subtree_depth(crossover_point2)
        parent_depth = get_node_depth(parent1, crossover_point1)
        if parent_depth + subtree_depth > max_depth:
            return copy.deepcopy(parent1), copy.deepcopy(parent2)
    
    # 执行交叉
    child1 = copy.deepcopy(parent1)
    child2 = copy.deepcopy(parent2)
    
    replace_subtree(child1, crossover_point1, copy.deepcopy(crossover_point2))
    replace_subtree(child2, crossover_point2, copy.deepcopy(crossover_point1))
    
    return child1, child2

可视化示例

Parent1:                        Parent2:
      +                               /
     / \                             / \
    ×   3                         sin   +
   / \                             |   / \
  x   2                           y   z   w

Crossover point1: ×              Crossover point2: sin

Child1:                          Child2:
      +                               /
     / \                             / \
   sin   3                          ×   +
   |                               / \   \
  y                              x   2  z   w

子树交叉的直观理解:把父代1的某个分支切下来,用父代2的对应分支替换。这样两个孩子各继承了一部分”基因”。

4.2.2 点交叉(Point Crossover)

保留大部分树结构,仅交换一个节点:

def point_crossover(parent1, parent2):
    node1 = random.choice(get_all_nodes(parent1))
    node2 = random.choice(get_all_nodes(parent2))
    
    child1 = copy.deepcopy(parent1)
    child2 = copy.deepcopy(parent2)
    
    replace_node(child1, node1, copy.deepcopy(node2))
    replace_node(child2, node2, copy.deepcopy(node1))
    
    return child1, child2

点交叉保守一些——只换了一个节点,改变较小。适合当你想微调解的时候用。

4.3 变异操作(Mutation)

变异相当于无性繁殖中的”基因突变”,给种群引入新东西,防止进化陷入死胡同。

4.3.1 子树变异

随机选择一棵子树,用新生成的随机树替换:

def subtree_mutation(tree, function_set, terminal_set, max_depth=6):
    child = copy.deepcopy(tree)
    mutation_point = random.choice(get_all_nodes(child))
    
    # 生成新的随机子树
    depth = random.randint(0, max_depth)
    new_subtree = grow_initialization(depth, function_set, terminal_set)
    
    replace_subtree(child, mutation_point, new_subtree)
    return child

子树变异是比较”激进”的变异,可能带来大的改变。

4.3.2 收缩变异(Shrink Mutation)

有时候一棵大子树其实没太大用,收缩变异就是把它”压缩”成一个简单的叶子:

def shrink_mutation(tree, terminal_set):
    child = copy.deepcopy(tree)
    mutation_point = random.choice(get_non_leaf_nodes(child))
    
    replace_node(child, mutation_point, random.choice(terminal_set))
    return child

比如把sin(x) * (x + y - z)中的(x + y - z)直接替换成x。这招对付”代码膨胀”特别有效。

4.3.3 Hoist变异

从父代树中随机选一个节点,把它”提拔”成新的根:

def hoist_mutation(tree):
    child = copy.deepcopy(tree)
    mutation_point = random.choice(get_all_nodes(child))
    
    # 用选中的节点创建新树
    return copy.deepcopy(mutation_point)

Hoist变异的好处是能快速”瘦身”,因为新树肯定比原来的小。这也是一种对付bloat的方法。

4.3.4 点变异(Point Mutation)

仅改变单个节点的值,最保守的变异:

def point_mutation(tree, function_set, terminal_set, rate=0.01):
    child = copy.deepcopy(tree)
    
    for node in get_all_nodes(child):
        if random.random() < rate:
            if node.is_function():
                node.value = random.choice(function_set)
            else:
                node.value = random.choice(terminal_set)
    
    return child

点变异只动一个节点,副作用小,适合细调。


五、适应度评估

5.1 适应度案例(Fitness Cases)

GP的适应度通过在多个测试案例上的表现来衡量。适应度函数的设计直接决定GP能学到什么。

class FitnessEvaluator:
    def __init__(self, fitness_cases):
        """
        fitness_cases: [(input1, input2, ..., expected_output), ...]
        """
        self.fitness_cases = fitness_cases
    
    def evaluate(self, tree):
        errors = []
        for inputs, expected in self.fitness_cases:
            try:
                result = tree.execute(*inputs)
                errors.append(abs(result - expected))
            except:
                errors.append(float('inf'))
        
        return {
            'raw_fitness': sum(errors),
            'mean_absolute_error': sum(errors) / len(errors),
            'hits': sum(1 for e in errors if e < 0.01),
        }

5.2 标准化适应度

对于符号回归问题,标准化均方误差(Normalized Mean Squared Error):

标准化版本

Warning

过拟合是GP面临的重要问题。需要使用独立的验证集来评估泛化能力。

5.3 适应度设计的经验

  1. 尽量用归一化误差:这样不同问题之间有可比性,也方便设定终止条件。

  2. 处理无效程序:GP经常产生除零、溢出、越界访问这类程序,要给它们一个很差的适应度(比如无穷大),而不是直接让程序崩溃。

  3. 多目标适应度:有时候我们不只想效果好,还想要代码简洁。这时候可以把多个目标加权组合,或者用Pareto多目标进化。


六、Koza的GP五大步骤

Koza在《Genetic Programming》中提出了GP的标准流程,这五步基本成了GP的”教科书流程”。

步骤1:确定终结符集

终结符集包括:

  • 输入变量:x, y, z等
  • 数值常量:0, 1, 2或随机生成的数
  • 零元函数:返回常数的函数

步骤2:确定函数集

根据问题领域选择函数,比如:

  • 数值问题:+, -, *, /, sin, cos, exp, log
  • 布尔问题:AND, OR, NOT, IF
  • 分类问题:决策阈值函数

关键是保证闭合性——函数之间能自由嵌套。

步骤3:确定适应度函数

适应度函数衡量个体好坏,通常是:

  • 误差函数(回归问题)
  • 正确率(分类问题)
  • 自定义评分(设计问题)

步骤4:确定控制参数

包括:

  • 种群大小(通常500-1000)
  • 交叉概率(通常80-90%)
  • 变异概率(通常10-20%)
  • 树深度限制(通常6-17)
  • 最大代数(通常50-500)

步骤5:确定终止条件和结果判定

终止条件通常是:

  • 达到最大代数
  • 找到完美解(适应度=0)
  • 适应度停滞

七、自动函数发现

7.1 ADF(Automatically Defined Functions)

Koza提出的ADF机制允许进化过程中自动定义和调用子函数:

class ADFTree:
    def __init__(self, main_tree, adf_definitions):
        self.main_tree = main_tree  # 主程序
        self.adf_definitions = adf_definitions  # {name: Tree}
    
    def execute(self, inputs):
        # 准备ADF参数
        adf_args = {f'ARG{i}': inputs[i] for i in range(len(inputs))}
        
        # 计算ADF
        for name, tree in self.adf_definitions.items():
            adf_args[name] = tree.execute(inputs)
        
        # 执行主程序
        return self.main_tree.execute(adf_args)

ADF树结构示例

ADF1(x, y):                       Main Program:
      +                              IF-ZERO
     / \                            /   |   \
    ×   z                          >   ADF1  ADF2
   / \                              |   |     |
  x   y                           x   y     z

ADF的价值在于:如果某个子问题在程序中反复出现,ADF可以让进化”发现”这个子问题,抽取成独立函数,避免重复表达。

7.2 模块化GP架构

class ModularGP:
    def __init__(self, n_adfs=3):
        self.n_adfs = n_adfs
        self.adf_functions = {}
    
    def evolve_with_adfs(self, fitness_cases):
        # 进化主程序和ADF定义
        pass
    
    def call_adf(self, adf_name, args):
        return self.adf_functions[adf_name].execute(args)

八、遗传编程的问题:代码膨胀(Bloat)

8.1 什么是Bloat?

Bloat,字面意思是”臃肿”,指的是GP进化过程中程序树可能无限增长,但适应度却不再提升。

举个例子:你的程序进化到某个程度,已经能95%解决问题了。剩下的5%很难通过”微调”解决,但树还在不断生长——因为交叉操作倾向于产生更大的树。

大树的”生存策略”:

  • 包含很多”中性”代码——删掉不影响功能,但留着也不碍事
  • 这些代码可以当作”跳板”,让后续进化更容易找到好的变种
  • 短视的适应度函数看不到大树的问题

Bloat会让GP变慢、占内存、最后解的可读性也很差。

8.2 深度限制策略

最简单的对策:限制树的最大深度。

class DepthLimitedGP:
    def __init__(self, max_depth=17):
        self.max_depth = max_depth
    
    def is_valid_tree(self, tree):
        return get_tree_depth(tree) <= self.max_depth
    
    def crossover_with_limit(self, parent1, parent2):
        child1, child2 = subtree_crossover(parent1, parent2)
        if not self.is_valid_tree(child1):
            return copy.deepcopy(parent1)
        if not self.is_valid_tree(child2):
            return copy.deepcopy(parent2)
        return child1, child2

限制深度的代价是可能限制进化的上限——有时候真正的解就需要很深的树。

8.3 惩罚机制

给”大块头”程序加点惩罚:

def size_penalized_fitness(raw_fitness, tree, penalty_coef=0.001):
    complexity = tree.count_nodes()
    return raw_fitness + penalty_coef * complexity

这样同样效果的程序,小的优于大的。

8.4 多目标优化

把适应度和大小拆成两个目标,用NSGA-II之类的多目标算法来进化:

from deap import tools
 
def evolve_multi_objective(population, fitness_cases):
    creator.create("FitnessMulti", base.Fitness, 
                    weights=(1.0, -1.0))  # 最大化适应度,最小化大小
    creator.create("Individual", list, fitness=creator.FitnessMulti)
    
    toolbox = base.Toolbox()
    # ... setup operators ...
    
    # NSGA-II多目标优化
    population = algorithms.eaMuPlusLambda(
        population, toolbox, mu=500, lambda_=500, 
        cxpb=0.9, mutpb=0.1, ngen=100
    )
    
    return population

多目标方法的好处是能得到一组Pareto最优解——从小到大各种复杂度的”够用解”都有。


九、经典应用案例

9.1 符号回归(Symbolic Regression)

符号回归是GP的招牌应用:给你一堆数据点,让你找到生成这些数据的数学公式。

def symbolic_regression_example():
    """
    目标函数: f(x) = sin(x) + x^2
    """
    # 生成训练数据
    x_train = np.linspace(-10, 10, 100)
    y_train = np.sin(x_train) + x_train**2
    
    # 定义函数集和终止符集
    function_set = ['+', '-', '*', '/', 'sin', 'cos', 'exp', 'log']
    terminal_set = ['x', *np.random.uniform(-1, 1, 10)]
    
    # 运行GP
    gp = GeneticProgramming(function_set, terminal_set, fitness_cases)
    best_solution = gp.evolve()
    
    print(f"Discovered: {best_solution.to_expression()}")
    # 输出类似: sin(x) + (x * x)

实战技巧

  1. 数据归一化:输入输出都做归一化,能大幅提升效果
  2. 函数集要”够用但不冗余”:sin和cos在很多场景都有用,但别加太多
  3. 保护性除法:防止除以零或者很小的数
  4. ephemeral随机常量:在(-1, 1)范围内生成,用evolution来微调

9.2 电路自动设计

Koza在《Genetic Programming III》中展示了GP自动设计电路的能力,包括:

  • 电路拓扑发现
  • 元件参数优化
  • 信号处理滤波器设计
class CircuitDesignGP:
    def __init__(self, circuit_constraints):
        self.constraints = circuit_constraints
        self.available_components = ['R', 'L', 'C', 'Diode', 'Transistor']
    
    def evaluate_circuit(self, circuit_tree):
        """
        评估电路性能指标
        - 频率响应
        - 功耗
        - 噪声特性
        """
        simulation_result = self.simulate(circuit_tree)
        return {
            'fitness': fitness_score(simulation_result),
            'validity': self.check_constraints(circuit_tree),
        }

GP设计的电路有多离谱?Koza团队用它设计出了比人类设计的更好的滤波器专利——当然,是在一堆人类已知的电路结构基础上做搜索,不是凭空创造。

9.3 游戏策略演化

class GameStrategyGP:
    def __init__(self, game_state_space):
        self.game_states = game_state_space
    
    def create_game_tree(self):
        function_set = ['IF-STATE', 'AND', 'OR', 'NOT']
        terminal_set = [f'S{i}' for i in range(len(self.game_states))]
        return grow_initialization(5, function_set, terminal_set)
    
    def evaluate_strategy(self, strategy_tree):
        wins = 0
        for game_state in self.game_states:
            action = strategy_tree.execute(game_state)
            if self.play_and_win(action, game_state):
                wins += 1
        return wins / len(self.game_states)

十、几何语义遗传编程

10.1 为什么需要语义GP?

传统的子树交叉是”语法层面”的操作——我们只看树的形状,不看树计算出来的值。但有时候两个语法差异很大的树,语义上却非常接近。

比如:

  • 树A:sin(x) * sin(x)
  • 树B:(x - sin(x) + x) * sin(x) / 2 - x * sin(x) / 2

从语法上看,这两棵树完全不同。但如果你把它们化简一下,会发现它们的语义其实很接近。

几何语义GP就是想解决这个问题——直接在”语义空间”里做操作。

10.2 几何语义交叉

Moraglio等人提出的几何语义交叉有个神奇的性质:保证子代在训练点上的表现不劣于任一父代

怎么做到的?核心思想是构造一个”混合树”:

def geometric_semantic_crossover(parent1, parent2, training_data):
    """
    几何语义交叉
    """
    # 构造混合树: T1 * R + T2 * (1 - R)
    # R是一个随机产生的"掩码树"
    r_tree = generate_random_tree()
    
    # 混合
    mixed_tree = create_mix_tree(parent1, parent2, r_tree)
    return mixed_tree

直觉上:这相当于在两个函数的输出之间做插值。如果两个父代都还行,插值出来的也不会太差。

10.3 几何语义变异

语义变异通过添加”随机方向”来探索语义空间:

def geometric_semantic_mutation(tree, mutation_step=0.5, training_data):
    """
    语义变异
    """
    # 生成随机目标点
    random_target = random.choice(training_data)
    
    # 生成随机子树
    random_subtree = generate_random_tree()
    
    # 构造变异树: T + MR
    mutated = create_mutation_tree(tree, random_subtree, mutation_step)
    return mutated

十一、模因计算(Meme Computation)

11.1 什么是模因计算?

模因计算是个有意思的概念——借鉴达尔文以前的”获得性遗传”思想。

在经典GP里,我们只在种群内做选择和繁殖。但模因计算允许个体从”文化知识库”里学习——把其他个体的高质量子树”借鉴”过来。

这类似于:

  • 经典进化:只能从父母那里继承基因
  • 模因计算:还能从”文化”中学习(类似人类的学习行为)

11.2 模因库的实现

class MemeLibrary:
    """
    模因库:存储进化过程中发现的好"想法"
    """
    def __init__(self, max_size=100):
        self.memes = []  # (subtree, fitness)
        self.max_size = max_size
    
    def add_meme(self, subtree, fitness):
        """添加一个模因"""
        self.memes.append((subtree, fitness))
        self.memes.sort(key=lambda x: -x[1])  # 按适应度排序
        self.memes = self.memes[:self.max_size]  # 淘汰差的
    
    def get_random_meme(self, tournament_size=3):
        """随机获取一个还算不错的模因"""
        candidates = random.sample(self.memes, 
                                   min(tournament_size, len(self.memes)))
        return max(candidates, key=lambda x: x[1])[0]
    
    def mutate_with_meme(self, tree, probability=0.1):
        """以一定概率用模因替换树的一部分"""
        if random.random() < probability:
            meme = self.get_random_meme()
            mutation_point = random.choice(get_all_nodes(tree))
            replace_subtree(tree, mutation_point, copy.deepcopy(meme))
        return tree

模因计算的效果取决于模因库的质量——如果库里有垃圾模因,可能会拖后腿。


十二、从零实现遗传编程:完整代码

好了,理论说够了,让我们来写一个能跑的GP实现。

import random
import copy
import operator
import math
 
class TreeNode:
    """表达式树节点"""
    def __init__(self, value):
        self.value = value
        self.children = []
    
    def __str__(self):
        if not self.children:
            return str(self.value)
        if len(self.children) == 1:
            return f"{self.value}({self.children[0]})"
        if len(self.children) == 2:
            return f"({self.children[0]} {self.value} {self.children[1]})"
        # 多参数函数
        args = ", ".join(str(c) for c in self.children)
        return f"{self.value}({args})"
    
    def evaluate(self, context=None):
        """计算树的值"""
        context = context or {}
        
        # 如果是叶子节点
        if not self.children:
            # 变量
            if isinstance(self.value, str) and self.value in context:
                return context[self.value]
            # 常量
            return self.value
        
        # 函数节点
        args = [child.evaluate(context) for child in self.children]
        
        # 内置函数
        if self.value == '+':
            return args[0] + args[1]
        elif self.value == '-':
            return args[0] - args[1]
        elif self.value == '*':
            return args[0] * args[1]
        elif self.value == '/':
            # 保护性除法
            if abs(args[1]) < 1e-10:
                return 1.0
            return args[0] / args[1]
        elif self.value == 'sin':
            return math.sin(args[0])
        elif self.value == 'cos':
            return math.cos(args[0])
        elif self.value == 'exp':
            return math.exp(args[0])
        elif self.value == 'log':
            return math.log(abs(args[0]) + 1e-10)
        elif self.value == 'sqrt':
            return math.sqrt(abs(args[0]))
        else:
            raise ValueError(f"Unknown function: {self.value}")
    
    def copy(self):
        """深拷贝"""
        new_node = TreeNode(self.value)
        new_node.children = [child.copy() for child in self.children]
        return new_node
    
    def get_all_nodes(self):
        """获取所有节点"""
        nodes = [self]
        for child in self.children:
            nodes.extend(child.get_all_nodes())
        return nodes
    
    def count_nodes(self):
        """节点数量"""
        return 1 + sum(child.count_nodes() for child in self.children)
    
    def depth(self):
        """树深度"""
        if not self.children:
            return 1
        return 1 + max(child.depth() for child in self.children)
 
 
class GeneticProgramming:
    """遗传编程主类"""
    
    def __init__(self, function_set, terminal_set, fitness_cases,
                 max_depth=8, pop_size=500, elite_size=10):
        self.function_set = function_set
        self.terminal_set = terminal_set
        self.fitness_cases = fitness_cases
        self.max_depth = max_depth
        self.pop_size = pop_size
        self.elite_size = elite_size
        
        # 函数元数
        self.arity = {
            '+': 2, '-': 2, '*': 2, '/': 2,
            'sin': 1, 'cos': 1, 'exp': 1, 'log': 1, 'sqrt': 1
        }
    
    def generate_random_tree(self, max_depth, grow=True):
        """生成随机树"""
        if max_depth == 0 or (grow and random.random() < 0.3):
            return TreeNode(random.choice(self.terminal_set))
        
        func = random.choice(self.function_set)
        node = TreeNode(func)
        n_children = self.arity.get(func, 0)
        
        for _ in range(n_children):
            child = self.generate_random_tree(max_depth - 1, grow)
            node.children.append(child)
        
        return node
    
    def generate_initial_population(self):
        """生成初始种群"""
        population = []
        for i in range(self.pop_size):
            # Ramped half-and-half
            depth = random.randint(2, self.max_depth)
            tree = self.generate_random_tree(depth, grow=(i % 2 == 0))
            population.append(tree)
        return population
    
    def evaluate_fitness(self, tree):
        """评估适应度(均方误差)"""
        total_error = 0
        
        for inputs, expected in self.fitness_cases:
            try:
                predicted = tree.evaluate(inputs)
                total_error += (predicted - expected) ** 2
            except:
                total_error += 1e10
        
        return total_error / len(self.fitness_cases)
    
    def tournament_selection(self, population, fitnesses, k=3):
        """锦标赛选择"""
        tournament = random.sample(range(len(population)), k)
        best_idx = min(tournament, key=lambda i: fitnesses[i])
        return population[best_idx].copy()
    
    def crossover(self, parent1, parent2):
        """子树交叉"""
        child1 = parent1.copy()
        child2 = parent2.copy()
        
        nodes1 = child1.get_all_nodes()
        nodes2 = child2.get_all_nodes()
        
        point1 = random.choice(nodes1)
        point2 = random.choice(nodes2)
        
        # 找到父节点以便替换
        def replace_in_tree(tree, old_node, new_node):
            for i, child in enumerate(tree.children):
                if child is old_node:
                    tree.children[i] = new_node.copy()
                    return True
                if replace_in_tree(child, old_node, new_node):
                    return True
            return False
        
        # 确保新子树不超过深度限制
        new_subtree = point2.copy()
        if new_subtree.depth() + self._get_node_depth(child1, point1) <= self.max_depth:
            replace_in_tree(child1, point1, new_subtree)
            replace_in_tree(child2, point2, point1.copy())
        
        return child1, child2
    
    def _get_node_depth(self, tree, target):
        """获取目标节点的深度"""
        def find_depth(node, current_depth):
            if node is target:
                return current_depth
            for child in node.children:
                depth = find_depth(child, current_depth + 1)
                if depth != -1:
                    return depth
            return -1
        return find_depth(tree, 1)
    
    def mutate(self, tree):
        """变异"""
        child = tree.copy()
        
        # 随机选择一个节点
        nodes = child.get_all_nodes()
        point = random.choice(nodes)
        
        # 用新树替换
        new_subtree = self.generate_random_tree(
            random.randint(1, self.max_depth // 2)
        )
        
        def replace_in_tree(tree, old_node, new_node):
            for i, child in enumerate(tree.children):
                if child is old_node:
                    tree.children[i] = new_node.copy()
                    return True
                if replace_in_tree(child, old_node, new_node):
                    return True
            return False
        
        replace_in_tree(child, point, new_subtree)
        return child
    
    def evolve(self, n_generations=100, crossover_rate=0.8, 
               mutation_rate=0.1, verbose=True):
        """主进化循环"""
        
        # 初始化
        population = self.generate_initial_population()
        best_ever = None
        best_ever_fitness = float('inf')
        
        for gen in range(n_generations):
            # 评估
            fitnesses = [self.evaluate_fitness(tree) for tree in population]
            
            # 找精英
            elite_indices = sorted(range(len(fitnesses)), 
                                  key=lambda i: fitnesses[i])[:self.elite_size]
            elites = [population[i].copy() for i in elite_indices]
            
            # 记录最佳
            best_idx = elite_indices[0]
            if fitnesses[best_idx] < best_ever_fitness:
                best_ever = population[best_idx].copy()
                best_ever_fitness = fitnesses[best_idx]
            
            if verbose and gen % 10 == 0:
                avg_fitness = sum(fitnesses) / len(fitnesses)
                print(f"Gen {gen}: Best={best_ever_fitness:.6f}, "
                      f"Avg={avg_fitness:.6f}, "
                      f"Best tree: {best_ever}")
            
            # 检查终止条件
            if best_ever_fitness < 1e-10:
                print(f"Perfect solution found at generation {gen}!")
                break
            
            # 生成新一代
            new_population = elites.copy()  # 保留精英
            
            while len(new_population) < self.pop_size:
                if random.random() < crossover_rate:
                    # 交叉
                    p1 = self.tournament_selection(population, fitnesses)
                    p2 = self.tournament_selection(population, fitnesses)
                    c1, c2 = self.crossover(p1, p2)
                    new_population.append(c1)
                    if len(new_population) < self.pop_size:
                        new_population.append(c2)
                else:
                    # 变异或复制
                    if random.random() < mutation_rate:
                        p = self.tournament_selection(population, fitnesses)
                        new_population.append(self.mutate(p))
                    else:
                        p = self.tournament_selection(population, fitnesses)
                        new_population.append(p.copy())
            
            population = new_population[:self.pop_size]
        
        return best_ever
 
 
def demo():
    """演示:用GP做符号回归"""
    
    # 目标函数:sin(x) + x^2
    def target(x):
        return math.sin(x) + x**2
    
    # 生成训练数据
    x_values = [x / 10 for x in range(-50, 51)]
    y_values = [target(x) for x in x_values]
    fitness_cases = list(zip([{'x': x} for x in x_values], y_values))
    
    # 定义函数集和终止符集
    function_set = ['+', '-', '*', '/', 'sin', 'cos']
    terminal_set = ['x', 0, 1, 2, 3, 
                   random.uniform(-1, 1) for _ in range(5)]
    
    # 跑GP
    gp = GeneticProgramming(
        function_set=function_set,
        terminal_set=terminal_set,
        fitness_cases=fitness_cases,
        max_depth=6,
        pop_size=300
    )
    
    best = gp.evolve(n_generations=50)
    print(f"\n最终结果: {best}")
    print(f"最终适应度: {gp.evaluate_fitness(best):.6f}")
    
    # 测试
    test_x = 2.5
    actual = target(test_x)
    predicted = best.evaluate({'x': test_x})
    print(f"\n测试 x={test_x}: 实际={actual:.4f}, 预测={predicted:.4f}")
 
 
if __name__ == "__main__":
    demo()

这段代码虽然简单,但覆盖了GP的核心要素:

  • 树结构表示
  • Ramped初始化
  • 锦标赛选择
  • 子树交叉
  • 变异
  • 精英保留

跑起来试试,你会发现它确实能发现sin(x) + x^2这样的公式——当然,也可能找到等价形式,比如x^2 + sin(x)sin(x) + x*x之类的。


十三、实用框架与工具

13.1 DEAP框架

虽然手写GP能帮你理解原理,但真想干活还得用成熟的框架。Python里最流行的GP库是DEAP。

from deap import base, creator, tools, algorithms, gp
import operator
import random
import math
 
def deap_gp_complete_example():
    """DEAP GP完整示例:符号回归"""
 
    # 定义原始集
    pset = gp.PrimitiveSetTyped("MAIN", [float], float)
    pset.renameArguments(ARG0='x')
 
    # 保护性除法
    def protectedDiv(left, right):
        try:
            return left / right
        except ZeroDivisionError:
            return 1
 
    # 添加算术运算
    pset.addPrimitive(operator.add, [float, float], float, name="add")
    pset.addPrimitive(operator.sub, [float, float], float, name="sub")
    pset.addPrimitive(operator.mul, [float, float], float, name="mul")
    pset.addPrimitive(protectedDiv, [float, float], float, name="div")
    pset.addPrimitive(operator.neg, [float], float, name="neg")
    pset.addPrimitive(math.cos, [float], float, name="cos")
    pset.addPrimitive(math.sin, [float], float, name="sin")
 
    # 添加常量
    pset.addEphemeralConstant("rand100", lambda: random.uniform(-100, 100), float)
    pset.addTerminal(1.0, float)
 
    # 创建适应度类(最小化)
    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_=2, max_=6)
    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):
        func = toolbox.compile(expr=individual)
        x = [random.uniform(-1, 1) for _ in range(100)]
        sqerrors = []
        for xi in x:
            try:
                y = xi**4 + xi**3 + xi**2 + xi  # 目标函数
                sqerrors.append((func(xi) - y) ** 2)
            except Exception:
                sqerrors.append(1000000)
        return (math.fsum(sqerrors) / len(x),)
 
    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)
 
    # 统计
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", math.mean)
    stats.register("min", min)
 
    # 运行
    pop = toolbox.population(n=300)
    algorithms.eaSimple(pop, toolbox, cxpb=0.9, mutpb=0.1, 
                       ngen=40, stats=stats, verbose=True)
 
    return pop, stats

十四、参数配置建议

GP的效果对参数设置挺敏感的,这里给点经验值:

参数推荐范围说明
种群大小500-1000小问题500够用,复杂问题可能需要1000+
最大深度6-17Koza早期用17,后来发现6-8就够了
交叉率80-90%GP中交叉比变异重要
变异率10-20%主要探索手段
锦标赛规模3-7大了选择压力强,小了多样性高
精英数1-10保留最好的解不被淘汰

实际调参的时候,先用小种群快速测试(省时间),找到靠谱的设置后再放大。


十五、常见问题与解决方案

问题1:过拟合

症状:训练集上效果很好,测试集上效果差。

原因:GP可能找到了”记住”训练数据的解,而不是真正的规律。

解决

  • 增加训练数据量
  • 用独立的验证集评估
  • 加复杂度惩罚
  • 减少最大深度

问题2:早熟收敛

症状:种群多样性快速消失,适应度停滞不前。

原因:选择压力太大,交叉产生的子代不如父代。

解决

  • 增加变异率
  • 用稳态进化而非代际进化
  • 引入小生境(niche)机制

问题3:无效程序

症状:很多程序产生NaN、Inf或除零错误。

原因:数值运算产生溢出或边界情况。

解决

  • 保护性除法
  • 限制指数函数的输入
  • 给异常程序很差的适应度

十六、进阶主题

16.1 强类型GP

普通GP里函数可以接受任何类型的参数,但强类型GP要求类型匹配。

# DEAP强类型示例
pset = gp.PrimitiveSetTyped("MAIN", [float, float], float)
pset.renameArguments(ARG0='x', ARG1='y')
 
# 只能接受两个浮点数的加法
pset.addPrimitive(operator.add, [float, float], float)
# 只能接受浮点数的非操作
pset.addPrimitive(operator.neg, [float], float)

强类型GP的好处是语法永远有效,减少了无效程序。

16.2 语法演化(Grammatical Evolution)

语法演化把”语法规则”和”进化搜索”分开:

# BNF语法
grammar = """
<expr> ::= <expr><op><expr> | (<expr><op><expr>) | <var>
<op> ::= + | - | * | /
<var> ::= x | y | 0 | 1 | 2
"""
 
# 进化搜索的是整数序列
# 用整数序列通过语法规则生成程序

这样可以用进化算法搜”程序编码”,而不是直接搜程序结构。

16.3 线性GP

不用树,用线性序列表示程序:

# 线性GP的"染色体"
chromosome = [
    ('LOAD', 'x'),      # 加载x
    ('SIN', None),      # sin操作
    ('LOAD', 'x'),      # 加载x
    ('LOAD', 'x'),      # 加载x
    ('MUL', None),      # x * x
    ('ADD', None),      # sin(x) + x*x
]

线性GP便于和遗传算法对接,也有利于并行执行。


十七、参考文献

  1. Koza, J. R. (1992). Genetic Programming: On the Programming of Computers by Means of Natural Selection. MIT Press.
  2. Koza, J. R. (1994). Genetic Programming II: Automatic Discovery of Reusable Programs. MIT Press.
  3. Koza, J. R., et al. (1999). Genetic Programming III: Darwinian Invention and Problem Solving. Morgan Kaufmann.
  4. Poli, R., Langdon, W. B., & McPhee, N. F. (2008). A Field Guide to Genetic Programming. Lulu.com.
  5. Banzhaf, W., et al. (1998). Genetic Programming: An Introduction. Morgan Kaufmann.
  6. Langdon, W. B., & Poli, R. (2002). Foundations of Genetic Programming. Springer.
  7. Moraglio, A., et al. (2012). “Geometric Semantic Genetic Programming.” PPSN XII, 21-31.
  8. Vanneschi, L., & Castelli, M. (2019). “A Survey of Semantic Genetic Programming.” Frontiers in Robotics and AI, 6, 37.
  9. Cramer, N. L. (1985). “A Representation for the Adaptive Generation of Simple Sequential Programs.” Proceedings of ICGA, 1985.
  10. Montana, D. J. (1995). “Strongly Typed Genetic Programming.” Evolutionary Computation, 3(2), 199-230.
  11. Ryan, C., & O’Neill, M. (1998). “Grammatical Evolution: Evolving Programs in an Arbitrary Language.” Proceedings of GP, 1998.

相关文档