2024-04-19

ハミルトニアンモンテカルロ法を今更理解する

趣旨

ハミルトニアンモンテカルロ法(HMC)を理解する(執筆当時ではまだ理解できてないけど)。

注意

物理ど素人なので完全には理解できてない。情報の正確性には注意。

理解してないなら「理解する」とか書くなよ

うるせーよ

では解説する。

HMCとは

Metropolis-Hastingsに代表されるこれまでのMCMCでは、あらかじめ与えられた初期値と初期提案分布からスタートし、棄却したりしなかったりしながらサンプリングを行うわけだが、提案分布が目的の分布と大きく外れていれば酷い出来になる。これを軽減するためにランダムウォークMetropolis-Hastingsがあるらしいのだが、あんまり探索的な(提案分布から推察すれば生起確率が低そうな)サンプリングを行うと高い確率で棄却されてしまう。要するに、HMCは強化学習における「探索と活用のトレードオフ」に似た問題を抱えている。これをハミルトン力学の観点から解決したのがHMC、らしい。つまり、ある程度探索的なサンプリングを確保しながら、同時に採択率を高く保てる方法がHMCなのである。

1. 初等物理量

ハミルトン「力学」とあるとおり、HMCの起源は物理学、より厳密には解析力学にある。よって、HMCの理解にはある程度の物理的知識が必要であるため、一旦物理のおさらいをしたい。ここで、時間を $\tau$ で表現することにする。

1-1. 位置、速度、加速度

空間上を動く物体を考える。この物体の時刻 $\tau$ における位置をベクトル $\boldsymbol{r} (\tau)$ で表現する。また、同時刻のおける物体の速度 $\boldsymbol{v} (\tau)$ は、速度が瞬間的な位置の変化を表現していることに注意すれば、

$\begin{align} \boldsymbol{v} (\tau) = \frac{d \boldsymbol{r} (\tau)}{d \tau} \end{align}$

と表現できるだろう。また、加速度 $\boldsymbol{a} (\tau)$ も同様に考えれば、

$\begin{align} \boldsymbol{a} (\tau) = \frac{d \boldsymbol{v} (\tau)}{d \tau} = \frac{d^2 \boldsymbol{r} (\tau)}{d \tau^2} \end{align}$

と書ける。

1-2. 質量、運動量、力

物体の質量を $m$ としたとき、時刻 $\tau$ における運動量 $\boldsymbol{p} (\tau)$ を

$\begin{align} \boldsymbol{p} (\tau) = m \boldsymbol{v} (\tau) \end{align}$

と定義する。これは誤解を恐れずに言うなら「物体の止めにくさ」を意味する。

例え高速で飛来する物体でも、それが野球ボールなら簡単に止められるだろう。一方で、ゆっくり動いていても、それが船ならまず人力では止められない。この辺りの違いを数値的に表現したのが運動量である。

また、運動量の変化を力 $\boldsymbol{F} (\tau)$ という。

$\begin{align} \boldsymbol{F} (\tau) = \frac{d \boldsymbol{p} (\tau)}{d \tau} \end{align}$

この表現と上の定義から、運動方程式と呼ばれる式が成立する。

$\begin{align} \boldsymbol{F} (\tau) = m \frac{d^2 \boldsymbol{r} (\tau)}{d \tau^2} = m \boldsymbol{a} (\tau) \end{align}$

この式は高校物理でもよく見るだろう。速度が光速に近いと成り立たないらしいが、そんなことは知らん。

1-3. 仕事

仕事とは、物体をある位置からある位置に動かすのに要するエネルギーのことである。具体的には仕事 $W$ は力の積分によって表される: 物体の移動経路を $C$ で表現すると、

$\begin{align} W = \int_C \boldsymbol{F} d \boldsymbol{r} \end{align}$

と表現される。ここで、積分は $\boldsymbol{r}$ について行われることに注意: 直感的には時間 $\tau$ を用いても良いような気がするが、仕事に置いて時間はさほど重要ではない。例えどれだけ時間をかけたとしても、結果として物体がほとんど動かなかったら、仕事としては小さい。極端な話、動かなかったら仕事は $0$ である。シビアやね。

上の式にあるように、仕事は経路 $C$ に依存するため、始点と終点以外の部分も考慮される: 両端点が同じでも、最短経路で動かした場合と、大回りしてたどり着いた場合では仕事の値が違う可能性がある。ここで、仕事が経路 $C$ に依存しない時、力 $\boldsymbol{F} (\tau)$ を保存力と呼ぶことにしよう。 $\boldsymbol{F} (\tau)$ が保存力であるという仮定のもとでは、先ほどの式は端点 $\boldsymbol{r}_s, \boldsymbol{r}_g$ のみを考慮した積分の式

$\begin{align} W = \int_{\boldsymbol{r}_s}^{\boldsymbol{r}_g} \boldsymbol{F} d \boldsymbol{r} \end{align}$

で表される。一般に仕事を評価する時、それは地点 $\boldsymbol{r} (0)$ から地点 $\boldsymbol{r} (T)$ まで、つまり時刻 $0$ から $T$ までの仕事を評価することがよくある。その場合は

$\begin{align} W(0, T) &= \int_{\boldsymbol{r} (0)}^{\boldsymbol{r} (T)} \boldsymbol{F} d \boldsymbol{r} \\\ &= \int_0^T \boldsymbol{F} (\tau) \frac{d \boldsymbol{r} (\tau)}{d \tau} d \tau \\\ &= \int_0^T \boldsymbol{F} (\tau) \boldsymbol{v} (\tau) d \tau \\\ \end{align}$

のように書ける。

2. 力学的エネルギー

運動する物体にはエネルギーを持っている。物体が存在する空間には摩擦や空気抵抗などが存在しないと仮定すると、物体が持つエネルギーはポテンシャルエネルギーと運動エネルギーの2種類に大別される。

2-1. ポテンシャルエネルギー

ポテンシャルエネルギーとは、物体が特定の位置に存在することによって潜在的に持つエネルギーで、位置エネルギーとも言う。保存力 $\boldsymbol{F}$ に対して、基準点 $\boldsymbol{r}_0$ を定めれば、以下のようなポテンシャルエネルギーを定義できる。

$\begin{align} U(\boldsymbol{r}) = -\int_{\boldsymbol{r}_0}^{\boldsymbol{r}} \boldsymbol{F} d \boldsymbol{r}\\\ \end{align}$

ここでもやはり、基準点を $\boldsymbol{r} (0)$ 、そして現在の位置を $\boldsymbol{r} (T)$ としてみると、

$\begin{align} U(\boldsymbol{r} (T)) = -\int_{\boldsymbol{r} (0)}^{\boldsymbol{r} (T)} \boldsymbol{F} d \boldsymbol{r} = -W(0, T) \\\ \end{align}$

と書ける。

例を挙げよう。今、鉛直上向きに $x$ 軸が伸びる座標系を考え、またこの系では鉛直下向きに質量に比例した力 $mg$ が与えられるものとする(要するに、高校力学と同じ環境)。この時、基準点を $0$ としたときの位置 $h$ にある物体のポテンシャルエネルギーは、力の向きに注意すると、

$\begin{align} U(h) = -\int_0^h (-mg) d x = mgh\\\ \end{align}$

と書ける。お馴染みの位置エネルギーである。

別の例として、 $x$ 軸だけが伸びる座標系において、基準点を $0$ としたとき位置 $h$ にある物体に力 $kh$ が与えられるものとする。この場合、

$\begin{align} U(h) = -\int_0^h kx d x = \frac{1}{2} k h^2\\\ \end{align}$

のようにポテンシャルエネルギーが与えられ、これはばねの弾性エネルギーと等価である。ここでいうポテンシャルエネルギーとは、高校物理における位置エネルギーより広い範囲のエネルギーであると言える。

2-2. 運動エネルギー

運動エネルギーは、以下の式で与えられるエネルギーである。

$\begin{align} K(\boldsymbol{p}) = \frac{1}{2} m {||\boldsymbol{v}||}^2 = \frac{1}{2m} {||\boldsymbol{p}||}^2 \\\ \end{align}$

さて、この運動エネルギーについての考察を得るため、両辺を $\tau$ で微分してみる。

$\begin{align} \frac{d K(\boldsymbol{p})}{d \tau} = m \boldsymbol{v} \cdot \frac{d \boldsymbol{v}}{d \tau} = m \boldsymbol{a} \cdot \frac{d \boldsymbol{r}}{d \tau} = \boldsymbol{F} (\tau) \boldsymbol{v} (\tau)\\\ \end{align}$

さらにこれを区間 $(0, T)$ で $\tau$ について積分してみると、

$\begin{align} K(\boldsymbol{p} (T)) - K(\boldsymbol{p} (0)) = \int_0^T \boldsymbol{F} (\tau) \boldsymbol{v} (\tau) d \tau = W(0, T) \\\ \end{align}$

が得られる。

2-3 エネルギー保存則

$\boldsymbol{F}$ が保存力である場合のポテンシャルエネルギー $U(\boldsymbol{r} (T))$ について、

$\begin{align} U(\boldsymbol{r} (T)) = -\int_{\boldsymbol{r} (0)}^{\boldsymbol{r} (T)} \boldsymbol{F} d \boldsymbol{r} = -W(0, T) \\\ \end{align}$

と書けることは上で確認した。ここで、この記法に従えば、 $U(\boldsymbol{r} (0)) = 0$ であることは容易に確認できるため、あえて

$\begin{align} U(\boldsymbol{r} (T)) - U(\boldsymbol{r} (0)) = -W(0, T) \\\ \end{align}$

と記載しておく。さて運動エネルギーの方はというと、

$\begin{align} K(\boldsymbol{p} (T)) - K(\boldsymbol{p} (0)) = W(0, T) \\\ \end{align}$

が成立しているのだった。2つの式の右辺を見てやると、

$\begin{align} U(\boldsymbol{r} (T)) - U(\boldsymbol{r} (0)) = -\lbrace K(\boldsymbol{p} (T)) - K(\boldsymbol{p} (0)) \rbrace \\\ \end{align}$

$\begin{align} \therefore U(\boldsymbol{r} (T)) + K(\boldsymbol{p} (T)) = U(\boldsymbol{r} (0)) + K(\boldsymbol{p} (0)) \\\ \end{align}$

が成立することがわかる。これをエネルギー保存則という。この法則は、ポテンシャルエネルギーと運動エネルギーの和は時間推移によって変化しないことを主張している。

2-4. ハミルトニアン、運動方程式

ここで、ハミルトニアン $H(\tau)$ を以下のように定義する。

$\begin{align} H(\tau) = U(\boldsymbol{r} (\tau)) + K(\boldsymbol{p} (\tau)) \end{align}$

このとき、 $H(\tau)$ は時間変化によらずに一定であると言える。

ここで、「物体」は空間上の曲面を転がっているとする:「物体」が存在する空間は位置に関わらず質量と加速度 $m, g$ に比例した力 $mg$ が与えられ、かつ「物体」の基準点からの高さ $h(\tau)$ は物体の(超)平面上の位置 $\boldsymbol{\theta} (\tau)$ に依存するものとする。この時のポテンシャルエネルギーは、

$\begin{align} U(h(\tau)) &= -\int_0^{h(\tau)} (-mg) d x = h(\tau) = mgh(\boldsymbol{\theta} (\tau)) \\\ \end{align}$

と書ける。ここで、 $U(h(\tau)) = U(\boldsymbol{\theta} (\tau))$ と表現しよう。さて、この時のハミルトニアンについて何が言えるだろうか。

曲面を転がる物体は、時間の変化によらずに、常にハミルトニアンが一定である。すなわち、

$\begin{align} \frac{d H(\tau)}{d \tau} &= 0 \\\ \end{align}$

