【機械学習】EMアルゴリズムをゼロから理解する

EMアルゴリズム(expectation maximization algorithm)は、潜在変数を有する確率モデルのパラメータの最尤推定値を求めるための近似的アルゴリズムです。

一般的に、潜在変数を仮定するモデルでは潜在変数は入手できず、観測値しか入手できませんが、そのような状況であっても近似的な手法を用いることで、モデルパラメータの最尤推定値を得ることができることから、応用上でも非常に多く使われています。特に、混合モデル、隠れマフコフモデル(HMM)、線形動的システム(Linear Dynamic System)などのモデルパラメータの学習に使われるほか、欠損値の多いデータでモデルを学習させたい時にも、EMアルゴリズムは利用されます。

EMアルゴリズムの全体像

この記事では、EMアルゴリズムの導出と解説を行なっていきますが、まず最初に、目指すべきアルゴリズムの全体像を見ていきましょう。

と、その前に、EMアルゴリズムの導出に入る前に、EMアルゴリズムについて学習していると度々登場することになる、完全データ集合と不完全データ集合という用語をまず導入します。これらの定義は、このようになっています。

完全データ集合と不完全データ集合

EMアルゴリズムのEステップとMステップ

EMアルゴリズムは、EステップとMステップの2ステップから構成されるアルゴリズムで、データから最適なパラメータを学習する手法です。

それぞれのステップで行う内容は、下記にまとめてあります。

\begin{equation}
\mathcal{Q} (\theta , \theta^{old}) = \sum_Z p(Z | X, \theta^{old}) ln p (X, Z | \theta)
\end{equation}

\theta^{new} = argmax_{\theta } \space \mathcal{L}

これらがEMアルゴリズムの全体像です。

とはいえ、これだけでは全く理解できないと思うので、以降でこの導出をやっていきます。

EMアルゴリズムの導出

それでは、EMアルゴリズムの導出に入っていきます。

今、観測値$\bm{X}$が得られているとします。

このとき、最尤推定の考え方だと、最尤推定値$\bm{\hat{\theta}}$は、対数尤度関数$ln p(\bm{X} | \theta)$を最小化するような$\bm{\theta}$を求めることになります。

\begin{equation}
\bm{\hat{\theta}} = \mathop{\arg\max}\limits_{\theta} 
ln p(\bm{X} | \theta)
\end{equation}

しかし、EMアルゴリズムが力を発揮するのは、対数尤度関数 $ln p(\bm{X} | \theta)$が複雑な形をしており、パラメータに関する偏微分などで最適解を求められないような場合です。

隠れ変数の導入

対数尤度関数 $ln p(\bm{X} | \theta)$が複雑な形をしている場合、どうすれば良いでしょうか。

EMアルゴリズムでは、対数尤度関数の形に潜在変数を導入することで、この問題にアプローチしている。隠れ変数を入れるとはどういうことか？と思う人も多いと思うが、今は天下り的にそういうものかと思ってください。

今、すべての観測変数$\bm{X} $と潜在変数$\bm{Z} $を示した確率モデルを考えます。

この完全データ同時分布を、$p(\bm{X}, \bm{Z} | \bm{\theta})$とすると、尤度関数$p(\bm{X} | \bm{\theta})$は、確率変数の周辺化を用いてこのように書くことができます。

\begin{equation}
\begin{split}
p(\bm{X} | \bm{\theta}) = \sum_Z p(\bm{X}, \bm{Z} | \bm{\theta})
\end{split}
\end{equation}

\begin{equation}
\begin{split}
ln p(\bm{X} | \bm{\theta}) = ln \sum_Z p(\bm{X}, \bm{Z} | \bm{\theta})
\end{split}
\end{equation}

(3)式により、最尤推定によって最大化したい対数尤度関数を、潜在変数を入れた形に書き直すことができました。

少し奇妙というか、どういう意図を持ってこのような式変形をしているのか、ということは、今後明らかになるので今は、天下り的にでも理解をしてください。

では、(3)式の右辺を最大化するにはどうしたら良いか、ということですが、ここでイェンゼンの不等式を利用することになります。

イェンゼンの不等式についてはこちらもご参考ください。

機械学習で登場するイェンゼンの不等式を学ぶ

