交換モンテカルロ法 exchange

exchange はレプリカ交換モンテカルロ法を用いてパラメータ探索を行う Algorithm です。

前準備

あらかじめ mpi4py をインストールしておく必要があります。:

python3 -m pip install mpi4py

入力パラメータ

サブセクション param と exchange を持ちます。

[param] セクション

探索空間を定義します。 mesh_path キーが存在する場合は離散空間を、そうでない場合は連続空間を探索します。

  • 連続空間

    • initial_list

      形式: 実数型のリスト。長さはdimensionの値と一致させます。

      説明: パラメータの初期値。 定義しなかった場合は一様ランダムに初期化されます。

    • unit_list

      形式: 実数型のリスト。長さはdimensionの値と一致させます。

      説明: 各パラメータの単位。

      探索アルゴリズム中では、各パラメータをそれぞれこれらの値で割ることで、 簡易的な無次元化・正規化を行います。 定義しなかった場合にはすべての次元で 1.0 となります。

    • min_list

      形式: 実数型のリスト。長さはdimensionの値と一致させます。

      説明: パラメータが取りうる最小値。

      モンテカルロ探索中にこの値を下回るパラメータが出現した場合、 ソルバーは評価されずに、値が無限大だとみなされます。

    • max_list

      形式: 実数型のリスト。長さはdimensionの値と一致させます。

      説明: パラメータが取りうる最大値。

      モンテカルロ探索中にこの値を上回るパラメータが出現した場合、 ソルバーは評価されずに、値が無限大だとみなされます。

  • 離散空間

    • mesh_path

      形式: ファイルパス

      説明: メッシュ定義ファイル。

    • neighborlist_path

      形式: ファイルパス

      説明: 近傍リスト定義ファイル。

[exchange] セクション

  • numsteps

    形式: 整数値。

    説明: モンテカルロ更新を行う回数。

  • numsteps_exchange

    形式: 整数値。

    説明: 「温度」のレプリカ交換を行う頻度。この回数だけモンテカルロ更新を行ったらレプリカ交換を実行します。

  • Tmin

    形式: 実数値。

    説明: 「温度」(\(T\))の最小値。

  • Tmax

    形式: 実数値。

    説明: 「温度」(\(T\))の最大値。

  • bmin

    形式: 実数値。

    説明: 「逆温度」(\(\beta = 1/T\))の最小値。 温度と逆温度はどちらか片方だけを指定する必要があります。

  • bmax

    形式: 実数値。

    説明: 「逆温度」(\(\beta = 1/T\))の最大値。 温度と逆温度はどちらか片方だけを指定する必要があります。

  • Tlogspace

    形式: 真偽値。 (default: true)

    説明: 「温度」を各レプリカに割り当てる際に、対数空間で等分割するか否かを指定します。

    true のときは対数空間で等分割します。

  • nreplica_per_proc

    形式: 整数。 (default: 1)

    説明: ひとつのMPI プロセスが担当するレプリカの数。

アルゴリズム補助ファイル

メッシュ定義ファイル

本ファイルで探索するグリッド空間を定義します。 1列目にメッシュのインデックス（実際には使用されません）、 2列目以降は探索空間の座標を指定します。

以下、サンプルを記載します。

1 6.000000 6.000000
2 6.000000 5.750000
3 6.000000 5.500000
4 6.000000 5.250000
5 6.000000 5.000000
6 6.000000 4.750000
7 6.000000 4.500000
8 6.000000 4.250000
9 6.000000 4.000000
...

近傍リスト定義ファイル

離散空間をモンテカルロ法で探索する場合、各点 \(i\) ごとに次に移動できる点 \(j\) を定めておく必要があります。 そのために必要なのが近傍リスト定義ファイルです。

1列目に始点の番号 \(i\) を記載し、 2列目以降に \(i\) から移動できる終点 \(j\) を列挙します。

近傍リスト定義ファイルをメッシュ定義ファイルから生成するツール py2dmat_neighborlist が提供されています。 詳細は 関連ツール を参照してください。

0 1 2 3
1 0 2 3 4
2 0 1 3 4 5
3 0 1 2 4 5 6 7
4 1 2 3 5 6 7 8
5 2 3 4 7 8 9
...

