EMアルゴリズムの手短な導出
観測されるデータをまとめて y と書くことにし, 次の対数尤度関数を最大にする \theta を求めたいとする.
\ell(\theta) = \log p(y|\theta)
ここで p(\cdot) は確率変数ごとに定義された密度を表す記号とした（つまり p(y) と p(z) が同じ関数とは限らない）.
適当な確率変数 z を用いて,
p(y|\theta) = \int p(y, z|\theta) \, dz
と書けるとしよう. ここでの積分記号は定義域全区間に渡る積分を意味する.
この記事ではそれ以外の範囲での積分は出てこないので混乱することは多分ないだろうと思う.
条件付き確率より,
p(z|y,\theta) = p(y, z|\theta) / p(y|\theta)
であるから, 両辺の対数の期待値を取り,
\log p(y|\theta) = E_z[\log p(y, z|\theta)] - E_z[\log p(z|y,\theta)] \tag{1}
も成り立つ. ここでの期待値は「yを所与としたときのzの条件付き期待値」,
E_z[f(z)] = \int f(z) p(z|y, \theta) \, dz
である.
さて, EM（Expectation-Maximization）アルゴリズムでは, (1)式の第1項を求め（期待値を取る操作なのでこれをEステップと呼ぶ）, 第1項が大きくなるようにパラメータを更新する（最大化を目指す操作なのでこれをMステップと呼ぶ）.
EMアルゴリズムでは(1)式の右辺の第2項は明示的に計算する必要がない.
なぜならば, ある\theta_0を固定して任意の\thetaを考えると,
\begin{aligned}
& E_z[\log p(z|y,\theta)] - E_z[\log p(z|y,\theta_0)]\\
&= \int \log \left(\frac{p(z|y,\theta)}{p(z|y,\theta_0)}\right) p(z|y,\theta_0) \, dz \\
& \le \log\left( \int \frac{p(z|y,\theta)}{ p(z|y,\theta_0)} p(z|y,\theta_0) \, dz \right) = \log 1 =0\\
\end{aligned}
が成り立つ. 3行目はイェンセンの不等式による.
よって, 「第1項が大きくなるようにパラメータを更新する」すなわち
E_z[\log p(y, z|\theta)] - E_z[\log p(y,z|\theta_0)] \ge 0
が満たされていれば,
\log p(y|\theta) -\log p(y|\theta_0) \ge 0
も成り立つことがわかる.
これがEMアルゴリズムの基礎である.
尤度と事前分布の積を最大化するMAP（Maximum A Posteriori）推定の場合も, 尤度に事前分布の密度をかけ合わせるだけで全く同様の議論ができる.
変分ベイズ法の手短な導出
上と同様, 観測されるデータをまとめてyと書くことにし, 潜在変数zと\thetaを用いて,
p(y,z,\theta) = p(y|z,\theta) p(z) p(\theta)
の形で書けるモデルを考える.
またここでの積分記号は定義域全区間に渡る積分を意味する. この記事内ではそれ以外の範囲での積分は出てこない.
いま, 事後分布として p(\theta, z|y) を求めたい. 一般には事後分布は解析的には求まらない.
そこで近似事後分布のクラスを
p(\theta, z|y) \approx q(\theta) q(z)
の形に制限し, q(\theta) と q(z) を求めることにする.
ここでq(\cdot)はやはり確率変数ごとに定義された密度を表す記号とした. p(\cdot)はモデルで測った確率を表していて, q(\cdot) は近似事後分布で測った確率を表している.
さて, 上では近似（\approx）の記号を曖昧に使ったが, これから次のKL（カルバック・ライブラ）情報量の意味での近似を考える.
D(q(z)q(\theta)\|p(\theta, z|y)) = \int\int \log \frac{q(z)q(\theta)}{p(\theta, z|y)}\, dz \, d\theta
\theta を固定したとき, q(\theta) による期待値を
E_{\theta} [f(\theta)] = \int f(\theta)q(\theta) \, d\theta
と表すことにすると
D(q(z)q(\theta)\|p(\theta, z|y)) = \int -E_{\theta}[\log p(\theta, z|y)] + E_{\theta}[\log q(z)] d\theta
と書け, これは
q(z) \propto \exp(E_{\theta}[\log p(\theta, z|y)]) \tag{2}
のとき最小になる.
z を固定したときも, q(z) による期待値を考えると,
q(\theta) \propto \exp(E_{z}[\log p(\theta, z|y)]) \tag{3}
のときKL情報量が最小になる.
変分ベイズ法では(2)式による q(z) の更新と, (3)式による q(\theta) の更新を交互に繰り返すことで事後分布を近似する q(\cdot) を求める.
特に指数型分布の混合分布モデルを考えるときなど, 変分ベイズ法はEMアルゴリズムとよく似たアルゴリズムが導かれるので, (2)式による q(z) の更新を「変分Eステップ」, (3)式による q(\theta) の更新を「変分Mステップ」と呼ぶこともある.
変分下限
さて, 変分ベイズ法についてほんの少しだけ違う視点から議論するために, 次の対数周辺尤度と呼ばれる量を考える:
\log p(y) = \log \int p(y, z,\theta)\, d z d\theta.
対数周辺尤度を次のように変形する:
\log p(Y) = \log \int q(z)q(\theta) \frac{p(Y, z, \theta)}{q(z)q(\theta)}\,dz d\theta.
イエンセンの不等式により、対数周辺尤度に対して以下が成り立つ。
\log p(Y) \ge \int \int q(z)q(\theta) \log \frac{p(y, z,\theta)}{q(z,\theta)}\,dz d\theta,
右辺を変分下限と呼ぶ. evidence lower bound (ELBO) と呼ばれることもある. 分野によって色々な呼ばれ方をしている.
変分下限は
\begin{aligned}
& \int \int q(z)q(\theta) \log \frac{p(y)p(y, z,\theta)}{p(y)q(z)q(\theta)}\,dz d\theta \\
& = \mathrm{const.} - \int \int q(z)q(\theta) \log \frac{q(z)q(\theta)}{p(z,\theta|y)}\,dz d\theta \tag{4}
\end{aligned}
と変形できる.
\mathrm{const.} は q(z)q(\theta) に依存しない項（定数）, 第2項は事後分布と q(z)q(\theta) のKL情報量である. すなわち, KL情報量の小さい事後分布の近似を求めることは, 変分下限を最大化することと等価であることがわかった.
仮にこの q(\theta) を1点まで退化した分布とすると（あるいは \theta を「定数」とすると）,
なので(4)式の第2項は
- \int q(z) \log \frac{q(z)}{p(z,\theta|y)}\,dz
である. これを(1)式のEMアルゴリズムの目的関数と比べると,
H_z = -\int q(z) \log q(z)\,dz
で定義される H_z の有無のみが異なっている. この H_z はエントロピーと呼ばれる.
VAE（Variational Auto-Encoder）などでは, 全体に共通のニューラルネットワークのウェイト（重み, 係数）パラメータの部分をMAP法で推定して, サンプルごとに異なる潜在変数を変分ベイズ法（ELBO最大化）で推定することが広く行われているようである. そのような場合には実装するコード上に明示的に現れるロス関数においても, エントロピー H_z があるかないかがEMアルゴリズムと変分ベイズ法の違いになるだろう.
（→実装編に続く…）
Discussion