イェンゼンの不等式について学んでいきます。イェンゼンの不等式は、凸関数において成り立つ不等式です。 なかなか抽象的な面もありますが、一般系な凸関数において成り立つので、汎用性が高く、機械学習や統計学の勉強をすると度々登場 […]

イェンゼンの不等式を利用するために、潜在変数を用いた(3)式で表した対数尤度関数を下記のように変形します。

ここでは、潜在変数$\bm{Z}$の確率分布$q(\bm{Z})$を用いて下記のように書き直すことができます。

\begin{equation}
\begin{split}
ln \sum_Z p(\bm{X}, \bm{Z} | \bm{\theta}) = ln \sum_Z q(\bm{Z}) 
\frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{q(\bm{Z}) }
\end{split}
\end{equation}

(4)式の変形は、分母と分子に$q(\bm{Z})$をかけただけの式変形です。

なぜこんなことをするのか？と面食らってしまうかもしれませんが、(4)式の右辺をみると、イェンセンの不等式を利用できる形になっていることがわかります。

こういうトリッキーな式変形は天下りに導入されると理解できずに、勉強する気が失せてしまいますよね。(4)式の変形する意味がいまいちピンとこない人は、イェンゼンの不等式の定義について理解してから、戻ってくることをお勧めします。

イェンゼンの不等式を用いると、(4)式は下限値を評価することができます。

\begin{equation}
\begin{split}
ln \sum_Z q(\bm{Z}) 
\frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{q(\bm{Z}) } 
&\geq 
\sum_Z q(\bm{Z}) ln
\frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{q(\bm{Z}) } \\
&= \mathcal{L}(q(\bm{Z}), \bm{\theta})
\end{split}
\end{equation}

ここで、下限値を$\mathcal{L}(q(\bm{Z}), \bm{\theta})$と書き直しました。

これは下限値の取り扱いを楽にするために、書き直しているだけなので、そのようなものだと思えば大丈夫です。

ここで、今、最大化したい対数尤度関数と、(5)式によって登場した下限値$\mathcal{L}(q(\bm{Z}), \bm{\theta})$はどれくらい離れているか見てみましょう。

\begin{equation}
\begin{split}
ln &p(\bm{X} | \bm{\theta}) -  \mathcal{L}(q(\bm{Z}), \bm{\theta})  \\
&= 
ln p(\bm{X} | \bm{\theta})  - 
\sum_Z q(\bm{Z}) ln
\frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{q(\bm{Z}) } \\
&= 
ln p(\bm{X} | \bm{\theta})  \sum_Z q(\bm{Z})  -  
\sum_Z q(\bm{Z}) ln
\frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{q(\bm{Z}) } \\
&= 
\sum_Z q(\bm{Z}) 
\biggr \{
ln p(\bm{X} | \bm{\theta}) - 
ln \frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{q(\bm{Z}) }
\biggr \} \\
&= 
\sum_Z q(\bm{Z}) 
\biggr \{
ln p(\bm{X} | \bm{\theta}) - 
ln 
\frac{p(\bm{Z} |\bm{X},  \bm{\theta})  p(\bm{X}| \bm{\theta}) }{q(\bm{Z}) }
\biggr \} \\
&= 
\sum_Z q(\bm{Z}) 
\biggr \{
ln p(\bm{X} | \bm{\theta}) - 
ln  p(\bm{Z} |\bm{X},  \bm{\theta}) - ln  p(\bm{X}| \bm{\theta}) + ln q(\bm{Z})
\biggr \} \\
&= 
- \sum_Z q(\bm{Z}) 
\biggr \{
\frac{p(\bm{Z} |\bm{X},  \bm{\theta}) }{q(\bm{Z}) }
\biggr \} \\
&= 
KL(q(\bm{Z}) || p(\bm{Z} |\bm{X},  \bm{\theta}) )
\end{split}
\end{equation}

となりました。

天下り的ですが、求めたい対数尤度関数と、そのイェンゼンの不等式から導出された下限値は、カルバックライブラー情報量を用いて、$KL(q(\bm{Z}) || p(\bm{Z} |\bm{X}, \bm{\theta}) )$のようになることがわかりました。

最後の式変形は、カルバックライブラー(KL情報量)の定義そのものを利用しています。

情報理論や確率変数のKLダイバージェンス(カルバック・ライブラー情報量)を解説

