Fire Engine

消防士→ITエンジニア→研究者

グラフィカルモデルに基づく因果探索手法の調査

最近,因果推論や因果探索に興味を持ち,勉強している.というのも最近,ゆううきさん と一緒に分散システムの異常の原因を即時に診断するための研究を進めている.原因を診断するためのアプローチとして,サーバやコンテナ等から取得できる様々なメトリック(CPU使用率やメモリ使用率など)を(グラフ理論における)ノードとして,因果グラフを構築することを考えている.メトリック同士の単なる「相関」ではなく,結果と原因の関係である「因果」を捉えようとするアプローチである.例えば,システムの障害が発生した場合,相関だけでは,AとBが関連がありそうというところまでしか言えないが,因果を特定できると理想的には,Aの原因はBであるといった議論ができるため,有用だと考えている.
実際に,前述のような因果グラフを構築して障害の原因を特定しようというアプローチは,以下の例に挙げるようにここ数年で増えている印象がある.

今回の記事は,因果探索手法の中でもグラフィカルモデルに基づくアプローチに絞って,まとめていく.以下のサーベイ論文がよくまとまっていたので全面的に参考にした.

Review of Causal Discovery Methods Based on Graphical Models(2019)

グラフィカルモデル

グラフィカルモデルとは,確率モデルをグラフを用いて記述したものである.ここでいうグラフとは, x軸, y軸があるグラフではなく,グラフ理論のことであり,ノード(頂点)とエッジ(枝)の集合で構成されるものである.グラフィカルモデルを用いることで,どの確率変数同士が直接的に関係しているかを視覚的に表現できる.
上記でエッジと書いた部分は,単なる線と矢印の二つの表現があり,それぞれを用いたグラフを無向グラフおよび有向グラフと呼ぶ.変数間の因果を表現する際に利用されるのは主に後者の有向グラフである.有向グラフの代表的な例としてベイジアンネットワークが挙げられる.

因果探索

因果探索とは,データから因果関係(因果グラフ)を推定する手法である.前述のグラフィカルモデルで考えると,考慮すべき確率変数(ノード)のデータは手元にある状態で,因果の構造,すなわちどのようにエッジが繋がっているかを推定するということになる.多くの分野において,データの根底にある因果関係を見つけ出し,それを利用することは非常に重要であるため,因果探索は,ゲノミクス,臨床医学,疫学など幅広い分野に応用されている. 次から,因果探索の具体的なアプローチについてまとめていく.

制約ベースの因果探索

制約ベースのアプローチでは,まず,観測変数にどのような条件付き独立性が成り立つかを推測する.次に推測された条件付き独立性を制約として,それを満たすような因果グラフを探索する.ここでは,代表的な制約ベースのアルゴリズムとして,PCアルゴリズム[1]を説明する.

PCアルゴリズム

PCアルゴリズムの処理の流れを以下の図を用いて説明する.前提として,PCアルゴリズムでは,未観測共通原因(潜在的な交絡因子)がないことを仮定している.

f:id:hirotsuru314:20201007095925p:plain
PCアルゴリズムによる因果グラフの推定

  1. 完全無向グラフを構築する.
  2. 条件付き独立性が成立するエッジをグラフから除外する.
    ここで,条件付き独立性を検索する際は,条件として与える変数の数を0から1ずつ増やしながら,条件付き独立性が成立しうる全ノードの組合せを対象に行われる.
  3. V-structureに基づきエッジの方向を決定する.
    V-structureでは,三つの変数A,B,Cが「A-B-C」でA-C間にエッジがなく,A,Cを条件付き独立とする条件ノード群にBが含まれていないならば,「A→B←C」の関係が成立するとして,無向グラフ中の一部の因果の方向を推定する.
  4. オリエンテーションルールに基づきエッジの方向を決定する.
    オリエンテーションルールは有向非巡回グラフ(DAG)の定義から導かれるグラフの構造に関するルールである.具体的なルールについては「Estimating High-Dimensional Directed Acyclic Graphs with the PC-Algorithm」にまとめてある.

PCアルゴリズムは,一般にノード数が増加すると上記2の条件付き独立性を検索する際のノードの組合せの数が増大し,計算時間がかかるという問題がある.これに対して,条件付き独立性の検定をGPUで並列化して高速に処理するといった提案[2]もあるようだ.また,PCアルゴリズムは未観測共通原因がある場合のFast causal inferece(FCI)アルゴリズム[3]などへ拡張されている.

スコアベースの因果探索

スコアベースのアプローチでは,同じ条件付き独立性を与える因果グラフの集合であるマルコフ同値類ごとに,モデルのよさを測るスコアをつける.ここでは,代表的なスコアベースのアルゴリズムとして,Greedy Equivalence Search(GES)アルゴリズム[4]を説明する.

GESアルゴリズム

GESアルゴリズムでは,PCアルゴリズムのように完全無向グラフから始めるのではなく,空のグラフから始まり,必要なエッジを追加したり,不要なエッジを削除していきながら因果グラフを構築する.アルゴリズムの各ステップでは,グラフに有向エッジを追加することで,ベイズ情報量規準(BIC)のようなスコアや統計的仮説検定のZスコアで計測される適合性が向上するかどうかが決定され,適合性を最も向上させるエッジが追加される.スコアをそれ以上改善できなくなった場合,GESアルゴリズムは,どのエッジを削除するとスコアが最も改善されるかをエッジごとに検証していく.このような手順でスコア(例えばBIC)を用いて因果グラフを推定していく.

構造方程式モデルに基づく因果探索(非ガウス性・非線形性の活用)

先に挙げた制約ベース・スコアベースの因果探索は,データの構造に対してできるだけ仮定をおかずに,データのみから因果グラフを推測するアプローチである.これは,後述する構造方程式モデルにおける関数形と誤差変数の分布に仮定をおかないノンパラメトリックアプローチに分類される.仮定をおかないため様々なデータに適用できるという利点があるが,因果の構造を正確に推定できないケースも多い.そこで,分析者の事前知識を仮定としてモデルに取り入れるアプローチがある.ここでは,仮定を取り入れる上で重要な概念である構造方程式モデルについて述べ,非ガウス性の仮定を利用したLiNGAMという手法について述べる.

構造方程式モデル

構造方程式モデルは,データの生成過程を記述するツールである.ここで, x yという二つの変数について,以下のような単純な因果グラフを考える.

f:id:hirotsuru314:20201007131606p:plain
因果グラフのサンプル

このような因果グラフにおける構造方程式を以下のように表す.

 \displaystyle
x=e_x \\
y=f_y(x, e_y) \\

ここで, x yは観測変数で, e_x e_yは未観測の変数であり,誤差変数という. e_y yの値を決定するために寄与しうる x以外の全ての変数をまとめて表した変数という意味をもつ.因果グラフをこのような方程式で表現した場合,分析者は以下の2点に仮定を導入することができる.

  1. 関数 f_yの仮定(例えば線形か非線形か)
  2. 誤差変数 e_x e_yの分布の仮定(例えば,ガウス分布か非ガウス分布か)

上記の二つに仮定を導入するものをパラメトリックアプローチ,関数には仮定をおく一方,誤差変数の分布には仮定をおかないものをセミパラメトリックアプローチと呼ぶ.ここでは,代表的なセミパラメトリックアプローチであるLiNGAMという手法を紹介する.

LiNGAM

