ポピュレーションアニーリングモンテカルロ法 `pamc`¶

pamc はポピュレーションアニーリングモンテカルロ法を用いてパラメータ探索を行う Algorithm です。

前準備¶

MPI 並列をする場合にはあらかじめ mpi4py をインストールしておく必要があります。:

python3 -m pip install mpi4py

入力パラメータ¶

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

[`param`] セクション¶

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

連続空間
- initial_list
  
  形式: 実数型のリスト。長さはdimensionの値と一致させます。
  
  説明: パラメータの初期値。定義しなかった場合は一様ランダムに初期化されます。
- unit_list
  
  形式: 実数型のリスト。長さはdimensionの値と一致させます。
  
  説明: 各パラメータの単位。
  探索アルゴリズム中では、各パラメータをそれぞれこれらの値で割ることで、簡易的な無次元化・正規化を行います。定義しなかった場合にはすべての次元で 1.0 となります。
- min_list
  
  形式: 実数型のリスト。長さはdimensionの値と一致させます。
  
  説明: パラメータが取りうる最小値。
  モンテカルロ探索中にこの値を下回るパラメータが出現した場合、ソルバーは評価されずに、値が無限大だとみなされます。
- max_list
  
  形式: 実数型のリスト。長さはdimensionの値と一致させます。
  
  説明: パラメータが取りうる最大値。
  モンテカルロ探索中にこの値を上回るパラメータが出現した場合、ソルバーは評価されずに、値が無限大だとみなされます。
離散空間
- mesh_path
  
  形式: ファイルパス
  
  説明: メッシュ定義ファイル。
- neighborlist_path
  
  形式: ファイルパス
  
  説明: 近傍リスト定義ファイル。

[`pamc`] セクション¶

numsteps

形式: 整数値。

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

形式: 整数値。

説明: 「温度」を下げる頻度。この回数モンテカルロ更新を行った後に温度が下がります。
numT

形式: 整数値。

説明: 「温度」点の数。
Tmin

形式: 実数値。

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

形式: 実数値。

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

形式: 実数値。

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

形式: 実数値。

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

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

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

形式: 整数。 (default: 1)

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

形式: 整数。 (default: 1)

説明: レプリカのリサンプリングを行う頻度。この回数だけ温度降下を行った後にリサンプリングが行われます。
fix_num_replicas

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

説明: リサンプリングの際、レプリカ数を固定するかどうか。

ステップ数について¶

numsteps, numsteps_annealing, numT の3つのうち、どれか2つを同時に指定してください。残りの1つは自動的に決定されます。

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

メッシュ定義ファイル¶

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

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

6.000000 6.000000
6.000000 5.750000
6.000000 5.500000
6.000000 5.250000
6.000000 5.000000
6.000000 4.750000
6.000000 4.500000
6.000000 4.250000
6.000000 4.000000
...

近傍リスト定義ファイル¶

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

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

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

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

出力ファイル¶

`RANK/trial_T#.txt`¶

