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

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

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




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

# PPCA
### 生成モデル
PPCA (Probabilistic PCA) では、観測データ $\mathbf{x} \in \mathbb{R}^D$ と潜在変数 $\mathbf{z} \in \mathbb{R}^d$ ($d < D$) の関係を以下のようにモデル化します。
1.潜在変数の生成
```math
p(\mathbf{z}) = \mathcal{N}(\mathbf{z} \mid \mathbf{0}, \mathbf{I}),
```
2.観測データの生成
```math
p(\mathbf{x} \mid \mathbf{z},\theta) = \mathcal{N}(\mathbf{x} \mid \mathbf{W}\mathbf{z} + \boldsymbol{\mu},\, \sigma^2 \mathbf{I})
```
ここで、$\mathbf{W} \in \mathbb{R}^{D \times d}$ は射影行列、$\boldsymbol{\mu} \in \mathbb{R}^D$ は平均、$\sigma^2$ はノイズ分散です。

この時、同時分布は以下のようになります。
```math
p(\mathbf{x}, \mathbf{z}|\theta) = p(\mathbf{z}) \, p(\mathbf{x} \mid \mathbf{z},\theta)
= \mathcal{N}(\mathbf{z} \mid \mathbf{0}, \mathbf{I}) \, \mathcal{N}(\mathbf{x} \mid \mathbf{W}\mathbf{z} + \boldsymbol{\mu},\, \sigma^2 \mathbf{I})
```


### Q関数
各データ点 $\mathbf{x}_n$ ($n=1,\dots,N$) について、潜在変数 $\mathbf{z}_n$ の事後分布 $p(\mathbf{z}_n \mid \mathbf{x}_n, \theta^{(t-1)})$ に対する期待値として、Q 関数は以下のように定義されます。


```math
Q(\theta; \theta^{(t)})  = \sum_{n=1}^{N} \int {p(\mathbf{z}_n \mid \mathbf{x}_n, \theta^{(t-1)})} \log p(\mathbf{x}_n, \mathbf{z}_n \mid \theta) d\mathbf{z}_n=\sum_{n=1}^{N} \int {p(\mathbf{z}_n \mid \mathbf{x}_n, \theta^{(t-1)})} \log \Biggl[\mathcal{N}(\mathbf{z}_n \mid \mathbf{0}, \mathbf{I}) \, \mathcal{N}(\mathbf{x}_n \mid \mathbf{W}\mathbf{z}_n + \boldsymbol{\mu},\, \sigma^2 \mathbf{I})
\Biggr] d\mathbf{z}_n
```
ここで、ガウス分布の部分を具体的に書き下します。
```math
\log\mathcal{N}(\mathbf{z}_n \mid \mathbf{0}, \mathbf{I})=
- \frac{1}{2}\|\mathbf{z}_n\|^2 - \frac{d}{2}\log(2\pi)
```

```math
\log\mathcal{N}(\mathbf{x}_n \mid \mathbf{W}\mathbf{z}_n + \boldsymbol{\mu},\, \sigma^2 \mathbf{I})=
- \frac{1}{2\sigma^2}\| \mathbf{x}_n - \mathbf{W}\mathbf{z}_n- \mu \|^2-\frac{D}{2}\log(2\pi \sigma^2)
```

これを代入すると以下の式が得られます。
```math
Q(\theta; \theta^{(t)})  = \sum_{n=1}^{N} \int r^{(t)}(\mathbf{z}_n) \Biggl[  - \frac{1}{2\sigma^2}\| \mathbf{x}_n - \mathbf{W}\mathbf{z}_n - \mu \|^2 - \frac{1}{2}\|\mathbf{z}_n\|^2 -\frac{D}{2}\log(2\pi \sigma^2) - \frac{d}{2}\log(2\pi) \Biggr] d\mathbf{z}_n
```
### 負担率
負担率$r^{(t)}(\mathbf{z}_n)$ですが、ベイズの定理を利用すると以下のような関係式が得られます。
```math
r^{(t)}(\mathbf{z}_n)=p(\mathbf{z}_n|\mathbf{x}_n,\theta^{(t)})=\frac{p(\mathbf{z}_n)p(x_n|z_n,\theta^{(t)})}{\int p(\mathbf{z}_n)p(\mathbf{x}_n|\mathbf{z}_n,\theta^{(t)})d\mathbf{z}_n}
=\frac{\mathcal{N}(\mathbf{z}_n \mid \mathbf{0}, \mathbf{I})\mathcal{N}(\mathbf{x}_n \mid \mathbf{W}\mathbf{z}_n + \boldsymbol{\mu},\, \sigma^2 \mathbf{I})}{\int \mathcal{N}(\mathbf{z}_n \mid \mathbf{0}, \mathbf{I}) \, \mathcal{N}(\mathbf{x}_n \mid \mathbf{W}\mathbf{z}_n + \boldsymbol{\mu},\, \sigma^2 \mathbf{I})d\mathbf{z}_n}
```
# GTM
GTM (Generative Topographic Mapping) では、低次元潜在空間上の一点から潜在変数$\{\mathbf{z}\}$が、一様な事前分布に従って生成されたと仮定します。
そして、観測データ $\mathbf{x} \in \mathbb{R}^D$ は、潜在変数を入力とする非線形マッピング $y(\mathbf{z}; \mathbf{W})$ を経由し、ガウスノイズが付加されて生成されると考えます。

