粒子フィルターを用いたサンプリング | About SMC Samplers

テンパリングを通じたもう一つの万能サンプラー

Particles
MCMC
Survey
Author

司馬博文

Published

12/14/2023

Modified

8/05/2024

概要
粒子フィルターは 30 年前に「万能」非線型フィルタリング手法として開発されたが，それは粒子系を輸送するメカニズムとしての万能性も意味するのであり，汎用サンプラーとしても「万能」であるのかもしれないのである．近年，最適化や最適輸送の理論と結びつき，その真の力がますます明らかになりつつある．本稿では現在までのサンプラーとしての SMC 手法に対する理解をまとめる．

1 テンパリングの歴史

SMC の文脈で，目標の分布 \(\pi_p\in\mathcal{P}(E)\) が複雑であるとき，これに至る \(\mathcal{P}(E)\) 上の道 \[ [p]\ni n\mapsto\pi_n\in\mathcal{P}(E) \] を通じて，より簡単な分布 \(\pi_1,\pi_2,\cdots\) から逐次的にサンプリングをする，というアイデアを 調温 (tempering) という（粒子フィルターの稿 も参照）．

この tempering という考え方は本質的に逐次的な発想を持っているが，元々は SMC の文脈とは全く独立に，MCMC を多峰性を持つ複雑な分布に対しても使えるように拡張する研究で提案された．これらの手法が自然と SMC へと接続する様子を population Monte Carlo (Iba, 2001b), (Ajay Jasra et al., 2007) というキーワードで理解されている．

まずはその歴史を概観する．いずれも，目標分布 \(\pi_p\) が多峰性をもち，MCMC がうまく峰の間を遷移できずに正しいサンプリングができない（収束が遅くなる）問題を解決する文脈の中で捉えられる．

1.1 模擬アニーリング (Kirkpartick et al., 1983)

これは MCMC とは関係がなく，もはやシミュレーション法でさえなく最適化手法であるが，「調温」の考え方を一気にポピュラーにした手法であった．1 汎用最適化手法として，半導体製造を通じて，電子工学・コンピュータ産業にも大きな影響を与えた手法である．

そもそも 焼きなまし (annealing) とは，物性物理の用語であり，鉄などの固体を極めて高音にして溶解させたのちに徐々に冷却することで，基底状態の構造を得るのに使われる技術であった．2

分布列を \(\pi_n\,\propto\,e^{-\frac{h(x)}{T_n}}\,dx\) \[ T_1>T_2>\cdots>T_n\searrow 0 \] と構成することで， \[ \pi_n\xrightarrow{n\to\infty}1_{\operatorname*{argmin}h}(x)\,dx \] であることを利用して，関数 \(h\) の最小値を見つけることができる．3

1.2 拡張アンサンブル法

MCMC を複数同時に実行する手法を 拡張アンサンブル法 という (Iba, 2001a)．これは正準集団などの物理的根拠のあるアンサンブルを用いるのではなく，人工のアンサンブルを導入してサンプリング効率を向上させると捉えられるために呼ぶ．4

multilevel sampling とも呼ばれる．5

一方で，次節 1.3 で扱う相関粒子法も含めて，複数のサンプルを用いる手法はとして population-based method とも呼ばれる (Iba, 2001b), (Ajay Jasra et al., 2007)．

1.2.1 傘サンプル法 (Torrie and Valleau, 1977)

(Torrie and Valleau, 1977) では系のポテンシャルに傘ポテンシャルと呼ばれる追加項を足すことで，本来なら到達できない状態からもサンプリングすることを可能にするアイデアであり，拡張アンサンブル法の最初の萌芽と捉えられる．

この傘ポテンシャルとして，上述の意味でのテンパリング分布をとることも提案されており，後述の種々のテンパリング法の先駆けともみなせるのである．6

1.2.2 MC3 / 並行テンパリング / 交換モンテカルロ (Geyer, 1991), (Hukushima and Nemoto, 1996)