$\begin{align} \frac{d U(\boldsymbol{\theta} (\tau))}{d \tau} + \frac{d K(\boldsymbol{p} (\tau))}{d \tau} &= 0 \end{align}$

これを変形する。

$\begin{align} \frac{d K(\boldsymbol{p} (\tau))}{d \tau} &= - \frac{d U(\boldsymbol{\theta} (\tau))}{d \tau} \\\ \frac{d K(\boldsymbol{p} (\tau))}{d \boldsymbol{p} (\tau)} \frac{d \boldsymbol{p} (\tau)}{d \tau} &= - \frac{d U(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta} (\tau)} \frac{d \boldsymbol{\theta} (\tau)}{d \tau} \\\ \end{align}$

で、なんか知らんが(!?)ここから以下のような方程式が言えるらしい(というか、下の方程式を満たすとき、上の方程式を満たすような条件が見つかるのだろう)。

$\begin{align} \frac{d \boldsymbol{p} (\tau)}{d \tau} &= - \frac{d U(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta} (\tau)} \\\ \frac{d \boldsymbol{\theta} (\tau)}{d \tau} &= \frac{d K(\boldsymbol{p} (\tau))}{d \boldsymbol{p} (\tau)} \\\ \end{align}$

これをハミルトンの運動方程式と言い、物体の運動を予測する微分方程式である。

2-5. オイラー法とリープフロッグ法

上で定義したハミルトンの運動方程式を用いて、ハミルトニアンの値を変化させないまま運動量 $\boldsymbol{p}$ や位置 $\boldsymbol{\theta}$ を変化させる、つまり物体の移動を表現することを考える。ここで、計算を簡略化するためにいささか恣意的な仮定を置く: 物体の質量と重力加速度について $m=1, g=1$ とする。これにより、

$\begin{align} \frac{d \boldsymbol{p} (\tau)}{d \tau} &= - \frac{d h(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta} (\tau)} \\\ \frac{d \boldsymbol{\theta} (\tau)}{d \tau} &= \boldsymbol{p} (\tau) \\\ \end{align}$

のように簡潔に書ける。

さて、当然解析的には解けるはずがないのだが、数値的には解けそうである。オイラー法と呼ばれる方法は、上に対する素直な解法を提供してくれている。

$\begin{align} \boldsymbol{p} (\tau+ \varepsilon) &= p(\tau) + \varepsilon \frac{d \boldsymbol{p} (\tau)}{d \tau} = p(\tau) - \varepsilon \frac{d h(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta}} \\\ \boldsymbol{\theta} (\tau+ \varepsilon) &= \boldsymbol{\theta} (\tau) + \varepsilon \frac{d \boldsymbol{\theta} (\tau)}{d \tau} = \boldsymbol{\theta} (\tau) + \varepsilon \boldsymbol{p} (\tau) \\\ \end{align}$

$\cdots$ まあ要するに最急降下法なのだが、あんまり精度が良くないらしい。なので、以下のような更新式によって精度を高めているらしい。

$\begin{align} \boldsymbol{p} (\tau+\frac{ \varepsilon}{2}) &= p(\tau) - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta}} \\\ \boldsymbol{\theta} (\tau+ \varepsilon) &= \boldsymbol{\theta} (\tau) + \varepsilon \boldsymbol{p} \left( \tau + \frac{ \varepsilon}{2} \right) \\\ \boldsymbol{p} (\tau+ \varepsilon) &= \boldsymbol{p} \left( \tau+\frac{ \varepsilon}{2} \right) - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta} (\tau+ \varepsilon))}{d \boldsymbol{\theta}} \\\ \end{align}$

これをリープフロッグ法と呼ぶ。おそらく位置変数 $\boldsymbol{\theta}$ と運動量変数 $\boldsymbol{p}$ を時刻 $\varepsilon/2$ だけずらしながら交互に更新する様をカエルに見立てたのだろう。

オイラー法でもリープフロッグ法でも、結局積分を数値積分で近似しているだけなので、誤差が生じうる。これのハンドリングは後述。

2-6. 位相空間

物理的な空間内にある物体の動きを考える際、時間 $\tau$ について考えることは重要であるが、この後統計学へ応用する際においてはさほど重要な概念ではない。この視点のもとでのエネルギーは、 $m=1, g=1$ に注意すると、

$\begin{align} U(\boldsymbol{\theta}) &= h (\boldsymbol{\theta}) \\\ K(\boldsymbol{p}) &= \frac{1}{2} {|| \boldsymbol{p} ||}^2 \end{align}$

のように表現できる。このときのハミルトニアンは

$\begin{align} H(\boldsymbol{\theta}, \boldsymbol{p}) = U(\boldsymbol{\theta}) + K(\boldsymbol{p}) = h (\boldsymbol{\theta}) + \frac{1}{2} {|| \boldsymbol{p} ||}^2 \end{align}$

のように位置と運動量の関数として表現される。

位置と運動量を座標とする空間を位相空間(phase space)と言う。

数学における位相空間(topological space)とはまた違う。いい加減にしろよ。

で、ハミルトニアンを固定した時、 $(\boldsymbol{\theta}, \boldsymbol{p})$ は特定の軌跡(これを特に等高線と呼ぶようだ)を描く。実空間内の物体は、常にこの等高線に沿って位置と運動量を変化させることになる。

さて、位相空間では2つの重要な性質がある。

可逆: 運動している物体に対して、任意の時間で運動を停止させ、これまでに受けてきた力を逆向きに全く同じ大きさで(つまり $-1$ 倍の運動量を)与えると、物体はこれまでたどってきた軌跡を正確に辿る
体積保存: 位相空間上の領域(点の集合)は、時間発展によって(各点が好き勝手動いて)形を変えることはあっても、その面積(体積とよく表記される)は不変である

可逆については直感的だろう。体積保存の性質は、リウビルの定理(Liouville's Theorem)と呼ばれる定理がその成立を保証してくれている。なぜ成立するかは他のサイトに譲るとして、上の2つの性質によって以下のような重要な式が導かれる。

$\begin{align} \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}, \boldsymbol{p}) = \text{Pr} (\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \end{align}$

これは、地点 $\boldsymbol{\theta}, \boldsymbol{p}$ から地点 $\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}$ に遷移する確率は、 $\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}$ から $\boldsymbol{\theta}, \boldsymbol{p}$ に遷移する確率と等しいことを意味している。

可逆性だけではこの式は導けない。なぜなら、物体に運動量 $\boldsymbol{p}$ を与える確率と、 $-\boldsymbol{p}$ を与える確率は等しいとは言えないからだ。リウビルの定理がその辺りを説き伏せてくれるのだろう。どこにも説明がない。マジでキレそう。

3. MCMCに求められる条件

ハミルトニアンの話を一旦脇に置いて、MCMCのおさらいをしよう。今我々が提案するマルコフ連鎖では、状態 $z$ の分布が目的の分布に収束する必要がある。

3-1. エルゴード性

まずは一意の分布に収束することについて考えよう。この収束先の分布のことを定常分布という。より具体的には、以下の数式を満たす分布 $p^* (z)$ のことを指す。

$\begin{align} p^* (z) = \int p(z | z') p^* (z') dz' \end{align}$

特定の状態にハマると抜け出せないようなマルコフ連鎖などは、この条件を満たさない。定常分布を持つための必要十分条件は、以下の3つの性質を全て満たすことである。

既約的: 任意の状態から他の任意の状態へ到達可能
非周期的: 周期性を持たない
正再帰的: 任意の状態は限りなく何度も訪問される

これら3つの性質を満たすとき、そのマルコフ連鎖はエルゴード性を持つという。

3-2. 詳細釣り合い条件

仮にマルコフ連鎖がエルゴード的であったとして、定常分布が目的の分布と全く違っては意味がない。ある確率分布 $p^* (z)$ が定常分布であるための必要十分条件は分かっていないが、十分条件は詳細釣り合い条件として知られる。

$\begin{align} p^* (z) p(z' | z) = p^* (z') p(z | z') \end{align}$

4. HMC

さて、ハミルトン力学の考え方をMCMCに導入することを考える。

4-1. ハミルトニアンの導入

今、目的の分布である事後分布 $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ と、それとは独立の標準正規分布 $\text{Pr}(\boldsymbol{p})$ が存在するとする。同時分布は

$\begin{align} \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) = \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) \text{Pr}(\boldsymbol{p}) \end{align}$

のように表現される。これの対数を取ると、

$\begin{align} \log \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) &= \log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) + \log \text{Pr}(\boldsymbol{p}) \\\ \end{align}$

と書ける。右辺の各項について考察を深めよう。

$\cdots$ と言っておいてアレだが、右辺第一項 $\log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ については仮定を置くだけである。ハミルトン運動方程式を導いた際に、物体は高さが位置に依存した関数 $h(\boldsymbol{\theta})$ で表現されるような曲面上を動くと書いた。ここでは、

$\begin{align} h(\boldsymbol{\theta}) = - \log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) \end{align}$

と置く。

右辺第二項 $\log \text{Pr}(\boldsymbol{p})$ については具体的に計算してみる。 $\text{Pr}(\boldsymbol{p})$ が標準正規分布であったことに注意すれば、

$\begin{align} \log \text{Pr}(\boldsymbol{p}) &= \log \left\lbrace \frac{1}{\sqrt{2 \pi}} \exp \left( -\frac{1}{2} {||\boldsymbol{p}||}^2 \right) \right\rbrace \\\ &= -\frac{1}{2} {||\boldsymbol{p}||}^2 + \text{Const.} \\\ \end{align}$

と書ける。これらを全部代入すると、各項が2章6節で定義したエネルギーと同じ形をしているため、

$\begin{align} \log \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) &= -h(\boldsymbol{\theta}) -\frac{1}{2} {||\boldsymbol{p}||}^2 + \text{Const.} \\\ &= -U(\boldsymbol{\theta}) -K(\boldsymbol{p}) + \text{Const.} \\\ &= -H(\boldsymbol{\theta}, \boldsymbol{p}) + \text{Const.} \end{align}$

のように、突然ハミルトニアンが出現する。要するに、

$\begin{align} \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) = \frac{1}{Z_H} \exp (-H(\boldsymbol{\theta}, \boldsymbol{p})) \end{align}$

のように書ける。ここで、定数項の影響は正規化定数 $1 / Z_H$ で吸収した。

4-2. 運動の解釈

曲面上の物体は、重力に従って高さ $h(\boldsymbol{\theta})$ が小さくなる方向に運動することは想像に難くないだろう。前節によれば、高さは

$\begin{align} h(\boldsymbol{\theta}) = - \log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) \end{align}$

のように定義されているので、高さ $h(\boldsymbol{\theta})$ が小さくなるということは、事後分布 $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ を最大化することに等しい。

ただし実際は物体は放置されるわけではない。前節で、標準正規分布に従って $\boldsymbol{p}$ がサンプリングされることを仮定した。 $\boldsymbol{p}$ とは運動量のことであるから、物体はただ重力に従って曲面上を転げ落ちるのではなく、ランダムな力を常に加えられながら転がることを意味する。

定期的に物体を止めて、ランダムな方向と力で物体を弾くとイメージすれば良い。

このランダムネスによって、有効なサンプリングが期待される。

4-3. 詳細釣り合い条件

ここで提案されているHMCは、要するに以下の処理を交互に行うマルコフ連鎖であると言える。

運動量 $\boldsymbol{p}$ の確率的な更新
2章4節で定義したハミルトン運動方程式に基づいた、力学系の更新

3章2節で、マルコフ連鎖は詳細釣り合い条件を満足することが求められると説明した。ハミルトン運動方程式に従って $\boldsymbol{\theta}, \boldsymbol{p}$ を変化させた場合、詳細釣り合い条件は満たされるのだろうか。ハミルトン力学における記法に則れば、詳細釣り合い条件は以下のように表現される。