注意：通常、GTMは潜在空間をＫ個からなるノードに離散化します。しかし本記事では、PPCAと比較するため離散化せずに扱います。
### 生成モデル
1.潜在変数の生成
潜在変数は一様分布から生成されます。
今、潜在空間の各次元の定義域を$[0,v]$として、空間全体の体積を$V_d$とすると以下の式になります。
```math
p(\mathbf{z}) = \mathcal{U}(\mathbf{z}\mid [0,v]^d)=\frac{1}{V_d}
```
2.観測データの生成
```math
p(\mathbf{x} \mid \mathbf{z},\theta) = \mathcal{N}(\mathbf{x} \mid y(\mathbf{z}; \mathbf{W}), \sigma^2 \mathbf{I})
```
この時、同時分布は以下のようになります。
```math
p(\mathbf{x}, \mathbf{z}|\theta) = p(\mathbf{z}) \, p(\mathbf{x} \mid \mathbf{z},\theta)
= \mathcal{U}(\mathbf{z}\mid [0,1]^d)\mathcal{N}(\mathbf{x} \mid y(\mathbf{z}; \mathbf{W}), \sigma^2 \mathbf{I})
```
### Q関数
Q関数は以下のように定義されます。
```math
Q(\theta; \theta^{(t)})  = \sum_{n=1}^{N} \int {p(\mathbf{z}_n \mid \mathbf{x}_n, \theta^{(t-1)})} \log p(\mathbf{x}_n, \mathbf{z}_n \mid \theta) d\mathbf{z}_n=\sum_{n=1}^{N} \int {p(\mathbf{z}_n \mid \mathbf{x}_n, \theta^{(t-1)})} \log \Biggl[\mathcal{U}(\mathbf{z}_n\mid [0,1]^d)\mathcal{N}(\mathbf{x}_n \mid y(\mathbf{z}_n; \mathbf{W}), \sigma^2 \mathbf{I})\Biggr] d\mathbf{z}_n
```
ここで、ガウス分布の部分を具体的に書き下します。
```math
\log \mathcal{N}(\mathbf{x} \mid y(\mathbf{z}; \mathbf{W}), \sigma^2 \mathbf{I})=\frac{D}{2}\log\frac{1}{2\pi\sigma^2} - \frac{1}{2\sigma^2}\| x_n - y(z_n; W) \|^2
```
これを代入すると以下の式が得られます。
```math
Q(\theta; \theta^{(t)})  = \sum_{n=1}^{N} \int r^{(t)}(\mathbf{z}_n) \left[- \frac{1}{2\sigma^2}\| \mathbf{x}_n - y(\mathbf{z}_n; W) \|^2  -\log V_d -\frac{D}{2}\log(2\pi \sigma^2)\right]d\mathbf{z}_n
```
### 負担率
```math
r^{(t)}(\mathbf{z}_n)=p(\mathbf{z}_n|\mathbf{x}_n,\theta^{(t-1)})=\frac{p(\mathbf{z}_n)p(x_n|z_n,\theta^{(t-1)})}{\int p(\mathbf{z}_n)p(\mathbf{x}_n|\mathbf{z}_n,\theta^{(t-1)})d\mathbf{z}_n}
=\frac{\mathcal{U}(\mathbf{z}_n \mid[0,v]^d)\mathcal{N}(\mathbf{x}_n \mid y(\mathbf{z}_n; \mathbf{W}^{(t-1)}), \sigma^2 \mathbf{I})}{\int \mathcal{U}(\mathbf{z}_n \mid[0,v]^d) \mathcal{N}(\mathbf{x} \mid y(\mathbf{z}_n; \mathbf{W}^{(t-1)}), \sigma^2 \mathbf{I})d\mathbf{z}_n}
```