LiNGAMはLinear Non-Gaussian Acyclic Modelの略であり,関数形に線形性を,誤差変数の分布に非ガウス分布を仮定するアプローチである.一見,二つの仮定をおいているようにも見えるが,非ガウス分布という仮定はガウス分布でさえなければどんな(連続)分布でも良いということからほぼ仮定をおいていないという解釈ができ,セミパラメトリックアプローチに分類されるようだ.LiNGAMについては,機械学習プロフェッショナルシリーズの「統計的因果探索」に詳しく解説されている.ここではさらっとアルゴリズムの概要をまとめる.
まずLiNGAMにおける仮定を整理する.

  • 関数形が線形
  • 誤差変数が非ガウス連続分布
  • 有向非巡回グラフ(DAG)
  • 誤差変数が独立性(未観測共通原因がないことを意味する)

p個の観測変数 x_1,x_2,...x_pに対するLiNGAMモデルは構造方程式モデルを用いて以下のように表される.

 \displaystyle
x_i=\sum_{i \neq j}^p b_{ij}x_i+e_i

これを行列・ベクトルで表現すると以下のようになる.

 \displaystyle
\mathbf{x}=\mathbf{B}\mathbf{x}+\mathbf{e}

ここで, p次元ベクトル \mathbf{x}は観測変数ベクトル, p× p行列 \mathbf{B}は係数行列, p次元ベクトル \mathbf{e}は誤差変数ベクトルである.そして, \mathbf{e}の成分 e_iは独立で,それぞれ非ガウス連続分布に従う. 係数行列 \mathbf{B}は,観測変数の分布 p(x)に基づいて一意に推定可能であり, \mathbf{B}の成分のゼロ・非ゼロパターンから因果グラフを描くことができる. 係数行列の推定には,独立成分分析(ICA)によるアプローチと回帰分析と独立性評価によるアプローチがある[5]. どちらのアプローチもcdt15/lingamで実装されている.OSSのコードも含めてLiNGAMはもう少し詳細に追いかけたい.

さいごに

今回は,サーベイ論文を参考にグラフィカルモデルを用いた因果探索の手法を整理した.ここでいくつかアルゴリズムを紹介したが,最近機械学習技術を活用した多くの手法が提案されている.例えば,深層強化学習[6],敵対的生成ネットワーク(GAN)[7]を活用したものなどがあり,まだまだ奥が深い.また,重要なイシューとして,確率変数が時系列データの場合にどのように時系列の情報を消失せずにモデリングするかというのも大変興味深い.今後も因果探索について着実に勉強していって,ブログにまとめていきたい.

参考

[1] An Algorithm for Fast Recovery of Sparse Causal Graphs
[2] cuPC: CUDA-Based Parallel PC Algorithm for Causal Structure Learning on GPU
[3] Causal inference in the presence of latent variables and selection bias
[4] Optimal structure identification with greedy search
[5] 統計的因果探索
[6] Causal Discovery with Reinforcement Learning
[7] SAM: Structural Agnostic Model, Causal Discovery and Penalized Adversarial Learning

経路積分量子モンテカルロ法を用いた量子アニーリングのシミュレータ

0または1の状態しかとれないビットを用いた古典的な計算機において量子力学の原理に基づいた量子アニーリングをシミュレートできるのだろうか.
量子アニーリングの理論について学んだ(前回の記事)のちに,このような疑問が生まれた.調べてみると,どうやら手元のPCでシミュレータを作るには経路積分量子モンテカルロ法という手法がよく用いられるようだった. そのため今回は,経路積分量子モンテカルロ法について学び,Pythonでの既存実装に触れてみることで,量子アニーリングについての理解をより深めていくことを目的とする.

目次

経路積分量子モンテカルロ法

今回,古典的な計算機を用いて量子アニーリングをシミュレートするために,経路積分量子モンテカルロ法という手法を用いる. 量子系をシミュレートするために,まず量子ゆらぎを含む横磁場イジングモデルの分配関数 Zを,対応する古典イジングモデルの分配関数で表す. そのとき,これから示す計算を追っていくと,横磁場イジングモデルの分配関数が1次元高い古典の分配関数に対応することがわかる.この結果を使って量子系の平衡状態における性質を古典系のマルコフ連鎖モンテカルロ法で調べる,というのが量子モンテカルロの基本的な考え方である.

以下,理論的な背景を少し詳細に見ていく.
スピン数が Nの場合,横磁場イジングモデルハミルトニアンは以下のように表される.

 \displaystyle
\begin{eqnarray}
\hat{H} & = & \hat{H_0}+\hat{H_q} \\
\hat{H_0} & = & -\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_i^z\hat{\sigma}_j^z-\sum_{i=1}^Nh_i\hat{\sigma}_i^z \\
\hat{H_q} & = & -\Gamma\sum_{i=1}^N\hat{\sigma}_i^x \\
\end{eqnarray}


ここで, \hat{\sigma_i^z}, \hat{\sigma_i^x}はパウリ行列の z,x成分であり,以下の行列で表される.

 \displaystyle
\hat{\sigma_i^z} = 
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix}
,
\hat{\sigma_i^x} = 
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix}


このハミルトニアンを用いて分配関数は以下のように表せる.

 \displaystyle
\begin{eqnarray}
Z & = & Tr[e^{-\beta\hat{H}}] \\
& = & Tr[e^{-\beta(\hat{H_0}+\hat{H_q})}] \\
& = & \sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid{e^{-\beta(\hat{H_0}+\hat{H_q})}}\mid\boldsymbol{\sigma_1}\rangle \\
\end{eqnarray}


ここで, \beta は温度の逆数(逆温度)である.
Suzuki・Trotter分解により,

 \displaystyle
e^{-\beta(\hat{H_0}+\hat{H_q})} = \lim_{M\rightarrow\infty}\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M


と表せることを用いると,

 \displaystyle
Z = \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M\mid\boldsymbol{\sigma_1}\rangle


ここでM個の積の間に以下の恒等演算子を挟み込むと,

 \displaystyle
\boldsymbol{1}  =  \sum_{\boldsymbol{\sigma_k}}\mid\sigma_k\rangle\langle\sigma_k\mid (k = 1, 2, ..., M)


分配関数は以下のように表される.

 \displaystyle
\begin{eqnarray}
Z & = & \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_k},\boldsymbol{\sigma_k'}}\langle\boldsymbol{\sigma_1}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_1'}\rangle\langle\boldsymbol{\sigma_1'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_2}\rangle \\
& \times & \langle\boldsymbol{\sigma_2}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_2'}\rangle\langle\boldsymbol{\sigma_2'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_3}\rangle \\
& \times & ・・・ \\
& \times & \langle\boldsymbol{\sigma_M}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_M'}\rangle\langle\boldsymbol{\sigma_M'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_1}\rangle
\end{eqnarray}


ここで, \mid\boldsymbol{\sigma_k}\rangle N 個のスピン系の直積空間を表す.

 \displaystyle
\mid\boldsymbol{\sigma_k}\rangle=\mid\boldsymbol{\sigma^z_{1,k}}\rangle\otimes\mid\boldsymbol{\sigma^z_{2,k}}\rangle\otimes・・・\otimes\mid\boldsymbol{\sigma^z_{N,k}}\rangle


分配関数の式から, k = 1, 2, ..., M は元の問題の持つ変数が定義された空間とは別の独立な変数のインデックスであり,トロッタ方向または虚時間方向と呼ばれる.
下図はSuzuki・Trotter分解を行なった際にトロッタ方向に次元が追加されている様子を示した概略図である.

