協方差自適應調整的進化策略（CMA-ES）

02-10

協方差自適應調整的進化策略(Covariance Matrix Adaptation Evolution Strategy, CMA-ES)是最有名，應用最多，性能最好的ES之一，在中等規模（變數個數大約在 3-300範圍內）的複雜優化問題上具有很好的效果。是Google Vizier: A Service for Black-Box Optimization 包含的四種用於調hyper-parameter 的方法之一。

按照慣例，先列舉幾個應用文章：

Path Integral Policy Improvement with Covariance Matrix Adaptation, ICML, 2012.

Optimizing Walking Controllers for Uncertain Inputs and Environments, ACM Transactions on Graphics 29. 4 (2010)

Flexible Muscle-Based Locomotion for Bipedal Creatures ， ACM Transactions on Graphics 32.6 (2013)

在強化學習和不確定環境下的行走控制取得了極好的效果。視頻2013年 Flexible Muscle-Based Locomotion for Bipedal Creatures

CMA-ES 演算法的基本特點有：

無梯度優化，不使用梯度信息,
局部搜索中無梯度演算法通常比梯度演算法慢，通常需要 $O(n)$ 倍的評估,
在複雜優化問題non-separable, ill-conditioned, or rugged/multi-modal 上表現良好。

1. CMA-ES 演算法步驟

CMA-ES的核心想法是通過對正態分布 $N(m_t,sigma_t^2C_t)$ 中協方差矩陣 $C$ 的調整來處理變數之間的依賴關係和scaling。演算法基本可以分成以下三步

採樣產生新解；
計算目標函數值；
更新分布參數 $m_t,~C_t,~sigma_t$

ES演算法設計的核心就是如何對這些參數進行調整，尤其是步長參數和協方差矩陣的調整，以達到儘可能好的搜索效果。對這些參數的調整在ES演算法的收斂速率方面有非常重要的影響。CMA-ES調整參數的基本思路是，調整參數使得產生好解的概率逐漸增大（沿好的搜索方向進行搜索的概率增大）。

CMA-ES 的搜索過程

1. 1 採樣產生新解

在CMA-ES中，每一次迭代從 $N(m_t,sigma_t^2C_t)$ 中產生 $lambda$ 個解(一般 $lambda$ 取值為log(n)的量級)

$x_i = m_t +sigma_t y_i, ~y_isim N(0,C_t)$

這裡的每一個 $y_iin R^n$ 是一個搜索方向。一般的， $y$ 可以通過協方差矩陣C的特徵分解 $C_t=BD^2B^T$ 或者Cholesky 分解由標準正態分布得到，即

$y_i=BDz_i, ~z_isim N(0,I).$

1. 2 計算目標函數值

對新產生的解計算對應的目標函數值 $f(x_i)$ ，並對這些目標函數值排序

$f(x_{1:lambda})leq f(x_{2:lambda}) leq cdotsleq f(x_{lambda:lambda})$

其中的下標表示在這些樣本中排第i位。對於截斷選擇來說，取前 $mu = lfloorfrac{lambda}{2} floor$ 個解用於更新分布參數。還存在其他的選擇方式，這裡為方便使用截斷選擇。值的注意的是，這裡的順序只依賴於目標函數值的比較，而不依賴於目標函數值本身，即comparison based，或者說 objective value-free。

1.3 分布參數更新

使用所選擇的解分別對分布參數 $m_t,~C_t,~sigma_t$ 分別進行獨立的更新。大體來說，對 $m_t, C_t$ 的更新使用的是最大似然估計 (ML-update)

1.3.1 均值

分布的均值即為所選擇的 $mu$ 個解的加權的最大似然估計

$m_{t+1}= sum_{i=1}^{mu}w_i x_{i:lambda}=m_t+ sum_{i=1}^{mu}w_i (x_{i:lambda}-m_t)$

此式表明，均值沿著平均搜索方向移動一步.

這裡通過組合所選擇的部分解得到下一代的均值的做法稱為multi-recombination，類似於其他進化演算法中的交叉(crossover operation)，即不同的解之間進行交換信息multi-recombination。
如果權重係數取相同的 $w_i = frac{1}{mu}$ , 那麼此式就是使用所選擇的解的最大似然估計。在實際演算法中，通常取 $w_1geq dots w_{mu}>0$ 以強調排在最前面的那些解。
更新項 $sum_{i=1}^{mu}w_i (x_{i:lambda}-m_t)$ 可以理解為對 $m$ 的自然梯度信息幾何優化，隨機優化, 與進化策略.

