（CMA-ES源码）协方差自适应进化策略（Covariance Matrix Adaptation Evolution Strategy，CMA-ES）——最好的单目标进化算法？

获取更多资讯，赶快关注公众号(名称：智能制造与智能调度，公众号：deeprlscheduler)吧！

关于此算法，我有点几点要说明：

1. 我没有花太多时间去整理这些公式的推导过程，这个还需要自己对数学基础有个基本的了解；
2. 源代码可关注公众号后，回复”CMA-ES“或”自适应协方差“获取链接；
3. 只要掌握几个关键的公式，然后再对照源码学习，就很快能理解算法了。

缩写

CMA：Covariance Matrix Adaptation，协方差自适应

EMNA：Estimation of Multivariate Normal Algorithm，多元正态估计算法

ES：Evolution Strategy，进化策略

( μ / μ { I , W } , λ ) − E S \left(\mu / \mu_{\{\mathrm{I}, \mathrm{W}\}}, \lambda\right)-\mathrm{ES} (μ/μ{I,W}​,λ)−ES：进化策略中存在$\mu$个父代，$\lambda$个子代，且所有$\mu$个父代都重组。

希腊符号

$\lambda \ge 2$：种群大小，样本大小，子代数量

$\mu \le \lambda$：父代数量，在种群中选择的搜索点的数量

${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}={\left({\sum }_{i=1}^{\mu }{w}_i^2\right)}^{-1}=1/\Vert \mathbf{w}{\Vert }^2$：the variance effective selection mass，方差有效选择质量

${\sigma }^{(g)}\in {\mathbb{R}}_{+}$：step-size，步长

拉丁符号

$\mathbf{B}\in {\mathbb{R}}^n$：正交矩阵，$\mathbf{B}$的列为$\mathbf{C}$的单位长度特征向量，对应于$\mathbf{D}$的对角元素

${\mathbf{C}}^{(g)}\in {\mathbb{R}}^{n\times n}$：第$g$代的协方差矩阵

$c_{ii}$：$\mathbf{C}$的对象元素

$c_{\mathrm{c}}\le 1$：协方差矩阵的秩一更新的累积学习率

$c_1\le 1-c_{\mu }$：协方差矩阵秩一更新的学习率

$c_{\mu }\le 1-c_1$：协方差矩阵秩$\mu$更新的学习率

$c_{\sigma }<1$：步长控制的累积学习率

$D\in {\mathbb{R}}^n$：对角矩阵，$D$的对焦元素为$C$的特征值的平方根，对应于$B$的各列

$d_i>0$：对角矩阵$D$的对角元素，$d_i^2$为$C$的特征值

$d_{\sigma }\approx 1$：步长更新的衰减参数

$E$：期望值

$f:{\mathbb{R}}^n\to \mathbb{R},\mathbf{x}\mapsto f(\mathbf{x})$：待最小化的目标函数（适应度函数）

$f_{\text{sphere }}:{\mathbb{R}}^n\to \mathbb{R},\mathbf{x}\mapsto \Vert \mathbf{x}{\Vert }^2={\mathbf{x}}^{\mathrm{T}}\mathbf{x}={\sum }_{i=1}^n{x}_i^2$

$g\in {\mathbb{N}}_0$：代数计数器，迭代次数

$\mathbf{I}\in {\mathbb{R}}^{n\times n}$：单位矩阵

${\mathbf{m}}^{(g)}\in {\mathbb{R}}^n$：第$g$代搜索分布的均值

$n\in \mathbb{N}$:搜索空间维度

$\mathcal{N}(0,\mathbf{I})$：0均值和单位协方差矩阵的多元正太分布，按照该分布生成的向量，其元素是独立且服从$\mathcal{N}(0,1)$标准正太分布

$\mathcal{N}(\mathbf{m},\mathbf{C})\sim \mathbf{m}+\mathcal{N}(0,\mathbf{C})$：均值为$\mathbf{m}\in {\mathbb{R}}^n$，协方差矩阵为$\mathbf{C}\in {\mathbb{R}}^{n\times n}$的多元正太分布，$\mathbf{C}$为对称正定矩阵

$\mathbf{p}\in {\mathbb{R}}^n$：进化路径，即策略经过多代所采取的一系列连续步

$w_i,$ where $i=1,\dots ,\mu$：重组权重

${\mathbf{x}}_k^{(g+1)}\in {\mathbb{R}}^n$：第$g+1$代中的第$k$个子代或个体，也将${\mathbf{x}}^{(g+1)}$称为搜索点或对象参数/变量，与候选解或设计变量同义

${\mathbf{x}}_{i:\lambda }^{(g+1)}$： ${\mathbf{x}}_1^{(g+1)},\dots ,{\mathbf{x}}_{\lambda }^{(g+1)}$中第$i$个最优个体，$i:\lambda$表示排名第$i$的个体索引，且$f\left({\mathbf{x}}_{1:\lambda }^{(g+1)}\right)\le f\left({\mathbf{x}}_{2:\lambda }^{(g+1)}\right)\le \cdots \le f\left({\mathbf{x}}_{\lambda :\lambda }^{(g+1)}\right)$

${\mathbf{y}}_k^{(g+1)}=\left({\mathbf{x}}_k^{(g+1)}-{\mathbf{m}}^{(g)}\right)/{\sigma }^{(g)}$：相当于${\mathbf{x}}_k=\mathbf{m}+\sigma {\mathbf{y}}_k$

CMA-ES是一种用于求解非线性非凸函数实值参数(连续域)优化的随机方法，相信很多读者在看了原作者的论文后会感到无助和绝望，有那么多公式，那么多符号，不知如何下手，不用慌，今天这篇文章将给你总结出最主要的核心点和算法思想，直接跳过前面那些让人头疼的数学基础，基础给出算法，然后再逐步分析算法中关键步骤。