積空間 \(\otimes_{n=1}^pE\) 上で \(\pi_1\otimes\cdots\otimes\pi_p\) を目標分布として MCMC を実行することを考えるのが MC3 (Metropolis-Coupled MCMC) (Geyer, 1991) である．

時折，不変分布を変えないような Metropolis 核による提案に従って，MCMC 鎖の位置を交換することで収束を加速する．

この手法は parallel tempering7 または exchange Monte Carlo (Hukushima and Nemoto, 1996) という名前による独立な提案に伴って 交換モンテカルロ または レプリカ交換法，8 さらには population-based MCMC9 とも呼ばれる．

特に，その分子動力学法版（REMD）(Sugita and Okamoto, 1999) が開発されてからは，分子シミュレーションの分野に広く受け入れられ，AMBER, CHARMM, GROMACS, NAMD などの汎用プログラムにも REMD が組み込まれた．(岡本祐幸, 2010)

マルチカノニカル法 1.2.5 や模擬テンパリング 1.2.4 では荷重を決定するために試行が必要であるが，並行テンパリングでは荷重は Boltzmann 因子であるため，このような予備試行は必要ない．10

しかしながら，全てのテンパリング手法に共通するように，交換の棄却率が高まりすぎないようにするためには隣り合う \(\pi_n,\pi_{n+1}\) を十分近く取る必要があり，すると必要な MCMC 鎖の数が極めて大きくなってしまうこともある．11

population-based (Iba, 2001b) というのは，\(p\) 個の粒子を展開して高温状態でも探索してもらい，定期的に粒子を交換することでその情報を互いに伝え合うメカニズムのように思えるために言う．12 この観点から見ると，「鎖の間の交換」とは，粒子の間の相互作用としては極めてナイーブなもので，粒子フィルターに見られるような遺伝的なアルゴリズムの導入でより効率化できるのではないか？という発想が出てくる．

1.2.3 進化的モンテカルロ

並行テンパリングに加えて，種々の population-based method が提案された．(Ajay Jasra et al., 2007) によるレビューも参照．

まずは Adaptive direction sampling (Gilks et al., 1994) がある．これは複数の粒子 \(\boldsymbol{x}:=\{x_t^n\}_{n=1}^p\) を，

  1. ある \(x_t^a\in\boldsymbol{x}\) を選んで，ここからアンカーポイント \(y\in E\) を何かしらの方法で定める．
  2. \(x_t^c\in\boldsymbol{x}\setminus\{x_t^a\}\) を選んで，1 で定めた \(y\in E\) の方向にランダムに動かす．

の繰り返しによって発展させていくことによりサンプリングする手法である．

このような手続きを，遺伝的アルゴリズムの考え方を取り入れてさらに推し進め，実際に MCMC としての収束レートを速めたのが 進化モンテカルロ (Liang and Wong, 2000), (Liang and Wong, 2001) である．

1.2.4 模擬テンパリング (Marinari and Parisi, 1992)

最適化手法である 焼きなまし法（または模擬アニーリング） (Kirkpartick et al., 1983) のサンプリングへの変形として提案されたのが 焼き戻し法，または 模擬テンパリング (simulated tempering) (Marinari and Parisi, 1992) である．13

模擬アニーリングでは温度は下がる一方であったのが，模擬テンパリングでは温度もある周辺分布に従って遷移する．模擬アニーリングは最終的にサンプルが最小値点の周りに集積して最適化問題を解くことが目的であったが，模擬テンパリングは高温状態においては多峰性分布が軟化され，峰の間を遷移しやすくなることを利用し，多峰性分布からの効率的なサンプリングを目指す．

模擬テンパリングは状態空間を \(E\times [p]\) に拡大して，その上でサンプリングを行うものともみなせる．14 \(E\times[p]\) 上の標的分布を \[ X|N=n\sim\pi_n \] を満たすようにし，\(N|X\) は適宜架橋分布 \(\{\pi_n\}\) を往来するよう設計することで，MC3 が \(p\) 本の MCMC を用いて実現していたことを，\(E\times [p]\) 上の MCMC 1つで効率的に実行する．

