粒子法の概観 – Hirofumi Shiba

概要

粒子法とは空間や分布を多数の粒子の集合として離散化して表現・計算する技術の総称である．シミュレーションからデータ同化まで幅広い応用を持つ．この記事ではこれらの技術を「粒子」という軸でひとつの記事にまとめることを試みる．

1 シミュレーションのための粒子法

1.1 はじめに

多くの物質は元々離散的であり，その最小単位は原子や分子である．

これをコンピュータ内に粒子として表現し，その現実世界での動きを模倣することを目指すシミュレーション法は分子動力学法 (MD: molecular dynamics) と呼ばれる．

アイデアは大変シンプルである．

原子論についての Feynman の名言

Feynman の講義録（のちに教科書になる）では，最も重要な科学知識の一つとして原子論を挙げている：

If, in some cataclysm, all of scientific knowledge were to be destroyed, and only one sentence passed on to the next generations of creatures, what statement would contain the most information in the fewest words? I believe it is the atomic hypothesis (or the atomic fact, or whatever you wish to call it) that all things are made of atoms—little particles that move around in perpetual motion, attracting each other when they are a little distance apart, but repelling upon being squeezed into one another. In that one sentence, you will see, there is an enormous amount of information about the world, if just a little imagination and thinking are applied. (Feynman et al., 1964)

1.2 分子動力学法

分子動力学法の文脈で 粒子法 (particle methods) と言ったとき，

古典的な質点には Newton 力学
剛体で体積を持つ場合は Euler 方程式
内部構造を持つ場合は Langrange 方程式
さらに一般の場合は Hamilton 方程式

に基づいて粒子の動きをシミュレーションすることを指す．

1.3 モンテカルロ法

分子動力学法は Monte Carlo 法と組み合わせて使われることが多い．

粒子の動きを上述の方程式に基づいて提案し，離散化誤差を Metropolis 法によって補正することで，系の物理量の分布を求めたり，平均値を求めたりする (Rapaport, 2004, p. 4), (Griebel et al., 2007, p. 17)．

そもそも MCMC の手法 (Metropolis et al., 1953) 自体も，Schrödinger 方程式に基づいた基底状態のシミュレーションを実行するために生まれたものである．

これが強力なのである．極めて強力であるが故に，計算科学の分野ではすでに，そして将来的には確実に，シミュレーションと数値実験が，実際の実験と似たような重要な役割を科学において担うことになると見られている．

1.4 推論のための粒子法との関係

この意味での粒子法と，統計学において粒子法と言った場合に想像されるであろう SMC や IPM (interacting particle methods) とは 実は数理的には同一物である ことを見たい．

一言で言えば，一般に粒子法はポテンシャル（ハミルトニアン）に従って粒子を伝播する方法であり，粒子フィルターは特にこれが Feynman-Kac 測度で与えられる場合に用いられる粒子法である．

2 推論のための粒子法

2.1 ベイズ逆問題

シミュレーションができるようになったならば，これを実際の実験データや観測データと付き合わせて比較したいということになる．

このようにデータに基づいてシミュレーションを校正する作業は 逆問題 (inverse problem) と呼ばれ，旧来は統計的推論との関係は希薄であった．

逆問題の例

逆散乱問題：反響定位 (echolocation) やソナー，地震探索 (seismic exploration)
断層映像法 (tomography)：海洋音響トモグラフィーや地震波トモグラフィーなど．
医用画像 (medical imaging)：X線撮影，CT, MRI など．
非破壊検査 (NDT: Non-Destructive Testing)：超音波探傷試験 (UT: Ultrasonic Testing / SAT: Scanning Acoustic Tomography)
認知科学：脳波 (EEG: Electroencephalogram)，脳磁図 (MEG: Magnetoencephalography) など．¹

逆問題の形式的な定義

逆問題は広く Banach 関数上の作用素 \(F:X\osup U\to Y\) の零点を求める問題として定式化される (定義1.2.1 Nakamura and Potthast, 2015, pp. 1–15)．

