混合正規分布モデル(GMM)とEMアルゴリズム
はじめに
混合正規分布モデル (GMM: Gaussian mixture model) は、凸凹した確率密度関数を複数の正規分布の和によって表現するもので、ここでは観測されたサンプル集合をもとに EM アルゴリズムで各正規分布のパラメータを推定する。
以下、太字はベクトルまたは行列を表し、行列 $\boldsymbol{A}$ の逆行列の転置を $\boldsymbol{A}^{-T}=\left(\boldsymbol{A}^{-1}\right)^T$ と略記する。
正規分布
$D$ 次元の正規分布は、平均ベクトル $\boldsymbol{\mu}\in\mathbb{R}^D$、非特異な分散共分散行列 $\boldsymbol{\Sigma}\in\mathbb{R}^{D\times D}$ をパラメータとして持ち、次のように表される。
\begin{eqnarray}
N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})
&=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right\} \label{Normal}
\end{eqnarray}
混合正規分布
混合正規分布を、正規分布 $N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})$ の線形結合として、次のように定義する。
\begin{eqnarray}
p(\boldsymbol{x}) &=& \sum_{k=1}^K \pi_k N(\boldsymbol{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \label{GMM}
\end{eqnarray}
ただし、結合係数 $\pi_k$ は確率に対応するもので、
\begin{eqnarray}
0\leq \pi_k\leq 1
\end{eqnarray}
かつ
\begin{eqnarray}
\sum_{k=1}^K \pi_k &=& 1 \label{sumpi1}
\end{eqnarray}
を満たすものとする。
混合正規分布モデルの対数尤度関数
観測されたサンプル $\boldsymbol{x}_n\in\mathbb{R}^D,\ n=1,2,3,\cdots N$ が式(\ref{GMM})の混合正規分布モデルから出てきたものと仮定し、その対数尤度関数
\begin{eqnarray}
\log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})
&=& \log \prod_{n=1}^N p(\boldsymbol{x}_n) \\
&=& \sum_{n=1}^N \log p(\boldsymbol{x}_n) \\
&=& \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)
\end{eqnarray}
を最大化 (正確には極大化) するパラメータ $\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}$ の組を推定する。
ここで $\boldsymbol{X}$ はサンプル集合 $\{\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots \boldsymbol{x}_N\}$ を行列で表したものである。
\begin{eqnarray}
\boldsymbol{X}
&=& \left(
\begin{array}{ccc}
& \boldsymbol{x}_1^T & \\
& \boldsymbol{x}_2^T & \\
& \vdots & \\
& \boldsymbol{x}_N^T & \\
\end{array}
\right) \in \mathbb{R}^{N\times D}
\end{eqnarray}
パラメータ推定の方針としては、対数尤度関数の極大点では、全ての $m=1,2,3,\cdots K$ で
\begin{eqnarray}
\left\{
\begin{array}{rcl}
\displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& 0\\
\displaystyle\frac{\partial}{\partial\boldsymbol{\Sigma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& 0\\
\displaystyle\frac{\partial}{\partial\pi_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& 0
\end{array}
\right.
\end{eqnarray}
が成り立つことを利用する。
負担率
$z_n \in \{1, 2, \cdots, K\}$ をサンプル $\boldsymbol{x}_n$ がどの正規分布から生成されたかを示す潜在変数とし、$P(z_n = m) = \pi_m$ とする。
サンプル $\boldsymbol{x}_n$ が $m$ 番目の正規分布から出てきた事後確率をベイズの定理で計算すると
\begin{eqnarray}
\gamma_{n,m}
&=& \frac{P(\boldsymbol{x}_n | z_n = m) P(z_n = m)}{\displaystyle\sum_{k=1}^K P(\boldsymbol{x}_n | z_n = k) P(z_n = k)} \nonumber\\
&=& \frac{N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \cdot \pi_m}{\displaystyle\sum_{k=1}^K N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \cdot \pi_k} \label{gamma}
\end{eqnarray}
となり、これを負担率と呼ぶことにする。
なお、サンプル $\boldsymbol{x}_n$ は $K$ 個ある正規分布のどれかから出てきたわけであるから、$m=1,2,3,\cdots K$ について $\gamma_{n,m}$ を足し合わせた全確率は当然のことながら 1 になる。
\begin{eqnarray}
\require{cancel}
\sum_{m=1}^K \gamma_{n,m}
&=& \sum_{m=1}^K \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{i=1}^K \pi_i N(\boldsymbol{x}_n|\boldsymbol{\mu}_i,\boldsymbol{\Sigma}_i)} \\
&=& \frac{\cancel{\displaystyle\sum_{m=1}^K \pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}}{\cancel{\displaystyle\sum_{m=1}^K \pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}} \\
&=& 1 \label{sumgamma}
\end{eqnarray}
対数尤度関数を極大化する $\boldsymbol{\mu}_m$ を求める
以下の導出では分母レイアウト記法を採用する。すなわち、スカラー関数をベクトルで微分した結果は列ベクトルとなる。
対数尤度関数を $\boldsymbol{\mu}_m$ で偏微分すると
\begin{eqnarray}
\frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})
&=& \frac{\partial}{\partial\boldsymbol{\mu}_m} \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \\
&=& \sum_{n=1}^N \frac{\partial}{\partial\boldsymbol{\mu}_m} \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \\
&& \text{一般に } \log'(f(\boldsymbol{x}))=\frac{f'(\boldsymbol{x})}{f(\boldsymbol{x})} \text{ であるから} \nonumber\\
&=& \sum_{n=1}^N \frac{\displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} \displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \\
&& k\neq m \text{ の場合、}\frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)=0 \text{ であるので} \nonumber\\
&=& \sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \label{dmu}
\end{eqnarray}
ここで式(\ref{Normal})から、正規分布 $N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)$ は
\begin{eqnarray}
N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)
&=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \nonumber \\
&& \boldsymbol{\mu}_m \text{ で偏微分することを考慮して } \boldsymbol{x} \text{ と } \boldsymbol{\mu} \text{ を入れ替え} \nonumber\\
&=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\}
\end{eqnarray}
これにより、式(\ref{dmu})の分子の $\pi_m$ の後ろは次のように書ける。
\begin{eqnarray}
&&\frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) && \nonumber\\
&&\hspace{3em}= \frac{\partial}{\partial\boldsymbol{\mu}_m} (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\
&&\hspace{3em}= (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}} \frac{\partial}{\partial\boldsymbol{\mu}_m} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\
&&\hspace{4em} \text{一般に } \exp'(f(\boldsymbol{x}))=\exp(f(\boldsymbol{x}))f'(\boldsymbol{x}) \text{ であるから} \nonumber\\
&&\hspace{3em}= (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\
&&\hspace{3em} \quad \times \frac{\partial}{\partial\boldsymbol{\mu}_m}
\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\
&&\hspace{3em}= N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)
\frac{\partial}{\partial\boldsymbol{\mu}_m}
\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\
&&\hspace{3em}= -\frac{1}{2} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)
\frac{\partial}{\partial\boldsymbol{\mu}_m}
\left\{(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \\
&&\hspace{4em} \text{一般に } \displaystyle\frac{\partial\boldsymbol{x}^T \boldsymbol{A} \boldsymbol{x}}{\partial\boldsymbol{x}}=(\boldsymbol{A}+\boldsymbol{A}^T)\boldsymbol{x} \text{ であるから}^* \nonumber\\
&&\hspace{3em}= -\frac{1}{2} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)
\left(\boldsymbol{\Sigma}_m^{-1}+\boldsymbol{\Sigma}_m^{-T}\right) (\boldsymbol{\mu}_m-\boldsymbol{x}_n) \\
&&\hspace{4em} \text{分散共分散行列 } \boldsymbol{\Sigma}_m \text{ は対称行列であるから } \boldsymbol{\Sigma}_m^{-1} \text{ も対称行列なので} \nonumber\\
&&\hspace{3em} \left(\boldsymbol{\Sigma}_m^{-1}+\boldsymbol{\Sigma}_m^{-T}\right)=2 \boldsymbol{\Sigma}_m^{-1} \nonumber\\
&&\hspace{3em}= - N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{\mu}_m-\boldsymbol{x}_n) \\
&&\hspace{3em}= N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{x}_n-\boldsymbol{\mu}_m)
\end{eqnarray}
* $\displaystyle\frac{\partial\boldsymbol{x}^T \boldsymbol{A} \boldsymbol{x}}{\partial\boldsymbol{x}}=(\boldsymbol{A}+\boldsymbol{A}^T)\boldsymbol{x}$ の 証明
よって式(\ref{dmu})は
\begin{eqnarray}
\frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})
&=& \sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \nonumber\\
&=& \sum_{n=1}^N \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{x}_n-\boldsymbol{\mu}_m)
\end{eqnarray}
\begin{eqnarray}
\frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})
&=& \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{x}_n-\boldsymbol{\mu}_m) \\
&& \text{極大点では勾配が } \boldsymbol{0} \text{ となるので} \nonumber\\
&=& \boldsymbol{0}
\end{eqnarray}
が成り立つには
\begin{eqnarray}
\sum_{n=1}^N \gamma_{n,m} \boldsymbol{\Sigma}_m^{-1} \boldsymbol{x}_n
&=& \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\Sigma}_m^{-1} \boldsymbol{\mu}_m
\end{eqnarray}
両辺に左から $\boldsymbol{\Sigma}_m$ を掛けて
\begin{eqnarray}
\sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n
&=& \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\mu}_m
\end{eqnarray}
よって
\begin{eqnarray}
\boldsymbol{\mu}_m
&=& \frac{\displaystyle\sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n}{\displaystyle\sum_{n=1}^N \gamma_{n,m}}
\end{eqnarray}
対数尤度関数を極大化する $\boldsymbol{\Sigma}_m$ を求める
分散共分散行列 $\boldsymbol{\Sigma}_m$ の逆行列を
\begin{eqnarray}
\boldsymbol{\Gamma}_m=\boldsymbol{\Sigma}_m^{-1} \label{G}
\end{eqnarray}
と置くと
\begin{eqnarray}
\frac{\partial}{\partial\boldsymbol{\Sigma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}_m)\ =\ \boldsymbol{0}
\end{eqnarray}
を満たす解は
\begin{eqnarray}
\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Gamma}_m^{-1})\ =\ \boldsymbol{0}
\end{eqnarray}
を解くことで計算できる。
\begin{eqnarray}
\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Gamma}_m^{-1})
&=& \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1}) \\
&=& \sum_{n=1}^N \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1}) \\
&& \text{一般に } \log'(f(\boldsymbol{x}))=\frac{f'(\boldsymbol{x})}{f(\boldsymbol{x})} \text{ であるから} \nonumber\\
&=& \sum_{n=1}^N \frac{\displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})} \\
&& k\neq m \text{ の場合、}\frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)=0 \text{ であるので} \nonumber\\
&=& \sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1})}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})} \label{dgamma}
\end{eqnarray}
式(\ref{Normal})から、正規分布 $N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})$ は
\begin{eqnarray}
N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})
&=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right\} \nonumber
\end{eqnarray}
であるから、式(\ref{dgamma})の分子の $\pi_m$ の後ろは次のように書ける。
\begin{eqnarray}
&& \frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1}) \nonumber\\
&=& (2\pi)^{-\frac{D}{2}} \frac{\partial}{\partial\boldsymbol{\Gamma}_m}|\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \label{dNdG}
\end{eqnarray}
一般に $\{f(\boldsymbol{x})g(\boldsymbol{x})\}'=f'(\boldsymbol{x})g(\boldsymbol{x})+f(\boldsymbol{x})g'(\boldsymbol{x})$ であることから、$\displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m}$ 以降は次のように書ける。
\begin{eqnarray}
&& \left(\frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}\right) \exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \nonumber\\
&& + |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \left[\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \right] \label{ddGm}
\end{eqnarray}
前半の ( ) 内は
\begin{eqnarray}
\frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}
&=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{-\frac{1}{2}} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m| \\
&& \text{一般に } \frac{\partial|\boldsymbol{X}|}{\partial\boldsymbol{X}}=|\boldsymbol{X}|\boldsymbol{X}^{-T} \text{ であるので}^* \nonumber\\
&=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{-\frac{1}{2}} |\boldsymbol{\boldsymbol{\Gamma}_m}|\boldsymbol{\boldsymbol{\Gamma}_m}^{-T} \\
&=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}\boldsymbol{\boldsymbol{\Gamma}_m}^{-T}
\end{eqnarray}
* $\displaystyle\frac{\partial|\boldsymbol{X}|}{\partial\boldsymbol{X}}=|\boldsymbol{X}|\boldsymbol{X}^{-T}$ の 証明
後半の [ ] 内は
\begin{eqnarray}
&&\hspace{-2em}\displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \\
&& \hspace{2em}\text{一般に } \exp'(f(\boldsymbol{x}))=\exp(f(\boldsymbol{x}))f'(\boldsymbol{x}) \text{ であるので} \nonumber\\
&& \hspace{2em}=\exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\}
\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \label{dexpdGm} \\
&& \hspace{2em}\text{一般に } \displaystyle\frac{\partial\boldsymbol{a}^T \boldsymbol{X} \boldsymbol{a}}{\partial\boldsymbol{X}}=\boldsymbol{a}\boldsymbol{a}^T \text{ であるから}^* \nonumber\\
&& \hspace{2em}=-\frac{1}{2}
\exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\}
(\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\end{eqnarray}
* $\displaystyle\frac{\partial\boldsymbol{a}^T \boldsymbol{X} \boldsymbol{a}}{\partial\boldsymbol{X}}=\boldsymbol{a}\boldsymbol{a}^T$ の 証明
よって式(\ref{ddGm})は
$\boldsymbol{E}_m = \exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\}$ と略記すると
\begin{eqnarray}
&& \hspace{-2em}\left(\frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}\right) \boldsymbol{E}_m
+ |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \left(\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \boldsymbol{E}_m \right) \\
&& \hspace{2em}=\frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \boldsymbol{\Gamma}_m^{-T} \boldsymbol{E}_m
- \frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \boldsymbol{E}_m (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T \\
&& \hspace{2em}\text{共通因子 } |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \boldsymbol{E}_m \text{ はスカラーなので括り出せて} \nonumber\\
&& \hspace{2em}=\frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \boldsymbol{E}_m
\left\{
\boldsymbol{\Gamma}_m^{-T}
- (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\right\}
\end{eqnarray}
式(\ref{dNdG})は、$(2\pi)^{-\frac{D}{2}}$ を式(\ref{ddGm})に掛けたものであるから
\begin{eqnarray}
&& \hspace{-2em}\frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1}) \nonumber\\
&& \hspace{2em}=(2\pi)^{-\frac{D}{2}} \frac{\partial}{\partial\boldsymbol{\Gamma}_m}|\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \nonumber\\
&& \hspace{2em}\text{式(\ref{ddGm})を代入して} \nonumber\\
&& \hspace{2em}=\frac{1}{2} (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Gamma}_m(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \nonumber\\
&& \hspace{2em}\quad \times \left\{
\boldsymbol{\Gamma}_m^{-T}
- (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\right\} \nonumber\\
&& \hspace{2em}\text{上記の } (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\{\cdots\} \text{ 部分は } N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \text{ なので} \nonumber\\
&& \hspace{2em}=\frac{1}{2} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)
\left\{
\boldsymbol{\Gamma}_m^{-T}
- (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\right\} \label{ddGN}
\end{eqnarray}
よって、式(\ref{dgamma})は
\begin{eqnarray}
&& \hspace{-2em}\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Gamma}_m^{-1}) \nonumber\\
&& \hspace{2em}=\sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1})}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})} \nonumber\\
&& \hspace{2em}\text{式(\ref{ddGN})を代入して} \nonumber\\
&& \hspace{2em}=\frac{1}{2}
\sum_{n=1}^N \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})}
\left\{
\boldsymbol{\Gamma}_m^{-T}
- (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\right\} \\
&& \hspace{2em}\boldsymbol{\Gamma}_m \text{ は } \boldsymbol{\Sigma}_m \text{ の逆行列なので } \boldsymbol{\Gamma}_m^{-T}=\boldsymbol{\Sigma}_m^T \nonumber\\
&& \hspace{2em}=\frac{1}{2}
\sum_{n=1}^N \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}
\left\{
\boldsymbol{\Sigma}_m^T
- (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\right\} \\
&& \hspace{2em}\text{式(\ref{gamma}) の } \gamma_{n,m} \text{ を使って書けば} \nonumber\\
&& \hspace{2em}=\frac{1}{2}
\sum_{n=1}^N \gamma_{n,m}
\left\{
\boldsymbol{\Sigma}_m^T
- (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\right\} \\
&& \hspace{2em}\boldsymbol{\Sigma}_m \text{ は対称行列なので } \boldsymbol{\Sigma}_m^T=\boldsymbol{\Sigma}_m \nonumber\\
&& \hspace{2em}=\frac{1}{2}
\sum_{n=1}^N \gamma_{n,m}
\left\{
\boldsymbol{\Sigma}_m
- (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\right\} \\
&& \hspace{2em}\text{極大点では勾配が } \boldsymbol{0} \text{ となるので} \nonumber\\
&& \hspace{2em}=\boldsymbol{0}
\end{eqnarray}
とすると
\begin{eqnarray}
\sum_{n=1}^N \gamma_{n,m}
\boldsymbol{\Sigma}_m
&=& \sum_{n=1}^N \gamma_{n,m}
(\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\end{eqnarray}
より、以下を得る。
\begin{eqnarray}
\boldsymbol{\Sigma}_m
&=& \frac{1}{\displaystyle\sum_{n=1}^N \gamma_{n,m}}
\sum_{n=1}^N \gamma_{n,m}
(\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T
\end{eqnarray}
対数尤度関数を極大化する $\pi_m$ を求める
混合係数 $\pi_m$ は $\displaystyle\sum_{m=1}^K \pi_m=1$ を満たさねばならないという強い制約があるため、Lagrange の未定乗数法により、対数尤度関数に $\lambda\left(\displaystyle\sum_{k=1}^K \pi_k -1\right)$ を加えたものを極大化することにする。
\begin{eqnarray}
&& \hspace{-2em}\frac{\partial}{\partial\pi_m} \left\{\log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) + \lambda\left(\sum_{k=1}^K \pi_k -1\right)\right\} \nonumber\\
&& \hspace{2em}=\left\{\frac{\partial}{\partial\pi_m} \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\right\}+\lambda \\
&& \hspace{2em}=\left\{\sum_{n=1}^N \frac{\partial}{\partial\pi_m} \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\right\}+\lambda \\
&& \hspace{2em}\text{一般に } \log'(f(\boldsymbol{x}))=\frac{f'(\boldsymbol{x})}{f(\boldsymbol{x})} \text{ であるから} \nonumber\\
&& \hspace{2em}=\sum_{n=1}^N \frac{\displaystyle\frac{\partial}{\partial\pi_m} \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} +\lambda \\
&& \hspace{2em}k\neq m \text{ の場合、}\frac{\partial}{\partial\pi_m} \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)=0 \text{ であるので} \nonumber\\
&& \hspace{2em}=\sum_{n=1}^N \frac{N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}+\lambda \\
&& \hspace{2em}\text{式(\ref{gamma}) の } \gamma_{n,m} \text{ を使って書けば} \nonumber\\
&& \hspace{2em}=\frac{1}{\pi_m} \sum_{n=1}^N \gamma_{n,m}+\lambda \\
&& \hspace{2em}=0
\end{eqnarray}
とすると、両辺に $\pi_m$ を掛けて
\begin{eqnarray}
\sum_{n=1}^N \gamma_{n,m} + \pi_m \lambda &=& 0
\end{eqnarray}
より
\begin{eqnarray}
\pi_m &=& -\frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m} \label{pim}
\end{eqnarray}
であるが、式(\ref{sumpi1})より
\begin{eqnarray}
\sum_{m=1}^K \pi_m &=& 1
\end{eqnarray}
であることを考えると
\begin{eqnarray}
\sum_{m=1}^K \pi_m &=& - \sum_{m=1}^K \frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m} &=& 1
\end{eqnarray}
でなければならず、変形してゆくと
\begin{eqnarray}
- \sum_{m=1}^K \frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m}
&=& - \frac{1}{\lambda} \sum_{m=1}^K \sum_{n=1}^N \gamma_{n,m} \\
&=& - \frac{1}{\lambda} \sum_{n=1}^N \sum_{m=1}^K \gamma_{n,m} \\
&& \text{式(\ref{sumgamma}) より } \sum_{m=1}^K \gamma_{n,m}=1 \text{ であるから} \nonumber\\
&=& - \frac{1}{\lambda} \sum_{n=1}^N 1 \\
&=& - \frac{N}{\lambda} \\
&=& 1
\end{eqnarray}
から
\begin{eqnarray}
\lambda &=& -N
\end{eqnarray}
従って式(\ref{pim})は
\begin{eqnarray}
\pi_m
&=& -\frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m} \nonumber\\
&=& \frac{1}{N} \sum_{n=1}^N \gamma_{n,m}
\end{eqnarray}
となる。
EMアルゴリズム
これまでに得た結果をまとめると
\begin{eqnarray}
\gamma_{n,m} &=& \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \nonumber
\end{eqnarray}
の下で対数尤度関数を極大化するパラメータは
\begin{eqnarray}
\left\{
\begin{array}{rcl}
\boldsymbol{\mu}_m
&=& \frac{\displaystyle\sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n}{\displaystyle\sum_{n=1}^N \gamma_{n,m}}, \\
\boldsymbol{\Sigma}_m
&=& \displaystyle\frac{1}{\displaystyle\sum_{n=1}^N \gamma_{n,m}} \displaystyle\sum_{n=1}^N \gamma_{n,m} (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T, \\
\pi_m
&=& \displaystyle\frac{1}{N} \sum_{n=1}^N \gamma_{n,m}
\end{array}
\right.\nonumber
\end{eqnarray}
であるが、$\displaystyle\sum_{n=1}^N \gamma_{n,m}$ は $m$ 番目の正規分布から出てきたサンプル数の推定値に相当するので、これを $N_m$ と書くことにすると、もう少し簡単になる。
\begin{eqnarray}
\left\{
\begin{array}{rcl}
\boldsymbol{\mu}_m
&=& \displaystyle\frac{1}{N_m} \sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n,\\
\boldsymbol{\Sigma}_m
&=& \displaystyle\frac{1}{N_m} \sum_{n=1}^N \gamma_{n,m} (\boldsymbol{x}_n-\boldsymbol{\mu}_m)(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T,\\
\pi_m
&=& \displaystyle\frac{N_m}{N}
\end{array}
\right.
\end{eqnarray}
しかし $\gamma_{n,m}$ の計算式の中に $\boldsymbol{\mu}_m, \boldsymbol{\Sigma}_m, \pi_m$ が、 $\boldsymbol{\mu}_m, \boldsymbol{\Sigma}_m, \pi_m$ の計算式の中に $\gamma_{n,m}$ があるため、この連立方程式は非線形でスパッと解くことができない。そこで次のように初期値から出発して、E-step と M-step を交互に反復する EM アルゴリズムにより、ジワジワと最適解に近づけてゆくようにする。
EMアルゴリズムの手順
- 初期値 $\boldsymbol{\mu}_m^{(0)}, \boldsymbol{\Sigma}_m^{(0)}, \pi_m^{(0)}$ を適当に設定し、$i=0$ とする
- $\boldsymbol{\mu}_m^{(i)}, \boldsymbol{\Sigma}_m^{(i)}, \pi_m^{(i)}$ を使い、$\gamma_{n,m}^{(i)}, N_m^{(i)}$ を計算する (E-step)
\begin{eqnarray}
\left\{
\begin{array}{rcl}
\gamma_{n,m}^{(i)} &=& \frac{\pi_m^{(i)} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m^{(i)},\boldsymbol{\Sigma}_m^{(i)})}{\displaystyle\sum_{k=1}^K \pi_k^{(i)} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k^{(i)},\boldsymbol{\Sigma}_k^{(i)})} \\
N_m^{(i)} &=& \displaystyle\sum_{n=1}^N \gamma_{n,m}^{(i)}
\end{array}
\right.
\end{eqnarray}
- $\gamma_{n,m}^{(i)}, N_m^{(i)}$ を使い、パラメータ $\boldsymbol{\mu}_m^{(i+1)}, \boldsymbol{\Sigma}_m^{(i+1)}, \pi_m^{(i+1)}$ を更新する (M-step)
\begin{eqnarray}
\left\{
\begin{array}{rcl}
\boldsymbol{\mu}_m^{(i+1)}
&=& \displaystyle\frac{1}{N_m^{(i)}} \sum_{n=1}^N \gamma_{n,m}^{(i)} \boldsymbol{x}_n, \\
\boldsymbol{\Sigma}_m^{(i+1)}
&=& \displaystyle\frac{1}{N_m^{(i)}} \sum_{n=1}^N \gamma_{n,m}^{(i)} (\boldsymbol{x}_n-\boldsymbol{\mu}_m^{(i)})(\boldsymbol{x}_n-\boldsymbol{\mu}_m^{(i)})^T, \\
\pi_m^{(i+1)}
&=& \displaystyle\frac{N_m^{(i)}}{N}
\end{array}
\right. \label{update}
\end{eqnarray}
- パラメータ $\boldsymbol{\mu}_m^{(i+1)}, \boldsymbol{\Sigma}_m^{(i+1)}, \pi_m^{(i+1)}$ で対数尤度関数
\begin{eqnarray}
\log p(\boldsymbol{X}|\boldsymbol{\pi}^{(i+1)},\boldsymbol{\mu}^{(i+1)},\boldsymbol{\Sigma}^{(i+1)})
&=& \sum_{n=1}^N \log \sum_{k=1}^K \pi_k^{(i+1)} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k^{(i+1)},\boldsymbol{\Sigma}_k^{(i+1)})
\end{eqnarray}
を計算し、十分収束したら終了する。さもなくば $i\leftarrow i+1$ としてステップ 2 に戻る
注意
EM アルゴリズムは、任意の初期値から出発して広域最適解に収束するとは限らない。
参考
対数尤度関数が十分に滑らかな場合は、式(\ref{update}) の $\boldsymbol{\Sigma}_m^{(i+1)}$ の式に使う $\boldsymbol{\mu}_m^{(i)}$ を更新後の $\boldsymbol{\mu}_m^{(i+1)}$ にすることもできる。
EMアルゴリズムの収束性
EMアルゴリズムの各反復で対数尤度が単調に増加（正確には非減少）することを示す。
対数尤度の下界
$\boldsymbol{\theta}=\{\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}\}$ をパラメータの集合、$\boldsymbol{Z}=\{z_1, z_2, \cdots, z_N\}$ を潜在変数の集合とする。
潜在変数 $\boldsymbol{Z}$ に関する任意の確率分布 $q(\boldsymbol{Z})$ に対して、対数尤度は次のように変形できる。
\begin{eqnarray}
\log p(\boldsymbol{X}|\boldsymbol{\theta})
&=& \log p(\boldsymbol{X}|\boldsymbol{\theta}) \underbrace{\sum_{\boldsymbol{Z}} q(\boldsymbol{Z})}_{=1} \\
&=& \sum_{\boldsymbol{Z}} q(\boldsymbol{Z}) \log p(\boldsymbol{X}|\boldsymbol{\theta}) \\
&& \text{条件付き確率の定義より } p(\boldsymbol{X},\boldsymbol{Z}|\boldsymbol{\theta}) = p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta}) p(\boldsymbol{X}|\boldsymbol{\theta}) \text{ であるから} \nonumber \\
&=& \sum_{\boldsymbol{Z}} q(\boldsymbol{Z}) \log \frac{p(\boldsymbol{X},\boldsymbol{Z}|\boldsymbol{\theta})}{p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta})} \\
&=& \sum_{\boldsymbol{Z}} q(\boldsymbol{Z}) \log \frac{p(\boldsymbol{X},\boldsymbol{Z}|\boldsymbol{\theta})}{q(\boldsymbol{Z})} \cdot \frac{q(\boldsymbol{Z})}{p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta})} \\
&=& \sum_{\boldsymbol{Z}} q(\boldsymbol{Z}) \log \frac{p(\boldsymbol{X},\boldsymbol{Z}|\boldsymbol{\theta})}{q(\boldsymbol{Z})}
+ \sum_{\boldsymbol{Z}} q(\boldsymbol{Z}) \log \frac{q(\boldsymbol{Z})}{p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta})} \\
&=& \mathcal{L}(q,\boldsymbol{\theta}) + \mathrm{KL}(q \| p)
\end{eqnarray}
ここで $\mathcal{L}(q,\boldsymbol{\theta})$ と $\mathrm{KL}(q \| p)$ は次のように定義した。
\begin{eqnarray}
\mathcal{L}(q,\boldsymbol{\theta})
&\triangleq& \sum_{\boldsymbol{Z}} q(\boldsymbol{Z}) \log \frac{p(\boldsymbol{X},\boldsymbol{Z}|\boldsymbol{\theta})}{q(\boldsymbol{Z})} \\
\mathrm{KL}(q \| p)
&\triangleq& \sum_{\boldsymbol{Z}} q(\boldsymbol{Z}) \log \frac{q(\boldsymbol{Z})}{p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta})}
\end{eqnarray}
$\mathcal{L}(q,\boldsymbol{\theta})$ は変分下界（ELBO: Evidence Lower BOund）または負の変分自由エネルギーと呼ばれる。
直接最大化が困難な対数尤度 $\log p(\boldsymbol{X}|\boldsymbol{\theta})$ の代わりに、この下界を最大化するのがEMアルゴリズムの基本的なアイデアである。
$\mathrm{KL}(q \| p)$ は $q(\boldsymbol{Z})$ と事後分布 $p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta})$ の間の Kullback-Leibler ダイバージェンス（KLダイバージェンス）であり、常に $\mathrm{KL}(q \| p) \geq 0$ が成り立つ（等号は $q(\boldsymbol{Z}) = p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta})$ のとき）。
したがって $\mathcal{L}(q,\boldsymbol{\theta})$ は対数尤度の下界となる。
\begin{eqnarray}
\log p(\boldsymbol{X}|\boldsymbol{\theta}) &\geq& \mathcal{L}(q,\boldsymbol{\theta})
\end{eqnarray}
この不等式は、$\mathcal{L}(q,\boldsymbol{\theta})$ が常に対数尤度以下であることを意味する。
対数尤度 $\log p(\boldsymbol{X}|\boldsymbol{\theta})$ を直接最大化するのは困難（$\log$ の中に $\sum$ があるため）だが、
下界 $\mathcal{L}(q,\boldsymbol{\theta})$ を最大化することで間接的に対数尤度を増加させることができる。
これがEMアルゴリズムの基本的なアイデアである。
E-stepの解釈
E-step では、現在のパラメータ $\boldsymbol{\theta}^{(i)}$ を固定して、下界 $\mathcal{L}(q,\boldsymbol{\theta}^{(i)})$ を最大化する $q(\boldsymbol{Z})$ を求める。
$\log p(\boldsymbol{X}|\boldsymbol{\theta}^{(i)}) = \mathcal{L}(q,\boldsymbol{\theta}^{(i)}) + \mathrm{KL}(q \| p)$ において、左辺は $q$ に依存しないので、$\mathcal{L}$ を最大化することは $\mathrm{KL}(q \| p)$ を最小化することと等価である。
$\mathrm{KL}(q \| p) \geq 0$ であり、等号成立は $q(\boldsymbol{Z}) = p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta}^{(i)})$ のときなので、E-step の最適解は
\begin{eqnarray}
q^{(i)}(\boldsymbol{Z}) &=& p(\boldsymbol{Z}|\boldsymbol{X},\boldsymbol{\theta}^{(i)})
\end{eqnarray}
である。このとき $\mathrm{KL}(q^{(i)} \| p) = 0$ となり、下界が対数尤度に一致する。
\begin{eqnarray}
\mathcal{L}(q^{(i)},\boldsymbol{\theta}^{(i)}) &=& \log p(\boldsymbol{X}|\boldsymbol{\theta}^{(i)})
\end{eqnarray}
M-stepの解釈
M-step では、E-step で求めた $q^{(i)}(\boldsymbol{Z})$ を固定して、下界 $\mathcal{L}(q^{(i)},\boldsymbol{\theta})$ を最大化するパラメータ $\boldsymbol{\theta}^{(i+1)}$ を求める。
\begin{eqnarray}
\boldsymbol{\theta}^{(i+1)} &=& \arg\max_{\boldsymbol{\theta}} \mathcal{L}(q^{(i)},\boldsymbol{\theta})
\end{eqnarray}
最大化により $\mathcal{L}(q^{(i)},\boldsymbol{\theta}^{(i+1)}) \geq \mathcal{L}(q^{(i)},\boldsymbol{\theta}^{(i)})$ が成り立つ。
対数尤度の単調増加性
以上より、各反復での対数尤度の変化を追跡すると
\begin{eqnarray}
\log p(\boldsymbol{X}|\boldsymbol{\theta}^{(i+1)})
&=& \mathcal{L}(q^{(i)},\boldsymbol{\theta}^{(i+1)}) + \mathrm{KL}(q^{(i)} \| p(\cdot|\boldsymbol{X},\boldsymbol{\theta}^{(i+1)})) \\
&\geq& \mathcal{L}(q^{(i)},\boldsymbol{\theta}^{(i+1)}) \\
&& \text{M-step で下界を最大化したので} \nonumber \\
&\geq& \mathcal{L}(q^{(i)},\boldsymbol{\theta}^{(i)}) \\
&& \text{E-step で } \mathrm{KL}=0 \text{ としたので} \nonumber \\
&=& \log p(\boldsymbol{X}|\boldsymbol{\theta}^{(i)})
\end{eqnarray}
したがって、EMアルゴリズムの各反復で対数尤度は単調に増加（非減少）する。
\begin{eqnarray}
\log p(\boldsymbol{X}|\boldsymbol{\theta}^{(i+1)}) &\geq& \log p(\boldsymbol{X}|\boldsymbol{\theta}^{(i)})
\end{eqnarray}
対数尤度は上に有界（確率なので $p(\boldsymbol{X}|\boldsymbol{\theta}) \leq 1$ より $\log p \leq 0$）であるため、単調増加数列は収束する。
ただし、収束先が大域的最適解である保証はなく、局所最適解や鞍点に収束する可能性がある。
参考文献・関連記事
- C.M. Bishop, "Pattern Recognition and Machine Learning", Springer, 2006
- 行列微分の公式集