また，MCMC の収束を大幅に加速する手法としても，遺伝学における複雑な事後分布からのサンプリングへの応用を念頭に独立に提案された (Geyer and Thompson, 1995)．

1.2.5 マルチカノニカル法 (Berg and Neuhaus, 1991)

マルチカノニカル法 (Berg and Neuhaus, 1991) もポテンシャルを人工的に変更する方法であり，この点で傘サンプリングの発展ともみなせ，Adaptive umberlla sampling とも呼ばれる (Iba, 2001a)．

物性物理学の分野から提案され，スピングラスの問題などでも大きな成果を挙げた．15

1.3 相関粒子系への発展

1.3.1 テンパリング遷移 (Neal, 1996)

tempered transitions では，架橋列 \(\{\pi_n\}\) をそれぞれの \(\pi_n\) を不変分布に持つ Markov 核を通じて１往復して探索し，その結果を元に \(\pi_p\) を効率的に探索するような MCMC の提案を構成する．16

この方法は混合モデルにおいて事後分布が多峰性を持つなどして Gibbs サンプラーがうまく収束しない場合でも，有効な MCMC サンプラーとなる (A. Jasra et al., 2005)．

また， \[ \pi_n(x)\,\propto\,\pi_0(x)e^{-\beta_nh(x)} \] と表せる際，架橋分布 \(\{\pi_n\}\) は温度比 \(\beta_n/\beta_{n+1}\) が一定になるように 幾何的に 取ることを提案しており，現在でも一般的な基準であるようである (Behrens et al., 2012)．

1.3.2 焼きなまし重点サンプリング (AIS: Annealed Importance Sampling) (Neal, 2001)

ここで初めて SMC の文脈にもテンパリングが輸入された．17 (Neal, 2001) は重点サンプリングによってあらゆる温度 \(\{\pi_n\}\) からの提案を効率的に採用する方法を模索した．

AIS は，各 \(\pi_i\) を不変分布とする MCMC 核 \(P_i\) について，\(\pi_0P_1P_2\cdots P_p\) を重点サンプリング法における提案分布に用いる方法である．

しかし，そのまま重点荷重を計算するのではなく，18 拡張された空間 \(E^{p+1}\) 上の目標分布 \[ \pi_p\otimes P_p^{-1}\otimes\cdots\otimes P_1^{-1} \] に対して \(P_p\otimes P_{p-1}\otimes\cdots\otimes P_1\otimes\pi_0\) を提案分布に用いたとして荷重荷重を計算する．19 実際には， \[ X_p\sim P_{p}(X_{p-1},-),\quad X_{p-1}\sim P_{p-1}(X_{p-2},-),\quad \cdots\quad X_1\sim P_1(X_0,-),\quad X_0\sim \pi_0 \] というように \(X_0\sim\pi_0\) を MCMC 核 \(P_1,\cdots,P_p\) で順に流し，最後にウェイト \[ w(X_{1:p}):=\frac{\pi_p(X_p)}{\pi_{p-1}(X_{p})}\frac{\pi_{p-1}(X_{p-1})}{\pi_{p-2}(X_{p-1})}\cdots\frac{\pi_2(X_2)}{\pi_1(X_2)}\frac{\pi_1(X_1)}{\pi_0(X_1)} \] を計算する．20

従って，本当は \(E^{p+1}\) 上で重点サンプリングを行っているが，\(x_p\) の成分のみに注目することで周辺分布では \(\pi_p\) に対する効率的な重点サンプリングが実現されている．

テンパリング遷移の後半のアルゴリズムを発展させた形とも見れる．

