CMA-ES

共分散行列適応進化戦略 (CMA-ES) は、数値最適化のための特別な種類の戦略です。進化戦略(ES)は、非線形または非凸連続最適化問題 の数値最適化のための確率的かつ導関数を使用しない方法です。これらは、進化アルゴリズムおよび進化計算のクラスに属します。進化アルゴリズムは、生物学的進化の原理、すなわち変異 (組み換えと突然変異による) と選択の繰り返される相互作用に広く基づいています。各世代 (反復) において、新しい個体 (候補解、 と表記) が、通常は確率的な方法で現在の親個体の変異によって生成されます。次に、いくつかの個体が、その適応度または目的関数値に基づいて、次の世代の親になるように選択されます。このように、世代シーケンスを通じて、 より良い -値を持つ個体が生成されます。 ${\displaystyle x}$${\displaystyle f(x)}$${\displaystyle f}$

進化戦略では、新しい候補ソリューションは通常、の多変量正規分布に従ってサンプリングされます。 再結合は、分布の新しい平均値を選択することです。 突然変異は、ランダムベクトル、つまり平均がゼロの摂動を追加することです。 分布内の変数間のペアワイズ依存関係は、共分散行列によって表されます。 共分散行列適応 (CMA) は、この分布の共分散行列を更新する方法です。 これは、関数が悪条件である場合に特に便利です。 ${\displaystyle \mathbb {R} ^{n}}$${\displaystyle f}$

共分散行列の適応は、古典的な最適化における準ニュートン法の逆ヘッセ行列の近似に似た、基礎となる目的関数の2次モデルを学習することに相当する。ほとんどの古典的な方法とは対照的に、基礎となる目的関数に関する仮定は少なくなる。候補ソリューションのランキング（または、同等のソート）のみが利用されるため、この方法では導関数も（明示的な）目的関数も必要ありません。たとえば、スイス式トーナメントでの候補ソリューション間のペアワイズ競争からランキングを導き出すことができます。

CMA-ES アルゴリズムでは、検索分布のパラメータを適応させるための 2 つの主な原則が活用されています。

まず、最大尤度原理は、成功する候補ソリューションと検索ステップの確率を高めるという考えに基づいています。分布の平均は、以前に成功した候補ソリューションの尤度が最大化されるように更新されます。分布の共分散行列は、以前に成功した検索ステップの尤度が高くなるように（増分的に）更新されます。両方の更新は、自然勾配降下法として解釈できます。また、結果として、CMA は、すべての主軸を保持しながら、成功した検索ステップの反復主成分分析を実行します。分布の推定アルゴリズムとクロスエントロピー法は非常によく似た考えに基づいていますが、成功する検索ステップではなく、成功するソリューションポイントの尤度を最大化することによって共分散行列を（非増分的に）推定します。

次に、戦略の分布平均の時間発展の 2 つのパスが記録されます。これらは検索パスまたは発展パスと呼ばれます。これらのパスには、連続するステップ間の相関に関する重要な情報が含まれています。具体的には、連続するステップが同様の方向に実行される場合、発展パスは長くなります。発展パスは 2 つの方法で利用されます。1 つのパスは、単一の成功した検索ステップの代わりに共分散行列適応手順に使用され、好ましい方向の分散の増加を可能な限り高速化します。もう 1 つのパスは、追加のステップ サイズ制御を実行するために使用されます。このステップ サイズ制御の目的は、分布平均の連続的な動きを期待どおりに直交させることです。ステップ サイズ制御は、早期収束を効果的に防止しながら、最適値への高速収束を可能にします。

以下では、最も一般的に使用される ( μ / μ _w、  λ )-CMA-ES の概要を示します。各反復ステップでは、λ 個の新しい候補ソリューションのうちμ個の最良のものの重み付けされた組み合わせを使用して、分布パラメータを更新します。メイン ループは、1) 新しいソリューションのサンプリング、2) 適合性に基づくサンプリングされたソリューションの並べ替え、3) 並べ替えられたサンプルに基づく内部状態変数の更新という 3 つの主要部分で構成されます。アルゴリズムの疑似コードは次のようになります。