0 $\left(\mu /{\mu }_{\mathrm{W}},\lambda \right)$-CMA-ES

其实对进化算法ES稍微有一些了解的童鞋，应该知道ES的大概思路，之前这篇文章（进化计算-进化策略（Evolutionary Strategies,ES）前世今生与代码共享）我们也有所介绍，那么本文介绍的CMA-ES其实是ES的自适应版本，它会记录一定迭代次数的种群历史，计算目标变量之间的协方差和方差信息，随后的搜索工作受到这些方差值的影响。只要记住一点，不论基因发生何种变化，产生的结果（性状）总遵循这零均值，某一方差的高斯分布。下面给出了CMA-ES算法的精简伪代码，是不是看起来没有那么复杂呢？

参数设置，包括$\lambda ,\mu ,w_{i=1\dots \mu },c_{\sigma },d_{\sigma },c_c,c_1,$ and $c_{\mu }$，这些参数含义见最开始的符号描述。

初始化，设置进化路径，协方差矩阵，当前代数，选择分布平均值和步长。

循环直到达到终止条件，从搜索点中采样新的种群，然后进行选择和重组，控制步长大小和协方差自适应。

那么接下来就分别看看这三步中具体是怎么实现的。

1 Set parameters参数设置

这一步中设置的参数主要是指策略参数，包括子代数量$\lambda$,父代数量$\mu$,重组权重$w_{i=1\dots \mu }$,步长控制的累积学习率$c_{\sigma }$,步长更新的衰减参数$d_{\sigma }$,协方差矩阵的秩一更新的累积学习率$c_{\mathrm{c}}$,协方差矩阵秩一更新的学习率$c_{1} $,协方差矩阵秩$\mu$更新的学习率$c_{\mu}$。这些参数具体干什么用这里暂时不介绍，后面用到的时候自然会解释。