$\begin{align} \text{Pr} (\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}, \boldsymbol{p}) = \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{x}) \text{Pr} (\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \end{align}$

2章6節で、位相空間においては以下の性質を満たすと書いたことを思い出して欲しい。

さらにここでハミルトニアンを保存しながら遷移することから、

$\begin{align} \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) = \text{Pr}(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{x}) = \frac{1}{Z_H} \exp (-H(\boldsymbol{\theta}, \boldsymbol{p})) \end{align}$

が成立する。以上の二つを組み合わせると詳細釣り合い条件が成立する。

4-4. 数値誤差と棄却

上の話は理想的な運動をシミュレートできればの話である。2章5節で紹介したようなハミルトンの運動方程式を解く方法は所詮は近似的なものであり、そこでは数値誤差、つまりハミルトニアンの変動が発生する。そのため、Metropolis-Hastingsアルゴリズムで登場した棄却サンプリングを適用する。

usapyoi.hatenablog.com

このアルゴリズムでは、分布 $q_k (\boldsymbol{z}^{(\tau)})$ からサンプリングしてきた $\boldsymbol{z}^{*}$ に対して、以下のような確率をもって採択しているのだった。

$\begin{align} A_k (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)}) = \min \left\lbrace 1, \frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \right\rbrace \end{align}$

今回の例では、目的の分布は $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ であり、状態遷移分布は $\text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)})$ のように表現される。後者の分布は対称的であることに注意すると、

$\begin{align} \frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} &\rightarrow \frac{\text{Pr}(\boldsymbol{\theta}^* | \boldsymbol{x}) \text{Pr} ( \boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)} | \boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) }{\text{Pr}(\boldsymbol{\theta}^{(\tau)} | \boldsymbol{x}) \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)})} \\\ &= \frac{\text{Pr}(\boldsymbol{\theta}^* | \boldsymbol{x}) }{\text{Pr}(\boldsymbol{\theta}^{(\tau)} | \boldsymbol{x})} \\\ &= \frac{ \text{Pr}(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{x}) }{ \text{Pr}(\boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)} | \boldsymbol{x})} \\\ &= \exp \left( H(\boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)}) - H(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \right) \end{align}$

と書ける。ただし、途中で $\text{Pr}(\boldsymbol{p})$ を分母分子にかけている。

結局、数値誤差によるズレは、以下のような確率で採択するMHアルゴリズム(より正確にはMetropolisアルゴリズム)でカバーするのであった。

$\begin{align} r^{(\tau)} = \min \left\lbrace 1, \exp \left( H(\boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)}) - H(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \right) \right\rbrace \end{align}$

$\exp$ の中身は実際はほぼ $0$ であり、すなわち採択率はほぼ $1$ になる。HMCの効率の良さはここから来てるんやね。ちなみにここでハミルトン力学とMetropolisアルゴリズムが融合するので、この手法はHybrid MCと呼ばれたりもする。どっちにしろHMCである。

4-4. 結局

RMC法におけるアルゴリズムは、以下のように表現される。

初期値 $\boldsymbol{\theta}^{(0)}$ 、ハイパーパラメタ $\varepsilon$ を定める。
ポテンシャルエネルギー $h(\boldsymbol{\theta})$ を目的分布 $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ の負の対数として定義する。
$\tau = 1, 2, \cdots, T$ において、以下のSTEP 4, 5, 6, 7を繰り返す。
$\boldsymbol{p}^{(\tau-1)}$ を標準正規分布 $\mathcal{N}(0,1)$ からサンプリングする。
後に示すリープフロッグ法を用いて $\boldsymbol{\theta}^{*}$ を計算する。
確率に従って採択する
- 採択された場合は、 $\boldsymbol{\theta}^{(\tau)} \leftarrow \boldsymbol{\theta}^{*}$ とする。
- 棄却された場合は、 $\boldsymbol{\theta}^{(\tau)} \leftarrow \boldsymbol{\theta}^{(\tau-1)}$ とする。つまりやり直し。
得られた $\boldsymbol{\theta}^{(\tau)}$ をサンプルとして保存する。

リープフロッグ法は以下の通り。

$\begin{align} \boldsymbol{p}^{(\tau-1/2)} &= \boldsymbol{p}^{(\tau-1)} - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta} ^{(\tau-1)})}{d \boldsymbol{\theta}} \\\ \boldsymbol{\theta}^{*} &= \boldsymbol{\theta}^{(\tau-1)} + \varepsilon \boldsymbol{p}^{( \tau - 1/2)} \\\ \boldsymbol{p}^{*} &= \boldsymbol{p}^{( \tau-1/2 )} - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta}^{*})}{d \boldsymbol{\theta}} \end{align}$

また採択率は以下の通りとなる。

$\begin{align} r^{(\tau)} = \min \left\lbrace 1, \exp \left( H(\boldsymbol{\theta}^{(\tau-1)}, \boldsymbol{p}^{(\tau-1)}) - H(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \right) \right\rbrace \end{align}$

まとめ

HMCについてまとめた。まとめられてないけど。多分一生腹落ちすることはないと思う。実装してみたい。

おまけ

かえるぴょこぴょこ。これまで変分ベイズの記事が最大文字数(1.5万文字ぐらい)だったが、この記事はそれを大きく上回り2.1万字となりました。

変分ベイズも見てやってね。

usapyoi.hatenablog.com

2024-04-13

ギブスサンプリングを今更理解する

趣旨

MCMCの一つであるギブスサンプリングを理解する。

MCMCって？

前の記事を参考にしてください。

usapyoi.hatenablog.com

では解説する。

ギブスサンプリングとは

MCMCの一つで、Metropolis-Hastingsアルゴリズムの特殊な例と言える。

問題設定

今、区間 $(0,1)$ で一様に分布する擬似乱数を発生させるアルゴリズムが得られている。これを元に、 $p(\boldsymbol{z})$ に従う $L$ 個のサンプル $\boldsymbol{z}^{(1)}, \boldsymbol{z}^{(2)}, \cdots, \boldsymbol{z}^{(L)}$ を抽出したい。ここで、 $\boldsymbol{z}$ は $D$ 次元変数 $\boldsymbol{z} = (z_1, z_2, \cdots, z_D)$ であるとする。

ギブスサンプリング

1. 手順

$\boldsymbol{z} = (z_1, z_2, \cdots, z_D)$ に初期値 $\boldsymbol{z}^{(0)} = (z_1^{(0)}, z_2^{(0)}, \cdots, z_D^{(0)})$ を与える。
各 $\tau = 1, 2, \cdots, T$ について、以下のSTEP 3, 4, 5を行う。
各 $d=1, 2, \cdots, D$ について、以下のSTEP 4, 5を行う。
$z_d^{(\tau)}$ を $p(z_d | \boldsymbol{z}_{\setminus d}^{(\tau-1)})$ からサンプリングする.。ただし $\boldsymbol{z}_{\setminus d}$ は $z_d$ 以外の $\boldsymbol{z}$ の要素。
得られた $z_d^{(\tau)}$ と $\boldsymbol{z}_{\setminus i}^{(\tau-1)}$ を組み合わせたもの $\boldsymbol{z}^*$ をサンプルとして保存する。

この手順によって得られる $\boldsymbol{z}$ は $DT$ なので、 $T = L/D$ とすれば最短のイテレーションで済む。ただし、初期値依存性を排除するために、STEP 5を行わない、いわゆる「burn-in」を設けることが必要。

2. なぜ上手くいくのか

ギブスサンプリングはMetropolis-Hastingsアルゴリズムの特殊な例である、と先に述べた。そうであることを示せば、あとは以前の記事で紹介した通りうまくいくことが示される。

2-1. 復習 Metropolis-Hastingsアルゴリズムの手順

現在の状態が $\boldsymbol{z}^{(\tau)}$ であるとき、分布 $q_k(\boldsymbol{z} | \boldsymbol{z}^{(\tau)})$ からサンプル $\boldsymbol{z}^*$ を抽出し、確率 $A_k (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})$ に従って採択する。ただし、

$\begin{align} A_k (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)}) = \min \left( 1, \frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \right) \end{align}$

また、状態 $\boldsymbol{z}^{(\tau)}$ の更新則は、

$\boldsymbol{z}^{(\tau + 1)} = \left\{ \begin{array}{cc} \boldsymbol{z}^* & \text{ if } \boldsymbol{z}^* \text{ is accepted} \\ \boldsymbol{z}^{(\tau)} & \text{ otherwise} \end{array} \right.$

で与えられる。提案分布はこのマルコフ連鎖が既約で非周期的であれば良い。十分条件として、任意の $\boldsymbol{z}_1, \boldsymbol{z}_2$ に対して $q(\boldsymbol{z}_2 | \boldsymbol{z}_1) > 0$ であることが挙げられる。

2-2. ギブスサンプリングとの対応関係

ギブスサンプリングでは、目的の(正規化していない)分布 $\tilde{p}(\boldsymbol{z})$ や提案分布 $q(\boldsymbol{z})$ は、どちらも目的分布 $p(\boldsymbol{z})$ になっている。また、ギブスサンプリングにおけるサンプル $\boldsymbol{z}^*$ の生成では、一つの要素 $z_k$ のみが新たに生成され、残る要素 $\boldsymbol{z}_{\setminus k}^*$ は直前の状態、つまり $\boldsymbol{z}^{(\tau)}$ を流用する。まとめると、

$\begin{align} \tilde{p}(\boldsymbol{z}^*) &= p(\boldsymbol{z}^*) = p(z_k^* | \boldsymbol{z}_{\setminus k}^{(\tau)}) p(\boldsymbol{z}_{\setminus k}^{(\tau)}) \\\ \tilde{p}(\boldsymbol{z}^{(\tau)}) &= p(\boldsymbol{z}^{(\tau)}) = p(z_k^{(\tau)} | \boldsymbol{z}_{\setminus k}^{(\tau)}) p(\boldsymbol{z}_{\setminus k}^{(\tau)}) \\\ q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)}) &= p (z_k^*| \boldsymbol{z}_{\setminus k}^{(\tau)}) \\\ q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*) &= p (z_k^{(\tau)}| \boldsymbol{z}_{\setminus k}^*) \\\ \end{align}$

この条件下でサンプリングを行った時、サンプル $\boldsymbol{z}^*$ の採択確率は

$\frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} = 1$

となる。よって、ギブスサンプリングは、Metropolis-Hastingsアルゴリズムの、サンプルが必ず採択されるような特殊な例であると言える。

まとめ

ギブスサンプリングについて勉強した。言うほど難しくはないっすね。PRMLではギブスサンプリングの欠点や変種などを紹介していたが、ここでは触れない。

おまけ

割と最近までギブソンサンプリングだと思っていたことを告白します。

2024-04-09

MCMCを今更理解する

趣旨

MCMC(マルコフ連鎖モンテカルロ)の理解を目的に、基本的なMetropolis-Hastingsアルゴリズムをまとめる。

なぜMCMC?

~~名前が可愛いから~~ 情報学徒として当然だよね。

では解説する。

MCMCとは

サンプリングを行うアルゴリズムの群。MCMCという特定のアルゴリズムがあるわけではなく、あくまで大雑把な指針。

背景

期待値を評価したい。

$\mathbb{E} \lbrack f \rbrack = \int f(\boldsymbol{z}) p(\boldsymbol{z}) d \boldsymbol{z}$

しかし直接計算できないとする。

サンプリング法

サンプリング法とは、標本平均で期待値を近似しようというアプローチ。分布 $p(\boldsymbol{z})$ に従う $L$ 個のサンプル $\boldsymbol{z}^{(1)}, \boldsymbol{z}^{(2)}, \cdots, \boldsymbol{z}^{(L)}$ を抽出して、