各温度点(#) ごと、モンテカルロサンプリングで提案されたパラメータと、対応する目的関数の値です。 1列目にステップ数、2列目にプロセス内のwalker 番号、3列目にレプリカの逆温度、4列目に目的関数の値、5列目からパラメータが記載されます。最後の2列はそれぞれレプリカの重み (Neal-Jarzynski weight) と祖先（計算開始時のレプリカ番号）です。

# step walker beta fx x1 weight ancestor
0 0.0 73.82799488298886 8.592321856342956 1.0 0
1 0.0 13.487174782058675 -3.672488908364282 1.0 1
2 0.0 39.96292704464803 -6.321623766458111 1.0 2
3 0.0 34.913851603463 -5.908794428939206 1.0 3
4 0.0 1.834671825646121 1.354500581633733 1.0 4
5 0.0 3.65151610695736 1.910894059585031 1.0 5
...

`RANK/trial.txt`¶

trial_T#.txt をすべてまとめたものです。

`RANK/result_T#.txt`¶

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

# step walker beta fx x1 weight ancestor
0 0.0 73.82799488298886 8.592321856342956 1.0 0
1 0.0 13.487174782058675 -3.672488908364282 1.0 1
2 0.0 39.96292704464803 -6.321623766458111 1.0 2
3 0.0 34.913851603463 -5.908794428939206 1.0 3
4 0.0 1.834671825646121 1.354500581633733 1.0 4
5 0.0 3.65151610695736 1.910894059585031 1.0 5
...

`RANK/result.txt`¶

result_T#.txt をすべてまとめたものです。

`best_result.txt`¶

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

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

`fx.txt`¶

各温度ごとに、全レプリカの情報をまとめたものです。 1列目は逆温度が、2列目と3列目には目的関数の期待値およびその標準誤差が、4列目にはレプリカの総数が、5列目には規格化因子（分配関数）の比の対数

\[\log\frac{Z}{Z_0} = \log\int \mathrm{d}x e^{-\beta f(x)} - \log\int \mathrm{d}x e^{-\beta_0 f(x)}\]

が、6列目にはモンテカルロ更新の採択率が出力されます。ここで \(\beta_0\) は計算している \(\beta\) の最小値です。

# $1: 1/T
# $2: mean of f(x)
# $3: standard error of f(x)
# $4: number of replicas
# $5: log(Z/Z0)
# $6: acceptance ratio
0.0 33.36426034198166 3.0193077565358273 100 0.0 0.9804
0.1 4.518006242920819 0.9535301415484388 100 -1.2134775491597027 0.9058
0.2 1.5919146358616842 0.2770369776964151 100 -1.538611313376179 0.9004
...

アルゴリズム解説¶

問題と目的¶

分布パラメータ \(\beta_i\) のもとでの配位 \(x\) の重みを \(f_i(x)\) と書くと（例えばボルツマン因子 \(f_i(x) = \exp\left[-\beta_i E(x)\right]\)）、 \(A\) の期待値は

\[\langle A\rangle_i = \frac{\int \mathrm{d}xA(x)f_i(x)}{\int \mathrm{d}x f_i(x)} = \frac{1}{Z}\int \mathrm{d}xA(x)f_i(x) = \int \mathrm{d}xA(x)\tilde{f}_i(x)\]

とかけます。ここで \(Z = \int \mathrm{d} x f_i(x)\) は規格化因子（分配関数）で、 \(\tilde{f}(x) = f(x)/Z\) は配位 \(x\) の確率密度です。

目的は複数の分布パラメータについてこの期待値および規格化因子（の比）を数値的に求めることです。

Annealed Importance Sampling (AIS) [1]¶

次の同時確率分布

\[\tilde{f}(x_0, x_1, \dots, x_n) = \tilde{f}_n(x_n) \tilde{T}_n(x_n, x_{n-1}) \tilde{T}_{n-1}(x_{n-1}, x_{n-2}) \cdots \tilde{T}_1(x_1, x_0)\]

を満たす点列 \(\{x_i\}\) を考えます。ここで

\[\tilde{T}_i(x_i, x_{i-1}) = T_i(x_{i-1}, x_i) \frac{\tilde{f}_i(x_{i-1})}{\tilde{f}_i(x_i)}\]

であり、 \(T_i(x, x')\) は \(\beta_i\) のもとでの配位 \(x\) から \(x'\) への遷移確率で、釣り合い条件

\[\int \mathrm{d}x \tilde{f}_i(x) T_i(x, x') = \tilde{f}_i(x')\]

を満たすようにとります（つまりは普通のMCMCにおける遷移確率行列）。

\[\int \mathrm{d} x_{i-1} \tilde{T}_i(x_i, x_{i-1}) = \int \mathrm{d} x_{i-1} \tilde{f}_i(x_{i-1}) T_i(x_{i-1}, x_i) / \tilde{f}_i(x_i) = 1\]

となるので、 \(\tilde{f}_n(x_n)\) は \(\tilde{f}(x_0, x_1, \dots, x_n)\) の周辺分布

\[\tilde{f}_n(x_n) = \int \prod_{i=0}^{n-1} \mathrm{d} x_i \tilde{f}(x_0, x_1, \dots, x_n)\]

です。これを利用すると、 \(\tilde{f}_n\) における平均値 \(\langle A \rangle_n\) は拡張した配位の重み付き平均として

\[\begin{split}\begin{split} \langle A \rangle_n &\equiv \int \mathrm{d} x_n A(x_n) \tilde{f}_n(x_n) \\ &= \int \prod_i \mathrm{d} x_i A(x_n) \tilde{f}(x_0, x_1, \dots, x_n) \end{split}\end{split}\]

と表せます。

さて、残念ながら \(\tilde{f}(x_0, x_1, \dots, x_n)\) に従うような点列を直接生成することは困難です。そこでもっと簡単に、

確率 \(\tilde{f}_0(x)\) に従う \(x_0\) を生成する
- 例えば MCMC を利用する
\(x_i\) から \(T_{i+1}(x_i, x_{i+1})\) によって \(x_{i+1}\) を生成する
- \(T_{i+1}\) は釣り合い条件を満たすような遷移確率行列なので、普通にMCMCを行えば良い

という流れに従って点列 \(\{x_i\}\) を生成すると、これは同時確率分布

\[\tilde{g}(x_0, x_1, \dots, x_n) = \tilde{f}_0(x_0) T_1(x_0, x_1) T_2(x_1, x_2) \dots T_n(x_{n-1}, x_n)\]

に従います。これを利用すると期待値 \(\langle A \rangle_n\) は

\[\begin{split}\begin{split} \langle A \rangle_n &= \int \prod_i \mathrm{d} x_i A(x_n) \tilde{f}(x_0, x_1, \dots, x_n) \\ &= \int \prod_i \mathrm{d} x_i A(x_n) \frac{\tilde{f}(x_0, x_1, \dots, x_n)}{\tilde{g}(x_0, x_1, \dots, x_n)} \tilde{g}(x_0, x_1, \dots, x_n) \\ &= \left\langle A\tilde{f}\big/\tilde{g} \right\rangle_{g, n} \end{split}\end{split}\]

と評価できます (reweighting method)。 \(\tilde{f}\) と \(\tilde{g}\) との比は、

\[\begin{split}\begin{split} \frac{\tilde{f}(x_0, \dots, x_n)}{\tilde{g}(x_0, \dots, x_n)} &= \frac{\tilde{f}_n(x_n)}{\tilde{f}_0(x_0)} \prod_{i=1}^n \frac{\tilde{T}_i(x_i, x_{i-1})}{T(x_{i-1}, x_i)} \\ &= \frac{\tilde{f}_n(x_n)}{\tilde{f}_0(x_0)} \prod_{i=1}^n \frac{\tilde{f}_i(x_{i-1})}{\tilde{f}_i(x_i)} \\ &= \frac{Z_0}{Z_n} \frac{f_n(x_n)}{f_0(x_0)} \prod_{i=1}^n \frac{f_i(x_{i-1})}{f_i(x_i)} \\ &= \frac{Z_0}{Z_n} \prod_{i=0}^{n-1} \frac{f_{i+1}(x_{i})}{f_i(x_i)} \\ &\equiv \frac{Z_0}{Z_n} w_n(x_0, x_1, \dots, x_n) \end{split}\end{split}\]

とかけるので、期待値は

\[\langle A \rangle_n = \left\langle A\tilde{f}\big/\tilde{g} \right\rangle_{g, n} = \frac{Z_0}{Z_n} \langle Aw_n \rangle_{g,n}\]

となります。規格化因子の比 \(Z_n/Z_0\) は \(\langle 1 \rangle_n = 1\) を用いると

\[\frac{Z_n}{Z_0} = \langle w_n \rangle_{g,n}\]

と評価できるので、 \(A\) の期待値は

\[\langle A \rangle_n = \frac{\langle Aw_n \rangle_{g,n}}{\langle w_n \rangle_{g,n}}\]

という、重み付き平均の形で評価できます。この重み \(w_n\) を Neal-Jarzynski 重みと呼びます。

population annealing (PA) [2]¶

AIS を使うと各 \(\beta\) に対する期待値を重み付き平均という形で計算できますが、 \(\beta\) の幅が大きくなると重み \(w\) の分散が大きくなってしまいます。そのため、適当な周期で確率 \(p^{(k)} = w^{(k)} / \sum_k w^{(k)}\) に従いレプリカをリサンプリングし、レプリカに割当られた重みをリセット \((w=1)\) します。

PAMC のアルゴリズムは次の擬似コードで示されます:

for k in range(K):
    w[0, k] = 1.0
    x[0, k] = draw_from(β[0])
for i in range(1, N):
    for k in range(K):
        w[i, k] = w[i-1, k] * ( f(x[i-1,k], β[i]) / f(x[i-1,k], β[i-1]) )
    if i % interval == 0:
        x[i, :] = resample(x[i, :], w[i, :])
        w[i, :] = 1.0
    for k in range(K):
        x[i, k] = transfer(x[i-1, k], β[i])
    a[i] = sum(A(x[i,:]) * w[i,:]) / sum(w[i,:])

リサンプリング手法として、レプリカ数を固定する方法[2]と固定しない方法[3]の2通りがあります。

参考文献¶

[1] R. M. Neal, Statistics and Computing 11, 125-139 (2001).

[2] K. Hukushima and Y. Iba, AIP Conf. Proc. 690, 200 (2003).

[3] J. Machta, PRE 82, 026704 (2010).

2DMAT 2.2.0 ドキュメント