同様の手法は自由エネルギーの推定の文脈で統計物理学で独立に提案されている (Jarzynski, 1997b), (Jarzynski, 1997a), (Crooks, 1998)．21

1.3.3 重点テンパリング (Gramacy et al., 2010)

こちらは模擬テンパリングを基にし，他の温度からの提案を保持しておく機構を提案している．

1.3.4 荷重を保つ模擬テンパリング (Tawn et al., 2020)

1.4 デノイジング拡散過程と最適架橋

簡単な分布からサンプリングをし，データの分布まで輸送するという発想は，生成モデリング，特に拡散過程のそれと同一である．

ここでは，近年の拡散過程とスコアマッチングの研究と SMC の交差について調べる．

1.4.1 SMC サンプラー (Del Moral et al., 2006)

1.5 その他の手法

1.5.1 多峰性の最適化に基づく対処

目標分布の峰を特定するタスクを MCMC から分離して，BFGS 法 に基づく最適化法によって先に解いてしまう手法が (Pompe and Łatuszyński, 2020) によって提案されている．

これにより探索した峰の全体を \(\mathcal{I}:=\{1,\cdots,I\}\) に格納し，拡大した状態空間 \(E\times\mathcal{I}\) 上で \(\widetilde{\pi}\) を対象とした MCMC を実行するが，この \(\widetilde{\pi}\) をさらに適応的に更新する Auxiliary Variable Adaptive MCMC を提案している．

2 テンパリングと最適化の関係

近年，(Chopin et al., 2023) で指摘されたように，テンパリングを通じた SMC サンプラーは，\(\mathcal{P}(E)\) 上での最適化としての解釈も持つことが理解されつつある．

References