f:id:hirotsuru314:20191210223523p:plain:w500
Suzuki-Trotter分解の概略図

これまでの結果から,分配関数を求めるには,以下の2つの値を求めればよいことがわかる.

 \displaystyle
\begin{eqnarray}
\langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle \tag{1} \\
\langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle \tag{2}
\end{eqnarray}


 (1) \hat{H_0}が対角行列であることから,

 \displaystyle
\begin{eqnarray}
\hat{\sigma}_i^z\mid\sigma_k\rangle & = & \hat{\sigma}_{i, k}\mid\sigma_k\rangle \\
\langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle & = & exp\biggl[\frac{\beta}{M}\biggl(\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}+\sum_{i=1}^Nh_i\hat{\sigma}_{i, k}\biggl)\biggl]\prod_{i=1}^N\delta(\hat{\sigma}_{i, k}, \hat{\sigma}_{i, k}')
\end{eqnarray}


となる.次に (2) は以下の式が成立する.

 \displaystyle
\begin{eqnarray}
\langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle = \biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\gamma}{M}\biggl)\biggl]^{\frac{N}{2}}exp\biggl[\frac{1}{2}\log \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1}\biggl)\biggl] \\
\end{eqnarray}


よって,分配関数 Zは以下のように表すことができる.

 \displaystyle
Z = \lim_{M\rightarrow\infty}\biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\Gamma}{M}\biggl)\biggl]^{\frac{NM}{2}}\sum_{\{\hat{\sigma}_{i, k}=\pm1\}}exp\biggl\{\sum_{k=1}^M\biggl[\sum_{i, j=1}^N\frac{\beta J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\frac{\beta h_{i}}{M}\hat{\sigma}_{i, k} \\
+\frac{1}{2}\ln \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}\hat{\sigma}_{i, k+1}\biggl)\biggl] \\


この結果から横磁場イジングモデルの分配関数は対応する磁場なしイジングモデルに1次元トロッタ方向に次元を追加したもので表されることがわかる.対応するハミルトニアンは以下の式で表される.

 \displaystyle
H_{eff} = -\sum_{i, j=1}^N\sum_{k=1}^M\frac{J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\sum_{k=1}^M\frac{h_{i}}{M}\hat{\sigma}_{i, k}-\frac{1}{2\beta}\ln \coth\biggl(\frac{\beta\Gamma}{M}\biggl)\sum_{i=1}^N\sum_{k=1}^M\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1} \\


シミュレーテッドアニーリングとの対比

これまでの計算から,横磁場イジングモデルの分配関数は対応する古典イジングモデルの分配関数に対して虚時間方向に次元を一つ追加したもので表せることがわかった. このことから,量子アニーリングのシミュレータは,一見単にトロッタ数 Mで並列化したシミュレーテッドアニーリングと同じもののように感じるかもしれないが,上記で計算したハミルトニアンを読み解いてみると,そうでないことがわかる.
ハミルトニアン H_{eff}の第1・2項を見ると,同じ層内( kが同一)での相互作用のみであるが,横磁場の効果を表す第三項では,層間の相互作用を持っている.( k k+1の相互作用)このことから時間発展とともに層同士が互いに独立ではなく,相互に作用し合うため,単にトロッタ数分の並列処理をしているわけではないようだ.
どのような相互作用でそれがどう基底状態探索に効いているのかなどについてはまた別途追っていきたい.

計算の流れ

計算の全体的な流れは以下のようになる.

  1. パラメータの初期化
    ・温度  T
    ・横磁場の強さ  \Gamma

  2. スピンの初期化  M層のそれぞれ N個のスピンをランダムに初期化する.

  3. ランダムに層と層中のスピンを選択し,スピンの向きを反転させる.

  4. スピン反転前後のエネルギーの差分を \Delta E とすると,確率 min(1, e^\frac{\Delta E}{T}) でスピンの向きの反転を受け入れる(受け入れない場合は棄却して変更なし)

  5. 任意の回数3,4のステップを繰り返す(モンテカルロステップ)

  6. 横磁場の強さを  \Gamma小さくする.

  7. 任意の回数5,6のステップを繰り返す(アニーリングステップ)

  8. 各層のエネルギーを求め(横磁場の項は除く)一番エネルギーが低い層のスピン配位を解とする

シミュレータを触ってみる

GitHub上ではいくつかの既存実装が存在したが,今回は以下のものを触ってみることにした.

github.com

必要なライブラリ環境をそろえて,リポジトリのトップディレクトリで以下のように実行すると,ヘルプが見れる.ちなみにPython 2系で実装されている.

$ python solve.py qmc -h
usage: Ising Solvers qmc [-h] [-problem PROBLEM] [-steps STEPS] [-dump DUMP]
                         [-spins SPINS] [-P P] [-T T] [-G0 G0] [-Gf GF]
                         [-e0 E0] [-ef EF]

optional arguments:
  -h, --help        show this help message and exit
  -problem PROBLEM  a file containing problem data in three column format, see
                    README for details
  -steps STEPS      number of monte carlo steps (MCS) used, good values depend
                    on problem size
  -dump DUMP        a file to dump move acceptances to, if desired -- slows
                    down performance
  -spins SPINS      initial spin configuration, if desired
  -P P              number of trotter slices (replicas), 20 to 80 is typical
  -T T              ambient temperature in K, 0.015 is typical
  -G0 G0            initial transverse field strength
  -Gf GF            final tranverse field strength, typically near zero
  -e0 E0            initial coupling pre-factor
  -ef EF            final coupling pre-factor

多くの引数があるが,基本的には指定しなくてもデフォルト値が入るようになっていて,必要であれば実行者が上書きできるようにしてある.
マストで渡さないといけない引数は-problem であり,これは,以下のような形式のファイルのパスを渡す必要がある.(sample_dataディレクトリ配下にサンプルが置いてある)

$ head -5 sample_data/ising12.txt
random 12x12 spin ising model with solution energy of -18,972,276
1 2 -28420
2 3 17885
3 4 -146311
4 5 128784

2行以降は以下の形式で

i j J_ij

すなわち,スピン間の相互作用を示す J_{ij}を設定するためのファイルであることがわかる.今回の最適化では実用的な問題が設定してあるのではなく,単にイジングモデルのスピン間にファイルで設定した相互作用を与え,そのイジングモデル基底状態を探索するのみである. また, -dump 引数を指定しておくと,実行の過程がファイル出力される.-steps でアニーリングステップ数が指定できる.
以下のように実行すると,用いたパラメータ等の出力が返ってくる.

$ python solve.py qmc -problem sample_data/ising12.txt -dump output.txt -steps 1000

Ising Spin Glass
+----------------+-------------------------------------------------------------------+
| Field          | Value                                                             |
+----------------+-------------------------------------------------------------------+
| data file      | sample_data/ising12.txt                                           |
| description    | random 12x12 spin ising model with solution energy of -18,972,276 |
| initial config | 0xe696019ff2add3594c219fcbfd9aef9b0dd4                            |
| current config | 0x313fac59f68f16d3ec6a17d6aca964b2025d                            |
| initial energy | 1627456.0000000007                                                |
| current energy | -15899168.000000006                                               |
+----------------+-------------------------------------------------------------------+