$\mathbb{E} \lbrack f \rbrack \approx \frac{1}{L} \sum_{l=1}^L f(\boldsymbol{z}^{(l)})$

と近似する。このアプローチで最も重要なのは、どうやって分布 $p(\boldsymbol{z})$ から独立にサンプリングするのか、という点。

問題設定

今、区間 $(0,1)$ で一様に分布する擬似乱数を発生させるアルゴリズムが得られている。これを元に、 $p(\boldsymbol{z})$ に従う $L$ 個のサンプル $\boldsymbol{z}^{(1)}, \boldsymbol{z}^{(2)}, \cdots, \boldsymbol{z}^{(L)}$ を抽出したい。更に、任意の $\boldsymbol{z}$ に対して $p(\boldsymbol{z})$ は得られないが、ある未知の正規化定数 $Z_p$ が存在して $\tilde{p}(\boldsymbol{z}) = Z_p p(\boldsymbol{z})$ ならば容易に計算できるものとする。

少し恣意的な問題設定に見えるかもしれないが、要するに $\tilde{p}(\boldsymbol{z})$ は分かるけど $Z_p = \int \tilde{p}(\boldsymbol{z}) d\boldsymbol{z}$ は分からないのだと考えれば良い。積分、難しいもんな。

棄却サンプリング

MCMCの前に、棄却サンプリングについて簡単に説明する。

1. 手順

$p(\boldsymbol{z})$ からサンプリングできないので、サンプリングが容易な分布 $q(\boldsymbol{z})$ を用意する(提案分布という)
定数 $k$ を任意の $\boldsymbol{z}$ について $kq(\boldsymbol{z}) \geq \tilde p(\boldsymbol{z})$ となるように定める
$\boldsymbol{z}^*$ を $q(\boldsymbol{z})$ からサンプリングし、確率 $\tilde{p}(\boldsymbol{z}^*) / \lbrace kq(\boldsymbol{z}^*) \rbrace$ で採択する( $k$ の定義から、この値は必ず $1$ 以下になる)

2. なぜ上手くいくのか

所望の分布 $p(\boldsymbol{z})$ と区別するため、確率を示したい場合は $\text{Pr}(\boldsymbol{z})$ と表記する。変数 $a$ を、ステップ3において採択されれば $a=1$ 、棄却されれば $a=0$ となるように定義する。この時、採択された $\boldsymbol{z}$ の分布は $\text{Pr}( \boldsymbol{z} | a=1)$ と書ける。ベイズの定理より、

$\begin{align} \text{Pr}(\boldsymbol{z} | a=1) &= \frac{\text{Pr}(\boldsymbol{z}) \text{Pr}(a=1 | \boldsymbol{z})}{\text{Pr}(a=1)} \end{align}$

と書ける。右辺を見ていこう。

$\begin{align} \text{Pr}(\boldsymbol{z}) &= q(\boldsymbol{z}) \\\ \text{Pr}(a=1 | \boldsymbol{z}) &= \frac{\tilde{p}(\boldsymbol{z})}{kq(\boldsymbol{z})} \\\ &= \frac{Z_p p(\boldsymbol{z})}{kq(\boldsymbol{z})} \\\ \text{Pr}(a=1) &= \int \text{Pr}(a=1, \boldsymbol{z}) d \boldsymbol{z} \\\ &= \int \text{Pr}(\boldsymbol{z}) \text{Pr}(a=1 | \boldsymbol{z}) d \boldsymbol{z} \\\ &= \int q(\boldsymbol{z}) \frac{\tilde{p}(\boldsymbol{z})}{kq(\boldsymbol{z})} d \boldsymbol{z} \\\ &= \frac{Z_p}{k} \int p(\boldsymbol{z}) d \boldsymbol{z} \\\ &= \frac{Z_p}{k} \end{align}$

これらを代入すれば $\text{Pr}(\boldsymbol{z} | a=1) = p(\boldsymbol{z})$ が得られる。

3. 欠点

せっかくサンプリングしたのに高い確率で棄却してしまうのは非効率的なので、 $k$ はできるだけ小さくあって欲しい。そのためには、 $q(\boldsymbol{z})$ は $p(\boldsymbol{z})$ に限りなく近いものであって欲しい。そんなこと最初から出来たら苦労しない。

MCMC

MCMCでは、提案分布 $q(\boldsymbol{z})$ を現在の状態に依存させることでサンプリングをより効率的に行うことを目指す。この時、提案分布は現在の状態 $\boldsymbol{z}^{(\tau)}$ に基づいた分布 $q(\boldsymbol{z} | \boldsymbol{z}^{(\tau)})$ として表現される。

MCMCの例: Metropolis-Hastingsアルゴリズム

1. 手順

また、状態 $\boldsymbol{z}^{(\tau)}$ の更新則は、

$\boldsymbol{z}^{(\tau + 1)} = \left\{ \begin{array}{cc} \boldsymbol{z}^* & \text{ if } \boldsymbol{z}^* \text{ is accepted} \\ \boldsymbol{z}^{(\tau)} & \text{ otherwise} \end{array} \right.$

2. なぜ上手くいくのか

~~あんまりよく分かりません~~ 既約で非周期的なマルコフ連鎖が定常分布を持つなら、必ずその定常分布に収束する(収束定理)。問題はその定常分布が所望の分布 $p(\boldsymbol{z})$ かどうかということ。

分布 $p(\boldsymbol{z})$ が定常分布であることを保証するための十分条件は、以下のような詳細釣り合い条件が挙げられる。

$p(\boldsymbol{z}^{(\tau)}) \text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*) = p(\boldsymbol{z}^*) \text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})$

ここで、 $\text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)$ は状態が $\boldsymbol{z}^{(\tau)}$ から $\boldsymbol{z}^*$ に遷移する確率を指し、今回の設定では

$\begin{align} \text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*) &= q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)}) A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)}) \end{align}$

と表される。 $\text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})$ も同様。この時、詳細釣り合い条件を少し変形させたもの

$\begin{align} \frac{\text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)}{\text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})} &= \frac{p(\boldsymbol{z}^*)}{p(\boldsymbol{z}^{(\tau)})} \end{align}$

が本当に成立するのか確かめるために、左辺について考える。上で定義した式を代入すれば、

$\begin{align} \frac{\text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)}{\text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})} &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)}) A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*) A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} \\\ &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)} \frac{A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} \end{align}$

と書ける。ここで、最右辺の右側の分数について更に考察する。分母分子について、

$\begin{align} A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)}) &= \min \left( 1, \frac{\tilde{p}(\boldsymbol{z}^*) q (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \right) \\\ A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*) &= \min \left( 1, \frac{\tilde{p}(\boldsymbol{z}^{(\tau)}) q (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})}{\tilde{p}(\boldsymbol{z}^*) q (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)} \right) \\\ \end{align}$

なのだが、よく見ると $\min$ の第二項が逆数になっているだけである。つまり、どちらかが $1$ より大きく、どちらかが $1$ より小さくなることが保証される。どちらがどうなるにせよ、先ほどの分数について、

$\begin{align} \frac{A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} = \frac{\tilde{p}(\boldsymbol{z}^*) q (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \end{align}$

となることは簡単な計算で分かるだろう。これを代入すれば、

$\begin{align} \frac{\text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)}{\text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})} &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)} \frac{A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} \\\ &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)} \frac{\tilde{p}(\boldsymbol{z}^*) q (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \\\ &= \frac{\tilde{p}(\boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)})} = \frac{p(\boldsymbol{z}^*)}{p(\boldsymbol{z}^{(\tau)})} \end{align}$

はい示せた。

結局

MCMCではまず提案分布 $q(\boldsymbol{z})$ を置き、何回かのサンプリングでそれを $p(\boldsymbol{z})$ に近似する。十分に近似できたらサンプルとして収集すれば良い。蒸留酒の最初の部分捨てるみたいだね。

まとめ

MCMCのコンセプトを理解するために、まずは基本的なMetropolis-Hastingsアルゴリズムをまとめた。変種としてギブスサンプリングやスライスサンプリングなどがあるが、~~勘弁してください~~今後理解する予定です。

おまけ

えむしーえむしー $\leftarrow$ かわいい

2024-04-04

変分ベイズを今更理解する

趣旨

変分ベイズを理解する。

なぜ変分ベイズを?

この前EMアルゴリズムを解説した。

usapyoi.hatenablog.com

じゃあ次は変分ベイズだね。

変分ベイズはEMより圧倒的に数式を追うのが難しいので頑張りましょう。

表記法

$\boldsymbol{X}$ : 観測変数(即ち入力)
$\boldsymbol{Z}$ : 潜在変数

モデルパラメタ $\boldsymbol{\theta}$ は?

変分ベイズの枠組みでは、モデルパラメタ $\boldsymbol{\theta}$ は"未知の変数"という意味で潜在変数 $\boldsymbol{Z}$ に含まれる。ここがEMアルゴリズムと大きく違う点である。

問題設定

同時分布 $p(\boldsymbol{X},\boldsymbol{Z})$ から、事後分布 $p(\boldsymbol{Z} | \boldsymbol{X})$ や、モデルエビデンス $p(\boldsymbol{X})$ を計算したい。事後分布がわかればそれを最大にする潜在変数 $\boldsymbol{Z}$ も分かる。そして変分ベイズでは $\boldsymbol{Z}$ にモデルパラメタが含まれるので、"最適な"モデルパラメタが $\mathbb{E} \lbrack p(\boldsymbol{Z} | \boldsymbol{X}) \rbrack$ としてわかるのである。

EMアルゴリズムでも事後分布 $p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta}^{(i)})$ をEステップで用いている。ただしEMにおける潜在変数はモデルパラメタを含まないので、 $p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta}^{(i)})$ を知っているからといってモデルパラメタが直ちに求まるわけではない。
モデルエビデンスの解釈は他サイトに譲る。~~知らんし~~

しかし、EMアルゴリズムが想定する状況ほど簡単には得られないものとする。

さて早速解法に $\cdots$ と行きたいところだが、前提と仮定が入り組みすぎてて面倒なので一旦後回しにする。うるせぇ早く見せろって人は第5章へどうそ。

1. $q(\boldsymbol{Z})$ による周辺化

$p(\boldsymbol{X},\boldsymbol{Z})$ 、 $p(\boldsymbol{Z} | \boldsymbol{X})$ 、 $p(\boldsymbol{X})$ の関係性は、EMアルゴリズムの時のように任意の分布 $q(\boldsymbol{Z})$ を導入して、

$\begin{align} \log p(\boldsymbol{X}) &= \log \left\lbrace \int p(\boldsymbol{X}, \boldsymbol{Z}) d \boldsymbol{Z}\right\rbrace \\\ &= \log \left\lbrace \int q(\boldsymbol{Z}) \frac{p(\boldsymbol{X}, \boldsymbol{Z})}{q(\boldsymbol{Z})} d \boldsymbol{Z} \right\rbrace \\\ &= \int q(\boldsymbol{Z}) \log \left\lbrace \frac{p(\boldsymbol{X}, \boldsymbol{Z})}{q(\boldsymbol{Z})} \right\rbrace d \boldsymbol{Z} - \int q(\boldsymbol{Z}) \log \left\lbrace \frac{p(\boldsymbol{Z} | \boldsymbol{X})}{q(\boldsymbol{Z})} \right\rbrace d \boldsymbol{Z} \end{align}$

と表され、さらに最右辺第一項を汎関数 $\mathcal{L}(q)$ 、第二項をKLダイバージェンス $\text{KL}(q || p)$ で表現すると、