採樣，選擇和更新

1.3.2 進化路徑/搜索路徑(evolution path)

CMA-ES中，對協方差矩陣C 的更新包含rank-1 和rank- $mu$ 兩項，其中的rank-1 項使用的是歷史搜索信息——進化路徑，其構造方式為

$p_{t+1}= (1-c)p_t+ sqrt{c(2-c)} sqrt{mu_w}frac{m_{t+1} - m_t}{sigma_t}$

它描述了分布均值的移動，並且將每次迭代中移動方向 $frac{m_{t+1} - m_t}{sigma_t}$ 做加權平均，使得這些方向中相反的方向分量相互抵消，相同的分量則進行疊加。這類似於神經網路優化中常用的Momentum。在神經網路中momentum起什麼做用？因此，進化路徑代表了最好的搜索方向之一。

這裡面的係數因子是按照如下方式設計的：

因子 $mu_w=frac{1}{sum_{i=1}^{mu}w_i^2}$ 的設計是根據 $sqrt{mu_w}frac{m_{t+1} - m_t}{sigma_t}sim N(0,C_t)$ .這是因為 $sqrt{mu_w}frac{m_{t+1} - m_t}{sigma_t}= sqrt{mu_w}sum_{i=1}^{mu}w_i y_{i:lambda}$ ，因此可以看成是一個從上述分布採樣得到的隨機向量（確切的說，如果 $x_{i:lambda}$ 是隨機選擇的）。
因子 $sqrt{c(2-c)}$ 的設計原理是 $(1-c)^2+ (sqrt{c(2-c)} )^2=1$ .這兩條被稱為平穩性條件(stationarity condition)，使得 $p_{t+1}$ 本身看起來像一個從當前分布 $N(0, C_t)$ 產的搜索方向 $p_{t+1} sim N(0, C_t)$ 。所以 $p_{t+1}$ 可以像一個mutation一樣用來更新協方差矩陣。
變化率/學習率 $c$ 的設計原理是 $c^{-1} propto n$ ，即學習率與所調整的變數自由度（參數個數）成反比。

1.3.3 協方差矩陣

協方差矩陣的更新的原理是

$argmax ~p(p_{t+1}|m,C),~~argmax ~prod_{i=1}^{mu} pleft(frac{x_{i:lambda} - m_t}{sigma_t}Big|m,C ight)$

在此基礎上，協方差的更新方式為

$C_{t+1}= (1-c_1-c_{mu})C_t+ c_1p_{t+1}p_{t+1}^T+ c_{mu}sum_{i=1}^{mu} w_i left(frac{x_{i:lambda} - m_t}{sigma_t} ight)left(frac{x_{i:lambda} - m_t}{sigma_t} ight)^T$

或者更加簡潔的寫成(由於前面採樣中 $y_{i}$ 和 $x_i$ 的關係)

$C_{t+1}= (1-c_1-c_{mu})C_t+ c_1p_{t+1}p_{t+1}^T+ c_{mu}sum_{i=1}^{mu} w_i y_{i:lambda} y_{i:lambda} ^T$

其中：

$C$ 的更新原理是增大沿成功搜索方向的方差，即增大沿這些方向採樣的概率.
第一項 rank-1 update 可以理解成直接把進化路徑看作一個成功的搜索方向. 最初的CMA-ES只包含 rank-1 更新( $c_{mu}=0$ )，這對小規模種群是性能優異，例如 $lambda = O(log n)$ , 以及極限情形的每一次只產生一個新解的(1+1)-CMA-ES.
更新的第二項的秩為 $min(mu, n)$ ，由於通常 $mu =O(ln n)$ 所以這一項常被稱為rank- $mu$ update。這實際上是使用所選擇的 $mu$ 個解的加權的最大似然估計（去除步長參數）。這一項rank $mu$ update可以解釋為對 $C$ 的自然梯度信息幾何優化，隨機優化, 與進化策略，但是和前面對 $m$ 的自然梯度使用的步長不同。當種群規模較大的時候 $lambda = {n, n^2}$ ，這項更新起著更加重要的作用。
學習率 $c_1,c_{mu}$ 的設計原理和上面一樣，也就是 $c_1approxfrac{2}{n^2},~c_{mu}approxfrac{mu_w}{n^2}$ ，即學習率與所調整的變數自由度（參數個數）成反比。