Solver: Path-Integral Quantum Monte Carlo
+-----------+------------+
| Parameter | Value      |
+-----------+------------+
| G0        | 3          |
| Gf        | 1e-06      |
| P         | 40         |
| T         | 0.015      |
| dump      | output.txt |
| e0        | 1e-06      |
| ef        | 4          |
| steps     | 1000       |
+-----------+------------+
To replicate starting conditions, run with arguments:
    qmc -dump output.txt -G0 3 -ef 4 -P 40 -steps 1000 -T 0.015 -problem sample_data/ising12.txt -Gf 1e-06 -e0 1e-06 -spins 0xe696019ff2add3594c219fcbfd9aef9b0dd4

出力ファイルを見てみると以下のような形式で出力されている.

$ head output.txt
Path-Integral Quantum Monte Carlo
{'dump': 'output.txt', 'G0': 3, 'ef': 4, 'P': 40, 'steps': 1000, 'T': 0.015, 'Gf': 1e-06, 'e0': 1e-06}
<class 'ising.SpinGlass'>(data_file="sample_data/ising12.txt", spin_configuration="0xe696019ff2add3594c219fcbfd9aef9b0dd4")
3,16.27456,0xe696019ff2add3594c219fcbfd9aef9b0dd4,-1
3.0,15.64962,0xe696019ff2add3594c219fcbfd9aef9b0cd4,7
3.0,13.09112,0xe696019ff2add3594c219fcbfd9aef9b0ddc,28
3.0,13.02468,0xe69e019ff2add3594c219fcbfd9aef9b0dd4,12
3.0,16.31304,0xe696019ff2add3194c219fcbfd9aef9b0dd4,30
3.0,15.9128,0xe696019ff2add3594c219fcbfd9aefbb0dd4,0
3.0,14.19152,0xe697019ff2add3594c219fcbfd9aef9b0dd4,24

4行目以降の結果は,左から横磁場の強さ,全層の中から最も低いエネルギーの値,スピンの状態を16進数で表したもの,何層目のスピンを反転させたのか.
このライブラリの実装では,1アニーリングステップの中で,それぞれの層の中から1スピンずつ反転させている.すなわちモンテカルロステップ数はトロッタ数と同じである.そのモンテカルロステップごとに1行の出力がある.

ステップを重ねるにつれてイジングモデルが低エネルギー状態に向かっていることを確認するために,横軸にアニーリングステップ,縦軸はエネルギーでグラフを書くと以下のようになる.

f:id:hirotsuru314:20191216173856p:plain

アニーリングステップが進むにつれてエネルギーが低い状態に向かっていることから,基底状態を探索している様子が見て取れる.

まとめ

古典的な計算機で量子アニーリングをシミュレートする代表的な手法として,経路積分量子モンテカルロ法がある.経路積分を用いる過程でSuzuki・Trotter展開によって得られるハミルトニアンの分配関数の近似値から,マルコフ連鎖モンテカルロ法を用いることで,量子アニーリングの挙動をシミュレートすることができる. 動作としては,問題とするイジングモデルをトロッタ数分だけコピーして,それらが互いに相互作用しながら基底状態を探索していく. 今回用いたqmcのような既存実装を用いると,実用的な問題ではないものの手元のPCで簡単に量子アニーリングの挙動をシミュレートすることができる.

参考

量子アニーリングの原理について

量子アニーリングという言葉を初めて耳にした時から,ずっと気になっていた.なぜなら学生時代に材料工学を専攻していた私にとって,「量子」と「アニーリング」という言葉はどちらも聞き馴染みのある言葉だったからである. 最近になって,やっと勉強を始めることができ,これから理論的な背景から応用までしっかり学んでいきたいと思っている. せっかくなので勉強した時のメモを少しずつブログにしていく予定である.

目次

量子アニーリングとは

量子アニーリングとは,量子力学の法則を利用して,ある種の難しい問題を解くための枠組みである. その大きな特徴は,組み合わせ最適化問題に対するメタヒューリスティックな解法として利用可能な点である. 組み合わせ最適化問題とは,離散値をとる多くの変数があるときに,それらの目的関数を最小化ないし最大化する問題である. この組み合わせ最適化問題をイジング模型基底状態探索として表現し,量子力学的なゆらぎを利用して探索する方法が量子アニーリングである.

量子アニーリングは,シミュレーテッドアニーリング量子力学版ともとらえることができる.温度の代わりに量子ゆらぎを導入し,古典的な確率過程の代わりに量子力学的な状態の重ね合わせを利用している.
量子アニーリングは1998年に以下の論文で提案されたものである.

arxiv.org

量子アニーリングをハードウェアの動作原理として実装したものには,カナダのベンチャー企業D-Wave Systemsによって開発されたD-Waveマシンなどがある.

イジング模型

量子アニーリングの原理を理解する出発点となるのがイジング模型である. 量子アニーリングでは,最適化問題の目的関数を統計力学のイジング模型で表現し,イジング模型の基底状態(最低エネルギー状態)を探索することで最適解を得る.
イジング模型のハミルトニアンは以下のように表される.

 \displaystyle
\hat{H_0}=-\sum_{i\verb|<|j}J_{ij}\hat{\sigma_i^z}\hat{\sigma_j^z}-\sum_{i=1}^Nh_i\hat{\sigma_i^z}

 \hat{\sigma_i^z}パウリ行列のz成分である. また, iはイジング模型における格子点の番号, J_{ij}は相互作用の強さ, h_iは局所磁場の強さを表している.

イジング模型における格子点は,量子アニーリングマシンにおけるスピンいわゆる量子ビットに対応する.
ここで,量子ビットがもつ上向きスピン +1および下向きスピン -1という状態を以下のようにベクトルで表現する.

 \displaystyle
\mid\uparrow\rangle=
\begin{pmatrix}
1 \\
0 \\
\end{pmatrix}
,
\mid\downarrow\rangle=
\begin{pmatrix}
0 \\
1 \\
\end{pmatrix}

パウリ行列のz成分を上記の2つのベクトルに作用させると,

 \displaystyle
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix}
\begin{pmatrix}
1 \\
0 \\
\end{pmatrix}
=
\begin{pmatrix}
1 \\
0 \\
\end{pmatrix}
,
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix}
\begin{pmatrix}
0 \\
1 \\
\end{pmatrix}
=-
\begin{pmatrix}
0 \\
1 \\
\end{pmatrix}


となるから, \binom{1}{0} \hat{\sigma_i^z}固有値1の固有ベクトル \binom{0}{1}固有値-1の固有ベクトルである. したがって, \hat{\sigma_i^z}を上向きスピンと下向きスピンが同時に存在する量子力学的な重ね合わせの行列表現であると解釈できる.

一般の組み合わせ最適化問題は,多くの場合,前述のようなスピン間の相互作用をもつイジング模型の基底状態探索問題として表現できる. 逆に言うと,基底状態が最適解になるように,解きたい最適化問題に合わせて J_{ij} h_iを設定することができれば,量子アニーリング最適化問題を解くことができる.

量子ゆらぎの導入

量子アニーリングでは,イジング模型の相互作用のパラメータや局所磁場のパラメータのみでなく,量子ゆらぎと呼ばれるものを導入する.  \hat{H_0}基底状態を求めるために,量子力学的なゆらぎを利用して状態探索を行うのが量子アニーリングの基本戦略である.
量子ゆらぎを導入したハミルトニアンを以下の式で表す.

 \displaystyle
