EMアルゴリズム(Expectation-Maximization Algorithm)は,潜在変数(隠れ変数)を含む確率モデルにおいて,観測データの対数尤度を最大化するパラメータを反復的に推定するための汎用的な手法である.直接最適化が困難な問題に対し,E ステップと M ステップを交互に繰り返すことで単調増加を保証しながら局所最適解へ収束する.
観測データを $\mathbf{X} = \{x_1, x_2, \ldots, x_N\}$,潜在変数を $\mathbf{Z} = \{z_1, z_2, \ldots, z_N\}$,モデルパラメータを $\boldsymbol{\theta}$ とする.$(\mathbf{X}, \mathbf{Z})$ の組を完全データ,$\mathbf{X}$ のみを不完全データ(観測データ)と呼ぶ.
目標は,観測データの対数尤度\[\ell(\boldsymbol{\theta}) = \log p(\mathbf{X} \mid \boldsymbol{\theta}) = \log \int p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})\, d\mathbf{Z}\]を最大化するパラメータ $\boldsymbol{\theta}^* = \arg\max_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta})$ を求めることである.$\mathbf{Z}$ の積分(離散の場合は和)が含まれるため,$\ell(\boldsymbol{\theta})$ の直接最大化は一般に困難となる.
完全データの結合分布は,観測モデルと潜在変数の事前分布の積として\[p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}) = p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\theta})\, p(\mathbf{Z} \mid \boldsymbol{\theta})\]と分解される.ここで $p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})$ はパラメータ $\boldsymbol{\theta}$ のもとでの潜在変数の事後分布である.
対数尤度に直接アプローチする代わりに,任意の確率分布 $q(\mathbf{Z})$ を補助分布として導入する.対数の凹性(イェンセンの不等式 $\log \mathbb{E}[f] \geq \mathbb{E}[\log f]$)を適用すると,\[\log p(\mathbf{X} \mid \boldsymbol{\theta})= \log \int q(\mathbf{Z}) \frac{p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{q(\mathbf{Z})}\, d\mathbf{Z}\geq \int q(\mathbf{Z}) \log \frac{p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{q(\mathbf{Z})}\, d\mathbf{Z}\]が成り立つ.右辺をエビデンス下界(ELBO: Evidence Lower BOund)と呼び,\[\mathcal{L}(q, \boldsymbol{\theta})= \mathbb{E}_{q(\mathbf{Z})}\!\left[\log p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})\right]- \mathbb{E}_{q(\mathbf{Z})}\!\left[\log q(\mathbf{Z})\right]\]と定義する.第1項は完全データ対数尤度の $q$ に関する期待値,第2項は $q(\mathbf{Z})$ のエントロピーである.
対数尤度と ELBO の差は KL ダイバージェンスで正確に表される:\[\log p(\mathbf{X} \mid \boldsymbol{\theta})= \mathcal{L}(q, \boldsymbol{\theta}) + \mathrm{KL}\!\left(q(\mathbf{Z}) \,\|\, p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})\right)\]ここで\[\mathrm{KL}\!\left(q \,\|\, p\right)= \int q(\mathbf{Z}) \log \frac{q(\mathbf{Z})}{p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})}\, d\mathbf{Z} \geq 0\]であるから,$\mathcal{L}(q, \boldsymbol{\theta}) \leq \log p(\mathbf{X} \mid \boldsymbol{\theta})$ が常に成立する.等号は $q(\mathbf{Z}) = p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})$ のとき,すなわち $\mathrm{KL} = 0$ のときに達成される.
現在のパラメータ推定値を $\boldsymbol{\theta}^{(t)}$ とする.E ステップでは,$q(\mathbf{Z})$ を現在のパラメータのもとでの事後分布に設定する:\[q(\mathbf{Z}) \leftarrow p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)})\]この操作により $\mathrm{KL} = 0$ となり,ELBO が対数尤度 $\ell(\boldsymbol{\theta}^{(t)})$ に一致する.
E ステップの実質的な計算は,Q関数(完全データ対数尤度の事後期待値)の評価である:\[Q(\boldsymbol{\theta},\, \boldsymbol{\theta}^{(t)})= \mathbb{E}_{\mathbf{Z} \mid \mathbf{X},\, \boldsymbol{\theta}^{(t)}}\!\left[\log p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})\right]= \int p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)}) \log p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})\, d\mathbf{Z}\]$Q$ 関数は $\boldsymbol{\theta}$ の関数であり,$\boldsymbol{\theta}^{(t)}$ は積分の重みを決める固定パラメータである.
完全データ対数尤度が指数型分布族に属する場合,E ステップでは十分統計量の事後期待値のみを計算すれば十分である.例えばガウス混合モデル(GMM)では,各データ点が各混合成分に属する負担率(responsibility)\[r_{nk} = p(z_n = k \mid x_n, \boldsymbol{\theta}^{(t)})= \frac{\pi_k^{(t)}\, \mathcal{N}(x_n \mid \mu_k^{(t)}, \Sigma_k^{(t)})} {\sum_{j=1}^{K} \pi_j^{(t)}\, \mathcal{N}(x_n \mid \mu_j^{(t)}, \Sigma_j^{(t)})}\]を計算することが E ステップに相当する.
M ステップでは,E ステップで得られた Q 関数を $\boldsymbol{\theta}$ について最大化し,パラメータを更新する:\[\boldsymbol{\theta}^{(t+1)}= \arg\max_{\boldsymbol{\theta}}\, Q(\boldsymbol{\theta},\, \boldsymbol{\theta}^{(t)})\]潜在変数が周辺化された後の Q 関数は,$\boldsymbol{\theta}$ に関して閉形式の解を持つことが多く(特に指数型分布族),直接の対数尤度最大化よりも大幅に容易となる.
GMMの例では,M ステップにおける各パラメータの更新式は\[N_k = \sum_{n=1}^{N} r_{nk}, \qquad\pi_k^{(t+1)} = \frac{N_k}{N}\]\[\mu_k^{(t+1)} = \frac{1}{N_k} \sum_{n=1}^{N} r_{nk}\, x_n\]\[\Sigma_k^{(t+1)} = \frac{1}{N_k} \sum_{n=1}^{N} r_{nk}\, (x_n - \mu_k^{(t+1)})(x_n - \mu_k^{(t+1)})^\top\]と解析的に得られる.
EMアルゴリズムの最も重要な性質は,各反復で対数尤度が単調に増加(非減少)することである.
まず,E ステップにより $q^{(t)}(\mathbf{Z}) = p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)})$ とおいたとき,ELBO は対数尤度に一致する:\[\mathcal{L}(q^{(t)}, \boldsymbol{\theta}^{(t)}) = \ell(\boldsymbol{\theta}^{(t)})\]次に,M ステップにより $\boldsymbol{\theta}^{(t+1)}$ は Q 関数を最大化するため,\[Q(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) \geq Q(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)})\]が成り立つ.対数尤度の分解式を用いると,\[\ell(\boldsymbol{\theta}^{(t+1)})= Q(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) + H(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)})\]ここで\[H(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})= -\int p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)}) \log p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})\, d\mathbf{Z}\]は KL ダイバージェンスの非負性から\[H(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) \geq H(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)})\]が従う.以上を組み合わせると,\[\ell(\boldsymbol{\theta}^{(t+1)})\geq Q(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) + H(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)})\geq Q(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)}) + H(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)})= \ell(\boldsymbol{\theta}^{(t)})\]が示され,各反復で対数尤度は単調非減少となる.
単調増加性と対数尤度の上界有界性から,EMアルゴリズムは必ず収束する.しかし収束先は一般に局所最適解であり,大域最適解が保証されない点に注意が必要である.
収束の固定点条件は\[\nabla_{\boldsymbol{\theta}}\, Q(\boldsymbol{\theta}, \boldsymbol{\theta}^*)\Big|_{\boldsymbol{\theta} = \boldsymbol{\theta}^*} = \mathbf{0}\]であり,この条件は $\nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}^*) = \mathbf{0}$ と等価であることが示される.すなわち,EMの固定点は観測データ対数尤度の停留点と一致する.
収束速度は一般に線形(一次)収束であり,完全データと不完全データの情報量の比に依存する.具体的には,欠損情報の割合が高いほど収束が遅くなる(フィッシャーの欠損情報原理).収束判定の実用的な基準としては,\[\left|\ell(\boldsymbol{\theta}^{(t+1)}) - \ell(\boldsymbol{\theta}^{(t)})\right| < \varepsilon\quad \text{or} \quad\left\|\boldsymbol{\theta}^{(t+1)} - \boldsymbol{\theta}^{(t)}\right\| < \varepsilon\]が広く用いられる.局所最適解への依存を緩和するための実践的な手法として,異なる初期値からの多重起動(multiple restarts),確率的 EM,焼きなまし EM などが提案されている.
一般化 EM(GEM):M ステップで Q 関数の完全な最大化の代わりに増加のみを要求する:\[Q(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) \geq Q(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)})\]単調増加性は依然として保証される.
変分 EM(Variational EM):$q(\mathbf{Z})$ を扱いやすい分布族(例:平均場近似 $q(\mathbf{Z}) = \prod_i q_i(z_i)$)に制限し,真の事後分布が解析的に求まらない場合に適用する.E ステップでは KL を最小化し,M ステップでは Q 関数の近似を最大化する.
モンテカルロ EM(MCEM):E ステップの期待値をモンテカルロサンプリングにより近似する:\[Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) \approx \frac{1}{S} \sum_{s=1}^{S} \log p(\mathbf{X}, \mathbf{Z}^{(s)} \mid \boldsymbol{\theta}),\quad \mathbf{Z}^{(s)} \sim p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)})\]事後分布からのサンプリングが可能だが解析的期待値が求まらない場合に有効である.
オンライン EM(確率的 EM):大規模データに対し,各反復でミニバッチを用いて十分統計量を逐次更新する:\[\tilde{s}^{(t+1)} = (1 - \eta_t)\, \tilde{s}^{(t)} + \eta_t\, s(\mathbf{X}_{\text{batch}}, \boldsymbol{\theta}^{(t)})\]ここで $\eta_t$ はステップサイズ,$s$ は十分統計量である.
EMアルゴリズムは,潜在変数を含む確率モデルにおける最尤推定のための強力かつ汎用的な枠組みである.その本質は,直接最大化が困難な対数尤度 $\ell(\boldsymbol{\theta})$ に対し,ELBO という代理目的関数を交互最適化することにある.
単調増加性 $\ell(\boldsymbol{\theta}^{(t+1)}) \geq \ell(\boldsymbol{\theta}^{(t)})$ は理論的に保証されており,アルゴリズムは必ず対数尤度の停留点(局所最適解または鞍点)へ収束する.一方,大域最適性の保証がない点,収束速度が線形に留まる場合がある点は主要な限界である.ガウス混合モデル,隠れマルコフモデル,因子分析,混合回帰など多様なモデルへの適用実績を持ち,変分 EM,モンテカルロ EM,オンライン EM などの拡張を通じて,より複雑・大規模な現代的機械学習問題へも広く応用されている.
Mathematics is the language with which God has written the universe.