Behrens, G., Friel, N., and Hurn, M. (2012). Tuning tempered transitions. Statistics and Computing, 22, 65–78.
Berg, B. A., and Neuhaus, T. (1991). Multicanonical algorithms for first order phase transitions. Physics Letters B, 267(2), 249–253.
Bouchard-Côté, A., Sankararaman, S., and Jordan, M. I. (2012). Phylogenetic inference via sequential monte carlo. Systematic Biology, 61(4), 579–593.
Chopin, N., Crucinio, F. R., and Korba, A. (2023). A connection between tempering and entropic mirror descent.
Chopin, N., and Papaspiliopoulos, O. (2020). An introduction to sequential monte carlo. Springer Cham.
Crooks, G. E. (1998). Nonequilibrium Measurements of Free Energy Differences for Microscopically Reversible Markovian Systems. Journal of Statistical Physics, 90(5), 1481–1487.
Del Moral, P., Doucet, A., and Jasra, A. (2006). Sequential Monte Carlo Samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3), 411–436.
Doucet, A., Grathwohl, W. S., Matthews, A. G. D. G., and Strathmann, H. (2022). Score-based diffusion meets annealed importance sampling. In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in neural information processing systems.
Geman, S., and Geman, D. (1984). Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6(6), 721–741.
Geyer, C. J. (1991). Markov chain monte carlo maximum likelihood. In E. M. Keramidas, editor, Computing science and statistics: Proceedings of the 23rd symposium on the interface, pages 156–163. Interface Foundation.
Geyer, C. J., and Thompson, E. A. (1995). Annealing markov chain monte carlo with applications to ancestral inference. Journal of the American Statistical Association, 90(431), 909–920.
Gilks, W. R., Roberts, G. O., and George, E. I. (1994). Adaptive direction sampling. Journal of the Royal Statistical Society. Series D (The Statistician), 43(1), 179–189.
Gramacy, R., Samworth, R., and King, R. (2010). Importance tempering. Statistics and Computing, 20, 1–7.
Hukushima, K., and Nemoto, K. (1996). Exchange monte carlo method and application to spin glass simulations. Journal of the Physical Society of Japan, 65(6), 1604–1608.
Iba, Y. (2001a). Extended ensemble monte carlo. International Journal of Modern Physics C, 12(05), 623–656.
Iba, Y. (2001b). Population monte carlo algorithms. 人工知能学会論文誌, 16(2), 279–286.
Jarzynski, C. (1997a). Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach. Phys. Rev. E, 56, 5018–5035.
Jarzynski, C. (1997b). Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett., 78, 2690–2693.
Jasra, A., Holmes, C. C., and Stephens, D. A. (2005). Markov Chain Monte Carlo Methods and the Label Switching Problem in Bayesian Mixture Modeling. Statistical Science, 20(1), 50–67.
Jasra, Ajay, Stephens, D. A., and Holmes, C. C. (2007). On population-based simulation for static inference. Statistics and Computing, 17, 263–279.
Kirkpartick, S., Gelatt, C. D., and Vecchi, M. P. (1983). Optimization by simulated annealing. Science, 220(4598), 671–680.
Liang, F., and Wong, W. H. (2000). EVOLUTIONARY MONTE CARLO: APPLICATIONS TO c p MODEL SAMPLING AND CHANGE POINT PROBLEM. Statistica Sinica, 10(2), 317–342.
Liang, F., and Wong, W. H. (2001). Real-parameter evolutionary monte carlo with applications to bayesian mixture models. Journal of the American Statistical Association, 96(454), 653–666.
Liu, J. S. (2004). Monte carlo strategies in scientific computing. Springer New York.
Lyubartsev, A. P., Martsinovski, A. A., Shevkunov, S. V., and Vorontsov‐Velyaminov, P. N. (1992). New approach to Monte Carlo calculation of the free energy: Method of expanded ensembles. The Journal of Chemical Physics, 96(3), 1776–1783.
Marinari, E., and Parisi, G. (1992). Simulated tempering: A new monte carlo scheme. Europhysics Letters, 19(6), 451–458.
Neal, R. M. (1996). Sampling from Multimodal Distributions Using Tempered Transitions. Statistics and Computing, 6, 353–366.
Neal, R. M. (2001). Annealed Importance Sampling. Statistics and Computing, 11, 125–139.
Pompe, E., and Łatuszyński, C. H. K. (2020). A framework for adaptive MCMC targeting multimodal distributions. The Annals of Statistics, 48(5), 2930–2952.
Sugita, Y., and Okamoto, Y. (1999). Replica-exchange molecular dynamics method for protein folding. Chemical Physics Letters, 314(1), 141–151.
Swendsen, R. H., and Wang, J.-S. (1986). Replica monte carlo simulation of spin-glasses. Physical Review Letters, 57(21), 2607–2609.
Tawn, N. G., Roberts, G. O., and Rosenthal, J. S. (2020). Weight-preserving simulated tempering. Statistics and Computing, 30, 27–41.
Torrie, G. M., and Valleau, J. P. (1977). Nonphysical sampling distributions in monte carlo free-energy estimation: Umbrella sampling. Journal of Computational Physics, 23(2), 187–199.
岡本祐幸. (2010). 「拡張アンサンブル法」の特集にあたって. アンサンブル, 12(2), 2_6–2_7.
酒井佑士. (2017). マルコフ連鎖モンテカルロ法における詳細つり合い条件の破れの効果と応用 (PhD thesis). 東京大学. Retrieved from https://repository.dl.itc.u-tokyo.ac.jp/records/50422