\hat{H}(t)=\hat{H_0}+\Gamma(t)\hat{H_1}


 \Gammaは量子ゆらぎの制御変数である.
初期状態 t=0においては, \Gammaを十分大きく選ぶことで右辺第1項を無視でき,第2項の量子ゆらぎの項のみとなる.  \Gammaを時間とともにだんだん小さくしていって,終状態 t=Tにおいては, \Gamma(T)=0となり,最適化したい関数の項のみが残る.
ここで,初期ハミルトニアン \hat{H_1}基底状態が自明に見つかるように \hat{H_1}を選ぶことが重要である. よく用いられる量子ゆらぎとして,横磁場があり,横磁場の効果を示すハミルトニアンを以下のように表す.

 \displaystyle
\hat{H_1}=-\sum_{i=1}^N\hat{\sigma_i^x}

 \hat{\sigma_i^z}はパウリ行列のx成分である.  \hat{\sigma_i^z}固有ベクトルは,

 \displaystyle
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix}
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
=
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
,
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix}
\begin{pmatrix}
1 \\
-1 \\
\end{pmatrix}
=-
\begin{pmatrix}
1 \\
-1 \\
\end{pmatrix}


 \binom{1}{1} \binom{1}{-1}の二つである.
これらはいずれも \binom{1}{0} \binom{0}{1}を同じ絶対値を持つ係数(1と1あるいは1と-1)で加えたものであるから, \hat{\sigma_i^z}は上向きと下向きのスピンが等確率で重ね合わされた状態を表している.

 \hat{H_1}基底状態は, 2^N個の状態が重ね合わされて,量子力学の意味で同時に存在する. 当然ながら,元のイジング模型のハミルトニアン \hat{H_0}基底状態はあらかじめわかっていないため, 2^N個の全ての状態の同じ確率の重ね合わせを出発点として探索を始めることは理にかなっている.

 \hat{H_1}の自明な基底状態を初期状態として選び, \hat{H}(t)で記述される系を時間に依存したシュレディンガー方程式に従って時間発展させる.

 \displaystyle
i\frac{\partial}{\partial t}\mid\Psi(t)\rangle=\hat{H}(t)\mid\Psi(t)\rangle


 \mid\Psi(t)\rangleは,各時刻における系の状態である.
このとき, \Gamma(t)の時間変化が十分ゆっくりであれば,後述する量子力学の断熱定理が適用できる.

断熱変化

ハミルトニアン \hat{H}(t)は時間変数 tに依存するが,系が断熱的に(ゆっくりと)変化していると見なせる場合には,各時刻における系の状態 \mid\Psi(t)\rangleは, tを固定パラメータとみなしたときの定常状態のシュレディンガー方程式

 \displaystyle
\hat{H}(t)\mid\Phi_t\rangle=E_0(t)\mid\Phi_t\rangle


基底状態 \mid\Phi_t\rangleに十分に近い.実際の状態 \mid\Psi(t)\rangleと定常基底状態 \mid\Phi_t\rangleいずれも規格化されているとして,これらの内積が1に十分近いことを意味している.

 \displaystyle
|\langle\Phi_t\mid\Psi(t)\rangle|^2=1-\delta^2\quad(\delta\ll1)


これは各瞬間の基底状態をほぼ正確にたどっていくということを示している.
こうして,横磁場項 \Gamma(t)\hat{H_1}により引き起こされる量子ゆらぎを \Gamma(t)を通じて適切に制御することにより,自明な初期状態から出発して,非自明な状態である最適問題の解に到達する.以上が量子アニーリングの基本的な原理である.

応用

最後に少しだけ量子アニーリングの応用について触れる.最初に述べた通り,量子アニーリングの応用先として期待されているのは組み合わせ最適化問題である.(実際には機械学習のサンプリングにも有効であるそうだ)
一般的に組み合わせ最適化問題の例としてよく登場するものは,巡回セールスマン問題やナップサック問題などが挙げられる.一方,「量子アニーリングの基礎」の中で,クラスタリングが例に挙げられていたのが個人的には興味深かった. クラスタリング組み合わせ最適化問題として考えたことはなかったが,ある点がどのクラスに属するかというのは直感的にも組み合わせであるし,異なるクラスに属する点の間の距離を最大化するという問題に読みかえれば組み合わせ最適化問題の目的関数として定式化できる. この辺りはもっと調べたいので,以下の論文などを読んで学んでいきたい.

arxiv.org

参考

GMOペパボ株式会社を退職しました

2019年7月24日が最終出社日でした。ペパボでは、ロリポップやへテムルといったレンタルサーバサービスのインフラエンジニアとして働いており、約1年半在籍していました。
今回は、ペパボでの一年半の振り返りと、転職に至った経緯などについて書いていきたいと思います。

ペパボでの一年半の振り返り

私はペパボカレッジの枠でペパボに採用されました。ペパボカレッジとは第二新卒エンジニア向け研修のことで、中途採用でも現場に配属される前にしっかり研修を受けることができます。私はエンジニア歴も浅いし、そもそもインフラの経験がなかったため、このペパボカレッジの応募はまたとないチャンスだと感じました。実際に研修では約一ヶ月間で広大なインフラの分野の知識を爆速で学んでいきます。具体的な研修内容などは以前のブログに書きました。頑張りたい人を応援してくれる本当にいい仕組みです。

研修が終わって、配属後も周りのエンジニアの方々に多大なるサポートをいただき、未経験の私でもなんとかインフラエンジニアとして仕事をすることができました。ペパボの「いるだけで成長できる環境」とはこのことだったのかと身に染みて感じました。 どれだけ成長できたかをざっと振り返ってみると、ペパボに入った時にはサーバーに入ってもlscdくらいしか叩けずにSSHって何?といった状態だった私が、入社後半年ちょいでSSHのプロキシサーバを書くまでに至ったので、おそらくかなり成長させてもらったんだろうなーと思います。

また、働く環境も私にとっては最高でした。私のチームは完全フレックスだったのですが、名ばかりのフレックスではなく、本当にフレックスでした。私は子どもがまだ小さいため、早めに出社して早めに退社したり、趣味の筋トレをしてから遅めに出社したりと本当に自分のリズムで働かせていただきました。残業時間も月に20時間を超えたことはありませんでした。スーパーホワイトな会社で心からおすすめできる会社です。実際に私の高校時代の友人をペパボに誘って、入社しました。

転職に至った経緯

以上のように、学ぶ環境も整っていて、生活リズムに合わせて働けるとても良い環境ではありましたが、一方で、ずっと心の中にもどかしい感情がありました。

私はもともとペパボの採用面接のときから、「いずれはインフラの領域に機械学習を活用する研究がしたい」と明言していました。配属後も常にそのことを目標に仕事をしてきました。インフラエンジニアになりたての時はインフラ関連の勉強でいっぱいいっぱいでしたが、最近ではコンテナのリソース割り当ての最適化の論文を読んだり、サーバの負荷状態の把握に機械学習を用いる構想について登壇をしたりなど、やりたいことの妄想は膨らんでいました。しかしながら、なかなか妄想を実行に移すことができませんでした。それは、私自身の力不足はもちろんですが、ロリポップのように長くやっているサービスの運用は非常に大変であることも要因の一つとしてあり、結局一年半で一度も機械学習を絡めた仕事をすることができませんでした。

