この記事は #kosen10s Advent Calendar 2017 の 9日目の記事です．前回は marin72_com さんの「#kosen10s だより 2017」でした． kosen10s，自由気ままにだけど秩序持ってのんびりやってていいところだと思います．僕みたいなプログラミング雑魚勢でもちゃんと居場所があるので，少しでもいいなと思った2010年卒の高専生はslackチームへおいでなさいまし．

はじめに

今回は自分の研究テーマの1つである進化計算について，ちょっと紹介しようかと思います．まだまだ勉強中であるため，ところどころぼかしますがご容赦を．

遺伝的アルゴリズム(Genetic Algorithm; GA)や蟻コロニー最適化(Ant Colony Optimization; ACO)など進化計算にも様々な種類がありますが，対象とする問題は基本的に同じです．進化計算が取り扱う問題においては，最適値を求めたい関数，つまり目的関数について次の1や2のような制約が課されている場合が多いです．

目的関数を明に書き表すことができない（ただし，入力値に対応する出力値を得ることはできる）
目的関数の導関数が不明である

このような制約の中で行う最適化は「Black-box最適化」という問題に属します．つまり，進化計算はBlack-box最適化を行うための最適化アルゴリズムと言うことができます．ちなみに，Black-box最適化アルゴリズムはいくつか提案されており，例えばベイズ的最適化(Bayesian Optimization)もこれに該当します．

Covariance Matrix Adaptation Evolution Strategies

今回は，数あるBlack-box最適化アルゴリズムの中から，性能が比較的良好なCMA-ES(Covariance Matrix Adaptation Evolution Strategies)*1について紹介します．CMA-ESは現在でも改良が続けられており様々な手法*2 *3が提案されています．

CMA-ESとは，目的関数への入力値（以後，サンプルと呼びます）を平均ベクトルと共分散行列をパラメータに持つ正規分布からサンプリングし，得られた出力値を基に正規分布のパラメータを更新，更新された正規分布からまた新たにサンプルをサンプリングし,...という操作を繰り返すことによって，最適なサンプルを得るという手法です．CMA-ESにおける最適化の概要図を次のFig.1に示します．オレンジの点線の形状を決めているのが分散共分散行列で，平均ベクトルが中心を決めています．

f:id:satuma77:20170623230903p:plain

Fig.1 CMA-ESによる最適化過程の概要図*4

また，最適化をアニメーションにしたものを次のFig.2に示します．最適解は左上の箇所となります．

f:id:satuma77:20171208235245g:plain

Fig.2 2次元Rastrigin関数に対するCMA-ESの最適化過程*5

Fig.2からわかるように，探索の過程でサンプルの広がりが大きく広がり，最適解で収束するという現象が見られます．この現象は，CMA-ESの工夫のひとつであるCommurative Step-Size Adaptation(CSA)によるものです．CMA-ESにおける共分散行列の更新には，CSA, Rank-1 Update, Rank-μ Updateの3種類の更新規則が適用されています．それぞれ説明すると非常に長くなりますので，今回はRank-μ Updateに絞って説明します．また，平均ベクトルの更新もRank-μ Updateで行われています．そこで，次節ではRank-μ Updateにおける平均ベクトルと共分散行列の更新方法について説明します．

Rank-μ Updateと確率的勾配降下法

ある時点 $t$ において， $n$ 番目のサンプルを $x_{n}$ ，平均ベクトル $m^{(t)}$ ，共分散行列を $C^{(t)}$ としたとき，CMA-ESのRank-μ Updateは次式で定義されます*6．

$m^{(t+1)} = m^{(t)} + \eta_{m} \sum_{i=1}^{\lambda} w_{\mathrm{rk}}(x_{i}) (x_{i} - m^{(t)})$

$C^{(t+1)} = C^{(t)} + \eta_{C} \sum_{i=1}^{\lambda} w_{\mathrm{rk}}(x_{i}) \left( x_{i}-m^{(t)})(x_{i}-m^{(t)})^{\mathrm{T}}-C^{(t)} \right)$

$\eta_{m}$ および $\eta_{C}$ は学習率（ユーザが指定するパラメータ）で，一回の更新に使用するサンプル数を $\lambda$ としています．また， $w_{\mathrm{rk}}(x_{i})$ はサンプル $x_{i}$ が $\lambda$ 個のサンプルの中で何番目に良い出力値かという「ランキング」に基づいた値を意味します．CMA-ESについて，これらの更新式はハンドメイドによって作成されたものとなっています．