set  //  反復ごとのサンプル数、少なくとも 2、通常は 4 より大きい
。initialize、、、、、 //終了しない間に  状態変数を初期化する
。do   // で反復する
    。do //新しいソリューションを
          サンプリングして評価する。sample_multivariate_normal(mean 、 covariance_matrix )
         ← with // ソリューションをソートする
      。// 後で必要になる。and ← update_m   // 平均をより良いソリューションに移動 
    ← update_ps //  等方性進化パスを更新する
    ← update_pc   // 異方性進化パスを更新する
    ← update_C   // 共分散行列を更新
    する← update_sigma   // 等方性パスの長さを使用してステップ サイズを更新する
。return or${\displaystyle \lambda}$ ${\displaystyle m}$${\displaystyle \sigma}$${\displaystyle C=I}$${\displaystyle p_{\sigma }=0}$${\displaystyle p_{c}=0}$   ${\displaystyle i}$  ${\displaystyle \{1\ldots \lambda \}}$ ${\displaystyle \lambda}$${\displaystyle x_{i}={}}$${\displaystyle {}=m}$${\displaystyle {}=\sigma ^{2}C}$${\displaystyle f_{i}=\operatorname {適応度} (x_{i})}$
    ${\displaystyle x_{1\ldots \lambda}}$${\displaystyle x_{s(1)\ldots s(\lambda )}}$${\displaystyle s(i)=\operatorname {argsort} (f_{1\ldots \lambda },i)}$${\displaystyle m'=m}$${\displaystyle 'mm'}$${\displaystyle x_{i}-m'}$       
    ${\displaystyle m}$${\displaystyle (x_{1},\ldots,x_{\lambda})}$${\displaystyle p_{\sigma}}$${\displaystyle (p_{\sigma },\sigma ^{-1}C^{-1/2}(mm'))}$${\displaystyle p_{c}}$${\displaystyle (p_{c},\sigma ^{-1}(mm'),\|p_{\sigma }\|)}$${\displaystyle C}$${\displaystyle (C,p_{c},(x_{1}-m')/\sigma ,\ldots ,(x_{\lambda }-m')/\sigma )}$${\displaystyle \sigma}$${\displaystyle (\sigma ,\|p_{\sigma }\|)}$ ${\displaystyle m}$${\displaystyle x_{1}}$

5 つの更新割り当ての順序は重要です。を最初に更新し、 の前に更新し、 を最後に更新する必要があります。 5 つの状態変数の更新方程式は、次のように指定されます。 ${\displaystyle m}$${\displaystyle p_{\sigma}}$${\displaystyle p_{c}}$${\displaystyle C}$${\displaystyle \sigma}$

探索空間の次元と反復ステップが与えられている。5つの状態変数は ${\displaystyle n}$${\displaystyle k}$

${\displaystyle m_{k}\in \mathbb {R} ^{n}}$、分布平均と最適化問題の現在の最良解、

${\displaystyle \sigma_{k}>0}$、ステップサイズ、

${\displaystyle C_{k}}$、対称かつ正定値の 共分散行列であり、${\displaystyle n\times n}$ ${\displaystyle C_{0}=I}$

${\displaystyle p_{\sigma}\in \mathbb {R} ^{n},p_{c}\in \mathbb {R} ^{n}}$2 つの進化パスがあり、最初はゼロ ベクトルに設定されています。

反復は、多変量正規分布から候補解をサンプリングすることから始まります。つまり、${\displaystyle \lambda >1}$${\displaystyle x_{i}\in \mathbb {R} ^{n}}$ ${\displaystyle \textstyle {\mathcal {N}}(m_{k},\sigma _{k}^{2}C_{k})}$${\displaystyle i=1,\ldots ,\lambda }$