また仕事では、私自身が自分の強みだと思っている「深く思考すること」がなかなかできず、リアクティブに急ぎの対応するような仕事が占める割合が非常に多いのが現状でした。もちろん、自分が優秀であれば、与えられた仕事をこなしながらも、やりたいことができるかもしれませんが、私にはできませんでした。そんな時に消防士からエンジニアに夢を持って転職した自分の感情が蘇り、本当にこのままでは何者にもなれないまま終えてしまいそうにだという焦りが抑えきれなくなり、転職を決意しました。

研究者の道を選んだ

2019年8月1日からさくらインターネット研究所の研究員として働かせていただくこととなりました。
正直、このタイミングでエンジニアから研究者になることに関してはかなり悩みました。それは、私が固定観念で「エンジニアとして優秀じゃないと研究者として活躍できない」と思っていたからです。
この固定観念による不安に対する一つの答えを与えてくれたのは、@matsumotoryさんが書いた下のブログでした。研究をいずれやりたいと思っているなら今すぐ研究に飛び込んだ方がいいという考えに変わるきっかけとなりました。

hb.matsumoto-r.jp

実は松本さんには、もう一年以上前に松本さんがペパボに在籍されていた時から、エンジニアとしての自分の理想と現実が乖離していることや、いずれ研究をやりたいけどまだやっていける自信が持てないことなど色々相談させていただいていました。会議室で二時間近く二人でお話させていただいたこともありました。
その時に松本さんがおっしゃった言葉のいくつかは時間をおいてからやっとその意味が自分の中で消化できて、私の人生に大きく影響を与えた言葉となりました。松本さんはペパボ時代から私の学生時代の研究(専攻は材料工学でした)の業績をしっかり見て下さっていて「過去にこれだけ研究で活躍できていた人間が分野が違うからできないなんてことはありえない」と言ってくれていました。その言葉が最終的に私に研究に飛び込む勇気を与えてくれて、研究への挑戦を心に決めました。

そんな中、現在さくらインターネット研究所の研究員をされている@yuuk1tさんや、@kumagalliumさんとお会いする機会に恵まれ、色々と研究のお話をさせていただく中で、研究に対する熱量や考え方に共感し「私も一緒に研究をやりたい」という気持ちが高まって、さくらインターネット研究所に応募させていただきました。

さいごに

お世話になったペパボの皆さんには「辞めてすみませんでした」ではなく「教えてもらったことを活かしてこれからさらに飛躍します」とポジティブな気持ちでいます!
一年半の間本当にありがとうございました!

遺伝的アルゴリズムの並列化とgoroutineによる実装

先日、「遺伝的アルゴリズムをGoで実装してみた」という記事を書きました。
この内容で2019年7月13日(土)に開催されるGo Conference'19 in Fukuokaに登壇させていただくことになったので、開発中のeagoというパッケージをもっと作り込んで行きたいと思います。
今回は、遺伝的アルゴリズム(以下、GA)の計算処理をgoroutineでサクッと並列化した話です。

概要

GAは複数の個体を用いた多点探索のアルゴリズムであるため、本質的に並列化と親和性が高い手法であると言えます。実用的な観点からも短い計算時間で良好な解を得ることが望まれているため、これまでに数多くの並列化手法が提案されています。

今回はGo言語の強みの一つであるgoroutineを使ってGAの処理の一部を並列化してみます。また、擬似問題に対してベンチマークを取り、実装した並列処理が全体の計算時間短縮に有効であることも検証しました。

代表的な並列化手法

マスタースレーブ方式

GAを並列化するための最も単純な手法は適応度計算を並列化することです。マスタースレーブ方式では下図のように、全体の制御および遺伝的操作を行うマスターと、適応度計算を行う複数のスレーブから構成されます。

f:id:hirotsuru314:20190619141143p:plain
マスタースレーブ方式

この方式は、同時に計算できる処理を並列化しただけの極めて単純かつ一般的な手法であり、GA固有の並列化処理ではありません。
今回はこちらをgoroutineを使って実装してみます。

島モデル

島モデルはGA自体を並列化するというアプローチをとっています。集団を複数の部分集団に分割することでGAを同時に実行し、並列処理を実現しています。
その概要を下図に示します。

f:id:hirotsuru314:20190619142435p:plain
島モデル

それぞれの部分集合を初期化して、評価・選択・交叉・突然変異などをそれぞれで行います。ここで本手法の重要な特徴となるのが、ある一定の条件を元に部分集団から一部の個体を選んで交換を行う点です。このプロセスにより、島モデルは並列化手法であるという側面のほかに、個体の多様性を維持するための手法としての側面も併せもちます。

goroutineによる並列化

今回はマスタースレーブ方式を実装してみます。前述の通り、マスタースレーブ方式では適応度計算を並列化するため、適応度関数(目的関数)が非常に複雑で計算処理に時間がかかる場合に有効な並列化であると言えます。
個体群に属するそれぞれの個体の適応度計算は独立して行えるため、この適応度計算をするためのスレーブとしてgoroutineを立ち上げて、並列で計算を行うようにします。マスターについては、それぞれのスレーブの計算結果を受け取るプロセスがいるというよりは、それぞれのスレーブが以下のように定義したIndividual構造体のFitnessに計算結果を書き込むことで、以降の処理で全ての結果を参照できるようにしています。

type Individuals []Individual

type Individual struct {
    Chromosome Genome
    Fitness    float64
    Evaluated  bool
}

実装にはsync.WaitGroupを用いました。
WaitGroupを用いると、独立した複数のタスクを平行で処理して、それら全ての終了を待ち合わせる処理が非常にシンプルな記述で実装できます。

func (indis Individuals) Evaluate(parallel bool) {
    var wg sync.WaitGroup
    if !parallel {
        for i := range indis {
            indis[i].Evaluate()
        }
    } else {
        for i := range indis {
            wg.Add(1)
            go func(i int) {
                defer wg.Done()
                indis[i].Evaluate()
            }(i)
        }
        wg.Wait()
    }
}

上記の Evaluate() の部分が適応度計算の関数です。
WaitGroupはカウンターのようになっていて、Addを呼び出すと引数に渡された整数だけカウンターを増やし、Doneを呼び出すとカウンターを1つ減らします。そしてWaitを呼び出すとカウンターがゼロになるまでブロックします。

パフォーマンス検証

実行時間の測定には、testingパッケージに備わっているベンチマークの仕組みを用いました。
並列化が有効かどうかを検証するために、実用的ではないですが、適応度関数の中で1 ms sleepするようにして、擬似的に計算処理に時間がかかる適応度関数を作りました。

ベンチマークのための全体コードの以下のようになります。

package main

import (
    "github.com/tsurubee/eago"
    "log"
    "testing"
    "time"
)

type Vector []float64

func (V Vector) Initialization() eago.Genome {
    return Vector(eago.InitFloatVector(2, 32, -32))
}

func (V Vector) Fitness() float64 {
    time.Sleep(1 * time.Millisecond) // Heavy process
    var s float64
    for _, x := range V {
        s += x * x
    }
    return s
}

func (V Vector) Mutation() {
    eago.AddNormalFloat(V, 0.5)
}

func (V Vector) Crossover(X eago.Genome) eago.Genome {
    return Vector(eago.BLXalpha(V, X.(Vector), 0.3))
}

func BenchmarkParallelGA(b *testing.B) {
    for i := 0; i < b.N; i++ {
        runGA(true)
    }
}

func BenchmarkNonParallelGA(b *testing.B) {
    for i := 0; i < b.N; i++ {
        runGA(false)
    }
}

