关键词

进化策略CMA-ES协方差矩阵自然进化策略步长控制
(μ+λ)-ES(μ,λ)-ES重组操作累积步长调整成功比率
权重更新自然梯度信息几何等距采样排名选择
HansenBeyer演化多样性全局收敛超参数自适应

摘要

进化策略(Evolution Strategies, ES)是最早的进化计算范式之一,由Ingo Rechenberg和Hans-Paul Schwefel于1960年代在德国斯图加特大学开创。与遗传算法侧重于离散编码不同,进化策略专门针对连续实数参数的优化问题,采用正态分布采样进行解的搜索。其标志性成果——协方差矩阵自适应进化策略(CMA-ES)——被认为是连续优化领域最成功的算法之一,在各种基准测试和实际应用中展现出卓越的性能。


一、进化策略入门:从遗传算法说起

1.1 进化策略和遗传算法到底什么关系

很多人第一次接触进化策略的时候都会有点懵:这不就是遗传算法吗?确实,两者都是受自然界进化现象启发的优化算法,都用”种群”的概念,都靠选择、变异、重组这些操作来迭代搜索。但如果你仔细看,会发现它们其实有挺大的区别。

遗传算法最早是John Holland那帮人搞出来的,核心是把问题的解编码成一串”基因”——最经典的就是二进制串,比如01101001这种。你做变异就是在这些二进制位上翻翻转转,选择就是让适应度高的个体有更多机会传宗接代。GA特别适合离散组合优化问题,比如背包问题、调度问题、路由问题。

进化策略则是Rechenberg他们在优化喷嘴形状的时候发明的,他们发现流体机械的参数是实数,直接在实数空间里折腾比编码解码要方便得多。所以ES从一开始就走的是实数优化的路线,不太搞二进制编码那套。

打个比方的话,GA像是先用手摇计算机算出二进制答案,再转换成十进制;ES则是直接在十进制空间里操作,省掉了中间那层转换。

1.2 进化策略的核心思想

ES的基本思路其实挺直白的:

  1. 维护一个分布:不是维护一堆显式的个体,而是维护一个概率分布(通常是均值和方差已知的正态分布)

  2. 从这个分布里采样:从分布里掏出几个候选解

  3. 评估这些解:看看每个候选解的适应度怎么样

  4. 更新分布参数:根据评估结果调整分布的参数,让好解出现的地方概率更高

  5. 重复:不断重复上面这个过程,直到找到足够好的解

这么说可能还是有点抽象。让我换个说法——想象你在玩一个扔飞镖的游戏,你的目标是找到墙上哪个位置的分数最高。ES的做法是:

  • 先猜一个大概位置(均值)
  • 先扔几镖看看大概落在哪(采样)
  • 根据分数高低调整你的瞄准位置和扔的力度(更新参数)
  • 继续扔,继续调整

关键是,这个”瞄准位置”和”力度”不是随便调的,而是用统计学的道理来调整,所以能比瞎试高效很多。

1.3 进化策略的”策略参数”是什么

你可能会好奇,ES说的”策略参数”到底指什么。在GA里,基因就是解本身,没什么额外的策略参数。但在ES里,分布的参数(均值、方差、协方差)本身就是你要学习的东西的一部分。

比如说,最简单的ES只有均值和方差两个策略参数。更复杂一点的ES还有协方差矩阵,这个矩阵描述了各个维度之间的相关性。更高级的版本甚至把学习率、选择压力这些也当作策略参数来调整。

这就好比GA的”基因”只是问题的解,而ES的”策略”包含了搜索行为的配置。你可以手动调这些策略参数,也可以让算法自己学怎么调——后者就是CMA-ES这类自适应算法干的事情。


二、(μ+λ)和(μ,λ)策略:两种选择方式的核心区别

2.1 什么是μ和λ

在ES的 notation里,μ表示父代数量,λ表示子代数量。所以(μ+λ)-ES就是用μ个父代产生λ个后代,(μ,λ)-ES也是用μ个父代产生λ个后代。那它们的区别在哪?

2.2 Plus策略(μ+λ):父母子女一起竞争

(μ+λ)-ES的”+“号意思是加上父母代一起竞争。流程是这样的:

  1. 从当前种群(或者从分布)中选μ个父代
  2. 这μ个父代产生λ个子代
  3. 把μ个父代和λ个子代放在一起,总共μ+λ个解
  4. 从这μ+λ个解里选最好的μ个作为下一代的父代

关键点在于:父母代有资格和子女竞争,只要父母足够好,可以一直存活下去。

这种策略的好处是”精英保留”——一旦找到了好解,不会轻易丢掉。对于那些最优解相对稳定的场景,这种策略很有效。但坏处是可能会被困在局部最优,因为旧解会一直占据位置。

2.3 Comma策略(μ,λ):父母代完全退休

(μ,λ)-ES的”,“号意思是不要父母代。流程是这样的:

  1. 从当前种群中选μ个父代
  2. 这μ个父代产生λ个子代(通常λ > μ)
  3. 只看子代,从λ个子代里选最好的μ个
  4. 父母代全部淘汰,下一代的父代全部来自子代

关键点在于:父母代没有资格参与竞争,必须全部交棒给下一代。

