MCMC

MCMCを実装するライブラリはいくつかあるようです。

PyMC3:Theano
PyMC4:TensorFlow
Pyro:PyTorch
NumPyro:JAX

手元の環境で動いたのはPyroのみ。環境構築に関してはシビアなようです。

web上のサンプルを見てみますと、書き方はほぼPyMC3と同じでした。動かすのは簡単（traceplot がないのは残念。seaborn で書きました）。

が、写経し終わると理解できていない箇所が露わに。

 y_model = pyro.sample('y_model', dist.Normal(a*x + b, sd), obs=y)

ここの obs の計算内容がわかりません。ここがキモなのはわかっているのですが、マニュアルに詳細が書かれていませんし、ソースも追えません。尤度×事前確率にどのようにかかわっているのでしょう？
おそらく、理解している方々が必要として使われているので、学び方が逆なのでしょうね。

ま、他の部分は理解できていることが分かりました。メモしておきましょう。

 **************************************
20201123 MCMC メモ
 
データ駆動科学

• データを出発点と考える=因果律を遡る
• 従来：変形係数を確定(モデルを確定)→応力、ひずみがノイズを含み確率的分布（→安全率を含めた設計）
• データ駆動：応力、ひずみが確定→変形係数が確率的分布（逆解析的発想）
• データｘ,y（ひずみ、応力）が得られたとき、パラメータa（変形係数）が取り得る確率分布（事後確率）を求めたい。
• 必要なのは、①パラメータaを与えたときにyが取り得る確率分布(尤度)、②aの確立分布（事前確立）

取得データD[x,y]

• 確定x：ひずみx
• 確定y：応力y

予測分布モデル　（例）物理モデル+ノイズ
　物理モデル y=ax+b

• 未確定a:平均μa,標準偏差σaとして設定(事前確率P(a)=N(μa,σa2)を仮定)
• 未確定b:平均μb,標準偏差σbとして設定(事前確率P(b)=N(μb,σb2)を仮定)
• 未確定σn:一様分布(U([n])を仮定)
  ←いずれもデータDを見て取りうる範囲を仮定
  

　ノイズ

• 重畳ノイズを正規分布として設定（尤度P(y_model)=N(ax+b,σn2)を仮定
  

予測分布モデルから θ[a:変形係数、b：切片、sd:標準偏差]の事後確率 p(θ|D)を求める。

• p(θ|D) = p(D|θ)p(θ)/p(D) ∝ p(D|θ)p(θ)
  
• 事後確率 ∝ 尤度×事前確率
• 事前確率、尤度に正規分布を仮定した場合、事後分布を解析的に求めることが可能（共役事前分布）。
• しかし、一般的には計算が困難。→MCMC法などで帰納的に事後確率（最大）を求める。

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

• モンテカルロ法：乱数を使ったシミュレーション手法
• マルコフ連鎖：現在の状態が、前時刻の状態のみに依存するモデル→メトロポリス法などの利用

1. θ0を仮定
2. 乱数εを加え、θ1を作成。
3. 確率密度の比ｒを計算r=p(D|θ1)p(θ1)/ p(D|θ0)p(θ0)
   ノイズを正規分布と仮定した場合、p(D|θ) ∝ exp(-NE(θ)/σ2)
   あるθ0、θ1における exp(-NE(θ)/σ2)p(θ) を求める
   ←obs=D を利用 ※点での推定
   
4. 確率密度の比ｒを評価
   r>1でθ1を採用
   r<=1かつ乱数U（1,0）＜ｒでθ1を採用、U>=ｒでx0を採用
5. １～４の繰り返し←exp(-NE(θ)/σ2)p(θ)の分布形状を推定

 

他手法との違い

• 最小二乗法:残差の平方和 E(θ) を最小化
  
• 最尤推定:a,b,σについて導関数=0の位置（パラメータの値）を求める(global minimum とは限らない)
  
• ベイズ推定:パラメータを確率的に扱う→分布(形)を求める 

******************************************

20201129

MAP推定：パラメータ値の同定のみの場合（パラメータの分布を求めない）

• 作成した予測分布モデル＋観測データ（obs）に対し、MAP（maximum a posterior :最大事後確率）推定でパラメータを同定
• pyro.condition で同定パラメータをモデルに反映

 

20201202

変分推論：分布も推定

• svi = SVI(model, guide, optimizer, loss=Trace_ELBO())でモデル構造を作成
  
• step method でデータDをモデルに渡して調整←obsの役目
• predictive で予測分布を作成（モンテカルロ近似）
  

概ね、obs の役割について理解しました。
数式は理解できていません。残念ながら後回しです。

 

20201204

数式もOK（おそらく）。

これで MCMC（メトロポリス法利用）は理解できたと思います。

PyMC

PyMCを使いたくて苦闘。

MCMCを理解する過程で必要となったのですが、これが動きません。

