非线性规划中的自适应协方差矩阵调整演化策略（CMA-ES）基础题
题目描述
考虑一个无约束非线性优化问题：
\[\min_{x \in \mathbb{R}^n} f(x) \]
其中 \(f(x)\) 是一个连续、可微但可能非凸、多峰的函数，其解析梯度难以获取或计算成本高昂。我们要求使用自适应协方差矩阵调整演化策略（CMA-ES） 来求解该问题。本题目标是通过一个具体实例，详细解释CMA-ES的基本原理、关键步骤和参数更新机制，确保你能理解这种基于种群的随机优化方法如何通过自适应调整搜索分布来高效寻找全局最优解。
假设我们以经典的Rastrigin函数为例：
\[f(x) = 10n + \sum_{i=1}^{n} \left[ x_i^2 - 10 \cos(2\pi x_i) \right] \]
其中 \(n = 2\)，搜索范围设定为 \(x_i \in [-5, 5]\)。该函数具有多个局部极小点，全局最小值在 \(x = (0,0)\) 处，函数值为 \(0\)。我们将通过CMA-ES寻找该全局最小值。
解题过程循序渐进讲解
第一步：理解CMA-ES的核心思想
CMA-ES是一种进化策略算法，其核心是通过维护一个多元正态分布 \(\mathcal{N}(m, \sigma^2 C)\) 来指导搜索：
- \(m \in \mathbb{R}^n\)：分布均值，代表当前最优解的估计位置。
- \(\sigma > 0\)：全局步长，控制搜索范围的大小。
- \(C \in \mathbb{R}^{n \times n}\)：协方差矩阵，描述搜索分布的形态（如椭圆的方向和拉伸程度）。
算法迭代过程如下：
- 采样：从当前分布中生成一组候选解（种群）。
- 评估：计算每个候选解的目标函数值。
- 选择与更新：根据函数值选择表现较好的解，并用这些解的信息更新分布参数 \(m, \sigma, C\)，使分布逐渐向全局最优区域移动并调整形状以适应函数地形。
第二步：初始化参数
对于 \(n = 2\) 的Rastrigin函数，我们初始化：
- 均值 \(m^{(0)}\)：在搜索范围 \([-5,5]^2\) 内随机均匀采样，例如 \(m^{(0)} = (2, -2)\)。
- 步长 \(\sigma^{(0)}\)：通常设为搜索范围的一半左右，例如 \(\sigma^{(0)} = 2\)。
- 协方差矩阵 \(C^{(0)}\)：初始化为单位矩阵 \(I_2\)，表示各向同性搜索。
- 种群大小 \(\lambda\)：根据经验公式 \(\lambda = 4 + \lfloor 3 \ln n \rfloor \approx 7\)（实际应用中常取稍大值，这里为简化取 \(\lambda = 10\)）。
- 父代数量 \(\mu\)：选择种群中最好的 \(\mu\) 个解用于更新，通常 \(\mu = \lfloor \lambda/2 \rfloor = 5\)。
- 其他算法参数：如累积更新权重、学习率等，采用CMA-ES标准设置（具体见后续步骤）。
第三步：采样生成候选解
在第 \(g\) 代，我们从当前分布 \(\mathcal{N}(m^{(g)}, (\sigma^{(g)})^2 C^{(g)})\) 中采样 \(\lambda\) 个候选解：
\[x_k^{(g)} = m^{(g)} + \sigma^{(g)} \cdot y_k, \quad y_k \sim \mathcal{N}(0, C^{(g)}), \quad k = 1, \dots, \lambda \]
其中 \(y_k\) 是通过协方差矩阵 \(C^{(g)}\) 生成的多元正态随机向量。实际操作中，我们通过 Cholesky 分解 \(C = A A^T\) 来高效采样：\(y = A z\)，\(z \sim \mathcal{N}(0, I)\)。
第四步：评估与排序
计算每个候选解的目标函数值 \(f(x_k^{(g)})\)，并按值从小到大排序（最小化问题）。选择前 \(\mu\) 个最优解作为“父代”解，记为 \(x_{i:\lambda}^{(g)}\)（\(i=1,\dots,\mu\)）。
第五步：更新均值 \(m\)
新的均值通过加权平均得到：
\[m^{(g+1)} = \sum_{i=1}^{\mu} w_i \cdot x_{i:\lambda}^{(g)} \]
其中权重 \(w_i\) 通常取正且和为1，赋予更好解更高权重，例如 \(w_i = \frac{\ln(\mu+0.5) - \ln i}{\sum_{j=1}^{\mu} (\ln(\mu+0.5) - \ln j)}\)。
第六步：更新演化路径和步长 \(\sigma\)
CMA-ES使用演化路径（evolution path） 来记录均值更新的连续方向，以自适应调整步长。
- 各向同性演化路径 \(p_\sigma\)：
\[ p_\sigma^{(g+1)} = (1 - c_\sigma) p_\sigma^{(g)} + \sqrt{c_\sigma (2 - c_\sigma) \mu_{\text{eff}}} \cdot \frac{m^{(g+1)} - m^{(g)}}{\sigma^{(g)}} \]
其中：
- \(c_\sigma\) 是学习率（通常 \(\approx 4/n\)）。
- \(\mu_{\text{eff}} = 1 / \sum_{i=1}^{\mu} w_i^2\) 是有效父代数量。
- 初始 \(p_\sigma^{(0)} = 0\)。
- 步长更新：
\[ \sigma^{(g+1)} = \sigma^{(g)} \cdot \exp\left( \frac{c_\sigma}{d_\sigma} \left( \frac{\| p_\sigma^{(g+1)} \|}{\mathbb{E} \| \mathcal{N}(0,I) \|} - 1 \right) \right) \]
其中 \(d_\sigma\) 是阻尼系数（通常 \(\approx 1\)）。当演化路径长度比随机游动预期长度长时，步长增加；反之减少。
第七步：更新协方差矩阵 \(C\)
协方差矩阵更新也依赖于演化路径，以学习搜索分布的形状。
- 各向异性演化路径 \(p_c\)：
\[ p_c^{(g+1)} = (1 - c_c) p_c^{(g)} + h_\sigma \sqrt{c_c (2 - c_c) \mu_{\text{eff}}} \cdot \frac{m^{(g+1)} - m^{(g)}}{\sigma^{(g)}} \]
其中 \(c_c\) 是学习率（通常 \(\approx 4/n\)），初始 \(p_c^{(0)} = 0\)。\(h_\sigma\) 是一个指示函数，当 \(\| p_\sigma \|\) 未过大时取1，防止在步长迅速增加时破坏 \(p_c\)。
2. 协方差矩阵更新：
\[ C^{(g+1)} = (1 - c_1 - c_\mu) C^{(g)} + c_1 \cdot p_c^{(g+1)} (p_c^{(g+1)})^T + c_\mu \sum_{i=1}^{\mu} w_i \cdot y_{i:\lambda}^{(g)} (y_{i:\lambda}^{(g)})^T \]
其中：
- \(c_1\) 是秩-one更新的学习率（利用演化路径方向）。
- \(c_\mu\) 是秩-μ更新的学习率（利用当前父代解的形状）。
- \(y_{i:\lambda}^{(g)} = (x_{i:\lambda}^{(g)} - m^{(g)}) / \sigma^{(g)}\) 是标准化的搜索步长。
更新使 \(C\) 逐渐对齐目标函数的等高线方向。
第八步：迭代与终止
重复步骤三至七，直到满足终止条件，例如：
- 步长 \(\sigma\) 小于阈值（如 \(10^{-8}\)）。
- 函数值变化很小或达到最大迭代次数。
示例演示（简化迭代）
假设经过几代迭代：
- 均值 \(m\) 逐渐向 \((0,0)\) 靠近。
- 步长 \(\sigma\) 自适应减小，实现精细搜索。
- 协方差矩阵 \(C\) 可能变为对角优势（对于Rastrigin函数，其变量分离，最优解各向同性，但CMA-ES仍会学习局部地形）。
最终，算法将收敛到全局最优解 \(x \approx (0,0)\)，函数值 \(f \approx 0\)。
总结
CMA-ES通过自适应调整正态分布的均值、步长和协方差矩阵，有效地在复杂、多峰的非线性函数中进行全局优化。其关键在于利用演化路径积累信息，平衡探索（大步长、各向同性）与利用（小步长、对齐地形）。对于梯度不可用或计算昂贵的黑箱优化问题，CMA-ES是一种强大且鲁棒的方法。