${\displaystyle {\begin{aligned}x_{i}\ &\sim \ {\mathcal {N}}(m_{k},\sigma _{k}^{2}C_{k})\\&\sim \ m_{k}+\sigma _{k}\times {\mathcal {N}}(0,C_{k})\end{aligned}}}$

2行目は、現在の最良解ベクトル（分布平均ベクトル）の不偏摂動（突然変異）として解釈することを示唆している。候補解は、最小化されるべき目的関数に基づいて評価される。ソートされた候補解を次のように 表す。${\displaystyle m_{k}}$${\displaystyle x_{i}}$${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} }$${\displaystyle f}$

${\displaystyle \{x_{i:\lambda }\mid i=1\dots \lambda \}=\{x_{i}\mid i=1\dots \lambda \}{\text{ および }}f(x_{1:\lambda })\leq \dots \leq f(x_{\mu :\lambda })\leq f(x_{\mu +1:\lambda })\leq \cdots ,}$

新しい平均値は次のように計算される。

${\displaystyle {\begin{aligned}m_{k+1}&=\sum _{i=1}^{\mu }w_{i}\,x_{i:\lambda }\\&=m_{k}+\sum _{i=1}^{\mu }w_{i}\,(x_{i:\lambda }-m_{k})\end{aligned}}}$

ここで、正の（再結合）重みの合計は1 になります。通常、重みは となるように選択されます。ここで、および以下で目的関数から使用される唯一のフィードバックは、インデックス によるサンプルされた候補ソリューションの順序付けです。 ${\displaystyle w_{1}\geq w_{2}\geq \dots \geq w_{\mu }>0}$${\displaystyle \mu \leq \lambda /2}$${\displaystyle \textstyle \mu _{w}:=1/\sum _{i=1}^{\mu }w_{i}^{2}\approx \lambda /4}$${\displaystyle i:\lambda }$

ステップ サイズは、累積ステップ サイズ適応(CSA)を使用して更新されます。これは、パス長制御とも呼ばれます。進化パス (または検索パス)が最初に更新されます。 ${\displaystyle \sigma _{k}}$${\displaystyle p_{\sigma }}$

${\displaystyle p_{\sigma }\gets \underbrace {(1-c_{\sigma })} _{\!\!\!\!\!{\text{discount factor}}\!\!\!\!\!}\,p_{\sigma }+\overbrace {\sqrt {1-(1-c_{\sigma })^{2}}} ^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\text{complements for discounted variance}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}\underbrace {{\sqrt {\mu _{w}}}\,C_{k}^{\;-1/2}\,{\frac {\overbrace {m_{k+1}-m_{k}} ^{\!\!\!{\text{displacement of }}m\!\!\!}}{\sigma _{k}}}} _{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\text{distributed as }}{\mathcal {N}}(0,I){\text{ under neutral selection}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}}$

${\displaystyle \sigma _{k+1}=\sigma _{k}\times \exp {\bigg (}{\frac {c_{\sigma }}{d_{\sigma }}}\underbrace {\left({\frac {\|p_{\sigma }\|}{\operatorname {E} \|{\mathcal {N}}(0,I)\|}}-1\right)} _{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\text{unbiased about 0 under neutral selection}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}{\bigg )}}$

どこ

${\displaystyle c_{\sigma }^{-1}\approx n/3}$は進化経路の逆方向の時間範囲であり、1より大きい（は指数関数的減衰定数を連想させる。ここで は関連する寿命であり、半減期である）。${\displaystyle p_{\sigma }}$${\displaystyle c_{\sigma }\ll 1}$${\displaystyle (1-c_{\sigma })^{k}\approx \exp(-c_{\sigma }k)}$${\displaystyle c_{\sigma }^{-1}}$${\displaystyle c_{\sigma }^{-1}\ln(2)\approx 0.7c_{\sigma }^{-1}}$

${\displaystyle \mu _{w}=\left(\sum _{i=1}^{\mu }w_{i}^{2}\right)^{-1}}$は分散有効選択質量であり、の定義により、${\displaystyle 1\leq \mu _{w}\leq \mu }$${\displaystyle w_{i}}$