PyMC3 が動いたと教えていただいたライブラリVer.の組み合わせを真似し、仮想環境を構築したのですが、ダメでした。PIP を混ぜてPyMC4にもしてみましたが、ダメ。Anaconda を最新版に入れなおしてもダメ、gcc 関連のパスを通してもダメ。お手上げでした。

よく引っかかったのは Theano。これ、Deep Learning にも流用できる数値計算ライブラリとして有名だったのですが、既に開発が終わって時が経っています。
ちなみに、PyMC4 は TF を利用。TFも２年ほど前はよく利用されていましたが、今は PyTorch に首位の座を奪われているように感じます。

バックエンドが変更になるたびにコードが動かなくなるのは避けたいですね。Docker しかないかな？

DO CONCURRENT with GPUs

Accelerating Fortran DO CONCURRENT with GPUs and the NVIDIA HPC SDK
https://developer.nvidia.com/blog/accelerating-fortran-do-concurrent-with-gpus-and-the-nvidia-hpc-sdk/

昨日のブログです。
OpenACC のデータ転送が面倒で困っていたところ、この記事を見かけました。

 All data movement between host memory and GPU device memory is performed implicitly and automatically under the control of CUDA Unified Memory.

これでACCと同程度のパフォーマンスが出たらラッキー。頻繁にメモリへのアクセスがあるので過剰な期待はしていないのですが、いくらか早くなれば嬉しい。早く試したい！

Docker はまだ 20.9 (SDK 単体も)。公開まであと1か月くらいでしょうか？

**************************************

20201217追記

先日、Docker イメージが公開されました。
早速試してみたところ、-stdpar=multicore では OpenMP と同じ結果、速度でした。
が、-stdpar=gpu では発散。うまく計算しません。
コンパイル情報を見てみると、gpu では DO CONCURRENT 内の DO LOOP も自動並列化で GPU に載せられていました。これかな？

NVIDIA HPC SDK

OpenACC を利用しようと考えていました。

ACC の実装は比較的手軽だそうですが、場合によっては CUDA に並ぶことがあるとのこと（以前、チャレンジしたときは失敗しました）。
私はプロではないので、C & CUDA によるゴリゴリチューニングはできません。で、OpenMP からの ACC で、あわよくば狙いです。

OpenACC を gcc がサポートしているようなので使おうかな？などと考えながら、いつもの PGI Community Edition をチェック。以前のVer.と変わらないなあ、と製品版の価格をチェック。すると、以下の表示。

PGI Professional 製品の販売は終了しました。（詳細）

( ﾟДﾟ)？！
3月で販売終了していたとのこと。知りませんでした。

PGI Compilers＆Tools は NVIDIA ブランドの NVIDIA HPC SDK ソフトウェアとして、新しく生まれ変わりました。

NVIDIA HPC SDK
https://developer.nvidia.com/hpc-sdk

ついに PGI の冠が取れたのね。
プロファイラーも以前のまま含まれています。しかも、GPU Cloud で配布されています。
https://ngc.nvidia.com/catalog/containers/nvidia:nvhpc

はい、素晴らしい。さすが NVIDIA。
もうこれは、Linux ＋ Docker ＋ HPC SDK 一択でしょう。

最新のアーキテクチャに対応したコンパイラーが無償で配布されるのも素晴らしいのですが、これ、Win版の Fortran コンパイラーとしても使われだすのではないでしょうか？

俄然、やる気になりました。使ってみましょう。

**************************************
20201113追記

使ってみました。
dockerイメージを走らせただけ（nvidia-docker を使わなくなっていました）。その上で何の問題もなくコンパイルできました。
計算中に GPU を使っています。kernels 構文を差し込んだだけですので、遅いのですが。
チューニングには時間を取られそうですが、制御は拍子抜けするくらい簡単でした。ありがたいですね。

SPH 基本5：長所と短所

SPHスキームの長所と短所

圧縮性（SPH法）
長所：陽解法のため逐次代入計算により時間更新を行う。計算が簡単。
短所：粒子の重なりや抜けが生じやすい。結果として数値的不安定性を招きやすい。
→安定させるためにタイムステップを短くする必要がある。（計算負荷増大）
→人口粘性の導入が必要。ただし、適度に設定しないと訛りすぎる。

弱圧縮性（WSPH法）
長所：陽解法。疑似非圧縮性流体の計算が可能。
短所：計算の安定性のために0.5%程度の圧縮性を許容することが必要。
→体積保存性に欠陥

非圧縮性（ISPH法）
・SPH法の微分演算子とMPS法の計算アルゴリズムの組み合わせ。逆も可。半陰解法。
長所：非物理的な人工粘性を導入する必要がなく安定した計算が実行可能。
短所：圧力場の解に激しい攪乱を伴う
→インターバル平均で攪乱は消失。物理ベースCGなどで利用可能。
→防波堤に作用する衝撃波圧など瞬時の値を扱う問題には使えない。
※陽解法、半陰解法のいずれにおいても、安定性向上に有効なのが粒子再配列（高精度粒子法として分類される）
・XSPH
・Particle Shifting：J＝-v∇Ｃ
・Optimized Particle Shifting

 