这种策略有个很形象的比喻叫”忘记机制”——你必须把过去忘掉,即使父母再好也得退休。这听起来很残酷,但对于需要持续探索的问题特别有用。在非稳态优化问题里(目标函数本身可能随时间变化),Comma策略能帮助算法快速适应变化。

2.4 什么时候用哪个

这个问题没有标准答案,但有一些经验之谈:

  • 如果你的问题比较稳定,最优解不会跑,应该用 (μ+λ)-ES,能更好地利用已有的好解
  • 如果你的问题是非凸的,有很多局部最优,或者目标可能变化,应该用 (μ,λ)-ES,强制探索新区域
  • 一般建议 λ/μ 的比例在 5-10 之间,太小了探索不足,太大了浪费评估次数
  • (μ,λ)-ES还有一个要求:λ要比μ大不少,否则子代数量不够选择压力会太大

2.5 一道思考题

假设你在优化一个随时间缓慢变化的目标函数,今天的最优点在x=5,明天可能就飘到x=7了。如果用(μ+λ)-ES会发生什么?

答案是算法可能一直守着昨天的最优点x=5,不愿意挪窝。但如果用(μ,λ)-ES,父母代全换新一代,那明天的探索就会更加积极。这就是为什么动态环境问题通常偏好Comma策略。


三、从(1+1)-ES到现代ES:一步步进化的故事

3.1 一切从最简单的开始

(1+1)-ES是ES的最简形态,只有一个父代,每次产生一个子代。流程极度简单:

  1. 当前解产生一个变异子代
  2. 如果子代比父代好,就替换;否则父代继续
  3. 重复

这不就是最朴素的爬山法吗?没错!(1+1)-ES本质上就是带随机因素的爬山法。但Rechenberg给它加了一个精妙的设计——1/5成功规则

3.2 1/5成功规则:自适应步长的雏形

在(1+1)-ES里,变异是用正态分布实现的,分布的σ(标准差)就是”步长”。步长太大,会在最优点附近跳来跳去错过最优;步长太小,爬山爬得太慢。

Rechenberg发现,在二维和三维空间里,如果成功率保持在1/5左右,算法就能高效收敛。具体来说:

  • 如果成功率 > 1/5,说明步长太小了,需要调大
  • 如果成功率 < 1/5,说明步长太大了,需要调小

调整的幅度大约是乘以0.85或除以0.85(这个数字来自理论推导)。

这个规则虽然简单,但它是现代ES自适应机制的思想源头——让算法自己学该用多大的步长,而不是靠人调。

3.3 (1+1)-ES的完整实现

class OnePlusOneES:
    def __init__(self, objective, bounds, n_iterations=1000):
        self.objective = objective
        self.bounds = bounds
        self.n_iterations = n_iterations
        self.dim = len(bounds[0])
    
    def mutate(self, parent, sigma):
        """正态分布变异"""
        child = parent + np.random.normal(0, sigma, self.dim)
        return np.clip(child, self.bounds[0], self.bounds[1])
    
    def run(self):
        # 初始化
        parent = np.random.uniform(self.bounds[0], self.bounds[1])
        parent_fitness = self.objective(parent)
        
        # 1/5成功规则参数
        sigma = 0.1
        success_count = 0
        target_success_rate = 0.2
        
        for t in range(self.n_iterations):
            # 生成子代
            child = self.mutate(parent, sigma)
            child_fitness = self.objective(child)
            
            if child_fitness < parent_fitness:  # 最小化问题
                parent = child
                parent_fitness = child_fitness
                success_count += 1
            
            # 1/5成功规则调整步长
            if t % (10 * self.dim) == 0:
                success_rate = success_count / (10 * self.dim)
                if success_rate > target_success_rate:
                    sigma *= 1 / 0.85  # 增加步长
                else:
                    sigma *= 0.85  # 减小步长
                success_count = 0
        
        return parent, parent_fitness

这个代码虽然简单,但展示了ES的核心循环:生成-评估-更新步长-重复。


四、协方差矩阵自适应CMA-ES:为什么需要自适应?

4.1 固定分布的局限性

早期的ES用的是固定均值和固定方差的正态分布。均值决定搜索的中心位置,方差决定搜索的范围。这些参数需要人工设定,但问题来了:

如果我设定的方差不对怎么办?

设大了,搜索范围太大,效率低;设小了,又容易困在局部最优。而且更麻烦的是,对于复杂的非球形目标函数,固定方差的正态分布根本不够用——你想往某个特定方向多探索一些,但固定分布做不到。

4.2 想象你在黑暗中找一个最高点

让我用一个场景来说明CMA-ES的核心思想:

想象你被蒙住眼睛,放在一座连绵起伏的山地上。你的任务是找到海拔最高的点。你只能靠脚踩、手摸来估计坡度,一步步往上爬。

一开始你对地形一无所知,你试探性地走几步,看看哪里海拔升高了。一段时间后,你对周围的地形有了大致的印象:东边好像坡度比较陡,西边比较平缓。于是你调整策略——往东走得步子小一点,往北走得步子大一点。

CMA-ES做的事情本质上就是这样:它通过观察已探索的点,逐步学习目标函数的局部几何结构(用协方差矩阵描述),然后调整采样分布,让搜索更有效率。