$\begin{align} \log p(\boldsymbol{X}) &= \mathcal{L}(q) + \text{KL}(q || p) \end{align}$

のように非常に簡素に書ける。

汎関数とは、ざっくり言えば「関数を入力とする関数」のこと。今回の例では $\mathcal{L}(q)$ は確率分布 $q(\boldsymbol{Z})$ を入力としている
KLダイバージェンスは、確率分布同士の違い(乖離度)を示す指標である。今回だと $q(\boldsymbol{Z})$ と $p(\boldsymbol{Z} | \boldsymbol{X})$ を比較しており、 $q(\boldsymbol{Z}) = p(\boldsymbol{Z} | \boldsymbol{X})$ の時最小値 $0$ を取る。調べましょう。
KLダイバージェンスは $0$ 以上であることが保証されている。よって、上の式から以下のように表現でき、これは以前Jensenの不等式だのなんだので導出した不等式と等価であることがわかる

$\begin{align} \log p(\boldsymbol{X}) &\geq \mathcal{L}(q) \end{align}$

さて、左辺は定数であるため、 $\mathcal{L}(q)$ が最大の時 $\text{KL}(q || p)=0$ 、つまり $q(\boldsymbol{Z})=p(\boldsymbol{Z} | \boldsymbol{X})$ が成立するのだが、EMアルゴリズムの時と違い $p(\boldsymbol{Z} | \boldsymbol{X})$ が計算できない状況である。

2. 計算できないなら近似するしかない-平均場近似

変分ベイズでは、 $q(\boldsymbol{Z})$ を以下のような式で拘束し、拘束条件下で $p(\boldsymbol{Z} | \boldsymbol{X})$ に近づけることを考える。

$\begin{align} q(\boldsymbol{Z}) = \prod_{i=1}^M q_i(\boldsymbol{Z}_i) \end{align}$

つまり、潜在変数 $\boldsymbol{Z}$ を $M$ 個のグループ $\boldsymbol{Z}_1, \boldsymbol{Z}_2, \cdots, \boldsymbol{Z}_M$ に分割して、グループ間は統計的に独立であると仮定している。

例えば、GMMでは混合割合 $\boldsymbol{\pi}$ のグループや、平均 $\boldsymbol{\mu}$ のグループ、また分散共分散行列 $\boldsymbol{\Sigma}$ のグループがあった。グループ内はともかく、グループ間でこれらが依存しあっているとは考えにくい。

この仮定は平均場近似と呼ばれる。

それの何が嬉しいの?

平均場近似によって、 $q(\boldsymbol{Z})$ 上の $\mathcal{L}(q)$ の最大化から各 $i$ についての $q_i(\boldsymbol{Z}_i)$ 上の $\mathcal{L}(q)$ の最大化へと分けることができる。以下で、 $\mathcal{L}(q)$ の中の $q_j(\boldsymbol{Z}_j)$ に依存する項を探してみよう(頑張りましょう)。

3. $\mathcal{L}(q)$ の中の $q_j(\boldsymbol{Z}_j)$ に依存する項を探す

3-1. 重積分への変形

まず、潜在変数 $\boldsymbol{Z}$ 全体の積分から、独立関係にある潜在変数のグループ $\boldsymbol{Z}_m$ それぞれの積分の繰り返し、即ち重積分で表現する(あとついでに $q(\boldsymbol{Z})$ も変換する)。

$\begin{align} \mathcal{L}(q) &= \int q(\boldsymbol{Z}) \log \left\lbrace \frac{p(\boldsymbol{X}, \boldsymbol{Z})}{q(\boldsymbol{Z})} \right\rbrace d \boldsymbol{Z} \\\ &= \int \int \cdots \int \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log \left\lbrace \frac{p(\boldsymbol{X}, \boldsymbol{Z})}{\prod_{i=1}^M q_i (\boldsymbol{Z}_i)} \right\rbrace d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ \end{align}$

3-2. $\log$ に注目して分解する

$\log$ の性質より2つに分解する。見るのも嫌になるがやっていることはまだ単純。