func runGA(parallel bool) {
    var v Vector
    ga := eago.NewGA(eago.GAConfig{
        PopulationSize: 30,
        NGenerations:   20,
        CrossoverRate:  0.8,
        MutationRate:   0.01,
        ParallelEval:   parallel,
    })
    ga.PrintCallBack = func() {} // Do not print messages while running
    if err := ga.Minimize(v); err != nil {
        log.Fatal(err)
    }
}

このファイルを *_test.goという形式で保存して、ファイルが存在するディレクトリで以下のようなコマンドを実行することでベンチマークが取得できます。

$ go test -bench .
goos: darwin
goarch: amd64
pkg: github.com/tsurubee/eago/examples/ga
BenchmarkParallelGA-4             50      37193820 ns/op
BenchmarkNonParallelGA-4           2     821147753 ns/op
PASS
ok      github.com/tsurubee/eago/examples/ga    4.365s

適応度関数内でのsleepの時間を1・5・10 msにして、並列化あり(parallel)並列化なし(non-parallel)における処理時間をグラフにすると以下のようになりました。

f:id:hirotsuru314:20190623114700p:plain

今回のベンチマークの条件では、個体数30、世代数20であるため、適応度計算は600回程度行われ、sleepが10 msのときに処理を並列化しない場合は、適応度計算だけでも6000 msかかることになります。
このことを考えると、適応度計算が重い場合に、今回のようなマスタースレーブ方式の並列化が、全体の処理時間の短縮に有効であることがわかります。

今後

Go Conference'19 in Fukuokaまでに以下の2点について取り組んで行きたいと思います。

  • 他の並列化手法の検証(島モデルなど)

  • 実用的な問題への適用

参考

www.morikita.co.jp

www.oreilly.co.jp

遺伝的アルゴリズムをGoで実装してみた

こんにちは!つるべーです!
最近は、進化計算と呼ばれるバイオミメティックな計算技法に興味を持っており、実装しながら勉強しています。
前回の記事では粒子群最適化(Particle Swarm Optimization: PSO)を実装しました
今回は、遺伝的アルゴリズムという生物の進化の過程を模倣して作られた最適解探索アルゴリズムをGoでスクラッチで書いてみました。

遺伝的アルゴリズム(Genetic Algorithm: GA)

GAとは何か?という問いに対する答えとして、以下の一文が端的にその特徴を表していると思います。

生物進化における遺伝と適者生存による自然淘汰の仕組みをソフトウェア的に模すことで複雑な問題に対する最適解を探索する手法
引用:「遺伝的アルゴリズム(Genetic Algorithm)を始めよう!」

GAでは、問題に対する解を個体の染色体、解の構成要素を遺伝子として表現し、複数の個体を進化させて、より環境に適応できる個体を見つけることで最適解を導いていきます。その過程をフローチャートとしてまとめたものが以下になります。

f:id:hirotsuru314:20190503224038p:plain

コード

github.com

eago(イーゴ)という進化計算パッケージを鋭意開発中です!
ちなみに、Pythonで書かれたdeapという素晴らしいライブラリがあり、遺伝的アルゴリズムを始めとした様々な進化計算の手法が実装されています。

最適化問題を解いてみる

今回はAckley functionと呼ばれる関数の最小化問題を解いてみたいと思います。数式およびグラフは以下の通りです。グラフは左側が3Dに描いたもので、右側が2次元に写像したもの(左のグラフを上から見たもの)です。

f:id:hirotsuru314:20190504145941p:plain

今回は n=2の条件を採用しています。
Ackley functionは多峰性関数であり、多くの局所解をもつため、最適化アルゴリズムベンチマーク関数として用いられることもあるようです。下のサイトに様々なベンチマーク関数が載っていて眺めるだけでも面白かったです。

Test functions for optimization - Wikipedia

グラフからわかるように、この関数は (x_1,x_2)=(0,0)で最小値 0を取るため、今回の探索の目的は、 (x_1,x_2)=(0,0)を見つけることになります。

サンプルコードは以下の通りです。

package main

import (
    "github.com/tsurubee/eago"
    "log"
    "math"
)

type Variables []float64

func (V Variables) Initialization() eago.Genome {
    return Variables(eago.InitFloatVector(2, 32, -32))
}

func (V Variables) Fitness() float64 {
    return Ackley(V)
}

func (V Variables) Mutation() {
    eago.AddNormalFloat(V, 0.5)
}

func (V Variables) Crossover(X eago.Genome) eago.Genome {
    return Variables(eago.BLXalpha(V, X.(Variables), 0.3))
}

func Ackley(X []float64) float64 {
    a, b, c, d := 20.0, 0.2, 2*math.Pi, float64(len(X))
    var s1, s2 float64
    for _, x := range X {
        s1 += x * x
        s2 += math.Cos(c * x)
    }
    return -a*math.Exp(-b*math.Sqrt(s1/d)) - math.Exp(s2/d) + a + math.Exp(1)
}

func main() {
    var v Variables
    ga := eago.NewGA(eago.GAConfig{
        PopulationSize: 30,
        NGenerations:   20,
        CrossoverRate:  0.8,
        MutationRate:   0.01,
    })
    if err := ga.Minimize(v); err != nil {
        log.Fatal(err)
    }
}

コードの詳細については最後に補足します。
上のコードをmain.goとして以下のように実行します。

$ go run main.go
Generation   1: Fitness=11.208 Solution=[-4.501 0.087]
Generation   2: Fitness=10.001 Solution=[-3.592 -0.833]
Generation   3: Fitness=9.165 Solution=[-3.841 -0.960]
Generation   4: Fitness=5.747 Solution=[-1.426 -0.235]
Generation   5: Fitness=4.258 Solution=[-0.824 0.380]
Generation   6: Fitness=2.797 Solution=[-0.843 -0.011]
Generation   7: Fitness=2.605 Solution=[-0.955 -0.030]
Generation   8: Fitness=2.605 Solution=[-0.955 -0.030]
Generation   9: Fitness=2.604 Solution=[-0.925 0.016]
Generation  10: Fitness=2.604 Solution=[-0.925 0.016]
Generation  11: Fitness=2.479 Solution=[-0.346 -0.054]
Generation  12: Fitness=1.216 Solution=[-0.184 -0.006]
Generation  13: Fitness=0.615 Solution=[-0.111 -0.001]
Generation  14: Fitness=0.501 Solution=[-0.096 -0.007]
Generation  15: Fitness=0.267 Solution=[-0.055 0.024]
Generation  16: Fitness=0.056 Solution=[0.010 -0.014]
Generation  17: Fitness=0.056 Solution=[0.010 -0.014]
Generation  18: Fitness=0.053 Solution=[-0.005 -0.015]
Generation  19: Fitness=0.031 Solution=[0.010 0.003]
Generation  20: Fitness=0.013 Solution=[-0.000 0.004]

出力結果を見てみると、世代数が大きくなるにつれて、Fitnessが 0、パラメータの解が (x_1,x_2)=(0,0)に収束しようとしているため、うまく探索できていることがわかります。

最後に、今回用いたコードとアルゴリズムの補足をしておきます。

実装とアルゴリズムの補足

Genomeインターフェース

実装では、下記のようなGenomeインターフェースを定義しており、上記のサンプルコードのVariablesはGenomeインターフェースを満たすように各メソッドを実装しています。

type Genome interface {
    Initialization() Genome
    Fitness() float64
    Mutation()
    Crossover(Genome) Genome
}