这些参数往往都有自己的默认值，文中也给出了它们的计算公式，其中${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{1}{{\sum }_{i=1}^{\mu }{w}_i^2}\ge 1$ and ${\sum }_{i=1}^{\mu }{w}_i=1$。

选择和重组：

$$\lambda =4+\lfloor 3\ln n\rfloor ,\mu =\lfloor {\mu }^{\prime }\rfloor ,{\mu }^{\prime }=\frac{\lambda }{2}\text{\hspace{1em}(44)}$$

其中$\lfloor \rfloor$表示floor操作，即不大于某值得最大整数。由于$n$为搜索空间维度，是由问题决定的且已知，所以$\lambda$可求，然后就可求得${\mu }^{\prime }$，进而求得$\mu$。

$$w_i=\frac{w_i^{\prime }}{{\sum }_{j=1}^{\mu }{w}_j^{\prime }},w_i^{\prime }=\ln \left({\mu }^{\prime }+0.5\right)-\ln i\text{ for }i=1,\dots ,\mu \text{\hspace{1em}(45)}$$

有了${\mu }^{\prime }$和$\mu$就可求得$w_j^{\prime }$,进而求得$w_i$。

步长控制：

$$c_{\sigma }=\frac{{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}+2}{n+{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}+5},d_{\sigma }=1+2\max \left(0,\sqrt{\frac{{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}-1}{n+1}}-1\right)+c_{\sigma }\text{\hspace{1em}(46)}$$

有了$w_i$后，就可以计算出$\mu_{\mathrm{eff}}=\frac{1}{\sum_{i=1}^{\mu} w_{i}^{2}} $，进而得到$c_{\sigma}$和$d_{\sigma}$。

协方差自适应：

$$c_{\mathrm{c}}=\frac{4+{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}/n}{n+4+2{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}/n}\text{\hspace{1em}(47)}$$

$$c_1=\frac{2}{(n+1.3{)}^2+{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}}\text{\hspace{1em}(48)}$$

$$c_{\mu }=\min \left(1-c_1,{\alpha }_{\mu }\frac{{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}-2+1/{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}}{(n+2{)}^2+{\alpha }_{\mu }{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}/2}\right)\text{ with }{\alpha }_{\mu }=2\text{\hspace{1em}(49)}$$

根据${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}$和$n$可以计算出$c_{\mathrm{c}}$,$c_1$和$c_{\mu }$。

说明：

式(46)~(49)所计算的默认参数通常是比较鲁棒的，因此根据经验，可以广泛地适用于不同的优化问题，这些参数通常不必进行修改。而子代数量$\lambda$,父代数量$\mu$,重组权重$w_{i=1\dots \mu }$相对而言不是很关键，可以在较大范围内进行选择而不影响自适应，一般推荐选择$\mu \le \lambda /2$，以及根据$w_1\ge \cdots \ge w_{\mu }$选择重组权重，种群大小$\lambda$可以适当增大，如果$\mu$和$w_i$使用了与$\lambda$有关的默认值，种群大小$\lambda$会对全局搜索性能有显著影响，增加$\lambda$通常会提高CMA-ES的全局搜索能力和鲁棒性，但会降低其收敛速度，收敛速度最多随$\lambda$线性下降。

2 Initialization初始化

设置进化路径$p_{\sigma }=0,p_c=0$，相当于清除了路径，设置协方差矩阵$C=\mathbf{I}$，当前迭代次数$g=0$。

选择分布均值$\mathbf{m}\in {\mathbb{R}}^n$和步长大小$\sigma \in {\mathbb{R}}_{+}$，这两个值其实是和问题相关的，最优解大致是应该位于$\mathbf{m}\pm 3\sigma (1,\dots ,1{)}^{\mathrm{T}}$内的，如果最优解期望位于初始搜索间隔$[a,b{]}^n$，我们可以在$[a,b{]}^n$内均匀随机选择初始搜索点$\mathbf{m}$，且$\sigma =0.3(b-a)$。针对不同变量的不同搜索间隔$\Delta s_i$可以通过$C$的不同初始化来体现，其中$C$的对角元素满足$c_{ii}={\left(\Delta s_i\right)}^2$。注意$\Delta s_i$不能在数量级上差距很大，否则就需要进行缩放。

3 迭代优化过程

3.1 采样新种群

在CMA-ES中，新搜索点（个体，子代）的种群是通过多元正态分布采样（回顾：给定方差或协方差，正太分布的交叉熵在所有分布中是最大的）生成的，在第$g$代采样搜索点的公式如下：
$${\mathbf{x}}_k^{(g+1)}\sim {\mathbf{m}}^{(g)}+{\sigma }^{(g)}\mathcal{N}\left(0,{\mathbf{C}}^{(g)}\right)\text{ for }k=1,\dots ,\lambda \text{\hspace{1em}(5)}$$

其中：

$\sim$表示左右分布相同。

$\mathcal{N}\left(0,{\mathbf{C}}^{(g)}\right)$：是均值为0，协方差为${\mathbf{C}}^{(g)}$的多元正太分布，其满足${\mathbf{m}}^{(g)}+{\sigma }^{(g)}\mathcal{N}\left(0,{\mathbf{C}}^{(g)}\right)\sim \mathcal{N}\left({\mathbf{m}}^{(g)},{\left({\sigma }^{(g)}\right)}^2{\mathbf{C}}^{(g)}\right)$。

${\mathbf{x}}_k^{(g+1)}\in {\mathbb{R}}^n$：第$g+1$代中的第$k$个子代（个体，搜索点）。

${\mathbf{m}}^{(g)}\in {\mathbb{R}}^n$：第$g$代搜索分布的均值。

${\sigma }^{(g)}\in {\mathbb{R}}_{+}$：第$g$代整体标准差，步长。

$C^{(g)}\in {\mathbb{R}}^{n\times n}$：第$g$代的协方差矩阵

$\lambda \ge 2$：种群大小，采样大小，子代数量。

为了定义完整的迭代步骤，剩下的问题就是如何计算第$g+1$代的${\mathbf{m}}^{(g+1)},$ ${\mathbf{C}}^{(g+1)}$ 和 ${\sigma }^{(g+1)}$

3.2 选择和重组：均值移动

搜索分布新的均值${\mathbf{m}}^{(g+1)}$就是从样本${\mathbf{x}}_1^{(g+1)},\dots ,{\mathbf{x}}_{\lambda }^{(g+1)}$中选择的$\mu$个点的加权平均值：

$${\mathbf{m}}^{(g+1)}=\sum _{i=1}^{\mu }{w}_i{\mathbf{x}}_{i:\lambda }^{(g+1)}\text{\hspace{1em}(6)}$$

$$\sum _{i=1}^{\mu }{w}_i=1,w_1\ge w_2\ge \cdots \ge w_{\mu }>0\text{\hspace{1em}(7)}$$

其中

• $\mu \le \lambda$为父代种群大小，即选择点的个数。
• $w_{i=1\dots \mu }\in {\mathbb{R}}_{+}$：用于重组的正权重系数。如果$w_{i=1\dots \mu }=1/\mu$，等式(6)相当于计算了$\mu$个选择点的均值。
• ${\mathbf{x}}_{i:\lambda }^{(g+1)}$：根据式(5)计算得到的${\mathbf{x}}_1^{(g+1)},\dots ,{\mathbf{x}}_{\lambda }^{(g+1)}$中第$i$个最优个体，$i:\lambda$表示排名第$i$的个体索引，且$f\left({\mathbf{x}}_{1:\lambda }^{(g+1)}\right)\le f\left({\mathbf{x}}_{2:\lambda }^{(g+1)}\right)\le \cdots \le f\left({\mathbf{x}}_{\lambda :\lambda }^{(g+1)}\right)$。

等式(6)通过从$\lambda$个子代中选择$\mu <\lambda$个实现截断选择(truncation selection),其实分配不同的权重$w_i$也相当于一种选择机制。等式(6)实现了加权中间重组(weighted intermediate recombination)，考虑了$\mu >1$个个体以获得一个加权平均值。

$${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}={\left(\frac{\Vert \mathbf{w}{\Vert }_1}{\Vert \mathbf{w}{\Vert }_2}\right)}^2=\frac{\Vert \mathbf{w}{\Vert }_1^2}{\Vert \mathbf{w}{\Vert }_2^2}=\frac{1}{\Vert \mathbf{w}{\Vert }_2^2}={\left(\sum _{i=1}^{\mu }{w}_i^2\right)}^{-1}\text{\hspace{1em}(8)}$$

上式会在后面重复使用到，可以解释为方差有效选择质量( variance effective selection mass)。根据(7)中$w_i$的定义，可以推导出$1\le {\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\le \mu$，且当重组权重相等即对于所有$i=1\dots \mu$都有$w_i=1/\mu$，则${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}=\mu$。通常，${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\approx \lambda /4$表明$w_i$的设置比较合理，典型的设置是$w_i\propto \mu -i+1,$ and $\mu \approx \lambda /2$。

3.3 自适应协方差矩阵

这一节将推导协方差矩阵的更新过程。将从某一代的单个总体开始估计协方差矩阵(3.3.1)，对于小种群，这种估计是不可靠的，必须提出一种自适应方法(3.3.2的秩$\mu$更新)，在极限情况下，每一代只能使用一个单点来更新(调整)协方差矩阵(3.3.3的秩1更新)，适应性可以通过利用连续步骤之间的累积依赖性来增强(3.3.3.1)，最后结合秩$\mu$和秩1更新方法(3.3.4)。

3.3.1 估计协方差矩阵

现在假设种群包含足够的信息，可以可靠地估计协方差矩阵。为了方便起见，在本节中假设${\sigma }^{(g)}=1$，对于${\sigma }^{(g)}\ne 1$时，除了常数因子外，这个公式是仍然成立。

为了根据$\mathcal{N}(0,\mathbf{I})$分布的样本重新估计协方差矩阵且保证$\mathbf{C}$的条件数$\mathrm{cond}(\mathbf{C})<10$，应该保证采样大小$\lambda \ge 4n$。

我们可以使用根据(5)采样的种群重新评估原始的协方差矩阵${\mathbf{C}}^{(g)}$，通过下面的经验协方差矩阵：
$${\mathbf{C}}_{\mathrm{e}\mathrm{m}\mathrm{p}}^{(g+1)}=\frac{1}{\lambda -1}\sum _{i=1}^{\lambda }\left({\mathbf{x}}_i^{(g+1)}-\frac{1}{\lambda }\sum _{j=1}^{\lambda }{\mathbf{x}}_j^{(g+1)}\right){\left({\mathbf{x}}_i^{(g+1)}-\frac{1}{\lambda }\sum _{j=1}^{\lambda }{\mathbf{x}}_j^{(g+1)}\right)}^{\mathrm{T}}\text{\hspace{1em}(9)}$$

经验协方差矩阵${\mathbf{C}}_{\mathrm{e}\mathrm{m}\mathrm{p}}^{(g+1)}$是${\mathbf{C}}^{(g)}$的无偏估计：假设${\mathbf{x}}_i^{(g+1)},i=1\dots \lambda$为随机变量(而不是已完成的采样)，那么就有$\mathrm{E}\left[\begin{array}{ll}{\mathbf{C}}_{\mathrm{e}\mathrm{m}\mathrm{p}}^{(g+1)}& {\mathbf{C}}^{(g)}\end{array}\right]={\mathbf{C}}^{(g)}$。现在考虑一种稍微不同的获取${\mathbf{C}}^{(g)}$估计量的方法。

$${\mathbf{C}}_{\lambda }^{(g+1)}=\frac{1}{\lambda }\sum _{i=1}^{\lambda }\left({\mathbf{x}}_i^{(g+1)}-{\mathbf{m}}^{(g)}\right){\left({\mathbf{x}}_i^{(g+1)}-{\mathbf{m}}^{(g)}\right)}^{\mathrm{T}}\text{\hspace{1em}(10)}$$

同样${\mathbf{C}}_{\lambda }^{(g+1)}$也是${\mathbf{C}}^{(g)}$的无偏估计，式(9)和(10)之间的明显不同在于参考均值不同，${\mathbf{C}}_{\mathrm{e}\mathrm{m}\mathrm{p}}^{(g+1)}$中使用的是实际采样的均值，而${\mathbf{C}}_{\lambda }^{(g+1)}$使用的则是使用采样分布的真实均值${\mathbf{m}}^{(g)}$，因此估计${\mathbf{C}}_{\mathrm{e}\mathrm{m}\mathrm{p}}^{(g+1)}$和${\mathbf{C}}_{\lambda }^{(g+1)}$可以有不同的理解：前者估计了采样点内的分布方差，后者估计了采样步${\mathbf{x}}_i^{(g+1)}-{\mathbf{m}}^{(g)}$的方差。

公式(10)对原始协方差矩阵进行评估，为了评估更好的协方差矩阵，对(10)进行了改进，同样采用(6)中的“加权选择机制”：

$${\mathbf{C}}_{\mu }^{(g+1)}=\sum _{i=1}^{\mu }{w}_i\left({\mathbf{x}}_{i:\lambda }^{(g+1)}-{\mathbf{m}}^{(g)}\right){\left({\mathbf{x}}_{i:\lambda }^{(g+1)}-{\mathbf{m}}^{(g)}\right)}^{\mathrm{T}}\text{\hspace{1em}(11)}$$

${\mathbf{C}}_{\mu }^{(g+1)}$是对选择步上分布的估计，如同${\mathbf{C}}_{\lambda }^{(g+1)}$为选择之前所有步上的原始分布的估计。从${\mathbf{C}}_{\mu }^{(g+1)}$中采样往往会再生出选择的即成功的步，从而给出一个“更好的”协方差矩阵的解释。

为了保证使用(5)、(6)h和(11)时${\mathbf{C}}_{\mu }^{(g+1)}$是可靠的估计量，方差有效选择质量${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}$必须足够大：${\mathbf{C}}_{\mu }^{(g)}$在$f_{\text{sphere }}(\mathbf{x})={\sum }_{i=1}^n{x}_i^2$的条件数小于10就需要${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\approx 10n$。下一步就来说明如何规避对${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}$的这种限制。

3.3.2 秩$\mu$更新

为了实现快速搜索，种群大小$\lambda$必须尽量小，由于${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\approx \lambda /4$，所以${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}$肯定也会很小，例如假设${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\le 1+\ln n$，那么根据式(11)不可能得到好的协方差矩阵的可靠估计，作为补救措施，将额外使用前几代的信息，例如经过足够多代之后，估计协方差矩阵在所有代上的均值为：
$${\mathbf{C}}^{(g+1)}=\frac{1}{g+1}\sum _{i=0}^g\frac{1}{{\sigma }^{(i{)}^2}}{\mathbf{C}}_{\mu }^{(i+1)}\text{\hspace{1em}(13)}$$

使用改均值就可以得到选择步上的可靠估计，为了使得不同代之间的${\mathbf{C}}_{\mu }^{(g)}$可比较，结合了不同的${\sigma }^{(i)}$.

在(13)中，所有的代的步上都有相同的权重，为了分配最近几代更高的权重，引入了指数平滑。取${\mathbf{C}}^{(0)}=\mathbf{I}$且学习率$0<c_{\mu} \leq 1 $，那么$\boldsymbol{C}^{(g+1)}$取值为：
$$\begin{array}{rl}{\mathbf{C}}^{(g+1)}& =\left(1-c_{\mu }\right){\mathbf{C}}^{(g)}+c_{\mu }\frac{1}{{\sigma }^{(g{)}^2}}{\mathbf{C}}_{\mu }^{(g+1)}\\ & =\left(1-c_{\mu }\right){\mathbf{C}}^{(g)}+c_{\mu }\sum _{i=1}^{\mu }{w}_i{\mathbf{y}}_{i:\lambda }^{(g+1)}{\mathbf{y}}_{i:\lambda }^{(g+1{)}^{\mathrm{T}}}\\ & ={\mathbf{C}}^{(g{)}^{1/2}}\left(\mathbf{I}+c_{\mu }\sum _{i=1}^{\mu }{w}_i\left({\mathbf{z}}_{i:\lambda }^{(g+1)}{\mathbf{z}}_{i:\lambda }^{(g+1{)}^{\mathrm{T}}}-\mathbf{I}\right)\right){\mathbf{C}}^{(g{)}^{1/2}}\end{array}\text{\hspace{1em}(14)}$$

其中

• $c_{\mu }\le 1$:协方差矩阵更新的学习率，对于$c_{\mu }=1$,则没有保留先验信息，${\mathbf{C}}^{(g+1)}=\frac{1}{{\sigma }^{(g{)}^2}}{\mathbf{C}}_{\mu }^{(g+1)}$，而对于$c_{\mu }=0$，则没有进行学习，${\mathbf{C}}^{(g+1)}={\mathbf{C}}_0$。因此$c_{\mu }\approx \min \left(1,{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}/n^2\right)$式合理的选择。
• ${\mathbf{y}}_{i:\lambda }^{(g+1)}=\left({\mathbf{x}}_{i:\lambda }^{(g+1)}-{\mathbf{m}}^{(g)}\right)/{\sigma }^{(g)}$
• ${\mathbf{z}}_{i:\lambda }^{(g+1)}={\mathbf{C}}^{(g{)}^{-1/2}}{\mathbf{y}}_{i:\lambda }^{(g+1)}$为变异向量。

这个协方差矩阵更新称为秩$\mu$更新，因为(14)中外积和的秩为$min(\mu ,n)$。其实下面几段可以不用关心的！！！

$1/c_{\mu }$是backward time horizon，约占总体权重的63%。由于式(14)可以扩展为加权和：
$${\mathbf{C}}^{(g+1)}={\left(1-c_{\mu }\right)}^{g+1}{\mathbf{C}}^{(0)}+c_{\mu }\sum _{i=0}^g{\left(1-c_{\mu }\right)}^{g-i}\frac{1}{{\sigma }^{(i{)}^2}}{\mathbf{C}}_{\mu }^{(i+1)}\text{\hspace{1em}(15)}$$

backward time horizon $\Delta g$约占总体权重的63%，定义为：
$$c_{\mu }\sum _{i=g+1-\Delta g}^g{\left(1-c_{\mu }\right)}^{g-i}\approx 0.63\approx 1-\frac{1}{\mathrm{e}}\text{\hspace{1em}(16)}$$
分解求和得到
$${\left(1-c_{\mu }\right)}^{\Delta g}\approx \frac{1}{\mathrm{e}}\text{\hspace{1em}(17)}$$
然后使用泰勒近似分解得到$\Delta g$：
$$\Delta g\approx \frac{1}{c_{\mu }}\text{\hspace{1em}(18)}$$

$c_{\mu }$的选择很关键，小值导致学习慢，大值导致失败，因为协方差矩阵退化。幸运的是，一个好的设置似乎在很大程度上独立于要优化的函数。一个好的一阶近似是$c_{\mu }\approx {\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}/n^2$，因此(14)的特征时域约等于$n^{2}/\mu_{\mathrm{eff}} $.

即使学习率$c_{\mu }=1$,适应协方差矩阵也不能在一代之内完成，原始样本分布的影响直到足够的代数才会消失。假设搜索代价固定(函数评估的次数)，较小的种群规模$\lambda$允许较大的代数，因此通常会导致协方差矩阵更快的自适应。

3.3.3 秩1更新

在3.1节，使用从单代中所有选择的步从零开始估计了完整的协方差矩阵，现在采取恰恰相反的视角，只使用一个选定的步骤来重复更新代序列中的协方差矩阵，即使上一代最好的样本出现的几率更高。首先，这一视角将给出适应规则(14)的正当性，其次引入所谓的进化路径，最终用于协方差矩阵的秩1更新。

3.3.3.1 不同的视角

考虑一种特殊的方法来产生均值为零的n维正态分布，令向量${\mathbf{y}}_1,\dots ,{\mathbf{y}}_{g_0}\in {\mathbb{R}}^n,g_0\ge n,$ span ${\mathbb{R}}^n$，且$\mathcal{N}(0,1)$表示服从独立(0,1)正太分布的随机数，那么
$$\mathcal{N}(0,1){\mathbf{y}}_1+\cdots +\mathcal{N}(0,1){\mathbf{y}}_{g_0}\sim \mathcal{N}\left(0,\sum _{i=1}^{g_0}{\mathbf{y}}_i{\mathbf{y}}_i^{\mathrm{T}}\right)\text{\hspace{1em}(19)}$$

是均值为0协方差${\sum }_{i=1}^{g_0}{\mathbf{y}}_i{\mathbf{y}}_i^{\mathrm{T}}$的正太分布随机向量。该随机向量是通过线性分布$\mathcal{N}(0,1){\mathbf{y}}_i$相加得到的。奇异分布$\mathcal{N}(0,1){\mathbf{y}}_i\sim \mathcal{N}\left(0,{\mathbf{y}}_i{\mathbf{y}}_i^{\mathrm{T}}\right)$考虑了所有均值为0的正太分布，生成最大似然估计向量${\mathbf{y}}_i$。

考虑(19)和(14)的轻微简化，试图深入了解协方差矩阵的适应规则，令(14)中的和仅包含一个被加项(例如$\mu =1$)，且令${\mathbf{y}}_{g+1}=\frac{{\mathbf{x}}_{1:\lambda }^{(g+1)}-{\mathbf{m}}^{(g)}}{{\sigma }^{(g)}}$，那么协方差矩阵秩1更新写为：
$${\mathbf{C}}^{(g+1)}=\left(1-c_1\right){\mathbf{C}}^{(g)}+c_1{\mathbf{y}}_{g+1}{\mathbf{y}}_{g+1}^{\mathrm{T}}\text{\hspace{1em}(20)}$$

右边的被加项秩为1，再在协方差矩阵${\mathbf{C}}^{(g)}$加上${\mathbf{y}}_{g+1}$的最大似然项，因此在下一代中生成${\mathbf{y}}_{g+1}$的概率就增加了。

3.3.3.2 累计：利用进化路径

我们已经使用选择的步骤${\mathbf{y}}_{i:\lambda }^{(g+1)}=\left({\mathbf{x}}_{i:\lambda }^{(g+1)}-{\mathbf{m}}^{(g)}\right)/{\sigma }^{(g)}$来更新(14)和(20)中的协方差矩阵，因为$\mathbf{y}{\mathbf{y}}^{\mathrm{T}}=-\mathbf{y}(-\mathbf{y}{)}^{\mathrm{T}}$，所以步骤的符号与协方差矩阵的更新无关，也就是说，符号信息不用于计算${\mathbf{C}}^{(g+1)}$。为了利用符号信息，引入了所谓的进化路径。

我们将策略在许多代上采取的连续步序列称之为一个进化路径。进化路径可以用一系列连续的步骤总和来表示，这种总和称为累积。为了构建一个进化路径，不考虑步长$\sigma$。例如，分布均值为$\mathbf{m}$的三步进化路径为：
$$\frac{m^{(g+1)}-m^{(g)}}{{\sigma }^{(g)}}+\frac{m^{(g)}-m^{(g-1)}}{{\sigma }^{(g-1)}}+\frac{m^{(g-1)}-m^{(g-2)}}{{\sigma }^{(g-2)}}\text{\hspace{1em}(21)}$$
实际上，为了构造进化路径，使用了指数平滑，且${\mathbf{p}}_{\mathrm{c}}^{(0)}=0$：
$${\mathbf{p}}_{\mathrm{c}}^{(g+1)}=\left(1-c_{\mathrm{c}}\right){\mathbf{p}}_{\mathrm{c}}^{(g)}+\sqrt{c_{\mathrm{c}}\left(2-c_{\mathrm{c}}\right){\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}}\frac{{\mathbf{m}}^{(g+1)}-{\mathbf{m}}^{(g)}}{{\sigma }^{(g)}}\text{\hspace{1em}(22)}$$
其中

• ${\mathbf{p}}_{\mathrm{c}}^{(g)}\in {\mathbb{R}}^n$表示第$g$代的进化路径
• $c_{\mathrm{c}}\le 1$，同样，$1/c_{\mathrm{c}}$为进化路径${\mathbf{p}}_{\mathrm{c}}$的backward time horizon 。

因子$\sqrt{c_{\mathrm{c}}\left(2-c_{\mathrm{c}}\right){\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}}$是${\mathbf{p}}_{\mathrm{c}}$的归一化常数，对于$c_{\mathrm{c}}=1$和${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}=1$，该因子变为1，且${\mathbf{p}}_{\mathrm{c}}^{(g+1)}=\left({\mathbf{x}}_{1:\lambda }^{(g+1)}-{\mathbf{m}}^{(g)}\right)/{\sigma }^{(g)}$。

考虑进化路径${\mathbf{p}}_{\mathrm{c}}^{(g+1)}$的协方差矩阵$C^{(g)}$秩1更新写作：
$${\mathbf{C}}^{(g+1)}=\left(1-c_1\right){\mathbf{C}}^{(g)}+c_1{\mathbf{p}}_{\mathbf{c}}^{(g+1)}{\mathbf{p}}_{\mathbf{c}}^{(g+1{)}^{\mathrm{T}}}\text{\hspace{1em}(26)}$$
(26)中一种经验有效的学习率设置为$c_1\approx 2/n^2$，如果$c_1=1$且$\mu =1$，那么等式(26)、(20)和(14)就是完全相同的。

使用进化路径更新$\mathbf{C}$是当${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}$较小时对(14)的显著改进，因为利用了连续步骤之间的相关性。

3.4 组合秩$\mu$更新和累积

最终的协方差矩阵CMA更新就是组合(14)和(26)：
$$\begin{array}{rl}{\mathbf{C}}^{(g+1)}=& \left(1-c_1-c_{\mu }\right){\mathbf{C}}^{(g)}\\ & +c_1\underset{\text{rank-one update }}{\underbrace{{\mathbf{p}}_{\mathrm{c}}^{(g+1)}{\mathbf{p}}_{\mathrm{c}}^{(g+1)}}}+c_{\mu }\underset{\text{rank-}\mu \text{ update}}{\underbrace{\sum _{i=1}^{\mu }{w}_i{\mathbf{y}}_{i:\lambda }^{(g+1)}{\left({\mathbf{y}}_{i:\lambda }^{(g+1)}\right)}^{\mathrm{T}}}}\end{array}\text{\hspace{1em}(27)}$$

其中

• $c_1\approx 2/n^2$
• $c_{\mu }\approx \min \left({\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}/n^2,1-c_1\right)$
• ${\mathbf{y}}_{i:\lambda }^{(g+1)}=\left({\mathbf{x}}_{i:\lambda }^{(g+1)}-{\mathbf{m}}^{(g)}\right)/{\sigma }^{(g)}$

如果$c_1=0$则等式(27)退化为(14)，如果$c_{\mu }=0$则退化为(26)。(27)融合了(14)和(26)的优势，一方面，秩$\mu$更新有效地利用了一代种群内的信息，另一方面，秩1更新的进化路径利用了代之间的相关性信息。前者在大群体中很重要，后者在小群体中作用明显。

3.4 步长控制

在上一节中介绍的协方差矩阵自适应并没有明确地控制分布的“总体尺度”，步长。协方差矩阵自适应对于每一步选择只在一个方向上增加尺度，通过因子1 - c1 - c隐式地减少旧信息来降低尺度。除了$C^{(g)}$的自适应规则以外，有两点原因来引入步长控制：

• 最优整体步长不能通过(27)很好地近似，特别是当${\mu }_{eff}$大于1时。
• (27)中协方差矩阵更新的最大可靠学习率太慢，以至于无法实现整个步长的竞争性变化率。

为了控制步长${\sigma }^{(g)}$，我们利用到了进化路径，即连续步总和。该方法可独立于协方差矩阵更新而应用，可表示为累积路径长度控制(cumulative path length control)、累积步长控制(cumulative step-size control)或累积步长自适应( cumulative step
length adaptation,CSA)。基于以下推理，充分利用了进化路径的长度：

• 只要进化路径很短，单个步骤就会相互抵消，笼统地讲，它们是反相关的。如果步长互相湮灭，步长应该减小。
• 只要进化的路径很长，单个的步骤都指向相似的方向，笼统地讲，它们是相关的。因为步骤是相似的，同样的距离可以通过更少但更长的步骤进入相同的方向。在极限情况下，连续的步长有相同的方向，它们可以被一个放大的单步长代替。因此，应该增加步长。
• 在理想情况下，期望步之间是(近似)垂直的，因此也就不相关。

为了确定进化路径是“长”还是“短”，我们将路径的长度与随机选择下的期望长度进行比较。在随机选择下，连续的步骤是独立的，因此是不相关的。如果选择导致进化路径比期望值长，则应该增加$\sigma$，如果选择导致进化路径比期望值短，则应该降低$\sigma$。在理想情况下，选择不会对进化路径的长度产生偏差，在随机选择下，进化路径的长度等于其期望长度。

实际上，为了构建进化路径${\mathbf{p}}_{\sigma }$，采用了(22)中相同的技术。但和(22)不同的是，这里构建的是共轭进化路径，因为(22)中的进化路径${\mathbf{p}}_c$的期望长度取决于其方向。初始化${\mathbf{p}}_{\sigma }^{(0)}=0$，共轭进化路径写为：
$${\mathbf{p}}_{\sigma }^{(g+1)}=\left(1-c_{\sigma }\right){\mathbf{p}}_{\sigma }^{(g)}+\sqrt{c_{\sigma }\left(2-c_{\sigma }\right){\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}}{\mathbf{C}}^{(g{)}^{-\frac{1}{2}}}\frac{{\mathbf{m}}^{(g+1)}-{\mathbf{m}}^{(g)}}{{\sigma }^{(g)}}\text{\hspace{1em}(28)}$$

其中

• ${\mathbf{p}}_{\sigma }^{(g)}\in {\mathbb{R}}^n$为第$g$的代的共轭进化路径；
• $c_{\sigma }<1$，$1/c_{\sigma }$为进化路径的backward time horizon；
• $\sqrt{c_{\sigma }\left(2-c_{\sigma }\right){\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}}$为归一化常数，见(22)；
• ${\mathbf{C}}^{(g{)}^{-\frac{1}{2}}}\stackrel{\text{ def }}{=}{\mathbf{B}}^{(g)}{\mathbf{D}}^{(g{)}^{-1}}{\mathbf{B}}^{(g{)}^{\mathrm{T}}}$，其中${\mathbf{C}}^{(g)}={\mathbf{B}}^{(g)}{\left({\mathbf{D}}^{(g)}\right)}^2{\mathbf{B}}^{(g{)}^{\mathrm{T}}}$是${\mathbf{C}}^{(g)}$的特征分解，${\mathbf{B}}^{(g)}$是特征向量的标准正交基，且对角矩阵${\mathbf{D}}^{(g)}$的对角元素为对应正特征值的平方根。

对于${\mathbf{C}}^{(g)}=\mathbf{I}$，则有${\mathbf{C}}^{(g{)}^{-\frac{1}{2}}}=\mathbf{I}$，那么(28)就等价于(22)。转换${\mathbf{C}}^{(g{)}^{-\frac{1}{2}}}$就是在由${\mathbf{B}}^{(g)}$给定的坐标系内对${\mathbf{m}}^{(g+1)}-{\mathbf{m}}^{(g)}$进行重新缩放。

因此转换${\mathbf{C}}^{(g{)}^{-\frac{1}{2}}}$使得${\mathbf{p}}_{\sigma }^{(g+1)}$的期望长度与其方向无关，对于任何已实现的协方差矩阵序列$C_{g=0,1,2,\dots }^{(g)}$，给定${\mathbf{p}}_{\sigma }^{(0)}\sim \mathcal{N}(0,\mathbf{I})$，在随机选择下都有${\mathbf{p}}_{\sigma }^{(g+1)}\sim \mathcal{N}(0,\mathbf{I})$。

为了更新${\sigma }^{(g)}$，比较了$\Vert {\mathbf{p}}_{\sigma }^{(g+1)}\Vert$与其期望长度$\mathrm{E}\Vert \mathcal{N}(0,\mathbf{I})$，即：
$$\ln {\sigma }^{(g+1)}=\ln {\sigma }^{(g)}+\frac{c_{\sigma }}{d_{\sigma }\mathrm{E}\Vert \mathcal{N}(0,\mathbf{I})\Vert }\left(\Vert {\mathbf{p}}_{\sigma }^{(g+1)}\Vert -\mathrm{E}\Vert \mathcal{N}(0,\mathbf{I})\Vert \right)\text{\hspace{1em}(29)}$$

其中

• $d_{\sigma }\approx 1$为阻尼参数，缩放$\ln {\sigma }^{(g)}$的改变幅度；
• $\mathrm{E}\Vert \mathcal{N}(0,\mathbf{I})\Vert =\sqrt{2}\Gamma \left(\frac{n+1}{2}\right)/\Gamma \left(\frac{n}{2}\right)\approx \sqrt{n}+\mathcal{O}(1/n)$表示服从$\mathcal{N}(0,\mathbf{I})$分布的随机向量的欧几里得范数的期望。

由于${\sigma }^{(g)}>0$(29)等价于
$${\sigma }^{(g+1)}={\sigma }^{(g)}\exp \left(\frac{c_{\sigma }}{d_{\sigma }}\left(\frac{\Vert {\mathbf{p}}_{\sigma }^{(g+1)}\Vert }{\mathrm{E}\Vert \mathcal{N}(0,\mathbf{I})\Vert }-1\right)\right)\text{\hspace{1em}(34)}$$

4 算法总结

别看卡面给出了那么多公式，但实际上算法优化过程只用到了**(5)(6)(22)(27)(28)(34)**。

各符号说明如下：

• ${\mathbf{y}}_k\sim \mathcal{N}(0,\mathbf{C}),$ for $k=1,\dots ,\lambda$为均值为0协方差矩阵为$\mathbf{C}$的多元正太分布实现；
• $B,D$来自于协方差矩阵$\mathbf{C}$的特征分解，即 $\mathbf{C}=\mathbf{B}{\mathbf{D}}^2{\mathbf{B}}^T={\mathbf{B}\mathbf{D}\mathbf{D}\mathbf{B}}^{\mathrm{T}}$ 。$\mathbf{B}$的列为特征向量的正交基，对角矩阵$\mathbf{D}$的对角元素为对应正特征值的平方根。
• ${\mathbf{x}}_k\in {\mathbb{R}}^n,$ for $k=1,\dots ,\lambda$，$\lambda$个搜索点的样本；
• $⟨\mathbf{y}{⟩}_{\mathrm{w}}={\sum }_{i=1}^{\mu }{w}_i{\mathbf{y}}_{i:\lambda }$，不考虑步长$\sigma$时的分布均值步；
• ${\mathbf{y}}_{i:\lambda }=\left({\mathbf{x}}_{i:\lambda }-\mathbf{m}\right)/\sigma$，如下；
• ${\mathbf{x}}_{i:\lambda }\in {\mathbb{R}}^n$，根据式(37)得到的${\mathbf{x}}_1,\dots ,{\mathbf{x}}_{\lambda }$中的$i$个最优点；
• ${\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}={\left({\sum }_{i=1}^{\mu }{w}_i^2\right)}^{-1}$方差有效选择质量，见式(8)，其满足$1\le {\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\le \mu$；
• $C^{-\frac{1}{2}}\stackrel{\text{ def }}{=}B{D}^{-1}{B}^{\mathrm{T}}$；
• $\mathrm{E}\Vert \mathcal{N}(0,\mathbf{I})\Vert =\sqrt{2}\Gamma \left(\frac{n+1}{2}\right)/\Gamma \left(\frac{n}{2}\right)\approx \sqrt{n}\left(1-\frac{1}{4n}+\frac{1}{21n^2}\right)$
• $h_{\sigma }=\{\begin{array}{ll}1& \text{ if }\frac{\Vert p_{\sigma }\Vert }{\sqrt{1-{\left(1-c_{\sigma }\right)}^{2(g+1)}}}<\left(1.4+\frac{2}{n+1}\right)\mathrm{E}\Vert \mathcal{N}(0,\mathbf{I})\Vert \\ 0& \text{ otherwise }\end{array}$，其中$g$为代数。如果$\Vert {\mathbf{p}}_{\sigma }\Vert$很大，亥维赛函数$h_{\sigma }$停滞对(42中)${\mathbf{p}}_{\mathrm{c}}$的更新，这可以防止在线性环境中即当步长太小时C轴的过快增长。当初始步长选择得过小或目标函数随时间变化时，这是有用的。
• $\delta \left(h_{\sigma }\right)=\left(1-h_{\sigma }\right)c_{\mathrm{c}}\left(2-c_{\mathrm{c}}\right)\le 1$无关紧要。