EMアルゴリズムをもとに潜在変数モデル群の違いを整理する（2/4回目）

本記事では4回に渡り、筆者が親しみのある潜在変数モデル群を題材に、共通の計算手法を用いてそれらのモデルを比較し、その性質の違いをまとめます。それにより、モデル間の性質や違いを理解する助けとなることを目指します。

第2回目は、離散値をとる潜在変数モデルのQ関数の形を整理します。
※本記事は[第1回目](https://qiita.com/PingpongMachine/items/c5e5d7b4e69556a2b686)の続編です。第1回目は、EMアルゴリズムと本記事で扱う潜在変数モデル群を説明しています。


# 本記事の内容
本記事では、K-meansとGMMのモデルのQ関数を形を示します。ここで、GMMは確率モデルであるため、確率モデルという節を設けて、確率分布$p(x,z|\theta)$の形を示します。また、$t$回目の計算の時点でQ関数を得るためには、パラメータ$\theta$と観測データ$x$をから潜在変数$z$の情報の計算値、負担率を事前に得る必要があります。そのため、負担率という節を設けて、計算を示します。

これらの潜在変数は離散値を取ります。これと、潜在変数が連続値を取るモデルと対応をとるため、カテゴリカル変数としての整数値として扱います。他の表現としてワンホット表現があり、この表現で説明する資料も多いと思います。表現の違いにご留意ください。

# K-means（ハード割り当て極限として）

### Q関数
K-means は厳密な意味での EM アルゴリズムではありませんが、GMM のハード割り当て極限として扱うことができます。
その場合、Q 関数は各データ点のクラスタ割り当てに基づく平方誤差の和となります。
```math
Q(\theta; \theta^{(t)}) = -\sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk}^{(t)} \| x_n - \mu_k \|^2
```
ここで$\theta= \\{\mu_k \\}$であり$\mu_k$は$k$番目のクラスタの中心ベクトルとなります。ここで、$r_{nk}^{(t)}$は、パラメータ$\theta$と観測データ$x$をから潜在変数$z$の情報の計算値のため、負担率を示しています。
### 負担率
負担率$r_{nk}^{(t)}$はクロネッカーのデルタ関数により以下の式になります。
```math
r_{nk}^{(t)}=\delta_{z_n^{(t)},k}
```
潜在変数$z_n^{(t)}\in \\{1,..,K\\}$は$k$個のクラスタのいずれかから得られると$r_{nk}^{(t)}=1$になることを意味します。
さて、負担率$r_{nk}^{(t)}$の計算のためには、$z_{n}^{(t)}$の値が必要です。K-means は以下の形で計算します。
```math
z_{n}^{(t)} = \arg \min_{k'} ||x_n-\mu_{k'}^{(t-1)}||^2 \\
```
上式は、データ$x_n$と$\mu^{(t)}_k$の２乗誤差が最小となるクラスタ$k$を$z_n$の値とするということを意味します。


<!--
したがって、Q関数の$n$番目の部分に関しては、$K$に対する和の計算が残りますが、潜在変数$z_n$の値が$k$と決まれば、$r_{nk}^{(t)}$より、$k$番目の計算のみになります。
-->







# GMM
### 生成モデル

1.潜在変数の生成
```math
p(z_n=k) = \pi_k,　\sum_k^K\pi_k= 1
```
2.観測データの生成
```math
p(x_n|z_n=k,\theta) = \mathcal{N}(x_n \mid \mu_k, \Sigma_k)
```
ここで$\theta=\\{(\mu_k,\Sigma_k)\\}$です。$z_n$の値が$k$に決まると、パラメータ$\theta$も$k$番目の変数に絞り込まれ、$k$番目のクラスタから$\mathcal{N}(x_n \mid \mu_k, \Sigma_k)$のガウス分布に従って観測データ$x_n$が生成されるということを意味します。

$z_n$の総和をとることで$x_n$の周辺分布が得られます。
```math
p(x_n|\theta) =\sum_k^K p(z_n=k)p(x_n|z_n=k,\theta) =\sum_k^K \pi_k \mathcal{N}(x_n \mid \mu_k, \Sigma_k)
```

### Q関数
対数尤度に基づき、Q 関数は次のようになります。
```math
Q(\theta; \theta^{(t)}) =\sum_{n=1}^{N} \sum_{k=1}^{K} p(z_n=k|x_n,\theta^{(t-1)}) \log \left[p(z_n=k)p(x_n|z_n=k,\theta)\right] =\sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk}^{(t)} \log \left[ \pi_k\mathcal{N}(x_n \mid \mu_k, \Sigma_k) \right]
```
ここで$p(z_n=k|x_n,\theta^{(t)})$はパラメータ$\theta$と観測データ$x$をから潜在変数$z$の情報の計算値のため、負担率を示しており、K-meansと対応を取るため、以下のように表記しました。
```math
r_{nk}^{(t)}=p(z_n=k|x_n,\theta^{(t-1)})
```

ここで、ガウス分布の部分を具体的に書き下します。
```math
 \begin{aligned} \log \mathcal{N}(x_n \mid \mu_k, \Sigma_k) &= \log \left[ \frac{1}{(2\pi)^{D/2} |\Sigma_k|^{1/2}} \exp\!\left( -\frac{1}{2}(x_n - \mu_k)^T \Sigma_k^{-1} (x_n - \mu_k) \right) \right] \\ &= -\frac{D}{2}\log(2\pi) - \frac{1}{2}\log |\Sigma_k| -\frac{1}{2}(x_n - \mu_k)^T \Sigma_k^{-1} (x_n - \mu_k). \end{aligned}
```
これを代入すると以下の式が得られます。
```math
Q(\theta; \theta^{(t)}) = \sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk}^{(t)} \Biggl[ \log \pi_k -\frac{D}{2}\log(2\pi) - \frac{1}{2}\log |\Sigma_k| -\frac{1}{2}(x_n - \mu_k)^T \Sigma_k^{-1} (x_n - \mu_k) \Biggr]
```
### 負担率
負担率$r_{nk}^{(t)}$ですが、ベイズの定理を利用すると以下のような関係式が得られます。
```math
r_{nk}^{(t)}=p(z_n=k|x_n,\theta^{(t-1)})=\frac{p(z_n=k)p(x_n|z_n=k,\theta^{(t-1)})}{\sum_{k=1}^{K}p(z_n=k)p(x_n|z_n=k,\theta^{(t-1)})}
=\frac{\pi_k\mathcal{N}(x_n \mid \mu_k^{(t-1)}, \Sigma_k^{(t-1)})}{\sum_{k=1}^K \pi_k \mathcal{N}(x_n \mid \mu_k^{(t-1)}, \Sigma_k^{(t-1)})}
```
負担率$r_{nk}^{(t)}$ですが、$t$回目に使用する右辺の変数はすべて既知のため、計算によって値が得られます。
# 比較
どちらもQ関数と負担率という情報を持つためそれぞれを比較します。
ここで、Q関数ですが、$x_n$を含む項はデータに合わせようとする部分、$x_n$を明示的に含まない項は、パラメータを調整することで、確率モデルとしての制約から逸脱しないようにする、正規化の部分と解釈できます。モデルが変わるとパラメータも変わるので、正規化の部分の比較は難しそうです。そこで、$x_n$を含む項のみに着目して比較を試みます。

### Q関数
K-meansのQ関数を再掲します。
```math
Q(\theta; \theta^{(t)}) = -\sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk}^{(t)} \| x_n - \mu_k \|^2
```
次に、GMMのQ関数のうち$x_n$を含む項を再喝します。
```math
Q(\theta; \theta^{(t)}) = \sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk}^{(t)} \Biggl[ -\frac{1}{2}(x_n - \mu_k)^T \Sigma_k^{-1} (x_n - \mu_k) \Biggr]
```
２つの式を比較します。
まず、2乗誤差の計算部分に分散共分散行列の逆行列$\Sigma_k^{-1}$が入っている点に違いがありますね。
そのため、K-meansは$k$番目のクラスタ$\mu_k$から等距離円上に観測データが分布する様子を表現することに対し、GMMは楕円上に観測データが分布する様子を表現します。

以上より、GMMはK-meansよりも表現の自由度が高い、より一般的なモデルと確認できました。

### 負担率
K-meansの負担率を再掲します。
```math
r_{nk}^{(t)}=\delta_{z_n^{(t)},k},z_{n}^{(t)} = \arg \min_{k'} ||x_n-\mu_{k'}^{(t-1)}||^2 
```
GMMの負担率と対応を取るためこの式を見直します。
GMMの負担率に対し共分散が等方でかつクラスタ間で等しいと考えます。その時、以下の式が得られます。
```math
r_{nk}^{(t)}=\frac{\pi_k\mathcal{N}(x_n \mid \mu_k^{(t-1)}, \sigma_{t-1}^2\mathbf{I})}{\sum_{k=1}^K \pi_k \mathcal{N}(x_n \mid \mu_k^{(t-1)}, \sigma_{t-1}^2\mathbf{I})}
```
この時、ガウス分布は以下のようになります。
```math
\mathcal{N}(x_n \mid \mu_k^{(t-1)}, \sigma_{t-1}^2\mathbf{I})=\frac{1}{(2\pi\sigma_{t-1}^2)^{D/2}}\exp\Biggl(-\frac{1}{2\sigma_{t-1}^2}||x_n - \mu_k^{(t-1)}||^2\Biggr)
```
$\sigma^2\rightarrow0$の極限を考えると、この場合、指数関数は非常に鋭くピークを持つため、分子は$x_n$に対し$||x_n-\mu_{k}^{(t-1)}||^2 $が最も小さい$k=k^{\*}$のみ支配的な値をとります。この時、潜在変数を用いて$z_n=k^{\*}$とします。$z_n=k^*$の時、分母により規格化されるので$r_{nk}^{(t)}=1$となります。

以上より、これまでと同等の式が成り立ちます。
```math
r_{nk} =
\begin{cases}
1, & \text{if } z_n = \arg\min_{j} \|x_n - \mu_j\|^2, \\
0, & \text{otherwise,}
\end{cases}
```
次に、GMMの負担率を再掲します。
```math
r_{nk}^{(t)}=\frac{\pi_k\mathcal{N}(x_n \mid \mu_k^{(t-1)}, \Sigma_k^{(t-1)})}{\sum_{k=1}^K \pi_k \mathcal{N}(x_n \mid \mu_k^{(t-1)}, \Sigma_k^{(t-1)})}
```
上記のことから、K-meansの負担率はGMM
の負担率の極限に相当し、潜在変数$z_n^{(t)}\in \\{1,..,K\\}$は$k$個のクラスタのいずれか１つから得られると$r_{nk}^{(t)}=1$になることを意味します。つまり、K-meansでは$k$番目のクラスタから生成されたか否かを明確に決定します。その意味で、ハードな割り当て方でクラスタ推定しています。

これに対し、GMMの負担率ですが他のクラスタから生成される確率を考慮しています。その意味で、ソフトな割り当て方でクラスタ推定しているといえます。

以上より、GMMはK-meansよりも、不確かさも考慮している点で、より一般的なモデルと確認できました。


# おわりに
最後まで読んでいただきありがとうございました。

潜在変数が離散値を取るモデルはワンホット表現を利用することが多いと思います。しかしこれを使うと、潜在変数が冪指数として登場する形となり、連続値を取るモデルとの対応が取りにくくなります。そのため、本記事ではカテゴリ値を取るように表現を工夫しました。またQ関数による比較のため、観測データに合わせようとする部分に着目するよう工夫しました。

Q関数で比較する記事を見たことがなく、試行錯誤の中の執筆になりました。少なくとも自分自身で納得いくレベルには仕上がったと感じてます。実装があるとより視覚的に読みやすくなると考えてますが、そちらは今後の課題とします。

知識至らずで曖昧な点があると思います。ご指摘やご質問などがございましたら、編集や回答いたします。コメントお待ちしております。