4.3 协方差矩阵到底表示什么

在二维情况下,协方差矩阵是一个2x2的对称矩阵:

  • 对角线上的元素 告诉你沿x1和x2方向的方差有多大
  • 非对角线元素 告诉你x1和x2之间的相关性

把这个协方差矩阵做特征值分解,你会得到两个特征向量(定义主轴方向)和两个特征值(定义沿主轴的缩放程度)。所以协方差矩阵本质上描述了一个椭圆形等高线的形状和朝向

如果目标函数的等高线是椭圆形而不是圆形,那么沿着长轴方向探索和沿着短轴方向探索应该用不同的步长。CMA-ES通过自适应学习这个椭圆的形状,把搜索引向正确的方向。

4.4 CMA-ES的两大自适应机制

CMA-ES有两个核心的自适应机制:

第一,协方差矩阵自适应。通过观察种群中的好个体,推断目标函数的局部几何。如果好个体都聚集在某个方向,协方差矩阵就会调整,让下一次采样更倾向于这个方向。

第二,步长自适应。通过跟踪演化路径(evolution path),判断当前的步长是否合适。如果种群一直在往某个方向移动,说明步长可能太小了;如果种群在原地打转,说明步长太大了。

这两个机制配合工作,让CMA-ES能在各种复杂地形上都表现良好。


五、从零实现CMA-ES:Python完整代码

下面是一个从零实现的CMA-ES,虽然不是工业级实现,但足以帮助你理解算法的每个细节。代码里用了大量注释解释每一步在做什么。

import numpy as np
 
