非負値行列因子分解(NMF)のアルゴリズム導出：ギブスサンプリング

## 目次
- [はじめに](#はじめに)
- [モデル](#モデル)
- [ギブスサンプリング](#ギブスサンプリング)
    - [潜在変数の条件付き分布の導出](#潜在変数の条件付き分布の導出)
    - [パラメータの条件付き分布の導出](#パラメータの条件付き分布の導出)
- [さいごに](#さいごに)
- [参考文献](#参考文献)

## はじめに
「ベイズ推論による機械学習入門」の学習ノートです。この記事は5.2節の非負値行列因子分解の内容です。「観測モデルをデルタ分布」、「事前分布をガンマ分布」とする非負値行列因子分解の事後分布をギブスサンプリングを用いて推論します。初学者による学習ノートであるため、解釈の違い等あればコメント頂ければ幸いです。

## モデル
各要素が非負値を持つ行列$\bf{X}\in\mathbb{N}^{N\times M}$を性の値を持つ二つの行列$\bf{W}\in\mathbb{N}^{N\times K}$と$\bf{H}\in\mathbb{N}^{K\times M}$に近似分解する。

``` math
\mathbf{X} \approx \mathbf{WH} \tag{1}
```

ここで$\bf{X}$を要素ごとに書き下すと次式のようになる。

``` math
\begin{eqnarray}
  X_{n,m} &=&\sum_{k=1}^{K}S_{n,k,m} \tag{2}\\
  &\approx& \sum_{k=1}^{K}W_{n,k}H_{k,m} \tag{3}
\end{eqnarray}
```

$S_{n,k,m}\approx \sum_{k=1}^{K}W_{n,k}H_{k,m}$となるような潜在変数$\bf{S}$は行列分解において必要ないが、導入することで効率的なアルゴリズムの導出が可能。

また潜在変数$\bf{S}$は次のようなポアソン分布に従って発生するようにモデル化する。

``` math
\begin{eqnarray}
  P(\mathbf{S}|\mathbf{W},\mathbf{H})
  &=& \prod_{n=1}^{N}\prod_{k=1}^{K}\prod_{m=1}^{M}
  P(S_{n,k,m}|W_{n,k},H_{k,m}) \\
  &=& \prod_{n=1}^{N}\prod_{k=1}^{K}\prod_{m=1}^{M} Poi(S_{n,k,m}|W_{n,k}H_{k,m}) \tag{4}
\end{eqnarray}
```
観測データ$\bf{X}$はデルタ分布から生成されると仮定する。

``` math
\begin{eqnarray}
  P(\mathbf{X}|\mathbf{S})
  &=& \prod_{n=1}^{N}\prod_{m=1}^{M}
  P(X_{n,m}|S_{n,:,m}) \\
  &=& \prod_{n=1}^{N}\prod_{m=1}^{M}
  Del(X_{n,m}|\sum_{k=1}^{K}S_{n,k,m}) \tag{5}
\end{eqnarray}
```
生成的な観点から言えば(5)式は、$\sum_{k=1}^{K}S_{n,k,m}$から得られた値がそのまま$X_{n,m}$の値として扱われるということ。

また$\bf{W,H}$の事前分布にはポアソン分布の共役事前分布であるガンマ分布を用いる。

``` math
\begin{eqnarray}
  P(\mathbf{W})
  &=& \prod_{n=1}^{N}\prod_{k=1}^{K}
  P(W_{n,k}) \\
  &=& \prod_{n=1}^{N}\prod_{k=1}^{K}
  Gam(W_{n,k}|a_W,b_W) \tag{6}
\end{eqnarray}
```

``` math
\begin{eqnarray}
  P(\mathbf{H})
  & = & \prod_{k=1}^{K}\prod_{m=1}^{M}
  P(H_{k,m}) \\
  & = & \prod_{k=1}^{K}\prod_{m=1}^{M}
  Gam(H_{k,m}|a_H,b_H) \tag{7}
\end{eqnarray}
```

改めてこれらの確率変数の関係性から同時分布を書き出すと次式のようになる。

``` math
\begin{eqnarray}
  P(\mathbf{X},\mathbf{S},\mathbf{W},\mathbf{H})
  & = & P(\mathbf{X}|\mathbf{S},\mathbf{W},\mathbf{H})P(\mathbf{S},\mathbf{W},\mathbf{H})  \\
  & = & P(\mathbf{X}|\mathbf{S})P(\mathbf{S}|\mathbf{W},\mathbf{H})P(\mathbf{W},\mathbf{H}) \\
  & = & P(\mathbf{X}|\mathbf{S})P(\mathbf{S}|\mathbf{W},\mathbf{H})P(\mathbf{W})P(\mathbf{H})\tag{8}
\end{eqnarray}
```
1行目から2行目の式変形において$\bf{X}$は$\bf{S}$から生成されるため、$\bf{W,H}$を除き、$\bf{S,W,H}$の同時分布をベイズの定理から分解した。2行目から3行目の式変形において$\bf{W,H}$は独立であるため分解した。


## ギブスサンプリング
先程構築したモデルから観測されていない変数$\bf{W,H}$及び$\bf{S}$の事後分布を推論するためのアルゴリズムを導出する。ここではギブスサンプリングを用いた推論を行う。

### 潜在変数の条件付き分布の導出
初めに潜在変数$\bf{S}$をサンプルするための条件付き分布($\bf{S}$の事後分布)$P(\mathbf{S}|\mathbf{X},\mathbf{W},\mathbf{H})$を求めていく。

観測データ$\bf{X}$に加えて$\bf{W,H}$が観測されたとした時、$\bf{S}$の事後分布は次のように書くことができる。

``` math
\begin{eqnarray}
    P(\mathbf{S}|\mathbf{X},\mathbf{W},\mathbf{H})
&=& \frac{P(\mathbf{S},\mathbf{X},\mathbf{W},\mathbf{H})}{P(\mathbf{X},\mathbf{W},\mathbf{H})} \\
&\propto& P(\mathbf{S},\mathbf{X},\mathbf{W},\mathbf{H}) \\
&=& P(\mathbf{X}|\mathbf{S})P(\mathbf{S}|\mathbf{W},\mathbf{H})P(\mathbf{W})P(\mathbf{H}) \\
&\propto& P(\mathbf{X}|\mathbf{S})P(\mathbf{S}|\mathbf{W},\mathbf{H})\tag{9}
\end{eqnarray}
```
モデル構築の際に確認した生成過程(各変数の依存関係)に従い、各項を分解した。また、適宜$\bf{S}$に関係ない項を省略し、比例関係に注目する。省略した項については最後に正規化することで対応できる。

(4),(5),(9)式から$\bf{S}$の事後分布は次のように書くことができる。

``` math
\begin{eqnarray}
  P(\mathbf{S}|\mathbf{X,W,H})\propto \prod_{n=1}^{N}\prod_{m=1}^{M}
  \biggl\{Del(X_{n,m}|\sum_{k=1}^{K}S_{n,k,m})\cdot\prod_{k=1}^{K}  Poi(S_{n,k,m}|W_{n,k}H_{k,m})\biggr\} \tag{10}
\end{eqnarray}
```

(10)式から$\bf{S}$の事後分布は対数を取るとn,mに関する和で分解でき,次のように書くことができる。

``` math
\begin{eqnarray}
  \ln P(\mathbf{S}|\mathbf{X,W,H}) = \sum_{n=1}^{N}\sum_{m=1}^{M}
  \biggl\{ \ln Del(X_{n,m}|\sum_{k=1}^{K}S_{n,k,m}) + \sum_{k=1}^{K}  
 \ln Poi(S_{n,k,m}|W_{n,k}H_{k,m})\biggr\}  \\ +const. \tag{11}
\end{eqnarray}
```

具体的に(11)式の二項目を計算すると次のようになる。
``` math
\begin{eqnarray}
  \ln Poi(S_{n,k,m}|W_{n,k}H_{k,m})
  &=&
  \ln \biggl\{\frac{{W_{n,k}H_{k,m}}^{S_{n,k,m}}}{S_{n,k,m}!}\exp(-W_{n,k}H_{k,m})\biggr\} \\
  &=& S_{n,k,m}(\ln W_{n,k}+ \ln H_{k,m}) - \ln S_{n,k,m}! + const. \tag{12}
\end{eqnarray}
```
適宜$\bf{S}$に関係ない項はconst.にまとめた。
また、(11)式の一項目は言い換えると$X_{n,m} = \sum_{k=1}^{K}S_{n,k,m}$なので$\bf{S}$の事後分布は試行回数を$X_{n,m}$とするような多項分布として表現できる。以上より$\bf{S}$は次のようにサンプルすることができる。

``` math
\begin{equation}
  \mathbf{S_{n,:,m}} \sim Mult(\mathbf{S_{n,:,m}}|{\mathbf{\hat{\pi}_{n,m}}},X_{n,m}) \tag{13}
\end{equation}
```
ただし、

``` math
\begin{eqnarray}
  {\hat{\pi}_{n,m}}^{(k)} \propto \exp(\ln W_{n,k}+ \ln H_{k,m})  \\\\
\biggl(s.t \sum_{k=1}^{K}{\hat{\pi}_{n,m}}^{(k)} = 1\biggr) \tag{14}
\end{eqnarray}
```

### パラメータの条件付き分布の導出

次に$\bf{W}$をサンプルするための条件付き分布$P(\bf{W}|\bf{X},\bf{S},\bf{H})$を求める。

$\bf{S}$の時と同様に$\bf{X,S,H}$は観測データと見なして事後分布を計算する。ただし、サンプルとして得られている変数はそれを用いることができる。

``` math
\begin{eqnarray}
    P(\mathbf{W}|\mathbf{X,S,H}) &=& \frac{P(\mathbf{X,S,W,H})}{P(\mathbf{X,S,H})} \\
&\propto& P(\mathbf{X,S,W,H}) \\
&=& P(\mathbf{X}|\mathbf{S,W,H})P(\mathbf{S}|\mathbf{W,H}) \\
&=& P(\mathbf{X}|\mathbf{S})P(\mathbf{S}|\mathbf{W,H})P(\mathbf{W})P(\mathbf{H}) \\
&\propto& P(\mathbf{S}|\mathbf{W,H})P(\mathbf{W}) \tag{15}
\end{eqnarray}
```
$\bf{S}$と同様、生成過程に従い各項を分解し、$\bf{W}$に関係ない項を省略した。
(4),(6),(14)式から$\bf{W}$の事後分布は次のように計算することができる。

``` math
\begin{eqnarray}
     P(\mathbf{W}|\mathbf{X,S,H}) \propto \prod_{n=1}^{N}\prod_{k=1}^{K}
  \biggl\{\prod_{m=1}^{M} Poi(S_{n,k,m}|W_{n,k}H_{k,m})
  \cdot Gam(W_{n,k}|a_W,b_W)\biggr\} \tag{16}
\end{eqnarray}
```
対数を取るとn,kの和について分解でき、次のように表現することができる。

``` math
\begin{eqnarray}
  \ln P(\mathbf{W}|\mathbf{X,S,H}) =
  \sum_{n=1}^{N}\sum_{k=1}^{K}\biggl\{
  \sum_{m=1}^{M}\ln Poi(S_{n,k,m}|W_{n,k}H_{k,m}) +
  \ln Gam(W_{n,k}|a_W,b_W) \biggr\} \tag{17}
\end{eqnarray}
```

(17)式一項目、二項目について$\bf{W}$に関係ない項をconst.にまとめると次のようになる。

``` math
\begin{eqnarray}
  \ln Poi(S_{n,k,m}|W_{n,k}H_{k,m}) &=&
  \ln\biggl\{\frac{{W_{n,k}H_{k,m}}^{S_{n,k,m}}}{S_{n,k,m}!}\exp(-W_{n,k}H_{k,m})\biggr\} \\
  &=& S_{n,k,m}\ln W_{n,k} -H_{k,m}W_{n,k} +const. \tag{18}
\end{eqnarray}
```

``` math
\begin{eqnarray}
  \ln Gam(W_{n,k}|a_W,b_W) &=&
  \ln \biggl\{\frac{b_W^{a_W}}{\Gamma(a_W)}W_{n,k}^{(a_W-1)}
  \exp(-W_{n,k}b_W)\biggr\} \\
  &=& (a_W-1) \ln W_{n,k} -b_W W_{n,k} + const. \tag{19}
\end{eqnarray}
```

よって(17)式は(19),(18)式がmについての和を取ることを考慮して次のように書くことができる。

``` math
\begin{eqnarray}
  \ln P(\mathbf{W}|\mathbf{X,S,H}) =
  \sum_{n=1}^{N}\sum_{k=1}^{K}\biggl\{
  \biggl(\sum_{m=1}^{M}S_{n,k,m}+a_W-1\biggr)\ln W_{n,k}-
  \biggl(\sum_{m=1}^{M}H_{k,m}+b_W\biggr)W_{n,k}\biggr\}  \\ +const. \tag{20}
\end{eqnarray}
```


(20)式から$\bf{W}$の事後分布はガンマ分布で表現することができる。以上より$\bf{W}$は次のようにサンプルすることができる。

``` math
\begin{equation}
  W_{n,k} \sim Gam(W_{n,k}|{\hat{a}^{(n,k)}_W},{\hat{b}^{(k)}_W})
\end{equation} \tag{21}
```
ただし、
``` math
\begin{eqnarray}
  {\hat{a}^{(n,k)}_W} &=& \sum_{m=1}^{M}S_{n,k,m} + a_W \tag{22} \\
  {\hat{b}^{(k)}_W} &=& \sum_{m=1}^{M}H_{k,m} + b_W \tag{23}
\end{eqnarray}
```

最後に$\bf{H}$をサンプルするための事後分布を計算する。ここで$P(\bf{H}|\bf{X,S,W})$を計算すると次のようになる。


``` math
\begin{eqnarray}
P(\mathbf{H}|\mathbf{X,S,W}) &=& \frac{P(\mathbf{X,S,W,H})}{P(\mathbf{X,S,W})} \\

&\propto& P(\mathbf{X,S,W,H}) \\

&=& P(\mathbf{X}|\mathbf{S,W,H})P(\mathbf{S}|\mathbf{W,H}) \\

&=& P(\mathbf{X}|\mathbf{S})P(\mathbf{S}|\mathbf{W,H})P(\mathbf{W})P(\mathbf{H}) \\

&\propto& P(\mathbf{S}|\mathbf{W,H})P(\mathbf{H}) \tag{24}
\end{eqnarray}
```
先程同様適宜$\bf{H}$に関係ない項は省略し比例関係に注目した。よって(4),(7),(24)式より事後分布を書くと次のようになる。

``` math
\begin{eqnarray}
  P(\mathbf{H}|\mathbf{X,S,W})\propto \prod_{k=1}^{K}\prod_{m=1}^{M}
  \biggl\{\prod_{n=1}^{N} Poi(S_{n,k,m}|W_{n,k}H_{k,m})
  \cdot Gam(H_{k,m}|a_H,b_H)\biggr\}　\tag{25}
\end{eqnarray}
```
ここで(25)式と$\bf{W}$の事後分布の式(16)式の間には対称性があるため、対数を取り具体的に計算したものは次のようになる。

``` math
\begin{eqnarray}
  \ln P(\mathbf{H}|\mathbf{X,S,W}) =
  \sum_{k=1}^{K}\sum_{m=1}^{M}\biggl\{
  \biggl(\sum_{n=1}^{N}S_{n,k,m}+a_H-1\biggr)\ln H_{k,m}-
  \biggl(\sum_{n=1}^{N}W_{n,k}+b_H\biggr)H_{k,m}\biggr\} 　\\ +const. \tag{26}
\end{eqnarray}
```

(26)式より、$\bf{H}$の事後分布もガンマ分布で表現することができ、次のようにサンプルすることができる。
``` math
\begin{equation}
  H_{k,m} \sim Gam(H_{k,m}|\hat{a}^{(k,m)}_H,\hat{b}^{(k)}_H)
\end{equation}
```

ただし、
``` math
\begin{eqnarray}
  \hat{a}^{(k,m)}_H &=& \sum_{n=1}^{N}S_{n,k,m} + a_H \\
  \hat{b}^{(k)}_H &=& \sum_{n=1}^{N}W_{n,k} + b_H
\end{eqnarray}

```

これでパラメータ$\bf{S,W,H}$をサンプルすることができる。


## さいごに
本記事は初学者による『ベイズ推論による機械学習入門』の学習時ノートです。『ベイズ推論による機械学習入門』ではNMFのアルゴリズム導出はベイズ推論による変分推論で行われていたため学習の一環として事後分布をギブスサンプリングによって求めました。解釈の違い等ございましたら、是非コメント頂ければと思います。

また、本記事作成にあたって参考にした書籍、ブログを参考文献に挙げさせて頂きます。

[次は](https://qiita.com/akamaru0108/items/ceb3339992adb4205dc6)実装した例を更新できればと思います。よろしくお願いいたします。

## 参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
- https://www.anarchive-beta.com/about