大きく分類すると，源泉問題 (source problem) と散乱問題 (scattering problem)，力学系逆問題 (dynamical systems inversion)，スペクトル逆問題 (spectral inverse problem) との４つに分けられる．

場の復元と形状の復元は散乱問題の特別な場合とみれる．スペクトル逆問題は (Kac, 1966) で有名である．

データ同化問題は力学系逆問題を逐次的に解く方法と位置付けられる (Nakamura and Potthast, 2015, pp. 1–20)．

しかし「系に対する，方程式などの形で与えられる事前知識をもとに，データから実際の構造に対する知識をアップデートしていく」という営みはまさにベイズ推論の一種である．

特にベイズによる取り扱いが統一的な数学的基礎を与えることに加え，従来からの無定義用語であった well-posedness に数学的に厳密な定義を与えることができることから，近年 ベイズ逆問題 (Bayesian inverse problem) と呼ばれる分野が進展を見せている (Stuart, 2010)．

さらに従来の逆問題のためのアルゴリズムは統計学における対応物を持つことが多い．例えば Tikhonov-Phillips 正則化と呼ばれる手続きは Ridge 回帰 (Hoerl and Kennard, 1970) に他ならない．データ同化においては変分推論が先駆けて古くから使われている (Sasaki, 1958)．

ベイズの取り扱いはこのようなアルゴリズムの開発に対しても統一的な基盤を与える上に，アルゴリズムの最適性の議論も可能になる (Nickl, 2017)．

2.2 データ同化

極めて現代的な逆問題の発展形が データ同化 (data assimilation) または 動的逆問題 (dynamic inverse problem) である．²

データ同化は地球科学，特に気象予測の分野で芽生えた考え方である．気象予測では大規模なシミュレーションが欠かせないが，同時に正確な予測のためには初期値の正確さが重要であることがその黎明期から認識されていた (Kalnay, 2002)．

しかし系が巨大すぎるあまり，観測データから得られる情報は極めて限定的であり，むしろ初期値を正確に知るためにもシミュレーションや構造的な知識が必要である．

このように順問題と逆問題とを逐次的に考えることがデータ同化の特徴である．もしかしたら data-simulation assimilation や data-model assimilation と呼ぶべきかもしれない．

現在データ同化は測定値とシミュレーションを同化させて精度の良い予測を行う営為全般に応用が進んでいる (Nakamura and Potthast, 2015)．

2.3 EnKF

データ同化のダイナミックな問題設定が制御理論に近く，Kalman フィルター (Kalman, 1960) が利用可能であることが 1980 年代に自覚された (Ghil et al., 1981)．

しかし数値気象予測では状態空間が極めて高次元であり，計算コストの削減が必要であった．

アンサンブルカルマンフィルター (EnKF: Ensemble Kalman Filter) (Evensen, 1994) はデータ同化の 2024 時点でもデファクトスタンダードである．

2.4 変分法

(Sasaki, 1958) など黎明期から変分法の応用が考えられていた．

この方法は (Lions, 1971) で最適制御の 随伴法 (adjoint method) (Pontryagin et al., 1962) と合流し，3D-VAR, 4D-VAR などの手法が提案されている (Talagrand and Courtier, 1987)．

随伴法は連続正規化流／ニューラル ODE のパラメータ推定でも用いられている．

2.5 粒子フィルターの位置付け

従来から粒子法が用いられていたが，粒子フィルターは今後の大規模データ世界におけるデータ同化の主力技術と目されている (Leeuwen et al., 2019)．

例えば気象予測の分野で，積雲対流 (cumulus convection) という現象は数 km のスケールで起こるが，大気循環の全体の記述に不可欠であることが知られている (浅井冨雄, 1967)．

これを観測する数 km 解像度の convective scale でのデータ同化は従来の \(10^4\) オーダーでの総観スケール (synoptic scale) でのデータ同化よりもはるかに非線型性が強いため，従来法の革新が迫れられており (Yano et al., 2018)，粒子フィルターが注目されている．

References

Evensen, G. (1994). Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics. Journal of Geophysical Research, 99(C5), 10143–10162.