1.3.4 步長更新

CMA-ES默認使用累積式步長調整(Cumulative step size adaptation， CSA) 。CSA是當前最成功、用的最多的步長調整方式。CSA的原理可以從方面來理解：相繼搜索的方嚮應該是共軛的。

如果相繼的搜索方向之間正相關（夾角小於pi/2），那麼表明步長太小，應該增大；
如果相繼搜索方向之間是負相關的，那麼表明步長太大，應該減小。

進化路徑

與前面的進化路徑相似，構造另一個進化路徑（有些文獻裡面稱為共軛路徑conjugate evolution path）

$s_{t+1}= (1-c_{sigma})s_t+ sqrt{c_{sigma}(2-c_{sigma})} sqrt{mu_w}Bsum_{i=1}^{mu} w_iz_{i:lambda}\ sigma_{t+1}=sigma_texpleft( frac{c_{sigma}}{d_{sigma}}left(frac{|s_{t+1}|}{E|N(0,I)|}-1 ight) ight)$

其中

更新項在很多論文里寫成 $Bsum_{i=1}^{mu} w_iz_{i:lambda}=C^{-frac{1}{2}}frac{m_{t+1}-m_t}{sigma_t}$ ，其中的 $C^{-frac{1}{2}}=BD^{-1}B^T$ 。因此，這個方向實際上是去掉尺度因子D之後的搜索方向。
在平穩性條件下有 $s_{t+1}sim N(0,I)$ ，即搜索路徑可以看成是一個n維標準正態分布的隨機向量。因此其模長服從卡方分布 $|s_{t+1}|sim chi(n)$ ，並且 $E(|s_{t+1}|) = sqrt{n}$ ，而 $E(|N(0,I)|)=sqrt{n}$ .因此，如果模長大於平均值，則指數上是正的，步長變大，否則指數上是負的，步長減小。這實現了前面的intuitive idea.
考慮另一種解釋. 在上述路徑中取 $c_{sigma}=1$ ，即不進行累積。這時候 $s_{t+1}=sqrt{mu_w}Bsum_{i=1}^{mu} w_iz_{i:lambda}$ 實際上是「平均搜索方向」，並且「大致」服從標準正態分布。從這個角度來說，步長調整可以理解為，如果所選擇的這些解的平均長度大於全體的平均長度，那麼長度就增加，反之如果較好的這部分的平均長度小於全體的平均，那麼長度應減小。這類似於xNES中步長調整的含義。在 $c_{sigma}<1$ 情況下的累積則代表通過歷史平均來消除/減小隨機性。
$c_{sigma}， d_{sigma}$ 是調整步長變化幅度的控制參數，通常設置為 $c_{sigma} propto frac{1}{n}, ~d_{sigma}>1$ .

註：在如果每次迭代步長几乎不變，大致有 $|s_{t+1}|sim sqrt{n}$ , 那麼會有如下近似

$(m_{t+1}-m_t)^TC_t^{-1}(m_{t+1}-m_t)approx 0$

即相繼的搜索方向關於協方差矩陣的逆 $C_t^{-1}$ 是共軛的，而在二次函數上，這個 $C_t^{-1}$ 收斂於Hessian矩陣（相差一個標量因子）。從這個角度來說， $s$ 被稱為共軛進化路徑。這個是很好的性質。

1.4 實驗：A Case Study on Rosenbrock

在12-d 的Rosenbrock 函數上的一次典型運行，橫軸為目標函數的評估次數(function evaluation)，縱軸分別為目標函數值，各個自變數的值，各軸向上的scaling和 $D_t/sigma_t$ . 可以看出，演算法逐漸逼近最優點的過程（右上，最優點的位置在(1,...,1)）就是演算法逐漸學得各個軸向和軸向上的scaling的過程（左下和右下）。

2. CMA-ES的一些性質