KLダイバージェンス(Kullback-Leibler divergence, KL情報量)は、2つの確率分布の距離を表す統計量として、統計学や機械学習分野で頻出の統計量となっています。 KLダイバージェンスを$D_{K […]

図にするとこのようなイメージです。

今回の目的は、対数尤度関数$p(\bm{X} | \bm{\theta})$を最大化することです。そして、今回変更することができるパラメータは、$q(\bm{Z})$と$\bm{\theta}) $の2つがあるわけです。

そのためには、KLダイバージェンスを最小化することと、下限値を引き上げることができます。

KLダイバージェンスの最小化について

まず、最初のステップとしてKLダイバージェンスの最小化に取り組んでいきます。KLダイバージェンスは、2つの分布の距離を示す量な統計量であり、2つの確率分布が一致するときに最小化されます。

このため、KLダイバージェンスを最小化するような$q^{new}(\bm{Z})$が次のように選べば良いことになります。

\begin{equation}
q^{new}(\bm{Z}) = p(\bm{Z} |\bm{X},  \bm{\theta})
\end{equation}

これにより、 新しい関数$q^{new}(\bm{Z})$は、潜在変数における、データ$X\bm{X}$とパラメータ$\bm{theta}$の事後分布に置き換えれば良いということになります。

下限値の最大化について

先ほどのKLダイバージェンスのところで、新しい関数$q^{new}(\bm{Z})$ は求まったので、続いては新しい関数$q^{new}(\bm{Z})$を利用して、下限値を最大化していきましょう。

下限値は、(5)式を用いるとこのようになります。

\begin{equation}
\begin{split}
\mathcal{L}(q^{new}(\bm{Z}), \bm{\theta}) &= 
\sum_Z q^{new}(\bm{Z}) ln
\frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{q^{new}(\bm{Z}) }  \\
&= 

\sum_Z p(\bm{Z} |\bm{X},  \bm{\theta}^{old}) ln
\frac{p(\bm{X}, \bm{Z} | \bm{\theta}) }{p(\bm{Z} |\bm{X},  \bm{\theta}^{old}) }  \\
&=
\sum_Z p(\bm{Z} |\bm{X},  \bm{\theta}^{old}) ln
p(\bm{X}, \bm{Z} | \bm{\theta} )+ const.  \\
&= 
\mathcal{Q}(\bm{\theta}, \bm{\theta}^{old}) + const.

\end{split}
\end{equation}

となります。(8)式を最終行は、$\mathcal{Q}(\bm{\theta}, \bm{\theta}^{old})$を導入しました。

$\mathcal{Q}(\bm{\theta}, \bm{\theta}^{old})$ は、EMアルゴリズムにおいて、重要となる式なので一度まとめておきます。

\begin{equation}
\mathcal{Q}(\bm{\theta}, \bm{\theta}^{old}) = 

\sum_Z p(\bm{Z} |\bm{X},  \bm{\theta}^{old}) ln
p(\bm{X}, \bm{Z} | \bm{\theta} )

\end{equation}

(9)式の$\mathcal{Q}$関数は、意味のわからない式になっているかもしれませんが、これは期待値の定義から、潜在変数$\bm{Z}$の事後分布における対数尤度関数の期待値になっています。

確率変数の期待値について徹底解説

確率変数の期待値は、機械学習の学習をする中で、非常によく登場する概念です。確率変数のモーメントとも繋がりがあり、統計学や機械学習の基本書を読んでも、最新の論文を読んでも登場するので、期待値の概念について正しく理解すること […]

よく論文等のEMアルゴリズムを利用した計算では、EMアルゴリズムの計算をする場面になると、突然潜在変数$\bm{Z}$の事後分布における対数尤度関数の期待値を計算し始めますが、それは最小化したい関数がまさにこの関数だったからなのです。

ここまでくると、EMアルゴリズムの全体像が結構見えてくるのではないのでしょうか？

EMアルゴリズムの計算では、結局のところ潜在変数$\bm{Z}$の事後分布における対数尤度関数の期待値を導出するEステップと、それをパラメタ$\bm{\theta}$で最適化するMステップなのですが、理論的背景はこの記事でまとめたのが理由となっています。

EMアルゴリズム自体は非常に機械学習の登竜門と言われるほど基本的でありながら理解するまで時間がかかると思うので、混合ガウス分布(GMM)のパラメタ推定などで何度も演習してみると、なんとなくやっていることが身についてくると思います。