${\displaystyle C_{k}^{\;-1/2}={\sqrt {C_{k}}}^{\;-1}={\sqrt {C_{k}^{\;-1}}}}$は の逆数の唯一の対称平方根であり、${\displaystyle C_{k}}$

${\displaystyle d_{\sigma }}$減衰パラメータは通常 1 に近くなります。またはの場合、ステップ サイズは変更されません。${\displaystyle d_{\sigma }=\infty }$${\displaystyle c_{\sigma }=0}$

が期待値より大きい場合にのみステップサイズが増加する。${\displaystyle \sigma _{k}}$${\displaystyle \|p_{\sigma }\|}$

${\displaystyle {\begin{aligned}\operatorname {E} \|{\mathcal {N}}(0,I)\|&={\sqrt {2}}\,\Gamma ((n+1)/2)/\Gamma (n/2)\\&\approx {\sqrt {n}}\,(1-1/(4\,n)+1/(21\,n^{2}))\end{aligned}}}$

小さい場合は減少する。このため、ステップサイズの更新は、適応が成功した後の連続ステップを-共役にする傾向がある。^[1]${\displaystyle C_{k}^{-1}}$${\displaystyle \textstyle \left({\frac {m_{k+2}-m_{k+1}}{\sigma _{k+1}}}\right)^{T}\!C_{k}^{-1}{\frac {m_{k+1}-m_{k}}{\sigma _{k}}}\approx 0}$

最後に、共分散行列が更新されますが、ここでもそれぞれの進化パスが最初に更新されます。

${\displaystyle p_{c}\gets \underbrace {(1-c_{c})} _{\!\!\!\!\!{\text{discount factor}}\!\!\!\!\!}\,p_{c}+\underbrace {\mathbf {1} _{[0,\alpha {\sqrt {n}}]}(\|p_{\sigma }\|)} _{\text{indicator function}}\overbrace {\sqrt {1-(1-c_{c})^{2}}} ^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\text{complements for discounted variance}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}\underbrace {{\sqrt {\mu _{w}}}\,{\frac {m_{k+1}-m_{k}}{\sigma _{k}}}} _{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\text{distributed as}}\;{\mathcal {N}}(0,C_{k})\;{\text{under neutral selection}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}}$

${\displaystyle C_{k+1}=\underbrace {(1-c_{1}-c_{\mu }+c_{s})} _{\!\!\!\!\!{\text{discount factor}}\!\!\!\!\!}\,C_{k}+c_{1}\underbrace {p_{c}p_{c}^{T}} _{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\text{rank one matrix}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}+\,c_{\mu }\underbrace {\sum _{i=1}^{\mu }w_{i}{\frac {x_{i:\lambda }-m_{k}}{\sigma _{k}}}\left({\frac {x_{i:\lambda }-m_{k}}{\sigma _{k}}}\right)^{T}} _{\operatorname {rank} \min(\mu ,n){\text{ matrix}}}}$

ここで は転置を表し、 ${\displaystyle T}$

${\displaystyle c_{c}^{-1}\approx n/4}$進化経路の逆方向の時間範囲であり、1より大きい。${\displaystyle p_{c}}$

${\displaystyle \alpha \approx 1.5}$そして、指示関数は の場合にのみ と評価され、言い換えれば となり、通常は となる。${\displaystyle \mathbf {1} _{[0,\alpha {\sqrt {n}}]}(\|p_{\sigma }\|)}$ ${\displaystyle \|p_{\sigma }\|\in [0,\alpha {\sqrt {n}}]}$${\displaystyle \|p_{\sigma }\|\leq \alpha {\sqrt {n}}}$

${\displaystyle c_{s}=(1-\mathbf {1} _{[0,\alpha {\sqrt {n}}]}(\|p_{\sigma }\|)^{2})\,c_{1}c_{c}(2-c_{c})}$指標がゼロの場合の小さな分散損失を部分的に補う。