在二次函數上，CMA-ES中的協方差矩陣收斂於Hessian矩陣的逆，並且收斂速率是(指數)線性的 $C(t) simeq H^{-1} e^{-gamma t}$ 。因此CMA-ES 是一種不使用梯度和函數值的variable-metric 方法。
相較於傳統數學優化方法，CMA-ES 的演算法行為在目標函數的單調變換下保持不變；相較於大部分進化演算法，則保持了variable metric的多種不變性，其中最重要的就是在搜索空間的正交變換（旋轉）下保持不變。
在演算法中使用一個獨立的步長參數是為了可以快速調整整個搜索區域的尺度。由於C的更新公式中的變化率因子都非常小（與變數個數的平方成反比），使用獨立的步長因子可以在不需要對C進行精細調整的情況下快速搜索到最優點附近。
從變換的角度來說，實際上是用 $C^{-frac{1}{2}}(x-m)$ 將目標函數變成儘可能可分的（separable），並且更加接近球函數（條件數減小）。經過變換之後的函數相比原有的目標函數更加容易優化。這個過程被進一步抽象為一種自適應編碼(Adaptive encoding, AE-CMA). 一個非常有意思的演算法是將坐標下降法和自適應編碼結合起來，得到一種自適應坐標下降 (adaptive coordinate descent).
很大程度上，ES演算法可以看成一個局部搜索。ES中實現全局搜索的主要方式是使用Restart 策略啟發式優化演算法中，如何使之避免陷入局部最優解？簡單來說，就是如果演算法在某個位置不再改進，那麼直接重新初始化換一個位置進行搜索。由於大多數多峰問題上使用較大的採樣規模（即種群規模）總是有益的，因此每次初始化通常會增大採樣規模（IPOP），或者交替使用較大的和較小的種群進行全局搜索和局部搜索（Bi-Pop）。

註：有些作者指出，這裡的 $sigma$ 不應該稱為步長，而應該稱為 scale factor 或者mutation strength，原因是有時候即使 $sigma o infty$ 演算法仍然能很好的收斂，只要 $det(C_t) o 0$ 並且更快。這裡為了方便我們就稱為步長了。

自適應坐標下降=坐標下降+自適應編碼

3. 關於ES理論方面

ES本身是進化計算領域相對來說理論比較完善的。早年間，隨機優化（包含模擬退火這類演算法）常常聲稱的優勢是全局收斂性，即在 $t o infty$ 時演算法以概率一收斂到全局最優。但是在ES領域則有不同的看法。

We believe that global convergence is rather meaningless in practice. The pure random search converge with probability one to the global optimum of functions belonging to a broad class, where the main assumption is that a neighborhood of the global optimum should be reachable by the search distribution with a positive probability. However, the convergence rate is sub-linear with degree 1/n, therefore, the running time is proportional to $(1/epsilon)^n$ . Even for moderate dimension, e.g., n=10, this is prohibitively slow in practice.

He (M. Powell) believed that local optimality is all one can get in the nonconvex case.

因此，ES領域的理論研究主要關注演算法的收斂速率方面。ES領域通常認為，全局收斂性實際上沒什麼意義，我們從來不可能真的讓 $t o infty$ ，而且這往往是通過犧牲局部搜索的效率換來的. 真正重要的是給定初始化之後，能夠快速地收斂到局部最優。如果這個局部最優仍然太差不可接受，那麼通過restart 重新搜索。

ES領域的理論分析（大部分情況下不包含協方差調整部分）主要包含以下三個方面（方法）：

Progress rate：這個是ES領域傳統的理論分析方法，H.-G. Beyer的整本書The Theory of Evolution Strategies都是基於此方法。Progress rate考慮的是演算法一步迭代中對目標函數的改進幅度的期望。在只考慮步長調整的情況下，通過對維度n取極限能極大的簡化問題的複雜性。目前這仍然是最重要的理論分析方法之一。
Convergence rate：收斂速率的定義方式基本類似於經典優化裡面收斂速率，除了搜索的點的序列是隨機的。近年來在這個方向有較大的進展Linear Convergence of Comparison-based Step-size Adaptive Randomized Search，最主要的結論之一就是證明了這種基於比較的演算法(comparison-based adaptive randomized search) 在一大類問題上能夠線性收斂，並且收斂速率 $gamma propto frac{1}{n}$ ，即與變數數目成反比(這是無梯度優化的普遍特點Some Comments on Gradient-Free Optimization)。
Runtime lower bound: 考慮的是為達到某個給定的精度所需要的迭代（function evaluation）的次數。由於不對維度取極限，能夠對較低維度的問題給出更加準確的bound，如 $(1,lambda)$ -ES在1/5法則的步長調整下runtime bound為 $nlambda ln(1/ epsilon)/sqrt{ln lambda}$ 。這方面的很多文章發在Theoretical Computer Science上。