これらの更新式は2010年に自然勾配と深い関連があることが明らかになっています*7．自然勾配*8は，確率分布モデルが形成する多様体が作る幾何的構造を分析する情報幾何学の中で発展した手法となります．

確率分布パラメータを「勾配」を用いて更新することを考えます．このとき，パラメータ間の擬距離にユークリッド距離ではない距離を採用した場合，ユークリッド空間における勾配は最急方向とは異なります．この場合の最急方向は，勾配にリーマン計量テンソルの逆行列をかけたものとなります．これが自然勾配と呼ばれるものです．情報幾何学において，確率分布モデルに対するパラメータ間の擬距離にはKL-Divergenceを採用する場合が多く，その場合のリーマン計量テンソルはフィッシャー情報行列(Fisher Information Matrix; FIM)となることが知られています．

話を戻してPure Rank-μ Updateの場合を考えます．多次元正規分布 $\pi(x; m, C)$ の対数尤度関数 $\ln \pi(x; m, C)$ に対して自然勾配を求めると

$\tilde{∇}_{m} \ln \pi(x; m, C) = (x_{n} - m)$

$\tilde{∇}_{C} \ln \pi(x; m, C) = (x_{i}-m)(x_{i}-m)^{\mathrm{T}}-C$

となります．この自然勾配を見るとPure Rank-μ Updateの更新式に類似した式がそのまま現れていることがわかります．

実際は，解くべき最適化問題を目的関数の期待値最小化問題 $\mathrm{argmin}_{x} \mathbb{E}_{\pi(x; m, C)}f(x)$ とした上で，自然勾配を用いるとPure Rank-μ Updateの式を導くことができます．詳しくは文献*9を参照していただければと思います．

追記(2017.12.09)

時間の許す限り，導出過程を書いていきます．目的関数の期待値最小化 $\mathrm{argmin}_{x} \mathbb{E}_{\pi(x; m, C)}f(x)$ を行うために， $\mathbb{E}_{\pi(x; m, C)}f(x)$ の勾配を計算します．文献*10に従って $\theta = \lbrack m^{\mathrm{T}}, \mathrm{vec}(C)^{\mathrm{T}} \rbrack ^{\mathrm{T}}$ とすることで，2つのパラメータに関する勾配をまとめて求めます． $\frac{\partial \log f(x)}{\partial x} = \frac{\partial f(x)}{\partial x} \frac{1}{f(x)}$ というテクニックを使えば，

$∇_{\theta} \mathbb{E}_{\pi(x; m, C)}f(x)$ $= ∇_{\theta} \int f(x) \pi(x; m, C) \mathrm{d}x$ $= \int f(x) ∇_{\theta} \pi(x; m, C) \mathrm{d}x$ $= \int f(x) \frac{∇_{\theta} \pi(x; m, C)}{\pi(x; m, C)} \pi(x; m, C) \mathrm{d}x$ $= \int \left( f(x) ∇_{\theta} \log \pi(x; m, C) \right) \pi(x; m, C) \mathrm{d}x$ $= \mathbb{E}_{\pi(x; m, C)} \left\lbrack f(x) ∇_{\theta} \log \pi(x; m, C) \right\rbrack$

とこのように勾配は期待値の形で書き表すことができます．今回 $f(x)$ が未知であるため，この期待値を代数的に求めることは難しいです．そこで，モンテカルロ法によって勾配を推定します． $\lambda$ 個のサンプルを使って勾配を近似すれば，

$∇_{\theta} \mathbb{E}_{\pi(x; m, C)}f(x) \approx \frac{1}{\lambda} \sum_{i=0}^{\lambda} f(x) ∇_{\theta} \log \pi(x; m, C)$

と求めることができます．このとき，FIMは次の式となります（文献*10）． $F = \lbrack \lbrack C^{-1}, 0 \rbrack, \lbrack 0, \frac{1}{2}(C^{-1} \bigotimes C^{-1}) \rbrack \rbrack$

FIMの逆行列を先程の勾配の式と掛け合わせると，自然勾配は

$\tilde{∇}_{m} \mathbb{E}_{\pi(x; m, C)}f(x) = \sum_{i=1}^{\lambda} f(x_{i}) (x_{i} - m^{(t)})$ $\tilde{∇}_{C} \mathbb{E}_{\pi(x; m, C)}f(x) = \sum_{i=1}^{\lambda} f(x_{i}) \left( x_{i}-m^{(t)})(x_{i}-m^{(t)})^{\mathrm{T}}-C^{(t)} \right)$