出力ファイル

RANK/trial.txt

各レプリカについて、モンテカルロサンプリングで提案されたパラメータと、対応する目的関数の値です。 1列目にステップ数、2列目にプロセス内のwalker 番号、3列目にレプリカの温度、4列目に目的関数の値、5列目以降にパラメータが記載されます。

# step walker T fx z1 z2
0 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
1 0 0.004999999999999999 0.0758494287185766 2.811346329442423 3.691101784194861
2 0 0.004999999999999999 0.08566823949124412 3.606664760390988 3.2093903670436497
3 0 0.004999999999999999 0.06273922648753057 4.330900869594549 4.311333132184154

RANK/result.txt

各レプリカについて、モンテカルロサンプリングで生成されたパラメータと、対応する目的関数の値です。 trial.txt と同一の書式です。

# step walker T fx z1 z2
0 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
1 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
2 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
3 0 0.004999999999999999 0.06273922648753057 4.330900869594549 4.311333132184154

best_result.txt

サンプリングされた全データのうち、目的関数の値が最小となったパラメータと、対応する目的関数の値です。

nprocs = 4
rank = 2
step = 65
fx = 0.008233957976993406
z1 = 4.221129370933539
z2 = 5.139591716517661

result_T#.txt

サンプリング結果を温度ごとにまとめ直したものです。 # は温度点の番号です。 ファイルの1 列目はステップ数、2列目は全体での walker 番号、3列目は目的関数の値、 4列目以降は探索パラメータの値です。

# T = 1.0
0 15 28.70157662892569 3.3139009347685118 -4.20946994566609
1 15 28.70157662892569 3.3139009347685118 -4.20946994566609
2 15 28.70157662892569 3.3139009347685118 -4.20946994566609
3 15 28.98676409223712 3.7442621319489637 -3.868754990884034

アルゴリズム解説

マルコフ連鎖モンテカルロ法

モンテカルロ法（モンテカルロサンプリング）では、パラメータ空間中を動き回る walker \(\vec{x}\) を重み \(W(\vec{x})\) に従って確率的に動かすことで目的関数の最適化を行います。 重み \(W(\vec{x})\) として、「温度」 \(T > 0\) を導入して \(W(\vec{x}) = e^{-f(\vec{x})/T}\) とすることが一般的です（ボルツマン重み）。 ほとんどの場合において、 \(W\) に基づいて直接サンプリングする (walker を生成する) のは不可能なので、 walker を確率的に少しずつ動かすことで、頻度分布が \(W\) に従うように時系列 \(\{\vec{x}_t\}\) を生成します (マルコフ連鎖モンテカルロ法, MCMC)。 \(\vec{x}\) から \(\vec{x}'\) へ遷移する確率を \(p(\vec{x}' | \vec{x})\) とすると、

\[W(\vec{x}') = \sum_{\vec{x}} p(\vec{x}' | \vec{x}) W(\vec{x})\]

となるように \(p\) を定めれば時系列 \(\{\vec{x}_t\}\) の頻度分布が \(W(\vec{x})\) に収束することが示されます（釣り合い条件） [1] 。 実際の計算では、より強い制約である詳細釣り合い条件