$\begin{align} \mathcal{L}(q) &= \int \int \cdots \int \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log \left\lbrace \frac{p(\boldsymbol{X}, \boldsymbol{Z})}{\prod_{i=1}^M q_i (\boldsymbol{Z}_i)} \right\rbrace d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ &= \int \int \cdots \int \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log p(\boldsymbol{X}, \boldsymbol{Z}) d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ & \quad - \int \int \cdots \int \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log \left\lbrace \prod_{i'=1}^M q_{i'} (\boldsymbol{Z}_{i'}) \right\rbrace d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ \end{align}$

3-3. それぞれ見ていく

ここからは、最右辺の2つの項のそれぞれについて式変形を試みる。

3-3-1. 第一項について

3-3-1-1. $q_j(\boldsymbol{Z}_j)$ の積分を一番外に持ってくる

積分の順序交換(数学的にそんなことして良いのかは触れない)によって、 $q_j(\boldsymbol{Z}_j)$ の積分を一番外側に持ってくる。

$\begin{align} \int \int \cdots \int &\left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log p(\boldsymbol{X}, \boldsymbol{Z}) d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ &= \int q_j (\boldsymbol{Z}_j) \left\lbrack \int \cdots \int \left\lbrace \prod_{i \neq q} q_i (\boldsymbol{Z}_i) \right\rbrace \log p(\boldsymbol{X}, \boldsymbol{Z}) d \boldsymbol{Z}_1 \cdots d \boldsymbol{Z}_M \right\rbrack d\boldsymbol{Z}_j \\\ \end{align}$

3-3-1-2. $\lbrack \cdot \rbrack$ の中身を簡単にする

上の式をよく見ると、 $\lbrack \cdot \rbrack$ の中身は期待値である(ここで、あくまで $\lbrace \boldsymbol{Z}_1, \cdots , \boldsymbol{Z}_{j-1}, \boldsymbol{Z}_{j+1}, \cdots , \boldsymbol{Z}_M \rbrace$ における期待値であって、 $\boldsymbol{Z}_j$ は含まれていないことに注意)。よって、

のように書ける。

3-3-1-3. 確率分布 $\tilde{p} (\boldsymbol{X}, \boldsymbol{Z}_j)$ を導入

さらに状況をわかりやすくするため、以下のような確率分布 $\tilde{p} (\boldsymbol{X}, \boldsymbol{Z}_j)$ を導入する。

$\begin{align} \tilde{p} (\boldsymbol{X}, \boldsymbol{Z}_j) = C \exp \left( \mathbb{E}_{i \neq j: \boldsymbol{Z}_i \sim q_i (\boldsymbol{Z}_i)} \lbrack \log p(\boldsymbol{X}, \boldsymbol{Z}) \rbrack \right) \end{align}$

ここで、 $C$ は規格化条件

$\begin{align} \int \int \tilde{p} (\boldsymbol{X}, \boldsymbol{Z}_j) d\boldsymbol{X} d\boldsymbol{Z}_j = 1 \end{align}$

を守るための定数だと考えれば良い。これの対数は

$\begin{align} \log \tilde{p}(\boldsymbol{X}, \boldsymbol{Z}_j) = \mathbb{E}_{i \neq j: \boldsymbol{Z}_i \sim q_i (\boldsymbol{Z}_i)} \lbrack \log p(\boldsymbol{X}, \boldsymbol{Z}) \rbrack + \text{Const.} \end{align}$

となる。これの導入により、最終的に第一項は

$\begin{align} \int \int \cdots \int &\left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log p(\boldsymbol{X}, \boldsymbol{Z}) d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ &= \int q_j (\boldsymbol{Z}_j) \hspace{1mm} \mathbb{E}_{i \neq j: \boldsymbol{Z}_i \sim q_i (\boldsymbol{Z}_i)} \lbrack \log p(\boldsymbol{X}, \boldsymbol{Z}) \rbrack d\boldsymbol{Z}_j \\\ &= \int q_j (\boldsymbol{Z}_j) \log \tilde{p}(\boldsymbol{X}, \boldsymbol{Z}_j) d\boldsymbol{Z}_j + \text{Const.} \\\ \end{align}$

と書けるのだった。

3-3-2. 第二項について

やることが今までとそこまで変わらないのでどんどん行きましょう。

3-3-2-1. $\log$ に注目して分解する

3-3-2-2. $q_j(\boldsymbol{Z}_j)$ の積分を一番外に持ってくる

$\begin{align} \int \int \cdots \int &\left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ &= \sum_{i'=1}^M \int \int \cdots \int \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log q_{i'} (\boldsymbol{Z}_{i'}) d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ &= \sum_{i'=1}^M \int q_j (\boldsymbol{Z}_j) \left\lbrack \int \cdots \int \left\lbrace \prod_{i \neq j} q_i (\boldsymbol{Z}_i) \right\rbrace \log q_{i'} (\boldsymbol{Z}_{i'}) d \boldsymbol{Z}_1 d \cdots d \boldsymbol{Z}_M \right\rbrack d \boldsymbol{Z}_j \\\ \end{align}$

3-3-2-3. 定数項をまとめる

今欲しいのは $q_j(\boldsymbol{Z}_j)$ に依存する項であってそれ以外は要らない。全部 $\text{Const.}$ に押し込んでしまおう。

$\begin{align} \int \int \cdots \int &\left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ &= \sum_{i'=1}^M \int q_j (\boldsymbol{Z}_j) \left\lbrack \int \cdots \int \left\lbrace \prod_{i \neq j} q_i (\boldsymbol{Z}_i) \right\rbrace \log q_{i'} (\boldsymbol{Z}_{i'}) d \boldsymbol{Z}_1 \cdots d \boldsymbol{Z}_M \right\rbrack d \boldsymbol{Z}_j \\\ &= \int q_j (\boldsymbol{Z}_j) \left\lbrack \int \cdots \int \left\lbrace \prod_{i \neq j} q_i (\boldsymbol{Z}_i) \right\rbrace \log q_j (\boldsymbol{Z}_j) d \boldsymbol{Z}_1 \cdots d \boldsymbol{Z}_M \right\rbrack d \boldsymbol{Z}_j + \text{Const.} \\\ \end{align}$

3-3-2-4. $\lbrack \cdot \rbrack$ の中身を簡単にする

今度の $\lbrack \cdot \rbrack$ の中身は $\log q_j (\boldsymbol{Z}_j)$ である。なぜなら $\log q_j (\boldsymbol{Z}_j)$ は $\lbrace \boldsymbol{Z}_1, \cdots , \boldsymbol{Z}_{j-1}, \boldsymbol{Z}_{j+1}, \cdots , \boldsymbol{Z}_M \rbrace$ に依存しないからである。よって、

$\begin{align} \int \int \cdots \int &\left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace d \boldsymbol{Z}_1 \cdots d \boldsymbol{Z}_M \\\ &= \int q_j (\boldsymbol{Z}_j) \left\lbrack \int \cdots \int \left\lbrace \prod_{i \neq j} q_i (\boldsymbol{Z}_i) \right\rbrace \log q_j (\boldsymbol{Z}_j) d \boldsymbol{Z}_1 \cdots d \boldsymbol{Z}_M \right\rbrack d \boldsymbol{Z}_j + \text{Const.} \\\ &= \int q_j (\boldsymbol{Z}_j) \log q_j (\boldsymbol{Z}_j) d \boldsymbol{Z}_j + \text{Const.} \\\ \end{align}$

お疲れ様でした。

3-4. 結局

上の式を全部代入すれば、

$\begin{align} \therefore \mathcal{L}(q) &= \int \int \cdots \int \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log p(\boldsymbol{X}, \boldsymbol{Z}) d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ & \quad - \int \int \cdots \int \left\lbrace \prod_{i=1}^M q_i (\boldsymbol{Z}_i) \right\rbrace \log \left\lbrace \prod_{i'=1}^M q_{i'} (\boldsymbol{Z}_{i'}) \right\rbrace d \boldsymbol{Z}_1 d \boldsymbol{Z}_2 \cdots d \boldsymbol{Z}_M \\\ &= \int q_j (\boldsymbol{Z}_j) \log \tilde{p}(\boldsymbol{X}, \boldsymbol{Z}_j) d\boldsymbol{Z}_j - \int q_j (\boldsymbol{Z}_j) \log q_j (\boldsymbol{Z}_j) d \boldsymbol{Z}_j + \text{Const.} \\\ &= -\text{KL}(q || \tilde{p})+ \text{Const.} \end{align}$

と書けるのだった。またKLダイバージェンスである。

4. で、どうすんの

第2章で述べたように、変分ベイズにおける最適化は各 $i$ についての $q_i(\boldsymbol{Z}_i)$ 上の $\mathcal{L}(q)$ の最大化と言える。そして、 $q_j(\boldsymbol{Z}_j)$ 上の $\mathcal{L}(q)$ の最大化は、上で得た結論

$\begin{align} \therefore \mathcal{L}(q) &= -\text{KL}(q || \tilde{p})+ \text{Const.} \end{align}$

の最大化、つまり $\text{KL}(q || \tilde{p})$ の最小化と解釈できる。KLダイバージェンスの性質から、最小化は

$\begin{align} \log q_j^{*} (\boldsymbol{Z}_j) &= \log \tilde{p}(\boldsymbol{X}, \boldsymbol{Z}_j) \\\ &= \mathbb{E}_{i \neq j: \boldsymbol{Z}_i \sim q_i (\boldsymbol{Z}_i)} \lbrack \log p(\boldsymbol{X}, \boldsymbol{Z}) \rbrack + \text{Const.} \end{align}$

の時成される。ただし、式を見れば明らかなように、 $\log q_j^{*} (\boldsymbol{Z}_j)$ の計算式に他の潜在変数グループ $\boldsymbol{Z}_i$ ( $i \neq j$ )が入ってしまっており、これは完全な解析解を示せていない。

$\log q_j^{*} (\boldsymbol{Z}_j)$ を計算したい $\rightarrow$ $\log q_i^{*} (\boldsymbol{Z}_i)$ が必要だから計算したい $\rightarrow$ $\log q_j^{*} (\boldsymbol{Z}_j)$ が必要だから計算したい $\rightarrow \cdots$ 以下ループ

これに対処するために、まず全ての $q_i (\boldsymbol{Z}_i)$ に初期値 $q_i^{(0)} (\boldsymbol{Z}_i)$ を与え、反復的に更新することで解を得る。これは収束することが保証されているらしい。

5. 解法

潜在変数 $\boldsymbol{Z}$ を $M$ 個の(統計的に独立と言えそうな)グループ $\boldsymbol{Z}_1, \boldsymbol{Z}_2, \cdots, \boldsymbol{Z}_M$ に分割し、それぞれに確率分布 $q_i (\boldsymbol{Z}_i)$ を定義して初期値 $q_i^{(0)} (\boldsymbol{Z}_i)$ を与える。
各 $j = 1, 2, \cdots, M$ について順番に、以下の式に従って更新する。

$\begin{align} \log q_j^{(t+1)} (\boldsymbol{Z}_j) &= \mathbb{E}_{i \neq j: \boldsymbol{Z}_i \sim q_i^{(t)} (\boldsymbol{Z}_i)} \lbrack \log p(\boldsymbol{X}, \boldsymbol{Z}) \rbrack + \text{Const.} \end{align}$

6. まとめ

変分ベイズをまとめた。変分ベイズとEMアルゴリズムの違いは以下の通り

EMアルゴリズムはモデルパラメタ $\boldsymbol{\theta}$ と潜在変数 $\boldsymbol{Z}$ は別という態度をとり、事前確率 $p(\boldsymbol{X} | \boldsymbol{\theta})$ を最大にするような $\boldsymbol{\theta}$ を求める。
変分ベイズはモデルパラメタ $\boldsymbol{\theta}$ は潜在変数 $\boldsymbol{Z}$ に含まれるという態度をとり、事後確率 $p(\boldsymbol{Z} | \boldsymbol{X})$ を近似して潜在変数 $\boldsymbol{Z}$ を求める。

特に、変分ベイズにおけるパラメタはEMアルゴリズムのように直接求まるわけではなく、あくまで事後分布の期待値、つまり $\mathbb{E} \lbrack p(\boldsymbol{Z} | \boldsymbol{X}) \rbrack$ として求まることに注意。

変分ベイズはGMM(またもや)などに使われる手法である。応用は頑張ってください。

7. おまけ

import numpy as np
from sklearn.mixture import BayesianGaussianMixture
model = BayesianGaussianMixture(n_components=2)
model.fit(X)

sklearn万歳

2024-04-03

EMアルゴリズムで今更GMMを解く

趣旨

前回紹介したEMアルゴリズムを用いて超有名なグラフ混合分布(GMM)を解く。煩雑な計算については極力述べず、統計的な観点から重要な点だけおさえる。

usapyoi.hatenablog.com

なんで今更?

~~うるさいなぁ~~理解を深めるため。

では解説する。

表記法

1. 入力とハイパーパラメタ

$D$ : 入力次元数
$N$ : 入力データの数
$K$ : 重ね合わせるガウス分布の数　未知だが推定不可なので決め打ちする(ハイパーパラメタ)
$\boldsymbol{X} = (\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots, \boldsymbol{x}_N) \in \mathbb{R}^{D \times N}$ : 総入力

2. モデルパラメタ

$\boldsymbol{\pi} = (\pi_1, \pi_2, \cdots, \pi_K) \in \mathbb{R}^{K}$ : モデルパラメタその1　各ガウス分布の混合係数
$\boldsymbol{\mu} = (\boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \cdots, \boldsymbol{\mu}_K) \in \mathbb{R}^{D \times K}$ : モデルパラメタその2　各ガウス分布の平均
$\boldsymbol{\Sigma} = (\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2, \cdots, \boldsymbol{\Sigma}_K) \in \mathbb{R}^{D \times D \times K}$ : モデルパラメタその3　各ガウス分布の分散共分散行列
$\boldsymbol{\theta}$ : 上の3つのモデルパラメタを総称したもの。分けて書く必要性が特にない時はこっちを使う。

3. 潜在変数

$\boldsymbol{Z} = (\boldsymbol{z}_1, \boldsymbol{z}_2, \cdots, \boldsymbol{z}_N) \in {\lbrace 0,1\rbrace}^{K \times N}$ : 潜在変数　対応する $\boldsymbol{x}_n$ がどのガウス分布に所属するかを示す(後述)

問題設定

以下のようなガウス分布の線形重ね合わせを考える。

$p(\boldsymbol{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \sum_{k=1}^K \pi_k \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$

ただし、 $\sum_{k=1}^K \pi_k = 1$ 。これを最大にするようなモデルパラメタ $\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma}$ を求めたい。

解法

Eステップ

負担率 $r_{nk}^{(i)}$ を計算する。

$r_{nk}^{(i)} = \frac{\pi_k^{(1)} \mathcal{N}(\boldsymbol{x}_n | \boldsymbol{\mu}_k^{(i)}, \boldsymbol{\Sigma}_k^{(i)})}{\sum_{\kappa=1}^K \pi_\kappa^{(i)} \mathcal{N}(\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)})}$

期待値の式は省略。

Mステップ

さっき省略した期待値を最大にするようにパラメタを更新する。ここで $N_k$ を以下のように定義する。

$N_k^{(i)} = \sum_{n=1}^{N} r_{nk}^{(i)}$

これを用いて、

$\begin{align} \pi_k^{(i+1)} &= \frac{N_k^{(i)}}{N} \\\ \boldsymbol{\mu}_k^{(i+1)} &= \frac{\sum_{n=1}^N r_{nk}^{(i)} \boldsymbol{x}_n}{N_k^{(i)}} \\\ \boldsymbol{\Sigma}_k^{(i+1)} &= \frac{\sum_{n=1}^N r_{nk}^{(i)} (\boldsymbol{x}_n - \boldsymbol{\mu}_k^{(i)}){(\boldsymbol{x}_n - \boldsymbol{\mu}_k^{(i)})}^T }{N_k^{(i)}} \\\ \end{align}$

なんでそうなるの?

EMアルゴリズムについては前の記事を参照のこと。

usapyoi.hatenablog.com

Eステップ

EMアルゴリズムによれば、Eステップでは期待値 $Q^{(i)}(\boldsymbol{\theta})$ の計算を行う。

さて、潜在変数 $\boldsymbol{Z}$ を導入し、まず $p(\boldsymbol{X}, \boldsymbol{Z} | \boldsymbol{\theta})$ 、次に $p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})$ を求めよう。

1. 潜在変数の導入

潜在変数 $\boldsymbol{z}_n$ を $\boldsymbol{x}_n$ がどのガウス分布に所属するかを示すOne-hotベクトルとして定義する。つまり、

$z_{nk} = \left\lbrace \begin{array}{cc} 1 & \text{if } \boldsymbol{x}_n \sim \mathcal{N} (\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \\ 0 & \text{otherwise} \end{array} \right.$

そして、 $\boldsymbol{z}_n$ を各 $n$ について並べたものを $\boldsymbol{Z}$ とする。

2. $p(\boldsymbol{X}, \boldsymbol{Z} | \boldsymbol{\theta})$ の計算

$\sum_{k=1}^K \pi_k = 1$ より、 $\boldsymbol{x}_n$ が $k$ 番目のガウス分布に所属する確率は $\pi_k$ と書ける。よって、

$\begin{align} p(\boldsymbol{z}_n | \boldsymbol{\pi}) &= \prod_{k=1}^K \pi_k^{z_{nk}} \\\ p(\boldsymbol{Z} | \boldsymbol{\pi}) &= \prod_{n=1}^N \prod_{k=1}^K \pi_k^{z_{nk}} \\\ \end{align}$

一方、 $\boldsymbol{z}_n$ を用いることで、以下のような関係性も成り立つ。

$\begin{align} p(\boldsymbol{x}_n | \boldsymbol{z}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= \prod_{k=1}^K {\mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}^{z_{nk}} \\\ p(\boldsymbol{X} | \boldsymbol{Z}, \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= \prod_{n=1}^N \prod_{k=1}^K {\mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}^{z_{nk}} \\\ \end{align}$

以上より、 $p(\boldsymbol{X}, \boldsymbol{Z} | \boldsymbol{\theta})$ について

$\begin{align} p(\boldsymbol{X}, \boldsymbol{Z} | \boldsymbol{\theta}) &= p(\boldsymbol{Z} | \boldsymbol{\pi}) p(\boldsymbol{X} | \boldsymbol{Z}, \boldsymbol{\mu}, \boldsymbol{\Sigma}) \\\ &= \prod_{n=1}^N \prod_{k=1}^K {\left\lbrace \pi_k \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right\rbrace}^{z_{nk}} \end{align}$

と書ける。

3. $p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})$ の計算

ベイズの定理より、

$\begin{align} p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)}) &= \frac{p(\boldsymbol{X}, \boldsymbol{Z} | \boldsymbol{\theta}^{(i)})}{p(\boldsymbol{X} | \boldsymbol{\theta}^{(i)})} \\\ &= \frac{\prod_{n=1}^N \prod_{k=1}^K {\left\lbrace \pi_k^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k^{(i)}, \boldsymbol{\Sigma}_k^{(i)}) \right\rbrace}^{z_{nk}}}{\prod_{n=1}^N \sum_{k=1}^K \pi_k^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k^{(i)}, \boldsymbol{\Sigma}_k^{(i)})} \\\ \end{align}$