${\displaystyle c_{1}\approx 2/n^{2}}$は共分散行列のランク1更新の学習率であり、

${\displaystyle c_{\mu }\approx \mu _{w}/n^{2}}$は共分散行列のランク更新の学習率であり、 を超えてはなりません。${\displaystyle \mu }$${\displaystyle 1-c_{1}}$

共分散行列の更新により、およびが からサンプリングされる尤度が増加する傾向があります。これで反復ステップが完了します。 ${\displaystyle p_{c}}$${\displaystyle (x_{i:\lambda }-m_{k})/\sigma _{k}}$${\displaystyle {\mathcal {N}}(0,C_{k+1})}$

反復ごとの候補サンプルの数 は事前に決定されておらず、広い範囲で変化する可能性があります。 などの値が小さいほど、検索動作がローカルになります。 などの値が大きくなると、検索はよりグローバルになります。 アルゴリズムは、再起動ごとに が係数 2 ずつ増加して繰り返し再起動されることがあります。[ ^2]設定（または、代わりに、たとえば が使用可能なプロセッサの数によって事前に決定されている場合）以外に、上記で紹介したパラメータは、特定の目的関数に固有のものではなく、したがってユーザーが変更することを意図したものではありません。 ${\displaystyle \lambda }$${\displaystyle \lambda =10}$${\displaystyle \lambda =10n}$${\displaystyle \mu _{w}\approx \lambda /4}$${\displaystyle \lambda }$${\displaystyle \lambda }$${\displaystyle \mu }$${\displaystyle \lambda }$

function  xmin = purecmaes % ( mu/mu_w, lambda ) - CMA - ES 
% -------------------- 初期化 --------------------------------   % ユーザー定義の入力パラメータ（編集する必要があります）strfitnessfct = 'frosenbrock' ; % 目的関数/適応度関数の名前N = 20 ; % 目的変数の数/問題次元xmean = rand ( N , 1 ); % 目的変数の初期点sigma = 0.3 ; % 座標ごとの標準偏差（ステップ サイズ）stopfitness = 1e-10 ; % 適応度 < stopfitness の場合は停止（最小化）stopeval = 1e3 * N ^ 2 ; % 関数評価の stopeval 回数後に停止% 戦略パラメータ設定: 選択  lambda = 4 + floor ( 3 * log ( N )); % 個体群のサイズ、子孫の数mu = lambda / 2 ; % 組み換えの親/ポイントの数weights = log ( mu + 1 / 2 ) - log ( 1 : mu ) ' ; % 重み付け組み換えの muXone 配列mu = floor ( mu ); weights = weights / sum ( weights ); % 組み換え重み配列を正規化mueff = sum ( weights ) ^ 2 / sum ( weights .^ 2 ); % 合計 w_i x_i の分散有効性  
  
      
                   
        
              
      
       
  
  
      
                   
     
            
         
   

  % 戦略パラメータ設定: 適応