除此之外，隨著對自然梯度和ES的關係認識的逐步加深，近年來有一些從自然梯度流（正態分布流形上的微分方程）的角度來分析演算法在極限情形下的收斂性質。這方面最主要的結論有：

$m(t) - m^* simeq H^{-1}(m_0 - m^*)e^{-gamma t}\ C(t)simeq H^{-1}e^{-gamma t}$

即梯度流方程是線性收斂的，並且收斂速率 $gamma propto frac{1}{sqrt{n}}$ 。這個速率遠大於前面提到的以及實際觀察到的 $gamma propto frac{1}{n}$ ，目前沒有很好的解釋. 這方面更多內容和資料可參考信息幾何優化，隨機優化, 與進化策略

4. CMA-ES 的一些改進

在前面的演算法步驟中，每一次迭代產生新解都需要進行矩陣分解。矩陣分解的複雜度與變數個數的三次方成正比，這使得演算法每次迭代需要耗費的計算非常多。同時注意到，在採樣中實際使用的並不是協方差矩陣本身，而是它的「平方根」，即滿足 $C=AA^T$ 的矩陣A。因此，一種更高效的方案是不直接更新和分解協方差矩陣，而是直接迭代更新矩陣A，從而將複雜度降到二次方。這就是Cholesky CMA-ES。值的注意的是，這裡並不一定要求A是上三角的。使用上三角矩陣的更新方式寫起來更加複雜，不便於向量化編程，但在採樣過程中的矩陣-向量乘法計算會少一半。

上面指出協方差矩陣的更新中，rank- $mu$ 部分可以解釋成自然梯度。那麼更直接的，如果考慮以平方根A為參數的正態分布，使用自然梯度下降就得到自然進化策略（natural evolution strategies, NES）. 這一思路下，一個更加有效、避免估計Fisher矩陣的方案是使用局部坐標系和指數映射的exponential NES。

為進一步降低演算法的內部複雜度，近年來有一些工作使用特殊形式的協方差矩陣，避免使用稠密矩陣造成的計算量。後面我們會單獨介紹這部分工作。

步長調整方式CSA具有極好的數學性質和性能，但是仍然有一些缺點。其一，CSA強烈依賴於標準正態分布，如果基礎不是標準的正態分布，CSA可能就會不成立，演算法性能會嚴重下降。其次，如果採樣規模（即種群規模）很大，CSA的性能會變差。近年來，開始有一些將傳統的1/5法則推廣到population的情形的方法，並且性能很好。不過目前這些方法普遍還缺少理論分析。

5. 程序

下面列舉一些CMA-ES的程序：

一個python的庫pypi-cma，是N. Hansen本人更新維護的，安裝和使用特別簡單.
CMA-ES 源代碼的主頁CMA-ES Source Code，包含各種語言寫的程序，尤其是多種C++的實現（包含並行libcmaes）和Matlab的版本CMA-ES Matlab。
Shark Machine Learning Library，這個包含著上面的網頁中，單獨拿出來說是因為這是一個廣泛的包含多種機器學習和進化計算方面的C++庫。這個庫在Windows上安裝略有不便，在Linux上一次編譯完成感覺太好了。

Wiki 主頁上有一個非常簡短的Matlab 程序CMA-ES - Wikipedia，只有最基本的演算法實現，可以用於了解演算法。如果使用的話，應該儘可能使用上面的完整的Code，包含了多種變體，對約束、雜訊的處理，以及restart 策略處理複雜優化問題。

部分參考文獻

[1] The CMA Evolution Strategy: A Tutorial

[2] Adaptive Coordinate Descent

[3] The Theory of Evolution Strategies

[4] Theory of Evolution Strategies: a New Perspective

[5] Linear Convergence of Comparison-based Step-size Adaptive Randomized Search

[6] Convergence Analysis of Evolutionary Algorithms That Are Based on the Paradigm of Information Geometry

[7] Impacts of invariance in search: When CMA-ES and PSO face ill-conditioned and non-separable problems

P.S. 寫的匆忙，肯定還有諸多謬誤，後面會逐步更新完善。