と書ける。ややこしくなってきたが大丈夫。

4. 期待値 $Q^{(i)}(\boldsymbol{\theta})$ の計算

以上の計算を踏まえ、いよいよ $Q^{(i)}(\boldsymbol{\theta})$ 、すなわち期待値の計算に移ろう。

$\begin{align} Q^{(i)}(\boldsymbol{\theta}) & = \mathbb{E}_{\boldsymbol{Z} \sim p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})} \left\lbrack \log p(\boldsymbol{X}, \boldsymbol{Z} | \boldsymbol{\theta}) \right\rbrack \\\ & = \mathbb{E}_{\boldsymbol{Z} \sim p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})} \left\lbrack \log \prod_{n=1}^N \prod_{k=1}^K {\left\lbrace \pi_k \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right\rbrace}^{z_{nk}} \right\rbrack \\\ &= \mathbb{E}_{\boldsymbol{Z} \sim p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})} \left\lbrack \sum_{n=1}^N \sum_{k=1}^K z_{nk} {\left\lbrace \log \pi_k + \log \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right\rbrace} \right\rbrack \\\ &= \sum_{n=1}^N \sum_{k=1}^K \mathbb{E}_{\boldsymbol{Z} \sim p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})} \left\lbrack z_{nk} \right\rbrack \left\lbrace \log \pi_k + \log \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right\rbrace \\\ \end{align}$

ここで負担率 $r_{nk}^{(i)}$ について

$\begin{align} r_{nk}^{(i)} &= \mathbb{E}_{\boldsymbol{Z} \sim p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})} \left\lbrack z_{nk} \right\rbrack \end{align}$

と定義すれば、

$\begin{align} Q^{(i)}(\boldsymbol{\theta}) &= \sum_{n=1}^N \sum_{k=1}^K \mathbb{E}_{\boldsymbol{Z} \sim p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})} \left\lbrack z_{nk} \right\rbrack \left\lbrace \log \pi_k + \log \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right\rbrace \\\ &= \sum_{n=1}^N \sum_{k=1}^K r_{nk}^{(i)} \left\lbrace \log \pi_k + \log \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right\rbrace \\\ \end{align}$

と多少は簡素に書ける。

4. $r_{nk}^{(i)}$ って具体的にはいくらなの?

$\begin{align} r_{nk}^{(i)} &= \mathbb{E}_{\boldsymbol{Z} \sim p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})} \left\lbrack z_{nk} \right\rbrack \\\ &= \sum_{\boldsymbol{Z}} z_{nk} \hspace{1mm} p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)}) \\\ \end{align}$

さて、前で示したような悍ましい式を $p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})$ に入れても良いのだが、その前にある重要な性質を確認しておく: $\boldsymbol{z}_n$ は統計的に独立している。つまり今回期待値を取りたいのは $z_{nk}$ なので、 $\nu \neq n$ における $\boldsymbol{z}_\nu$ などは考慮しなくても良いのである。よって先ほどの期待値は

$\begin{align} r_{nk}^{(i)} &= \sum_{\boldsymbol{z}_n} z_{nk} \hspace{1mm} p(\boldsymbol{z}_n | \boldsymbol{X}, \boldsymbol{\theta}^{(i)}) \\\ &= \sum_{\boldsymbol{z}_n} z_{nk} \hspace{1mm} \frac{\prod_{\kappa=1}^K {\left\lbrace \pi_\kappa^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)}) \right\rbrace}^{z_{n\kappa}}}{\sum_{\kappa=1}^K \pi_k^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)})} \\\ \end{align}$

と簡単化される(ここで、分母分子のインデックスを $k$ と区別するため $\kappa$ とした)。いやまだ簡単ではないだろと思うかもしれない。ここでさらにもう一つトリックを入れる: $\boldsymbol{z}_n$ はOne-hotベクトルである。つまり $z_{nk}$ の取りうる値は $\lbrace 0, 1\rbrace$ しかない。 $z_{nk} = 0$ の時は $\sum$ の中身も $0$ になるのだから、 $z_{nk}=1$ のときだけ考えれば良い。この時 $k' \neq k$ を満たす $k'$ について $z_{nk'}=0$ となることにも注意しよう。

$\begin{align} r_{nk}^{(i)} &= \sum_{\boldsymbol{z}_n} z_{nk} \hspace{1mm} \frac{\prod_{\kappa=1}^K {\left\lbrace \pi_\kappa^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)}) \right\rbrace}^{z_{n\kappa}}}{\sum_{\kappa=1}^K \pi_\kappa^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)})} \\\ &= {\left. \frac{\prod_{\kappa=1}^K {\left\lbrace \pi_\kappa^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)}) \right\rbrace}^{z_{n\kappa}}}{\sum_{\kappa=1}^K \pi_\kappa^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)})} \right|}_{z_{nk} = 1, z_{nk'}=0} \\\ &= \frac{ \pi_k^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k^{(i)}, \boldsymbol{\Sigma}_k^{(i)})}{\sum_{\kappa=1}^K \pi_\kappa^{(i)} \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_\kappa^{(i)}, \boldsymbol{\Sigma}_\kappa^{(i)})} \\\ \end{align}$

が成立する。なおこの負担率が何を意味するのかは触れない。ヨソ行け。

Mステップ

Mステップでは、Eステップで定義した期待値 $Q^{(i)}(\boldsymbol{\theta})$ :

$\begin{align} Q^{(i)}(\boldsymbol{\theta}) &= \sum_{n=1}^N \sum_{k=1}^K r_{nk}^{(i)} \left\lbrace \log \pi_k + \log \mathcal{N} (\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right\rbrace \\\ \end{align}$

を、 $\boldsymbol{\theta}$ (つまり、 $\boldsymbol{\pi}$ 、 $\boldsymbol{\mu}$ 、 $\boldsymbol{\Sigma}$ )について最大化する。偏微分の計算は~~面白くないので~~ ~~めんどくさいので~~ ~~もう疲れたので~~今ここで述べたい要旨からずれるので省略。計算的にはしんどいけど発想的にはそこまでしんどくないので多分大丈夫。なんなら他サイトあるし。しかし1点、拘束条件: $\sum_k \pi_k = 1$ にだけ注意。これの最適化はLagrangeの未定乗数法を用いる。

まとめ

EMアルゴリズムをGMMに適用するという情報学徒なら1度はやることを改めて勉強しなおした。Pythonの実行コードもあるので(そんなものいくらでもありそうだが)参考にどうぞ。

github.com

おまけ

K*gglerの皆さんは理論なんてどうでもいいだろうね、だってこれでいいし。

import numpy as np
from sklearn.mixture import GaussianMixture

model = GaussianMixture(n_components=2)
model.fit(X)

2024-04-02

EMアルゴリズムを今更理解する

趣旨

EMアルゴリズムを理解する。

なぜEMアルゴリズム?

~~なんか名前がかっこいいから~~機械学習を学ぶ人間として当然だから。

他にも解説記事いっぱいあるけど?

表記揺れ、定義の揺れが酷すぎて比較しづらい。例えば後述のEステップについて揺れがある。

では解説する。

表記法

$\boldsymbol{x}$ : 観測変数(即ち入力)
$\boldsymbol{z}$ : 潜在変数(観測されない変数　欠損データと呼ばれたりするが、実在性は特に重要ではない)
$\boldsymbol{\theta}$ : モデルパラメタ

問題設定

最尤推定がしたい。つまり、事前確率 $p(\boldsymbol{x} | \boldsymbol{\theta})$ を最大にするような $\boldsymbol{\theta}$ を求めたい。

$\boldsymbol{\theta}_{\text{ML}} \leftarrow \max_{\boldsymbol{\theta}} \log p(\boldsymbol{x} | \boldsymbol{\theta})$

しかし、直接偏微分して解くことが難しく、代わりに $p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})$ は陽に得られている状況を考える。

解法

$\boldsymbol{\theta}$ に適当な初期値 $\boldsymbol{\theta}^{(0)}$ を与えた後、以下のEステップ、Mステップを収束するまで繰り返す。

Eステップ

まず以下のような事後確率 $q(\boldsymbol{z})$ を計算する。

$q(\boldsymbol{z}) = p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta}^{(i)})$

次に先ほど計算した $q(\boldsymbol{z})$ を用いて、以下のような期待値を $Q^{(i)}(\boldsymbol{\theta})$ とする。

$\begin{align} Q^{(i)}(\boldsymbol{\theta}) & = \mathbb{E}_{\boldsymbol{z} \sim q(\boldsymbol{z})} \left\lbrack \log p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) \right\rbrack \\\ & = \int q(\boldsymbol{z}) \log p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) d \boldsymbol{z} \\\ \end{align}$

期待値をとってるからE(xpectation)ステップ。

Mステップ

先ほど計算した $Q^{(i)}(\boldsymbol{\theta})$ を $\boldsymbol{\theta}$ について最大化し解を $\boldsymbol{\theta}^{(i+1)}$ とする。

$\boldsymbol{\theta}^{(i+1)} \leftarrow \text{argmax}_{\boldsymbol{\theta}} Q^{(i)}(\boldsymbol{\theta})$

最大化しているからM(aximization)ステップ。

なんでうまくいくの?

細かい計算は飛ばして、ざっくり概要を述べる。~~めんどくさいし。~~

Eステップ

1. 潜在変数による周辺化

目的関数である $\log p(\boldsymbol{x} | \boldsymbol{\theta})$ を潜在変数 $\boldsymbol{z}$ についての周辺化によって表現する。

$\log p(\boldsymbol{x} | \boldsymbol{\theta}) = \log \int p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) d \boldsymbol{z}$

2. $q(\boldsymbol{z})$ を導入

任意の確率分布 $q(\boldsymbol{z})$ を導入し、以下のように式変形する。

$\log \int p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) d \boldsymbol{z} = \log \int q(\boldsymbol{z}) \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} d \boldsymbol{z}$

3. Jensenの不等式を導入

3-1. Jensenの不等式って?

関数 $f$ が(下に)凸であるとき、

$\begin{align} u(x) \geq 0, &\quad \int u(x) dx = 1 \\\ &\Rightarrow \int u(x) f(v(x)) dx \geq f(\int u(x) v(x) dx) \end{align}$

が成立し、等号成立は以下の2種類のいずれかを満足すると成立する。

$v(x) = m$ ( $x$ の値に関わらず一定)
$f(x)$ が原点を通る直線である(まず満足しない)

これをJensen(イェンゼン、イェンセン)の不等式と言う。なんでこれが成り立つのかは他のサイトに譲る。ともかくこれを適用することを試みる。

3-2. さっきの式に適用する

具体的には以下の右辺。

関数 $f$ が $\log$ であって欲しい。残念ながら $\log$ は凸は凸でも上に凸だ。しかし逆に考えれば $-\log$ とすれば下に凸になってくれる。よって、

$\begin{align} f(\cdot) &\rightarrow -\log (\cdot) \\\ u(x) &\rightarrow q(\boldsymbol{\theta}) \\\ v(x) &\rightarrow \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} \\\ \end{align}$

のような対応関係を結べば、

$f(\int u(x) v(x) dx) \rightarrow -\log \int q(\boldsymbol{z}) \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} d \boldsymbol{z}$

と分かり、不等式は

$\log \int q(\boldsymbol{z}) \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} d \boldsymbol{z} \geq \int q(\boldsymbol{z}) \log \left( \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} \right) d \boldsymbol{z}$

と書ける(不等号の向きに注意せよ)。

3-3. 等号成立条件

では等号成立はどの時だろうか。Jensenによれば、それは $v(x) = m$ の時である( $\log$ が直線な訳ないので)。即ち、

$\frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} = m$

$\therefore p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) = m q(\boldsymbol{z})$