となります．ここで悪スケールな関数に対してうまく最適化できるように，自然勾配を目的関数 $f(x)$ のスケールに依存しないように変更します．その変換方法として目的関数値 $f(x)$ のかわりに，解のランキングに基づいた値 $w_{\mathrm{rk}}(x_{i})$ を使用します．これで，Pure Rank-μ Updateが自然勾配によって更新するアルゴリズムということが確認されました．

Pure Rank-μ UpdateのPython実装

先程のPure Rank-μ Updateについては，私が現在作成しているPython用の進化計算ライブラリEvoltierにて実装済みですので，軽く紹介します．次にEvoltierにおけるPure Rank-μ Updateの実行方法を示します．（リポジトリにあるmain.pyを一部改編）

肝心のPure Rank-μ Updateの部分はevoltier.optimizers.GaussianNaturalGradientOptimizerにあります．実際に動かす際には，最適化したい関数などをevoltier.optimizers.updater.Updaterに入れてから最適化を行うように作成しています．Updaterの初期化におけるパラメータについて重要な部分のみ触れていきます．pop_sizeは $\lambda$ に相当し，thresholdは停止条件に関するパラメータです．具体的には目的関数の値がこの値を下回るサンプルを得られれた時点で停止するというものになります．max_iterも停止条件に関するパラメータであり，1試行における最大更新回数を設定するものとなります．

実行速度については100次元のSphre関数で $\lambda=1000$ としたとき，評価値が $10^{-30}$ となる値を見つけるまでに約40秒ほどでした．まだ，並列化には対応していないこともあり個体数が多いときはmultiprocessingモジュールの利用によって高速化も期待されます．

おわりに

今回はBlack-box最適化アルゴリズムの1つであるCMA-ESについて，説明しました．また，Pure Rank-μ Updateによる分布更新が自然勾配との関連があることも紹介しました．最後に現在作成中のライブラリEvoltierにおいて，Pure Rank-μ Updateによる最適化を実行する方法を示しました．CMA-ESについて，その派生アルゴリズムやRank-1 Updateなどの他の更新ルール，制約付き問題に対する適用方法など，もし時間があればご紹介できればと思います．また，Evoltierについてもゆっくりですが開発を進めていき，時期を見てPypiに上げて実際に使えるようにできたらと思います．とりあえずCuPyによるGPU対応も盛り込むつもりで作っています．もしよければStarください．

次回はNKMR6194さんの「なぜpthread_mutexは機能するのか？」です!

*1:N. Hansen and A. Ostermeier, “Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation,” Proceedings of IEEE International Conference on Evolutionary Computation, pp. 312–317, 1996.

*2:N. Hansen, “Benchmarking a BI-population CMA-ES on the BBOB-2009 noisy testbed,” Proceedings of the Genetic and Evolutionary Computation Conference - GECCO ’09, pp. 2389-2396, 2009.

*3:Y. Akimoto and N. Hansen, “Projection-Based Restricted Covariance Matrix Adaptation for High Dimension,” Proceedings of the Genetic and Evolutionary Computation Conference - GECCO ’16, pp. 197-204, 2016.

*4:https://en.wikipedia.org/wiki/CMA-ES

*5:http://blog.otoro.net/2017/10/29/visual-evolution-strategies/

*6:S. Shirakawa, Y. Akimoto, K. Ouchi, and K. Ohara, “Sample Reuse in the Covariance Matrix Adaptation Evolution Strategy Based on Importance Sampling,” in Proceedings of the 2015 on Genetic and Evolutionary Computation Conference - GECCO ’15, 2015, pp. 305–312.

*7:Y. Akimoto, Y. Nagata, I. Ono, and S. Kobayashi, “Bidirectional relation between CMA evolution strategies and natural evolution strategies,” Parallel Problem Solving from Nature (PPSN XI), pp. 154–163, 2010.

*8:S. Amari, “Natural Gradient Works Efficiently in Learning,” Neural Computation, vol. 10, no. 2, pp. 251–276, 1998.

*9:Y. Ollivier, L. Arnold, A. Auger, and N. Hansen, “Information-Geometric Optimization Algorithms: A Unifying Picture via Invariance Principles,” Journal of Machine Learning Research, vol. 18, no. 18, pp. 1–65, 2017.

*10:Y. Akimoto, “Analysis of a Natural Gradient Algorithm on Monotonic Convex-Quadratic-Composite Functions,” Proceedings of the fourteenth international conference on Genetic and evolutionary computation conference - GECCO '12, pp. 1293–1300, 2012.

satuma-portfolio

共分散行列適応進化戦略(CMA-ES)と自然勾配法