はじめに

• MCMCの一種
• 目標: ある分布 $\pi (x)$からのサンプリングを行いたい
  1. Metropolis-Hastingsアルゴリズム (MH)
  2. Hamiltonian Monte Carlo (HMC)
  3. Langevin Dynamics (Metropolis-adjusted Langevin Algorithm)
  4. Stochastic Gradient Langevin Dynamics (SGLD)
• の順に見ていくと理解しやすい

Metropolis-Hastings

• Metropolis-Hastingsについては既知のもとする
  • 提案分布 $q(z)$を元に判定関数を用いて受容・棄却を行うMCMC
  • cf. GibbsSamplingとHM
  • 提案分布が対称, すなわち $q(z|z_{\ast })=q(z_{\ast }|z)$ならば, 判定関数はただの $min(1,\frac{\pi (z_{\ast })}{\pi (z)})$となる
    • これを単にメトロポリス法と呼ぶ

Hamiltonian Monte Carlo

• Hamiltonian Monte Carlo (HMC)

  • ここでハミルトニアンが登場する
    $$H:=H(q,p;t)=\mathcal{K}(p)+\mathcal{U}(q)$$

  • HMCのお気持ち

    • 解析力学に則り, ある系における座標と運動量の時間に関する軌道を元に, サンプリングを行う
    • 例えば, $(z,p)$があるとして, 一定時間が経ったときの終点 $(z_{\ast },p_{\ast })$が得られればサンプル取れるよね
    • このとき, 判定関数がハミルトニアンの差で表せるので嬉しい
      • ハミルトニアン自体は系全体のエネルギーそのものなので, エネルギー保存則が成り立つ.
      • ゆえにハミルトニアンの差分は本来0なので嬉しい.

  • 簡単のため, 質量を $1$としておく ($\mathcal{K}$が下記のように表せるため)
    $$\mathcal{K}(p)=\frac{1}{2}{p}^{\mathrm{\top }}\mathrm{I}p$$

  • サンプリングしたい分布 $\pi (z)\propto \tilde{\pi }(z)$について, 補助変数 $p$を導入して, $\pi (z,p):=\pi (z)=\pi (z)\pi (p)$と拡張する.

  • また $\pi (p)=\mathcal{N}(p|0,I)$とし, $\ln \tilde{\pi }(z)=-\mathcal{U}(z)$とすると,
    $$\pi (p)=\frac{1}{{(2\pi )}^{k/2}}\exp (-\frac{{\mathbf{p}}^{\mathrm{\top }}\mathrm{I}\mathbf{p}}{2})$$

  • から,
    $$\begin{array}{rl}\pi (z,p)& =\pi (z)\pi (p)\\ & =\exp (\ln \pi (z)+\ln \pi (p))\\ & \propto \exp (\ln \tilde{\pi }(z)-\frac{p^{\mathrm{\top }}\mathrm{I}p}{2})\\ & =\exp (\ln \tilde{\pi }(z)-\mathcal{U}(p))\\ & =\exp (-\mathcal{K}(p)-\mathcal{U}(p))\\ & =\exp (-H(z,p))\end{array}$$

  • となる

  • ここで, 提案分布を $\mathcal{N}(0,I)$とすると, Metropolis-Hastings→メトロポリス法となるので,

    • 1. $p$を $\mathcal{N}(0,I)$からサンプリング
    • 2. Leapfrog法で $(z,p)$から候補点 $(z_{\ast },p_{\ast })$を取得 (ステップサイズ $ϵ$・ステップ数 $L$)

  • とすれば, 判定関数 $A(z_{\ast },p_{\ast })$は
    $$A(z_{\ast },p_{\ast })=min(1,\frac{\pi (z_{\ast },p_{\ast })}{\pi (z,p)})=\exp (-H(z_{\ast },p_{\ast })+H(z,p))$$

  • となる. だが, そもそもハミルトニアンは系全体のエネルギーそのものなので時間変化せず, 常に $H(z_{\ast },p_{\ast })$と $H(z,p)$は一致するので, 判定関数 $A(z_{\ast },p_{\ast })$は常に $1$になる

    • めでたしめでたし
    • (実際はLeapfrogの誤差の影響で棄却されることもある)

  • ここからはハイパラのお話

    • Leapfrogではステップサイズ $ϵ$・ステップ数 $L$というパラメタが存在する
    • $ϵ$が小さければ探索空間が小さく非効率だし, 大きすぎると棄却される
    • なので, 基本は
      • ① $L$固定で $ϵ$動かす
      • ② $ϵ$固定で $L$動かす
    • のどちらかだが, 後者は後者で計算コストが高すぎる
    • したがって, ここらへんは慎重に選択する必要がある