で、 $m$ っていくらなの? それは $\boldsymbol{z}$ について周辺化してやれば分かる。

$\begin{align} \int p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) d\boldsymbol{z} &= \int m q(\boldsymbol{z}) d\boldsymbol{z} \\\ p(\boldsymbol{x} | \boldsymbol{\theta}) &= m \\\ \end{align}$

と言うわけで、 $m=p(\boldsymbol{x} | \boldsymbol{\theta})$ であり、等号成立条件は

$q(\boldsymbol{z}) = \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{m} = \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{p(\boldsymbol{x} | \boldsymbol{\theta})} = p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta})$

と書き表せるのだった。

3-4. 等号が成立したら何が嬉しいの?

先ほどの不等式

において左辺は本来は $q(\boldsymbol{z})$ に依存しない値であるから、 $q(\boldsymbol{z})$ をいじって等号を成立させることは右辺(下限)の最大化を意味する。よって、目的関数の値は変わってないが、それの下限の値が最大化され、続くMステップにおいて目的関数の値を増加させるのに貢献するのである。

3-5. Eステップとやってること違わない?

右辺について $\log$ に着目し式展開を行うと、

$\begin{align} \int q(\boldsymbol{z}) \log &\left( \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} \right) d \boldsymbol{z} \\\ &= \int q(\boldsymbol{z}) \log p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) d \boldsymbol{z} + \text{Const.} \\\ &= \mathbb{E}_{\boldsymbol{z} \sim q(\boldsymbol{z})} \left\lbrack \log p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) \right\rbrack + \text{Const.} \\\ \end{align}$

と書ける。よって、

$q(\boldsymbol{z}) = p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta})$ とした上で $\mathbb{E}_{\boldsymbol{z} \sim q(\boldsymbol{z})} \left\lbrack \log p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) \right\rbrack$ を求めることは
$\int q(\boldsymbol{z}) \log \left( \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} \right) d \boldsymbol{z}$ の最大値を求めること、ひいては
$\log p(\boldsymbol{x} | \boldsymbol{\theta})$ の(あくまで)下限の $q(\boldsymbol{z})$ における最大値を求めることに等しい。

3-6. 厳密には

実際のEステップでは、 $\boldsymbol{\theta}$ の値として $\boldsymbol{\theta}^{(i)}$ が得られている。よって、先ほどまでの議論は実は $\log p(\boldsymbol{x} | \boldsymbol{\theta}^{(i)})$ の下限を最大化していたのだった。ではなぜそう書かないのか。おそらくだが、続くMステップで $\boldsymbol{\theta}$ の値をいじるため、定数 $\boldsymbol{\theta}^{(i)}$ ではなく変数 $\boldsymbol{\theta}$ として扱うことを強調したかったのだろう。

Mステップ

先のEステップにおいて、 $\log p(\boldsymbol{x} | \boldsymbol{\theta})$ の下限を $q(\boldsymbol{z})$ について最大化した。Mステップでは、その下限を今度は $\boldsymbol{\theta}$ について最大化している。

Jensenの不等式における右辺はEステップによって

$\begin{align} \int q(\boldsymbol{z}) \log &\left( \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta})}{q(\boldsymbol{z})} \right) d \boldsymbol{z} \\\ &= \mathbb{E}_{\boldsymbol{z} \sim q(\boldsymbol{z})} \left\lbrack \log p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}) \right\rbrack + \text{Const.} \\\ &= Q^{(i)}(\boldsymbol{\theta}) + \text{Const.} \\\ \end{align}$

のように表記できる。この式の最右辺を $\boldsymbol{\theta}$ について最大化することは、 $\log p(\boldsymbol{x} | \boldsymbol{\theta})$ の(あくまで)下限を最大化することと等価である。しかしここで問題が発生する。

等号成立の崩壊

Eステップで導入したJensenの不等式において、等号成立条件を

$q(\boldsymbol{z}) = p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta})$

と書いた。そして、実際は定数 $\boldsymbol{\theta}^{(i)}$ を用いて、

$q(\boldsymbol{z}) = p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta}^{(i)})$

のように計算した。しかし、 $\boldsymbol{\theta}$ について最大化、つまり $\boldsymbol{\theta}$ の値をいじって $\boldsymbol{\theta}^{(i+1)}$ としたとき、

$q(\boldsymbol{z}) \neq p(\boldsymbol{z} | \boldsymbol{x}, \boldsymbol{\theta}^{(i+1)})$

となってしまう。要するに、旧来の( $\boldsymbol{\theta}^{(i)}$ を用いた) $q(\boldsymbol{z})$ のままでは、不等式

$\log \int q(\boldsymbol{z}) \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}^{(i+1)})}{q(\boldsymbol{z})} d \boldsymbol{z} \geq \int q(\boldsymbol{z}) \log \left( \frac{p(\boldsymbol{x}, \boldsymbol{z} | \boldsymbol{\theta}^{(i+1)})}{q(\boldsymbol{z})} \right) d \boldsymbol{z}$

の等号成立条件を満たさなくなってしまう。これは何を意味するのか?

等号成立の崩壊は進歩を意味する

Mステップでやっていることは、等号が成立していた不等式

の右辺を最大化することで等号を不成立にすることである。最大化しているのだから、右辺はもちろん増加する。だが、等号が不成立なのだから、左辺はそれ以上に増加するのである。左辺というのはすなわち目的関数なので、目的関数の値が増加する結果となるのである。

Eステップに戻る

等号が不成立になっちゃったのなら、また $q(\boldsymbol{z})$ をいじって成立させればいいじゃない、そしてまた最大化して不成立にさせればいいじゃない。それを繰り返せばいつかは目的関数の値は上限に到達し、 $\boldsymbol{\theta}^{(i)}$ は最尤解 $\boldsymbol{\theta}_{\text{ML}}$ に限りなく近づくだろう。

結局

EMアルゴリズムは、 $\log p(\boldsymbol{x} | \boldsymbol{\theta})$ の下限を

$q(\boldsymbol{z})$ について最大化(目的関数 $=$ 下限)
$\boldsymbol{\theta}$ について最大化(目的関数 $>$ 下限になり結果的に目的関数の値が増加)

し続けるアルゴリズムなのである。終わり。

具体的にはどこで使われてるの?

GMM(混合ガウス分布)があまりにも有名だが、IRT(項目反応理論)でも用いられていたりする。予想しないタイミングでひょっこり顔を出すので注意。

GMMだとどういう更新則になるの?

~~知らん。その辺のサイトに転がっているので確認せよ。~~書いたよ。

usapyoi.hatenablog.com

まとめ

EMアルゴリズムを私なりに説明した。難解なことをやっているようで(実際しっかり理解しようとすると相当厄介なのだが)実際はそこまででもない。KLダイバージェンスの説明や行間の式変形、また理論的に解が得られる保証を省略したので、概要を掴んだ後はより詳しいサイトに行ってみよう。

おまけ

変分ベイズとかもかっこいいよね。響きが。

2024-02-16

Pythonライブラリ集

備忘録

地味に役立つPythonライブラリを紹介(随時追加予定)。派手に役立つもの(numpyとか)はいちいち紹介しない。

ydata_profiling

docs.profiling.ydata.ai

DataFrameの特徴をサクッと可視化するライブラリ。欠損値とか統計量とか散布図とかはこれで表現できるので、EDAに便利。

from ydata_profiling import ProfileReport
profile = ProfileReport(df, title="Pandas Profiling Report")

networkx + pyvis

networkx.org

pyvis.readthedocs.io

グラフネットワークの処理をするnetworkxと、それをインタラクティブに可視化するpyvisの組み合わせ。colab環境の場合コンソール上で可視化することはできず、一旦htmlに出力してダウンロードしなければならないのがちょっと不便か。

import networkx as nx
from pyvis.network import Network
G = nx.Graph()
G.add_nodes_from([1, 2, 3, 4, 5, 6])
G.add_edges_from([(1, 4), (1, 5), (1, 6), (3, 5), (3, 6)])
net = Network(notebook=True, cdn_resources='remote')
net.from_nx(G)
net.show("example.html")

polars

pola.rs

Pandas互換。Pandasと比べて謎の仕様が少なく、またRustベースで書かれていて動作も高速らしい。qiitaの記事が有用そう。

qiita.com

import polars as pl
df = pl.read_csv("train.csv")

趣旨

注意

理解してないなら「理解する」とか書くなよ

HMCとは

1. 初等物理量

1-1. 位置、速度、加速度

1-2. 質量、運動量、力

1-3. 仕事

2. 力学的エネルギー

2-1. ポテンシャルエネルギー

2-2. 運動エネルギー

2-3 エネルギー保存則

2-4. ハミルトニアン、運動方程式

2-5. オイラー法とリープフロッグ法

2-6. 位相空間

3. MCMCに求められる条件

3-1. エルゴード性

3-2. 詳細釣り合い条件

4. HMC

4-1. ハミルトニアンの導入

4-2. 運動の解釈

4-3. 詳細釣り合い条件

4-4. 数値誤差と棄却

4-4. 結局

まとめ

おまけ

趣旨

MCMCって？

ギブスサンプリングとは

問題設定

ギブスサンプリング

1. 手順

2. なぜ上手くいくのか

2-1. 復習 Metropolis-Hastingsアルゴリズムの手順

2-2. ギブスサンプリングとの対応関係

まとめ

おまけ

趣旨

なぜMCMC?

MCMCとは

背景

サンプリング法

問題設定

棄却サンプリング

1. 手順

2. なぜ上手くいくのか

3. 欠点

MCMCの例: Metropolis-Hastingsアルゴリズム

1. 手順

2. なぜ上手くいくのか

結局

まとめ

おまけ

趣旨

なぜ変分ベイズを?

表記法

モデルパラメタは?

問題設定

1. による周辺化

2. 計算できないなら近似するしかない-平均場近似

それの何が嬉しいの?

3. の中のに依存する項を探す

3-1. 重積分への変形

3-2. に注目して分解する

3-3. それぞれ見ていく

3-3-1. 第一項について

3-3-1-1. の積分を一番外に持ってくる

3-3-1-2. の中身を簡単にする

3-3-1-3. 確率分布を導入

3-3-2. 第二項について

3-3-2-1. に注目して分解する

3-3-2-2. の積分を一番外に持ってくる

3-3-2-3. 定数項をまとめる

3-3-2-4. の中身を簡単にする

3-4. 結局

4. で、どうすんの

5. 解法

6. まとめ

7. おまけ

趣旨

モデルパラメタ $\boldsymbol{\theta}$ は?

1. $q(\boldsymbol{Z})$ による周辺化

3. $\mathcal{L}(q)$ の中の $q_j(\boldsymbol{Z}_j)$ に依存する項を探す

3-2. $\log$ に注目して分解する

3-3-1-1. $q_j(\boldsymbol{Z}_j)$ の積分を一番外に持ってくる

3-3-1-2. $\lbrack \cdot \rbrack$ の中身を簡単にする

3-3-1-3. 確率分布 $\tilde{p} (\boldsymbol{X}, \boldsymbol{Z}_j)$ を導入

3-3-2-1. $\log$ に注目して分解する

3-3-2-2. $q_j(\boldsymbol{Z}_j)$ の積分を一番外に持ってくる

3-3-2-4. $\lbrack \cdot \rbrack$ の中身を簡単にする

2. $p(\boldsymbol{X}, \boldsymbol{Z} | \boldsymbol{\theta})$ の計算

3. $p(\boldsymbol{Z} | \boldsymbol{X}, \boldsymbol{\theta}^{(i)})$ の計算

4. 期待値 $Q^{(i)}(\boldsymbol{\theta})$ の計算

4. $r_{nk}^{(i)}$ って具体的にはいくらなの?

2. $q(\boldsymbol{z})$ を導入