# 比較
どちらもQ関数と負担率という情報を持つためそれぞれを比較します。
ここで、Q関数ですが、$x_n$を含む項はデータに合わせようとする部分、$x_n$を明示的に含まない項は、パラメータを調整することで、確率モデルとしての制約から逸脱しないようにする、正規化の部分と解釈できます。モデルが変わるとパラメータも変わるので、正規化の部分の比較は難しそうです。そこで、$x_n$を含む項のみに着目して比較を試みます。
### Q関数
PPCAのQ関数のうち、$\mathbf{x}_n$を含む項を再掲します。
```math
Q(\theta; \theta^{(t)})  = \sum_{n=1}^{N} \int r^{(t)}(\mathbf{z}_n) \Biggl[  - \frac{1}{2\sigma^2}\| \mathbf{x}_n - \mathbf{W}\mathbf{z}_n - \mu \|^2  \Biggr] d\mathbf{z}_n
```
GTMのQ関数のうち、$\mathbf{x}_n$を含む項を再掲します。
```math
Q(\theta; \theta^{(t)})  = \sum_{n=1}^{N} \int r^{(t)}(\mathbf{z}_n) \left[- \frac{1}{2\sigma^2}\| \mathbf{x}_n - y(\mathbf{z}_n; W) \|^2  \right]d\mathbf{z}_n
```
２つの式を比較します。
まず、観測データ$\mathbf{x}_n$との潜在変数$\mathbf{z}_n$
からの線型写像の誤差をとるか非線形写像の誤差を取るかの違いがありますね。

以上より、GTMはPPCAよりも表現の自由度が高い、より一般的なモデルと言えます。
### 負担率
PPCAの負担率の分子の部分を再掲します。
```math
r^{(t)}(\mathbf{z}_n)
\propto\mathcal{N}(\mathbf{z}_n \mid \mathbf{0}, \mathbf{I})\mathcal{N}(\mathbf{x}_n \mid \mathbf{W}^{(t-1)}\mathbf{z}_n + \boldsymbol{\mu}^{(t-1)},\, \sigma^2 \mathbf{I})
```
次に、GTMの負担率の分子の部分を再掲します。
```math
r^{(t)}(\mathbf{z}_n)
\propto \mathcal{U}(\mathbf{z}_n \mid[0,v]^d) \mathcal{N}(\mathbf{x} \mid y(\mathbf{z}_n; \mathbf{W}^{(t-1)}), \sigma^2 \mathbf{I})
```
さて、PPCAの負担率n分子は$z$に対して2次形式です。そのため、負担率もガウス分布となります。

```math
r^{(t)}(\mathbf{z}_n) = \mathcal{N}\Bigl(\mathbf{z}\,\Bigl|\, M^{-1} \mathbf{W}^T (\mathbf{x} - \boldsymbol{\mu}),\, \sigma^2 M^{-1}\Bigr),
```
ここで、
```math
M = \mathbf{W}^T \mathbf{W} + \sigma^2 \mathbf{I}.
```
を利用しました。

一方、GTMの負担率の分子ですが、$\mathbf{z}_n$と$y(\mathbf{z}_n; \mathbf{W}^{(t-1)})$は非線形な関係です。例えば$\mathbf{z}_n$多項式などが考えられます。この場合は$\mathbf{z}_n$に対して2次形式でないため、負担率はガウス分布のような形にならず、非線形関数の複雑さに応じて複雑な分布になります。

以上より、GTMはPPCAよりも負担率の形状が複雑です。
負担率は、潜在変数の割り当て方の情報ですから、潜在変数の分布を知りたい時には、負担率の解釈のしやすさも指標として大切です。その点いうとGTMはPPCAよりも解釈が難しいと確認できました。

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

潜在変数が連続値を取るモデルとしてPPCAとGTMを説明しました。Q関数としては、PPCAよりもGTMの方がモデルの表現の自由度が高い一方で、負担率の複雑となることを説明しました。負担率の形状は、解釈のしやさに関わる点で大切だ、というメッセージのために、負担率の形を分析するところを工夫しました。

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

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