cc = ( 4 + mueff / N ) / ( N + 4 + 2 * mueff / N ); % C の累積の時定数cs = ( mueff + 2 ) / ( N + mueff + 5 ); % シグマ制御の累積の t 定数c1 = 2 / (( N + 1.3 ) ^ 2 + mueff ); % C のランク 1 更新の学習率cmu = min ( 1 - c1 , 2 * ( mueff - 2 + 1 / mueff ) / (( N + 2 ) ^ 2 + mueff )); % ランク mu 更新の場合dumps = 1 + 2 * max ( 0 , sqrt (( mueff - 1 ) / ( N + 1 )) - 1 ) + cs ; % シグマの減衰% 通常は 1 に近い% 動的 (内部) 戦略パラメータと定数を初期化pc = zeros ( N , 1 ); ps = zeros ( N , 1 ); % C とシグマの進化パスB = eye ( N , N ); % B は座標系を定義しますD = ones ( N , 1 ); % 対角線 D はスケーリングを定義しますC = B * diag ( D .^ 2 ) * B ' ; % 共分散行列 C invsqrtC = B * diag ( D .^- 1 ) * B '          
        
          
           
          
                                                      
  
          
                           
                          
                    
        ; % C^-1/2 eigeneval = 0 ; % B と D の更新を追跡chiN = N ^ 0.5 * ( 1 - 1 / ( 4 * N ) + 1 / ( 21 * N ^ 2 )); % 期待値% ||N(0,I)|| == norm(randn(N,1)) % -------------------- 生成ループ -------------------------------- counteval = 0 ; % 次の 40 行に 20 行の興味深いコードが含まれていますwhile counteval < stopeval % ラムダ子孫を生成して評価しますfor k = 1 : lambda arx (:, k ) = xmean + sigma * B * ( D .* randn ( N , 1 )); % m + sig * Normal(0,C) arfitness ( k ) = feval ( strfitnessfct , arx (:, k )); % 目的関数呼び出しcounteval = counteval + 1 ; end % 適応度でソートし、加重平均を計算して xmean に格納します[ arfitness , arindex ] = sort ( arfitness ); % 最小化xold = xmean ; xmean = arx (:, arindex ( 1 : mu )) * weights ; % 再結合、新しい平均値% 累積: 進化パスを更新しますps = ( 1 - cs ) * ps ... + sqrt ( cs * ( 2 - cs ) * mueff ) * invsqrtC * ( xmean - xold ) / sigma ; hsig = norm ( ps ) / sqrt (    
                          
    
                                      
  
      
     
    
      
       
                     
              
            
      
    
      
          
        
           
    
      
          
                    
        1 - ( 1 - cs ) ^ ( 2 * counteval / lambda )) / chiN < 1.4 + 2 / ( N + 1 ); pc = ( 1 - cc ) * pc ... + hsig * sqrt ( cc * ( 2 - cc ) * mueff ) * ( xmean - xold ) / sigma ;    
         
                   

      % 共分散行列を適応 C 
artmp = ( 1 / sigma ) * ( arx (:, arindex ( 1 : mu )) - repmat ( xold , 1 , mu )); C = ( 1 - c1 - cmu ) * C ...                   % 古い行列を考慮  + c1 * ( pc * pc ' ...                  % ランク 1 更新+ ( 1 - hsig ) * cc * ( 2 - cc ) * C ) ...  % hsig==0 の場合はマイナー修正+ cmu * artmp * diag ( weights ) * artmp ' ; % ランク mu 更新          
           
               
                         
                   

      % ステップ サイズ sigma を適応します。
sigma = sigma * exp (( cs / dumps ) * ( norm ( ps ) / chiN - 1 )); % counteval - eigeneval > lambda / ( c1 + cmu ) / N / 10の場合、C を B*diag(D.^2)*B' (対角化) に分解して、O(N^2) を達成します。eigeneval = counteval ; C = triu ( C ) + triu ( C , 1 ) ' ; % 対称性を強制します。[ B , D ] = eig ( C ); % 固有分解、B== 正規化された固有ベクトルD = sqrt ( diag ( D )); % D は標準偏差のベクトルです。invsqrtC = B * diag ( D .^- 1 ) * B ' ; end % ブレーク、適応度が十分に高いか条件が 1e14 を超える場合、より良い終了方法が推奨されますif arfitness ( 1 ) <= stopfitness || max ( D ) > 1e7 * min ( D ) break ; end             
    
      
             
            
               
                       
                    
                
      
    
      
               
          
      

  end % while、生成ループの終了 

  xmin = arx (:, arindex ( 1 )); % 最後の反復のベストポイントを返します。% xmean はさらに % 良くなると予想されることに注意してください。end % ---------------------------------------------------------------   function f = frosenbrock ( x ) if size ( x , 1 ) < 2 error ( 'dimension must be greater one' ); end f = 100 * sum (( x ( 1 : end - 1 ) .^ 2 - x ( 2 : end )) .^ 2 ) + sum (( x ( 1 : end - 1 ) - 1 ) .^ 2 ); end