Feynman, R. P., Leighton, R. B., and Sands, M. (1964). The feynman lectures of physics. In,Vol. I. Addison-Wesley.

Ghil, M., Cohn, S., Tavantzis, J., Bube, K., and Isaacson, E. (1981). Applications of estimation theory to numerical weather prediction. In L. Bengtsson, M. Ghil, and E. Källén, editors, Dynamic meteorology: Data assimilation methods, pages 139–224. New York, NY: Springer New York.

Griebel, M., Zumbusch, G., and Knapek, S. (2007). Numerical simulation in molecular dynamics: Numerics, algorithms, parallelization, applications,Vol. 5. Springer Belin, Heidelberg.

Hoerl, A. E., and Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55–67.

Kac, M. (1966). Can one hear the shape of a drum? The American Mathematical Monthly, 73(4P2), 1–23.

Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35–45.

Kalnay, E. (2002). Atmospheric modeling, data assimilation and predictability. Cambridge University Press.

Leeuwen, P. J. van, Künsch, H. R., Nerger, L., Potthast, R., and Reich, S. (2019). Particle filters for high-dimensional geoscience applications: A review. Quarterly Journal of the Royal Meteorological Society, 145(723), 2335–2365.

Lions, J. L. (1971). Optimal control of systems governed by partial differential equations,Vol. 170. Springer Berlin, Heidelberg.

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953). Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6), 1087–1092.

Nakamura, G., and Potthast, R. (2015). Inverse modeling: An introduction to the theory and methods of inverse problems and data assimilation. IOP Publishing.

Nickl, R. (2017). On Bayesian Inference for Some Statistical Inverse Problems with Partial Differential Equations. In Bernoulli news,Vol. 24. Bernoulli Society.

Pontryagin, L. S., Boltyanskii, V. G., Gamkrelidze, R. V., and Mishchenko, E. F. (1962). The mathematical theory of optimal processes. John Wiley & Sons.

Rapaport, D. C. (2004). The art of molecular dynamics simulation. Cambridge University Press.

Sasaki, Y. (1958). An ObJective analysis based on the variational method. 氣象集誌. 第2輯, 36(3), 77–88.

Stuart, A. M. (2010). Inverse problems: A bayesian perspective. Acta Numerica, 19, 451–559.

Talagrand, O., and Courtier, P. (1987). Variational assimilation of meteorological observations with the adjoint vorticity equation. I: theory. Quarterly Journal of the Royal Meteorological Society, 113(478), 1311–1328.

Ter-Pogossian, M. M., Phelps, M. E., Hoffman, E. J., and Mullani, N. A. (1975). A positron-emission transaxial tomograph for nuclear imaging (PETT). Radiology, 114(1), 89–98.

Yano, J.-I., Ziemiański, M. Z., Cullen, M., Termonia, P., Onvlee, J., Bengtsson, L., … Wyszogrodzki, A. A. (2018). Scientific challenges of convective-scale numerical weather prediction. Bulletin of the American Meteorological Society, 99(4), 699–710.

岡本良夫, and 本間生夫. (2002). 脳波とその逆問題解析. 電気学会論文誌. C, 電子・情報・システム部門誌, 122(9), 1417–1425.

浅井冨雄. (1967). 細胞状積雲対流の特性について. Journal of the Meteorological Society of Japan. Ser. II, 45(3), 251–260.

Footnotes

「MRI や PET では脳内における血流量や糖代謝等の空間分布が mm オーダーの分解能で直接に計測できるのに対し、EEG や MEG では頭表面上で計測される電位分布や磁場分布から脳内の起電力分布を推定するという問題、いわゆる脳電位逆問題や脳磁場逆問題を解決しなければならない。」(岡本良夫 and 本間生夫, 2002)．fMRI や陽電子断層撮像法 (PET: Positron Emission Tomography) (Ter-Pogossian et al., 1975) は深いところまで見れるが時間分解能が低いという欠点もある．↩︎
データ同化に対するレトロニムとして，逆問題を 古典的逆問題 (classical inverse problem) と呼ぶこともある (Nakamura and Potthast, 2015, pp. 1–5)．↩︎