\[p(\vec{x} | \vec{x}') W(\vec{x}') = W(\vec{x})p(\vec{x}' | \vec{x})\]

を課すことがほとんどです。 両辺で \(vec{x}\) についての和を取ると釣り合い条件に帰着します。

\(p\) を求めるアルゴリズムはいくつか提案されていますが、 2DMAT では Metropolis-Hasting 法 (MH法) を用います。 MH 法では、遷移プロセスを提案プロセスと採択プロセスとに分割します。

  1. 提案確率 \(P(\vec{x} | \vec{x}_t)\) で候補点 \(\vec{x}\) を生成します

    • 提案確率 \(P\) としては \(\vec{x}_t\) を中心とした一様分布やガウス関数などの扱いやすいものを利用します

  2. 提案された候補点 \(\vec{x}\) を採択確率 \(Q(\vec{x}, | \vec{x}_t)\) で受け入れ、 \(\vec{x}_{t+1} = \vec{x}\) とします - 受け入れなかった場合は \(\vec{x}_{t+1} = \vec{x}_t\) とします

採択確率 \(Q(\vec{x} | \vec{x}_t)\) は

\[Q(\vec{x} | \vec{x}_t) = \min\left[1, \frac{W(\vec{x})P(\vec{x}_t | \vec{x}) }{W(\vec{x}_t) P(\vec{x} | \vec{x}_t)} \right]\]

とします。 この定義が詳細釣り合い条件を満たすことは、詳細釣り合いの式に代入することで簡単に確かめられます。 特に、重みとしてボルツマン因子を、 提案確率として対称なもの \(P(\vec{x} | \vec{x}_t) = P(\vec{x}_t | \vec{x})\) を用いたときには、

\[Q(\vec{x} | \vec{x}_t) = \min\left[1, \frac{W(\vec{x})}{W(\vec{x}_t)} \right] = \min\left[1, \exp\left(-\frac{f(\vec{x}) - f(\vec{x}_t)}{T}\right) \right]\]

という更に簡単な形になります。

\(\Delta f = f(\vec{x}) - f(\vec{x}_t)\) とおいて、 \(\Delta f \le 0\) のときに \(Q = 1\) となることを踏まえると、 MH 法によるMCMC は次のようになります。

  1. 現在地点の近くからランダムに次の座標の候補を選び、目的関数 \(f\) の値を調べる

  2. \(\Delta f \le 0\) ならば（山を下る方向ならば）移動する

  3. \(\Delta f > 0\) ならば採択確率 \(Q = e^{-\Delta f / T}\) で移動する

  4. 1-3 を適当な回数繰り返す

得られた時系列のうち、目的関数の値が一番小さいものを最適解とします。 3 番のプロセスのおかげで、 \(\Delta f \sim T\) ぐらいの山を乗り越えられるので、局所最適解にトラップされた場合にも脱出可能です。

レプリカ交換モンテカルロ法

モンテカルロ法による最適化では、温度 \(T\) は非常に重要なハイパーパラメータとなっています。 モンテカルロ法では、温度 \(T\) 程度の山を乗り越えられますが、逆にそれ以上の深さの谷からは容易に脱出できません。 そのため、局所解へのトラップを防ぐためには温度を上げる必要があります。 一方で、 \(T\) よりも小さい谷は谷として見えなくなるため、得られる \(\min f(\vec{x})\) の精度も \(T\) 程度になり、精度を上げるためには温度を下げる必要があります。 ここから、最適解を探すためには温度 \(T\) を注意深く決める必要があることがわかります。

この問題を解決する方法として、温度 \(T\) を固定せずに更新していくというものがあります。 たとえば、焼きなまし法 (simulated annealing) では、温度をステップごとに徐々に下げていきます。 焼戻し法 (simulated tempering) は、温度をハイパーパラメータではなく、サンプリングすべきパラメータとして扱い、（詳細）釣り合い条件を満たすように更新することで、加熱と冷却を実現します。温度を下げることで谷の詳細を調べ、温度を上げることで谷から脱出します。 レプリカ交換モンテカルロ法 (replica exchange Monte Carlo) は焼戻し法を更に発展させた手法で、並列焼戻し法 (parallel tempering) とも呼ばれます。 レプリカ交換モンテカルロ法では、レプリカと呼ばれる複数の系を、それぞれ異なる温度で並列にモンテカルロシミュレーションします。 そして、ある一定間隔で、（詳細）釣り合い条件を満たすように他のレプリカと温度を交換します。 焼戻し法と同様に、温度を上下することで谷を調べたり脱出したりするのですが、各温度点について、かならずレプリカのどれかが対応しているため、全体として特定の温度に偏ることがなくなります。 また、複数の MPI プロセスを用意してそれぞれレプリカを担当させることで簡単に並列化可能です。 数多くのレプリカを用意することで温度間隔が狭まると、温度交換の採択率も上がるため、大規模並列計算に特に向いたアルゴリズムです。 有限温度由来の「ぼやけ」がどうしても生まれるので、モンテカルロ法の結果を初期値として minsearch をするのがおすすめです。

脚注