Leapfrog

• 一旦ここでLeapfrogのおさらいをする
  • ハミルトニアンの正準方程式をLeapfrogで解こうとすると,
    $$\begin{array}{rl}{p}_i(t+\frac{ϵ}{2})& =p_i(t)-\frac{ϵ}{2}\frac{d\mathcal{U}}{d{z}_i}\\ z_i(t+ϵ)& =z_i(t)+ϵp_i(t+\frac{ϵ}{2})\end{array}$$
  • 一つの式にまとめると
    $$z_i(t+ϵ)=z_i(t)-\frac{{ϵ}^2}{2}\frac{d\mathcal{U}}{d{z}_i}+ϵp_i(t)$$

Langevin Dynamics

• Langevin Dynamics (LD)

  • LDは先程のHMCの特殊例である. 具体的には $L=1$としたときのHMCである

  • LDにおいて, あるモデルのパラメタ $\theta$に関して, 確率分布 $\pi (\theta )$からのサンプリングを考えてみる

  • 上のLeapfrogの式より
    $$\begin{array}{rl}\mathrm{\Delta }\theta :={\theta }_{t+1}-{\theta }_t& =-\frac{{ϵ}^2}{2}\mathrm{\nabla }\mathcal{U}+ϵp_i(t)\\ & =\frac{{ϵ}^2}{2}\mathrm{\nabla }\ln \pi (\theta )+ϵp_i(t)(\therefore \ln \tilde{\pi }(\theta )=-\mathcal{U}(\theta ))\end{array}$$

  • ここで, データ $(X,Y)$について, $\pi (\theta )\leftarrow \pi (\theta |Y,X)$であるとき, ベイズより
    $$\begin{array}{rl}\mathrm{\Delta }\theta & =\frac{{ϵ}^2}{2}\mathrm{\nabla }\ln \pi (\theta |Y,X)+ϵp_t\\ & =\frac{{ϵ}^2}{2}(\mathrm{\nabla }\ln \pi (Y|\theta ,X)+\mathrm{\nabla }\ln \pi (\theta ))+ϵp_t\\ & =\frac{{ϵ}^2}{2}(\sum _n\mathrm{\nabla }\ln \pi (y_n|\theta ,x_n)+\mathrm{\nabla }\ln \pi (\theta ))+ϵp_t\end{array}$$

  • ただし, $p_t$は先程の説明通り $\mathcal{N}(0,I)$からサンプリング

  • なので, これはノイズを加えた勾配法と同じようになる

Stochastic Gradient Langevin Dynamics

• Stochastic Gradient Langevin Dynamics (SGLD)
  • お待ちかねのSGLDだが, これは単純に上の式の第一項をサブサンプリングするだけ
  • なので更新幅は
    $$\mathrm{\Delta }\theta =\frac{{ϵ}^2}{2}(\frac{N}{M}\sum _{n=1}^M\mathrm{\nabla }\ln \pi (y_n|\theta ,x_n)+\mathrm{\nabla }\ln \pi (\theta ))+ϵp$$
  • これで, 棄却率がほぼ0 ∧ 計算量激減のMCMCが出来上がったことになる