Footnotes

  1. この前にも，Umbrella sampling (Torrie and Valleau, 1977) が本質的には密度の調温のアイデアを用いていた．(Liu, 2004, pp. 206–) Section 10.1 も参照．↩︎

  2. 分子動力学 (molecular dynamics) などの文脈では Metropolis 法はちょうど分子運動のシミュレーションになっていることを踏まえれば，これを simulated annealing と呼ぶのは極めて鮮やかなアナロジーとなっている．焼きなまし法自体も，シミュレーション可能になったのである．↩︎

  3. (Geman and Geman, 1984) によると，各 \(\pi_n\) における MCMC move の回数を \(N_n\) とした場合，\(O(\log(N_1+\cdots+N_n))\) のオーダーで \(T_n\) を（十分遅く）変化させれば，この手法はほとんど確実に \(\operatorname*{argmin}h\) 内に収束する．(Liu, 2004, pp. 209–) 10.2 節も参照．↩︎

  4. (岡本祐幸, 2010), (Iba, 2001a) など．↩︎

  5. (Liu, 2004, pp. 205–) Chapter 10. Multilevel Sampling and Optimization Methods も参照．↩︎

  6. (Liu, 2004, p. 207) も参照．↩︎

  7. (Chopin et al., 2023), (Liu, 2004, p. 4) でも (Geyer, 1991) を引用して PT と呼んでいる．一方で物理学の分野では (Hukushima and Nemoto, 1996) の exchange Monte Carlo や (Swendsen and Wang, 1986) などの文献もある．前者は (Liu, 2004, p. 4) が “is reminiscent of parallel tempering (Geyer, 1991)” と指摘しており，後者は (Bouchard-Côté et al., 2012) などが引用している．↩︎

  8. 最終講義 スピングラスと計算物性物理 p.34 も参照．温度の違う熱浴につけたレプリカをシミュレートして，時々交換する，という見方ができるためにこう呼ぶ．↩︎

  9. (Ajay Jasra et al., 2007) は (Geyer, 1991) を指して population-based MCMC と呼んでおり，SMC も含めて population-based simulation と呼んでいる．population-based という言葉自体は (Iba, 2001b) からとったという．“we define a population-based simulation method as one which, instead of sampling a single (independent/dependent) sample, generates a collection of samples in parallel” と定義しており，大きく MCMC によるものと逐次重点サンプリングベースのものの２流儀あるとしている．(Liu, 2004, pp. 225–) 第11章なども参照．↩︎

  10. (岡本祐幸, 2010) など．↩︎

  11. (Behrens et al., 2012, p. 66) も参照．↩︎

  12. (Iba, 2001a) が良い解説を与えていると (Ajay Jasra et al., 2007) でも言及されている．ただし，(Iba, 2001a) はこの並行テンパリングだけでなく，模擬テンパリング，multicanonical Monte Carlo (Berg and Neuhaus, 1991) / Adaptive Umbrella Sampling (Torrie and Valleau, 1977) を総称して 拡張アンサンブル法 (Extended Ensemble Monte Carlo) と呼んでサーベイしていることに注意．↩︎

  13. (Lyubartsev et al., 1992) が引用されることもある．(酒井佑士, 2017), (岡本祐幸, 2010) など．method of expanded ensemble とも呼ばれる (岡本祐幸, 2010), (Iba, 2001a)．↩︎

  14. 記法 \([p]=\{1,\cdots,p\}\) は 本サイトの数学記法一覧 を参照↩︎

  15. その後すぐに分子シミュレーションの分野にも導入された．(岡本祐幸, 2010) も参照．↩︎

  16. (Behrens et al., 2012) も参照．↩︎

  17. (Chopin and Papaspiliopoulos, 2020, p. 33) で，SMC を調温に初めて応用した論文として紹介されている．p.352 では “An early version of SMC tempering (without resampling)” としている．↩︎

  18. \(\pi_p(x_p)/\pi_0P_1P_2\cdots P_p\) は多くの場合計算不能である．↩︎

  19. ただし，\(P_i^{-1}\) とは，\[ P_i(x_{i-1},x_i)\pi_{i-1}(x_i-1)=\pi_i(x_i)P_i^{-1}(x_{i-1},x_i) \] で定まる確率核とした．\(\otimes\) の記法はこちらも参照．↩︎

  20. このウェイトの表示は，\(P_i^{-1}/P_i=\pi_{i-1}/\pi_i\) が成り立つことから直ちに従う．↩︎

  21. (Doucet et al., 2022) も参照．↩︎