このように実装した理由は、eagoを使う側で自由に独自に実装した選択や交叉などの手法を組み込めるようにするためです。できればメジャーな手法については全てeago側に実装して、特殊なケース以外は全てeagoに実装されたメソッドを呼び出すだけで事足りるようにしていきたいと思っています。

実数値遺伝的アルゴリズム(Real coded GA)

今回の最適化では、解表現に浮動小数点数を用いており、実数値GAと呼ばれる手法を実装しました。
GAにおいては、各遺伝子を0か1で表現したバイナリエンコーディングが最も一般的な解表現だと思います。このような解の表現の仕方により、探索過程の各ステップ(選択・交叉・突然変異)で用いる手法も大きく異なってきます。
それぞれのステップにおいて、今回用いた手法を説明します。

選択

今回用いている選択の手法は、トーナメント法と呼ばれるものです。
任意のn個の個体を組にして、それらの中で最も適応度が高いものを選択するという手法になります。
他の手法として、ルーレット法やランキング法などがあります。

交叉

交叉にはブレンド交叉(BLX-α)という実数値GA特有の手法を用いました。
BLX-αについては下のサイトがわかりやすく解説されていました。

実数型GAに於ける交叉法の改良

親となる2個体によってできる長方形を一定数倍引き伸ばした長方形で囲まれる範囲内でランダムに次世代の個体が生成されます。

突然変異

突然変異として、個体の各遺伝子に任意の確率で一様乱数を加算するようなものを実装しました。一般的に突然変異は、局所解に陥ることを防ぐために取り入れられます。

さいごに

GAは、解の表現の仕方や、選択・交叉・突然変異など各遺伝子操作に様々な手法が提案されていて、何を使うかは問題を解く側に委ねられており、以前実装したPSOと比べて、非常に自由度が高いという印象を受けました。なので、GAの理論や実装をまだまだ深掘りしながら楽しめそうです!
今後は、処理の並列化や多目的最適化への拡張などを中心に勉強していきたいと思っています。

粒子群最適化(Particle Swarm Optimization: PSO)をGoで実装してみた

粒子群最適化とは、群知能による最適化手法の一種です。この手のバイオミメティクス(生物模倣)によるアプローチは、私の学生時代の専門である材料工学でも非常に盛んに研究されていましたが、データサイエンスでも応用されているのを知って興味を持ったので、実装しながら勉強しました。

粒子群最適化(Particle Swarm Optimazation: PSO

粒子群最適化とは、鳥や魚などの群れに見られる社会的行動のシミュレーションを基にモデル化されたヒューリスティックな最適解探索アルゴリズムです。

全体的な最適化の流れを説明します。 まず最初に探索空間内で、決められた個数の粒子を全て初期化します。初期化とは空間内の位置をランダムに決めることを指します。

それぞれの粒子は、目的関数により適合度(どれだけその解が優れているか)を算出できます。全ての粒子をある移動ルールを元に動かしながら、適合度が高い位置を探していくのがPSOの探索方法です。

ここでPSOにおいて重要な移動ルールについて説明します。

粒子 P_{i}の時刻 tから t+1の移動速度 \vec{v}_{i}(t+1)は以下の式で表されます。

f:id:hirotsuru314:20190409205526p:plain:w500

ここで重要なのが、グローバルベスト \vec{g}(t)とパーソナルベスト \vec{p}_{i}(t)です。
グローバルベストとは全ての粒子群で過去の探索も含めて最も適合度が高い位置、パーソナルベストとは特定の1粒子におけるこれまでの探索の中で最も適合度が高い位置を示します。これは繰り返しの探索の中で随時更新されていく変数です。
 Iは慣性係数と呼ばれる定数で、前の移動速度がどれだけ維持されるかを表します。 A_{g}および A_{p}は加速係数と呼ばれる定数です。

この数式からわかるように粒子群の中ではグローバルベストの情報が共有されていて、その方向に向かうように移動速度を調整しながら、より良い解の探索を繰り返します。これはあたかも餌を見つけた鳥が仲間にその情報を伝えて、仲間が餌のある位置に集まってくるようです。

作ったもの

github.com

『Evolutionary Algorithm implemented in Go』の略で名前を付けました。 進化計算全般に興味があるので、ここに実装しながら勉強したいと思っています。

最適化問題を解いてみる

今回は下のシンプルな目的関数の最小値を最適解として、PSOでパラメータを探索してみました。

f:id:hirotsuru314:20190409212421p:plain:w120

 N=2としたときの関数のグラフを描くと以下のようになります。

f:id:hirotsuru314:20190409215708p:plain

左側が3Dに描いたもので、右側が2次元に写像したもの(左のグラフを上から見たもの)です。
グラフや数式から直感的にわかるように、この関数は (x_1,x_2)=(0,0)で最小値 0を取るため、今回の探索の目的は、 (x_1,x_2)=(0,0)を見つけることになります。

サンプルコードは以下の通りです。

package main

import (
    "github.com/tsurubee/eago"
    "log"
)

func objectiveFunc(x []float64) float64 {
    return x[0] * x[0] + x[1] * x[1]
}

func main() {
    pso := eago.NewDefaultPSO()
    pso.NParticle =  5
    pso.NStep = 20
    pso.Min = -20
    pso.Max = 10

    if err := pso.Minimize(objectiveFunc, 2); err != nil {
        log.Fatal(err)
    }
}

上のコードをmain.goとして以下のように実行します。

$ go run main.go
Step   0: Fitness=8.124 Position=[-2.168 1.851]
Step   1: Fitness=8.124 Position=[-2.168 1.851]
Step   2: Fitness=7.394 Position=[-1.618 2.186]
Step   3: Fitness=6.958 Position=[1.068 2.412]
Step   4: Fitness=6.958 Position=[1.068 2.412]
Step   5: Fitness=4.809 Position=[-2.106 0.611]
Step   6: Fitness=2.215 Position=[-1.052 1.053]
Step   7: Fitness=1.535 Position=[0.580 1.095]
Step   8: Fitness=0.163 Position=[0.023 0.403]
Step   9: Fitness=0.163 Position=[0.023 0.403]
Step  10: Fitness=0.055 Position=[-0.175 0.156]
Step  11: Fitness=0.055 Position=[-0.175 0.156]
Step  12: Fitness=0.055 Position=[-0.175 0.156]
Step  13: Fitness=0.055 Position=[-0.175 0.156]
Step  14: Fitness=0.055 Position=[-0.175 0.156]
Step  15: Fitness=0.055 Position=[-0.175 0.156]
Step  16: Fitness=0.055 Position=[-0.180 0.151]
Step  17: Fitness=0.044 Position=[0.019 0.210]
Step  18: Fitness=0.013 Position=[-0.111 0.023]
Step  19: Fitness=0.008 Position=[-0.083 -0.031]

出力結果を見てみると、Stepが進むにつれて、Fitnessが 0、パラメータの解が (x_1,x_2)=(0,0)に収束しようとしていることがわかります。
粒子が速度を調整しながら最適解に収束して行く様子は下のサイトのアニメーションが非常にわかりやすかったです。

Particle swarm optimization - Wikiwand

さいごに

群知能とか生物模倣的なアプローチはもう少し勉強を続けて深掘りしていきたい。
あと、もっと実践的なデータに対して適用してみて、自分の今の領域における応用についても考えていきたい。

参考

shop.ohmsha.co.jp