SPH 基本4：圧力場の振動

弱圧縮性SPH数値スキームでは、圧力場がスプリアス振動を示します。

状態方程式を用いた場合、密度の変化が非常に小さくても、圧力の変化は指数関数的に生じます。それらは、物理的にはありえないような巨大な圧力を生じさせます。

対策
Density filtering: Shepard Filter (0th order), Moving Least Squares (1st order)
人工拡散項の追加：δ-SPH, Fourtakas et al. (2019)
Riemann solver：Godunov SPH

δ-SPH
WSPH法の数値的不安定（引張不安定）の制御手法として、δ-SPH法が近年よく用いられる。δ-SPH法では密度勾配をTaylor級数の一次精度で算定することにより拡散項の評価精度を向上し、自由表面の跳躍を抑制している。
後藤仁志「粒子法　連続体・混相流・粒状体のための計算科学」より
Molteni, D. and A. Colagrossi. (2009)A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH.
S.Marrone et al.(2011) δ-SPH model for simulating violent impact flows

Fourtakas et al. (2019)
Local uniform stencil (LUST) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models

Godunov SPH
Riemann問題の解を適用し、必要最低限の粘性を自動的に導入する手法。あらたなパラメータ設定が不要。
Inutsuka et al. (2002) Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver

SPH 基本3：精度

SPH 空間補間が、平滑化長を h とすると、 𝑂(ℎ2) に等しい精度で収束。

1Dの場合、

\begin{align*}
\left\langle A\left(r\right)\right\rangle=\int_{\mathrm{\Omega}}{A\left(r^\prime\right)W\left(r-r^\prime\right)dr^\prime}\end{align*}
A(r’) を r 周りでテーラー展開すると
\begin{align*}
A\left(r^\prime\right)\cong A\left(r\right)+\frac{\partial A\left(r\right)}{\partial r}\left(r-r^\prime\right)+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r^2}\left(r-r^\prime\right)^2+・・・\end{align*}
SPHでは

\begin{align*}
\left\langle A\left(r^\prime\right)\right\rangle=A\left(r\right)\int_{\mathrm{\Omega}}{W\left(r-r^\prime\right)dr^\prime}\end{align*}

\begin{align*}+\frac{\partial A\left(r\right)}{\partial r}\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}

\begin{align*}+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r}\int_{\mathrm{\Omega}}{\left(r-r^\prime\right)^2W\left(r-r^\prime\right)}dr^\prime\end{align*}

\begin{align*}+・・・\end{align*}

 1．正規化：影響域で積分すると１になる。

\begin{align*}
\int_{\mathrm{\Omega}}{W\left(r-r^\prime\right)dr^\prime=1}\end{align*}
 2．偶関数：中心に対し点対称

\begin{align*}
W\left(r-r^\prime\right)=W\left(r^\prime-r\right)\end{align*}

\begin{align*}
\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}

\begin{align*}=\int_{\mathrm{\Omega}}{\left(r-r^\prime\right)W\left(r^\prime-r\right)dr^\prime}\end{align*}

\begin{align*}=-\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}

\begin{align*}=0\end{align*}

3．\begin{align*}q=\frac{r-r^\prime}{h}, w\left(r-r^\prime\right)=\frac{\alpha}{h}f\left(q\right), 𝑓: ℝ+ とすると\end{align*}

\begin{align*}\alpha\int_{\mathrm{\Omega}_q} f\left(q\right)dq=1\end{align*}

\begin{align*}dr=dq・h\end{align*}

\begin{align*}
\alpha\int_{\mathrm{\Omega}_q}{\left(r-r^\prime\right)^2f\left(q\right)dq}=\alpha h^2\int_{\mathrm{\Omega}_q}{q^2f\left(q\right)dq}\end{align*}
1, 2, 3より、

\begin{align*}
\left\langle A\left(r^\prime\right)\right\rangle=A\left(r\right)+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r}\alpha h^2\int_{\mathrm{\Omega}_q}{q^2f\left(q\right)dq}+・・・\end{align*}

\begin{align*}=A\left(r\right)+O\left(h^2\right)\end{align*}

SPH 基本2：カーネル

 カーネル関数の基本的な性質

• 正規化：影響域で積分すると１になる。
• compact support：影響域が有限。区間の端では関数値が０になる。
• 非負性：kernel は定義域内では非負である。
• 単調減少性：中心から離れるに従って単調に減少する。（近接している粒子ほど強く影響しあう）。
• Dirac のδ関数への収束性：smoothing lengthが０となる極限で Kernel が関数 f(x) 自体に一致。
• 偶関数：中心に対し点対称
• 平滑性：kernel とその微分が連続かつ滑らかな関数形状を有する。（計算が安定） 

後藤仁志「粒子法　連続体・混相流・粒状体のための計算科学」より