class SimpleCMAES:
    """
    简化版CMA-ES实现,保留核心逻辑
    适合学习用,生产环境请用 pycma 库
    """
    
    def __init__(self, objective, dim, x0=None, sigma=0.3, popsize=None):
        self.objective = objective
        self.dim = dim
        self.sigma = sigma  # 全局步长
        
        # 种群规模:维度越大,λ也应该越大
        if popsize is None:
            self.lambda_ = int(4 + 3 * np.log(dim))
        else:
            self.lambda_ = popsize
        
        # 初始搜索中心
        if x0 is None:
            self.xmean = np.random.randn(dim)
        else:
            self.xmean = np.array(x0, dtype=float)
        
        # 协方差矩阵初始化为单位阵
        # 初始时搜索是各向同性的(球形)
        self.C = np.eye(dim)
        
        # 演化路径(用于步长控制和协方差更新)
        self.pc = np.zeros(dim)  # 协方差矩阵的演化路径
        self.ps = np.zeros(dim)  # 步长的演化路径
        
        # 权重设置(半数法则)
        self.weights = np.array([max(self.lambda_/2 - i, 1) 
                                for i in range(self.lambda_//2)])
        self.weights = self.weights / sum(self.weights)
        
        # 有效个体数,用于计算学习率
        self.mu_eff = 1 / sum(self.weights**2)
        
        # 学习率(根据标准公式计算)
        self.cc = 4 / (dim + 4)
        self.c1 = 2 / ((dim + 1.3)**2 + self.mu_eff)
        self.cmu = min(1 - self.c1, 
                      2 * self.mu_eff / ((dim + 2)**2 + 2 * self.mu_eff))
        self.cs = (self.mu_eff + 2) / (dim + self.mu_eff + 5)
        self.damps = 1 + 2 * max(0, np.sqrt((self.mu_eff - 1)/(dim + 1)) - 1) + self.cs
        
        # 记录最优
        self.best_fitness = float('inf')
        self.best_x = self.xmean.copy()
        self.history = []
    
    def sample_population(self):
        """从多元正态分布中采样"""
        # Cholesky分解:确保采样的协方差正确
        # C = B @ B.T,其中B是下三角矩阵
        B = np.linalg.cholesky(self.C)
        
        # 采样λ个个体
        self.population = np.zeros((self.lambda_, self.dim))
        for i in range(self.lambda_):
            # z是标准正态分布的样本
            z = np.random.randn(self.dim)
            # x = mean + sigma * B @ z
            # 注意这里B @ z 等价于从N(0, C)中采样
            self.population[i] = self.xmean + self.sigma * B @ z
        
        # 评估所有个体
        self.fitness = np.array([self.objective(x) for x in self.population])
        
        # 记录全局最优
        best_idx = np.argmin(self.fitness)
        if self.fitness[best_idx] < self.best_fitness:
            self.best_fitness = self.fitness[best_idx]
            self.best_x = self.population[best_idx].copy()
    
    def update_distribution(self):
        """更新均值、协方差矩阵和步长"""
        # 按适应度排序
        sorted_idx = np.argsort(self.fitness)
        sorted_pop = self.population[sorted_idx]
        
        # 1. 更新均值
        old_mean = self.xmean.copy()
        self.xmean = sum(w * sorted_pop[i] for i, w in enumerate(self.weights))
        
        # 2. 更新步长(累积步长调整)
        # 演化路径ps记录了均值的累计移动方向
        y = (self.xmean - old_mean) / self.sigma
        self.ps = (1 - self.cs) * self.ps + np.sqrt(self.cs * (2 - self.cs)) * y
        
        # 如果ps的范数接近期望值,说明步长合适
        # 如果ps的范数大于期望值,说明种群一直在移动,步长太慢
        # 如果ps的范数小于期望值,说明种群在原地打转,步长太快
        expected_norm = np.sqrt(self.dim) * (1 - 1/(4*self.dim) + 1/(21*self.dim**2))
        self.sigma *= np.exp((self.cs / self.damps) * (np.linalg.norm(self.ps) / expected_norm - 1))
        
        # 3. 更新协方差矩阵的演化路径pc
        self.pc = (1 - self.cc) * self.pc + \
                  np.sqrt(self.cc * (2 - self.cc)) * (self.xmean - old_mean) / self.sigma
        
        # 4. 更新协方差矩阵
        # rank-1更新:用演化路径pc的信息
        rank1 = self.c1 * np.outer(self.pc, self.pc)
        
        # rank-μ更新:用好个体的信息
        # 对于每个好个体,计算它相对于旧均值的偏移
        rank_mu = np.zeros((self.dim, self.dim))
        for i, w in enumerate(self.weights):
            if w > 0:
                y = (sorted_pop[i] - old_mean) / self.sigma
                rank_mu += w * np.outer(y, y)
        
        # 组合两种更新
        self.C = (1 - self.c1 - self.cmu) * self.C + rank1 + rank_mu
        
        # 保持对称性(数值误差可能破坏对称性)
        self.C = (self.C + self.C.T) / 2
    
    def run(self, max_iterations=1000, tol=1e-10, verbose=True):
        """运行CMA-ES"""
        for iteration in range(max_iterations):
            # 采样和评估
            self.sample_population()
            
            # 记录历史
            self.history.append(self.best_fitness)
            
            # 更新分布
            self.update_distribution()
            
            # 输出进度
            if verbose and iteration % 50 == 0:
                print(f"迭代 {iteration:4d}: "
                      f"最优适应度 = {self.best_fitness:.6f}, "
                      f"步长σ = {self.sigma:.6f}")
            
            # 早停检查
            if self.sigma < tol:
                if verbose:
                    print(f"步长收敛,提前停止")
                break
        
        return self.best_x, self.best_fitness
 
 
# 使用示例
if __name__ == "__main__":
    # 找一个球的最低点
    def sphere(x):
        return np.sum(x**2)
    
    dim = 10
    cma = SimpleCMAES(
        objective=sphere,
        dim=dim,
        x0=np.random.randn(dim) * 5,  # 从远处开始
        sigma=2.0,
        popsize=20
    )
    
    best_x, best_f = cma.run(max_iterations=500)
    print(f"\n最终结果: f(x) = {best_f:.10f}")
    print(f"最优解 x = {best_x}")

这段代码实现了一个功能完整的CMA-ES,核心包括:

  1. 多元正态分布采样:通过Cholesky分解确保采样的协方差正确
  2. 均值更新:用加权平均,权重根据排名分配
  3. 步长自适应:通过演化路径的范数判断步长是否合适
  4. 协方差自适应:结合rank-1和rank-μ两种更新

六、CMA-ES的Python库:cma库实战使用

6.1 为什么要用现成的库

虽然上面的代码能让你理解CMA-ES的原理,但真要用在生产环境,还是推荐用 Hansen 大神的 cma 库(也叫 pycma)。这个库经过了十几年的打磨,考虑了大量边界情况和数值稳定性问题。

安装很简单:

pip install cma

6.2 最简单的用法

import numpy as np
import cma
 
def rosenbrock(x):
    """Rosenbrock香蕉函数,非凸优化的经典测试函数"""
    return sum(100 * (x[i+1] - x[i]**2)**2 + (x[i] - 1)**2 for i in range(len(x)-1))
 
# 10维Rosenbrock函数,从x0开始优化
x0 = np.zeros(10)  # 或者随机初始化:np.random.randn(10)
sigma0 = 0.5
 
# es是一个优化器对象
es = cma.CMAEvolutionStrategy(x0, sigma0)
 
# 优化循环
while not es.stop():
    # ask(): 生成候选解
    solutions = es.ask()
    
    # 评估这些解
    fitnesses = [rosenbrock(x) for x in solutions]
    
    # tell(): 把评估结果告诉优化器,更新分布
    es.tell(solutions, fitnesses)
    
    # 打印进度
    print(f"代数: {es.countiter}, "
          f"最优: {es.best.fitness:.6f}, "
          f"σ: {es.sigma:.6f}")
 
# 优化结束,获取结果
print(f"找到的最优解: f(x) = {es.best.fitness:.10f}")
print(f"最优x: {es.best.x}")

6.3 带边界约束的用法

很多实际问题是有边界的,比如神经网络的权重通常在[-1, 1]或者[-0.5, 0.5]之间。CMA-ES原生支持边界约束:

import cma
 
def objective(x):
    return sum(x**2)
 
# 定义每个维度的边界
bounds = [(-1, 1)] * 20  # 20维,每维都是[-1, 1]
 
es = cma.CMAEvolutionStrategy(
    x0=np.random.randn(20) * 0.1,
    sigma0=0.3,
    inopts={
        'bounds': bounds,  # 边界约束
        'popsize': 50,     # 种群规模
    }
)
 
while not es.stop():
    solutions = es.ask()
    
    # 方案1:让cma库自动裁剪到边界内
    # 内部会自动做反射或其他处理
    
    fitnesses = [objective(x) for x in solutions]
    es.tell(solutions, fitnesses)
 
print(f"最优适应度: {es.best.fitness}")

6.4 异步并行评估

如果你的目标函数计算很慢(比如要跑一个完整的神经网络训练),可以用异步模式加速:

import cma
import multiprocessing as mp
 
def objective(x):
    # 模拟一个耗时的评估
    import time
    time.sleep(0.1)
    return sum(x**2)
 
# 使用进程池进行并行评估
with mp.Pool(processes=4) as pool:
    es = cma.CMAEvolutionStrategy(np.random.randn(10), 0.5)
    
    while not es.stop():
        solutions = es.ask()
        
        # 并行评估
        fitnesses = pool.map(objective, solutions)
        
        es.tell(solutions, fitnesses)

6.5 常用参数说明

参数含义典型值
x0初始解随机或给定
sigma0初始步长0.3~1.0
popsize种群规模λ4+3*ln(n)
bounds边界约束[[low, high], …]
maxiter最大迭代次数1000+

七、进化策略 vs 梯度下降:什么时候ES更优?

7.1 梯度下降的优势

先说梯度下降的长处:

速度快。如果目标函数光滑、可导、有梯度,梯度下降通常收敛很快。随机梯度下降虽然有噪声,但在大规模问题上依然高效。

理论成熟。凸优化的理论非常完善,你知道什么情况下梯度下降能收敛到全局最优。

资源省。每次迭代只评估一个点(或者小批量),内存和计算都很省。

7.2 梯度下降的短板

但梯度下降也有明显的局限:

需要梯度。如果目标函数不可导,或者梯度难以计算(比如黑盒仿真器),梯度下降就用不了。

容易陷入局部最优。对于非凸问题,梯度下降只能保证找到局部最优。

对噪声敏感。如果梯度估计不准,收敛会出问题。

需要调学习率。手动调学习率是个技术活,虽然有Adam这样的自适应方法,但超参数选择依然麻烦。

7.3 进化策略的优势场景

ES正好能补上梯度下降的短板:

不需要梯度。ES只关心目标函数的输入输出,不管梯度存在与否。这对于仿真优化、强化学习中的策略优化特别有用。

全局搜索能力强。种群搜索机制让ES更容易跳出局部最优。虽然不能保证找到全局最优,但通常能找到还不错的结果。

自适应步长和协方差。不用手动调学习率和二阶信息。

鲁棒。对噪声不敏感,适合评估带噪声的目标函数。

7.4 什么时候选ES

下面这些场景,ES通常是更好的选择:

  1. 黑盒优化问题:目标函数的内部结构未知,只能黑盒调用评估
  2. 非凸、非光滑问题:存在多个局部最优,或者函数有突变
  3. 强化学习策略优化:策略的适应度评估天然带噪声,而且梯度难以获得
  4. 超参数优化:神经网络架构、学习率等超参数之间可能存在复杂交互
  5. 仿真优化:每次评估要跑一个仿真器,评估很慢但可以并行

7.5 什么时候选梯度下降

反过来说,下面这些场景梯度下降更合适:

  1. 凸优化问题:有理论保证全局收敛,没必要用更复杂的ES
  2. 需要快速收敛:对速度要求高,而且梯度可用
  3. 超大规模问题:维度上千上万,ES的种群开销太大
  4. 在线学习:数据流式输入,需要增量更新

7.6 一个直观的对比

特性梯度下降进化策略
梯度需求必须不需要
收敛速度(凸问题)
全局搜索能力
参数敏感性高(学习率)低(自适应)
对噪声的容忍度
每次迭代开销大(种群评估)
理论保证强(凸)

八、ES在超参数优化中的应用

8.1 超参数优化为什么难

训练过神经网络的人都知道,超参数调优是个大坑。学习率、批量大小、网络深度、卷积核大小、Dropout率、L2正则化系数……这些超参数之间可能存在复杂的交互,而且目标函数(验证集性能)对超参数通常是非凸、非光滑的。

更重要的是,每次评估一组超参数都要重新训练模型,成本很高。所以超参数优化是个评估代价昂贵的问题。

8.2 ES天然适合超参数优化

ES特别适合这个场景:

  1. 不需要梯度:验证集性能对超参数的梯度根本没法算
  2. 自适应:不同超参数的重要性不同,ES能自动学习
  3. 对噪声容忍:验证集性能天然有波动,ES不 care

8.3 实战:用ES优化一个简单模型的超参数

import numpy as np
import cma
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
from sklearn.datasets import make_classification
 
# 准备数据
X, y = make_classification(n_samples=1000, n_features=20, random_state=42)
 
def evaluate_hyperparams(params):
    """
    评估一组超参数
    params: [n_estimators, max_depth, min_samples_split, min_samples_leaf, max_features]
    """
    n_estimators = int(params[0])  # 树的数量
    max_depth = int(params[1]) if params[1] > 1 else None  # 最大深度
    min_samples_split = int(params[2])  # 分裂所需最小样本
    min_samples_leaf = int(params[3])  # 叶节点最小样本
    max_features = params[4]  # 分裂特征比例
    
    # 克隆模型
    clf = RandomForestClassifier(
        n_estimators=max(10, n_estimators),
        max_depth=max_depth,
        min_samples_split=max(2, min_samples_split),
        min_samples_leaf=max(1, min_samples_leaf),
        max_features=max_features if max_features > 0 else 'sqrt',
        random_state=42,
        n_jobs=-1
    )
    
    # 交叉验证评估
    scores = cross_val_score(clf, X, y, cv=5, scoring='accuracy')
    return -scores.mean()  # CMA-ES是最小化,所以取负
 
# 设定超参数搜索范围
# [n_estimators, max_depth, min_samples_split, min_samples_leaf, max_features]
lower_bounds = [10, 1, 2, 1, 0.1]
upper_bounds = [200, 20, 20, 10, 1.0]
bounds = list(zip(lower_bounds, upper_bounds))
 
# 初始化ES
x0 = [50, 10, 5, 2, 0.5]  # 初始猜测
es = cma.CMAEvolutionStrategy(x0, 0.5, inopts={'bounds': bounds, 'popsize': 20})
 
# 优化
n_evaluations = 0
while not es.stop() and n_evaluations < 100:
    solutions = es.ask()
    fitnesses = [evaluate_hyperparams(x) for x in solutions]
    es.tell(solutions, fitnesses)
    n_evaluations += len(solutions)
    
    print(f"评估 {n_evaluations}: 最优准确率 {-es.best.fitness:.4f}")
 
print(f"\n最优超参数: {es.best.x}")
print(f"最优准确率: {-es.best.fitness:.4f}")

这个例子展示了用CMA-ES优化随机森林超参数的全过程。虽然代码写起来比网格搜索长,但搜索效率高得多——ES能用更少的评估次数找到更好的超参数组合。

8.4 和贝叶斯优化的对比

超参数优化另一个常用方法是贝叶斯优化(用高斯过程或Tree-Parzen Estimator建模)。两者各有优劣:

  • 贝叶斯优化:利用历史评估结果建模,适合评估代价很高、迭代次数有限的情况
  • ES:并行性好,不容易陷入局部最优,适合评估代价可接受、需要全局搜索的情况

实践中可以结合两种方法:前期用ES快速探索,后期用贝叶斯精细调优。


九、实战:用ES优化神经网络权重

9.1 这个思路可行吗

理论上,用进化策略优化神经网络权重是完全可行的——神经网络权重就是一堆实数,ES能在实数空间里搜索。但实际上存在几个问题:

  1. 维度爆炸:现代神经网络动不动就几百万、几千万个参数,用ES的话每次评估都要前向传播整个种群,太慢了
  2. 评估代价高:每次评估要跑一遍前向传播(可能还要反向传播算其他东西)
  3. 存储开销大:要维护一个种群,每个个体都是完整的网络权重

不过在某些场景下这个思路是有价值的,比如:

  • 小规模网络:几KB~几MB的网络,可以接受ES的开销
  • 评估可并行:评估代价虽然高,但可以并行化
  • 需要跳出局部最优:对于某些强化学习任务,梯度方法容易陷入局部最优

9.2 一个简化示例:优化MLP权重

import numpy as np
import cma
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
 
# 生成测试数据
X, y = make_moons(n_samples=500, noise=0.2, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
 
# 简单的MLP实现
class SimpleMLP:
    def __init__(self, layer_sizes):
        self.weights = []
        self.biases = []
        for i in range(len(layer_sizes) - 1):
            # Xavier初始化
            w = np.random.randn(layer_sizes[i], layer_sizes[i+1]) * np.sqrt(2.0 / layer_sizes[i])
            b = np.zeros(layer_sizes[i+1])
            self.weights.append(w)
            self.biases.append(b)
    
    def set_weights_from_vector(self, vec):
        """从一维向量设置权重"""
        idx = 0
        for i in range(len(self.weights)):
            w_size = self.weights[i].size
            b_size = self.biases[i].size
            self.weights[i] = vec[idx:idx+w_size].reshape(self.weights[i].shape)
            self.biases[i] = vec[idx+w_size:idx+w_size+b_size]
            idx += w_size + b_size
    
    def predict_proba(self, X):
        h = X
        for i in range(len(self.weights) - 1):
            h = np.tanh(h @ self.weights[i] + self.biases[i])
        return h @ self.weights[-1] + self.biases[-1]  # 输出层无激活
    
    def accuracy(self, X, y):
        preds = (self.predict_proba(X) > 0).astype(int).ravel()
        return np.mean(preds == y)
 
# 网络结构
layer_sizes = [2, 16, 8, 1]
mlp = SimpleMLP(layer_sizes)
 
# 计算总参数数量
total_params = sum(w.size + b.size for w, b in zip(mlp.weights, mlp.biases))
print(f"总参数数量: {total_params}")
 
# 目标函数
def fitness(vector):
    mlp.set_weights_from_vector(vector)
    return -mlp.accuracy(X_train, y_train)  # 负准确率(最小化)
 
# 初始化
x0 = np.random.randn(total_params) * 0.1
 
# 用CMA-ES优化
es = cma.CMAEvolutionStrategy(x0, 0.1, inopts={'popsize': 50, 'maxiter': 200})
 
# 优化循环
for generation in range(200):
    solutions = es.ask()
    fitnesses = [fitness(x) for x in solutions]
    es.tell(solutions, fitnesses)
    
    if generation % 20 == 0:
        mlp.set_weights_from_vector(es.best.x)
        train_acc = mlp.accuracy(X_train, y_train)
        test_acc = mlp.accuracy(X_test, y_test)
        print(f"代 {generation:3d}: 训练准确率 {train_acc:.3f}, 测试准确率 {test_acc:.3f}")
 
# 最终结果
mlp.set_weights_from_vector(es.best.x)
print(f"\n最终训练准确率: {mlp.accuracy(X_train, y_train):.3f}")
print(f"最终测试准确率: {mlp.accuracy(X_test, y_test):.3f}")

这个例子用CMA-ES优化了一个小MLP的权重,展示了基本流程。虽然在小网络上效果还可以,但扩展到大网络就力不从心了。

9.3 神经进化领域的实际应用

在学术界和工业界,进化策略在神经进化领域有以下实际应用:

  1. NEAT(NeuroEvolution of Augmenting Topologies):用进化算法同时优化网络结构和权重,是神经进化的经典方法
  2. OpenAI的ES:Salimans等人在2017年用大规模ES在MuJoCo物理模拟任务上取得了和深度强化学习相当的结果
  3. CMA-NEAT:CMA-ES和NEAT的结合,用于演化更复杂的网络结构

这些方法的核心思想是:用进化策略优化权重或结构,而不是依赖梯度。


十、进化策略的前沿:OpenAI的ES与神经架构搜索

10.1 2017年的惊人结果

2017年,OpenAI的Salimans等人发表了一篇很有意思的论文”Evolution Strategies as a Scalable Alternative to Reinforcement Learning”。他们用进化策略(更准确地说是高斯扰动方法)在一个大规模的强化学习任务上取得了和PPO、SAC等主流算法相当的结果。

这在当时挺让人意外的——大家以为进化策略这种”老古董”方法在深度学习时代早就过时了,结果它居然还能打。

10.2 OpenAI ES的核心思想

OpenAI的ES并不是严格意义上的种群进化,而是采用了一种更简单的策略:对当前最优解加随机噪声来探索。

具体来说:

  1. 维护一个参数向量 θ(就是神经网络的权重)
  2. 生成N个噪声样本:θ + σ·ε,其中ε是从标准正态分布采样的
  3. 评估这N个扰动样本的适应度
  4. 用加权平均更新θ

这本质上是一种”无梯度”的优化方法,特别适合强化学习场景,因为强化学习的梯度本身就有高方差、难以估计的问题。

10.3 ES的优势:并行化简单

OpenAI的ES论文特别强调了并行化的优势。RL算法通常需要GPU来加速前向传播,而ES的每个噪声样本可以独立评估,所以可以用大量的CPU核心并行处理。

论文中,他们用1440个CPU核心进行并行评估,达到了和单GPU的PPO相当的样本效率。

10.4 神经架构搜索(NAS)

神经架构搜索是另一个ES大放异彩的领域。目标是自动搜索最优的网络结构(层的类型、连接方式、超参数等)。

由于架构空间通常是离散的、不可导的,梯度方法不好用。而ES——特别是结合了CMA-ES的方法——能够有效地探索这个空间。

一些经典的神经进化方法:

  • NEAT:用进化算法演化神经网络的拓扑结构和权重
  • AmoebaNet:Google用进化算法搜索到的网络架构,在ImageNet上取得了当时最好的准确率
  • LEO (Learning Evolution Embeddings):结合了贝叶斯优化和进化的方法

10.5 未来的发展方向

进化策略领域还有一些活跃的研究方向:

  1. 大规模并行ES:继续探索如何利用更多CPU/GPU资源
  2. 混合方法:把ES和梯度方法结合,比如用ES探索局部最优区域,再用梯度精细调优
  3. 自适应噪声:自动调整探索的程度,适应不同的问题
  4. 分布式ES:多节点、多机房的并行ES

十一、进阶主题:自然进化策略

11.1 什么是自然进化策略

自然进化策略(NES)是Wierstra等人在2014年提出的理论框架,它把进化策略和统计学中的自然梯度联系了起来。

在普通的梯度下降里,你沿着最陡的方向更新参数。但在参数空间里,“最陡”的方向和参数的度量有关——如果参数缩放不同,同样的方向实际上会有不同的步长效果。

自然梯度解决这个问题:它定义了一种”与参数无关”的梯度方向,基于Fisher信息矩阵对梯度进行重归一化。

11.2 从信息几何角度看ES

NES的核心洞察是:进化策略本质上是在优化一个概率分布的参数,让这个分布产生高适应度的样本。

这个概率分布空间有一个几何结构——它是一个黎曼流形,曲率由Fisher信息矩阵定义。自然梯度就是在黎曼流形上最陡下降的方向。

这个视角很优雅,它把”进化”和”统计推断”联系了起来。

11.3 xNES算法

xNES(指数自然进化策略)是一个实用的NES变体。它的特点是:

  1. 用指数映射来确保协方差矩阵始终正定
  2. 参数更新在切空间进行,然后映射回流形
  3. 实现相对简单,效果不错
class ExponentialNES:
    def __init__(self, dim, lr=0.01):
        self.dim = dim
        self.lr = lr
        self.mean = np.zeros(dim)
        self.A = np.eye(dim)  # 协方差矩阵的Cholesky因子
    
    def sample(self, n):
        """采样"""
        z = np.random.randn(n, self.dim)
        return self.mean + z @ self.A.T
    
    def fit(self, samples, fitnesses):
        """用适应度更新分布"""
        # 排名转换为utility
        ranks = np.argsort(np.argsort(fitnesses))
        utilities = np.array([max(0, np.log(n/2+1) - np.log(n-r)) 
                             for r in ranks])
        utilities /= sum(utilities)
        utilities -= 1.0 / len(utilities)  # 中心化
        
        # 计算梯度
        z = (samples - self.mean) @ np.linalg.inv(self.A.T)
        grad_mean = z.T @ utilities
        grad_A = utilities @ (np.eye(len(utilities)) - utilities[0] * z @ z.T)
        
        # 指数更新
        self.mean += self.lr * self.A @ grad_mean
        self.A = A @ expm(self.lr * grad_A)
    
    def run(self, objective, n_iter=1000):
        for _ in range(n_iter):
            samples = self.sample(100)
            fitnesses = [objective(x) for x in samples]
            self.fit(samples, fitnesses)

11.4 NES vs CMA-ES

NES和CMA-ES有很多相似之处,也有关键的区别:

方面NESCMA-ES
理论基础信息几何、自然梯度进化算法理论
更新方式梯度上升utility的期望协方差矩阵自适应
参数化协方差矩阵或Cholesky协方差矩阵
权重分配适应度排名→utility固定权重或排名权重
实践表现相当相当

实际上,两种方法在实践中表现相近,选择哪个更多是个人偏好。


十二、实战技巧与常见问题

12.1 如何设置初始步长σ

初始步长是CMA-ES最重要的超参数。经验法则:

  • 如果你知道最优解在x0附近,σ可以小一点(比如0.1)
  • 如果你不确定解在哪里,σ可以大一点(比如1.0)
  • 一个参考:把σ设为问题规模(参数空间的典型尺度)的1/10到1/3

如果σ设错了会怎样?设大了,算法会在初期浪费很多评估;设小了,算法可能被困在局部最优,或者演化路径不稳定。

12.2 如何设置种群规模λ

种群规模的典型值是:

其中n是维度。这个公式的经验来源是:在维度为n的问题上,λ应该足够大,以便能在n+1个点上拟合一个二次函数。

λ越大,探索能力越强,但每次迭代的计算开销也越大。对于困难问题,可以适当增大λ;对于简单问题,λ可以小一点。

12.3 算法不收敛怎么办

如果CMA-ES表现不好,可以检查以下几点:

  1. σ是不是太大或太小:观察步长的变化,如果一直在增大或减小,说明σ的初始值可能不对
  2. λ是不是太小:种群太小可能导致多样性不足
  3. 评估是否有噪声:如果目标函数有随机性,增加λ或重复评估可以降低噪声
  4. 是否陷入局部最优:可以尝试重启算法,从不同的初始点开始

12.4 边界约束的处理

CMA-ES原生不支持边界约束(虽然可以扩展)。常见的处理方法:

  1. 拒绝采样:采样到的解如果超出边界就拒绝,重新采样
  2. 反射:超出边界的部分反射回来
  3. 裁剪:直接裁剪到边界内(不推荐,会扭曲分布)
  4. 变量变换:对有界变量做对数变换,比如用log(x - low)和log(high - x)作为搜索空间

12.5 早停条件

CMA-ES有多种内置的早停条件:

  • tolx:步长变化太小
  • tolfun:适应度变化太小
  • tolstagnation:多代适应度没有提升
  • maxiter:达到最大迭代次数

可以根据需要设置这些阈值。


参考文献

  1. Rechenberg, I. (1973). Evolutionsstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution. Frommann-Holzberg.
  2. Schwefel, H.-P. (1981). Numerical Optimization of Computer Models. John Wiley & Sons.
  3. Hansen, N., & Ostermeier, A. (2001). Completely Derandomized Self-Adaptation in Evolution Strategies. Evolutionary Computation, 9(2), 159-195.
  4. Hansen, N., et al. (2015). CMA-ES: Evolution Strategies with Covariance Matrix Adaptation for Objective Functions with Time-dependency. Evolutionary Computation, 23(4), 643-669.
  5. Wierstra, D., et al. (2014). Natural Evolution Strategies. Journal of Machine Learning Research, 15, 949-980.
  6. Akimoto, Y., & Hansen, N. (2016). Block-Boundary Covariance Matrix Adaptation. Proceedings of GECCO, 2016.
  7. Salimans, T., et al. (2017). Evolution Strategies as a Scalable Alternative to Reinforcement Learning. arXiv preprint arXiv:1703.03864.

相关文档