【Python】非負値行列因子分解（NMF）における推論：ギブスサンプリング

・はじめに

「ベイズ推論による機械学習入門」の学習ノートです。この記事は5.2節の非負値行列因子分解の内容です。「観測モデルをデルタ分布」、「事前分布をガンマ分布」とする非負値行列因子分解に対するギブスサンプリングをPythonで実装します。

【数式読解編】

非負値行列因子分解(NMF)のアルゴリズム導出：ギブスサンプリング ベイズ推論による機械学習入門」の学習ノートです。この記事は5.2節の非負値行列因子分解の内容です。「観測モデルをデルタ分布」、「事前分布をガンマ分布」とする非負値行列因子分解の事後分布をギブスサンプリングを用いて推論します。

初学者による学習ノートであるため、解釈の違い等あればコメント頂ければ幸いです。

またコーディングについても初学者であるため、分かりにくい点が多いと思いますがご了承ください。

目次

• ・はじめに
• 【数式読解編】
• ・Pythonでやってみる
  • ・モデルの設定
  • ・事前分布、初期値の設定
  • ・ギブスサンプリング
    • ・コード全体
    • ・処理の確認
  • ・推論結果の確認
    • ・最後のパラメータの確認
    • ・潜在変数のパラメータの確認
• ・さいごに
• ・参考文献

・Pythonでやってみる

　人工的に生成したデータを用いて、ベイズ推論を行う。

　利用するライブラリを読み込む。

# 使用するライブラリ
import numpy as np
from scipy.stats import poisson, gamma, multinomial #　ポアソン分布、ガンマ分布、多項分布
import matplotlib.pyplot as plt

・モデルの設定

　まずは観測モデルを設定する。

ここでデータ$\bf{X}$は$N\times M$のデータであり、$N\times K$の行列$\bf{W}$と$K\times M$の行列$\bf{H}$に近似分解するモデルを考える。

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

この例では観測モデルをデルタ分布$Del(X_{n,m}|\sum_{k=1}^{K}S_{n,k,m})$とする。

また潜在変数$\bf{S}$は$Poi(S_{n,k,m}|W_{n,k}H_{k,m})$から生成されるとする。

そのため、初めに$\bf{S}$を生成するための$\bf{W,H}$をパラメータとして生成する。
$\bf{W,H}$の超パラメータを設定し、ガンマ分布に従う乱数として生成する。

# W, H の超パラメータを設定
a_w, b_w = 1, 1
a_h, b_h = 1, 1

# データのサイズを指定
N = 10
M = 10
K = 5

# W,Hを生成
np.random.seed(seed=0) # 乱数シードを指定
truth_W = np.random.gamma(shape=a_w, scale=1 / b_w, size = (N,K))
truth_H = np.random.gamma(shape=a_h, scale=1 / b_h, size = (K,M))

この$\bf{W,H}$から、ポアソン分布に従う$\bf{S}$を生成し、$\bf{X}$を生成する。

# ポアソン分布に従う乱数からSを生成
truth_S = np.zeros((N,K,M))
np.random.seed(seed=1)
for n in range(N):
    for k in range(K):
        for m in range(M):
            truth_S[n][k][m] = np.random.poisson(lam=truth_W[n][k]*truth_H[k][m])

# Sのkについて和を取りXとする
X = np.sum(truth_S, axis = 1)

　観測モデルを作図する。

# 観測モデルの作図
fig, ax = plt.subplots()
ax.tick_params(which='both', direction='out')
ax.set_title('X')
ax.set_xlabel('n')
ax.set_ylabel('m')

im = ax.imshow(X)
fig.colorbar(im, label='colorbar')

・事前分布、初期値の設定

　次に潜在変数のモデル$P(\bf{S}|\bf{W,H})$の共役事前分布を設定する。ポアソン分布のパラメータ$\bf{W,H}$に対するガンマ分布事前分布$P(\bf{W})$、$P(\bf{H})$として次のように設定する。

\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{3}
\end{eqnarray}

\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{4}
\end{eqnarray}

$\bf{W,H}$の事前分布のパラメータ(超パラメータ)を設定し、事前分布に従う$\bf{W,H}$を生成して初期値とする。

#  W,Hの事前分布のパラメータを指定
a_w, b_w = 5,5
a_h, b_h = 5,5

# W,Hを生成
W_nk = np.random.gamma(shape = a_w, scale=1 / b_w, size = (n,k))
H_km = np.random.gamma(shape = a_h, scale=1 / b_h, size = (k,m))

$\bf{W}$をW_nk、$\bf{H}$をH_kmとして生成した。これらの初期値の積をsyoki_Xとして表し、観測データの時と同様に作図する。

# shoki = WH として計算
syoki_X = np.dot(W_nk,H_km)

・ギブスサンプリング

　ギブスサンプリングによりパラメータと潜在変数を推定する。

・コード全体

pi_nkm,S_nkmについては添字を使って代入するため、受け皿を作成した。

import tqdm

# 試行回数を指定
MaxIter = 10000

# 受け皿を作成
pi_nkm = np.zeros((N,K,M))
S_nkm = np.zeros((N,K,M))

# ギブスサンプリング
for i in tqdm.tqdm(range(MaxIter)):
    
    # pi を計算するためのパラメータ
    ln_W_nk = np.log(W_nk) 
    ln_H_km = np.log(H_km)
    
    # S の事後分布のパラメータ pi を計算:式(6)
    for n in range(N):
        for k in range(K):
            for m in range(M):
                pi_nkm[n][k][m] = np.exp(ln_H_km[k][m]+ln_W_nk[n][k]) 
                
    # pi をkについて規格化
    pi_nkm /= np.sum(pi_nkm, axis=1, keepdims=True)
        
    # 規格化した pi を用いて S をサンプル:式(5)
    np.random.seed(seed=i+100)
    for n in range(N):
        for k in range(K):
            for m in range(M):
                S_nkm[n,:,m] = np.random.multinomial(n=X[n][m], pvals=pi_nkm[n,:,m])
    
    #  W の事後分布のパラメータを計算
    a_w_hat_nk = np.sum(S_nkm, axis=2) + a_w
    b_w_hat_k = np.sum(H_km, axis=1) + b_w

    # W をサンプル
    np.random.seed(seed=i+200)
    W_nk = np.random.gamma(shape=a_w_hat_nk, scale=1 / b_w_hat_k)

    #  H の事後分布のパラメータを計算
    a_h_hat_km = np.sum(S_nkm, axis=0) + a_h
    b_h_hat_k = np.sum(W_nk, axis=0) + b_h
    
    # W をサンプル
    np.random.seed(seed=i+300)
    H_km = (np.random.gamma(shape=a_h_hat_km.T, scale=1 / b_h_hat_k)).T

# W と H の積を計算
NMF_WH = np.dot(W_nk,H_km)

・処理の確認

続いて、ギブスサンプリングで行う処理の確認を行う。

・潜在変数$\bf{S}$のサンプリング
　$\bf{S}$の条件付き事後分布のパラメータを計算し、新たな$\bf{S}$を生成する。

# pi を計算するためのパラメータ
    ln_W_nk = np.log(W_nk) 
    ln_H_km = np.log(H_km)
    
    # S の事後分布のパラメータ pi を計算:式(6)
    for n in range(N):
        for k in range(K):
            for m in range(M):
                pi_nkm[n][k][m] = np.exp(ln_H_km[k][m]+ln_W_nk[n][k]) 
                
    # pi をkについて規格化
    pi_nkm /= np.sum(pi_nkm, axis=1, keepdims=True)
        
    # 規格化した pi を用いて S をサンプル:式(5)
    np.random.seed(seed=i+100)
    for n in range(N):
        for k in range(K):
            for m in range(M):
                S_nkm[n,:,m] = np.random.multinomial(n=X[n][m], pvals=pi_nkm[n,:,m])

　$\bf{S}$の条件付き事後分布は次のように書くことができる。

\begin{equation}
  \mathbf{S_{n,:,m}} \sim Mult(\mathbf{S_{n,:,m}}|{\mathbf{\hat{\pi}_{n,m}}},X_{n,m}) \tag{5}
\end{equation}

ただし、

\begin{eqnarray}
  \hat{\pi}^{(k)}_{n,m} \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{6}
\end{eqnarray}

　n,mを指定した際のパラメータ${\bf{\pi}}^k = (\pi_{n,1,m},\cdots,\pi_{n,K,m})$を計算し、kについて規格化する。その結果をpi_nkmとする。

　更新したパラメータpi_nkmを持ち、試行回数を$X_{n,m}$とする多項分布に従い、$\bf{S}$をサンプルする。更新した$\bf{S}$を用いて他の変数を更新する。

・潜在変数のパラメータ$\bf{W,H}$のサンプリング
　$\bf{W}$の条件付き事後分布のパラメータを計算し、新たな$\bf{W}$を生成する。

    # W の事後分布のパラメータを計算:式(8),(9)
    a_w_hat_nk = np.sum(S_nkm, axis=2) + a_w 
    b_w_hat_k = np.sum(H_km, axis=1) + b_w 

    # W をサンプル:式(7)
    np.random.seed(seed=i+200)
    W_nk = np.random.gamma(shape=a_w_hat_nk, scale=1 / b_w_hat_k)

$\bf{W}$の条件付き事後分布は次のように書くことができる。

\begin{equation}
  W_{n,k} \sim Gam(W_{n,k}|\hat{a}^{(n,k)}_W,\hat{b}^{(k)}_W)
\end{equation} \tag{7}

ただし、

\begin{eqnarray}
  {\hat{a}^{(n,k)}_W} &=& \sum_{m=1}^{M}S_{n,k,m} + a_W \tag{8} \\
  {\hat{b}^{(k)}_W} &=& \sum_{m=1}^{M}H_{k,m} + b_W \tag{9}
\end{eqnarray}

$\bf{W}$の条件付き事後分布のパラメータを(8),(9)式で計算し、a_w_hat_nk,b_w_hat_kとする。次に$\bf{W}$のサンプルを(7)式で行いW_nkとして更新した。

　最後に$\bf{H}$の条件付き事後分布のパラメータを計算し、新たな$\bf{H}$を生成する。

    # H の事後分布のパラメータを計算:(11)(12)式
    a_h_hat_km = np.sum(S_nkm, axis=0) + a_h
    b_h_hat_k = np.sum(W_nk, axis=0) + b_h
    
    # W をサンプル:(10)式
    np.random.seed(seed=i+300)
    H_km = (np.random.gamma(shape=a_h_hat_km.T, scale=1 / b_h_hat_k)).T

$\bf{H}$の条件付き事後分布は次のように書くことができる。

\begin{equation}
  H_{k,m} \sim Gam(H_{k,m}|\hat{a}^{(k,m)}_H,\hat{b}^{(k)}_H) \tag{10}
\end{equation} 

ただし、

\begin{eqnarray}
  \hat{a}^{(k,m)}_H &=& \sum_{n=1}^{N}S_{n,k,m} + a_H \tag{11} \\
  \hat{b}^{(k)}_H &=& \sum_{n=1}^{N}W_{n,k} + b_H \tag{12}
\end{eqnarray}

$\bf{H}$の条件付き事後分布のパラメータを(11),(12)式で計算し、a_h_hat_km,b_h_hat_kとする。次に$\bf{H}$のサンプルを(10)式で行いH_kmとして更新した。

　$\bf{W,H}$を更新することができたので、今度は$\bf{W,H}$のサンプルから潜在変数$\bf{S}$を更新する。この処理を指定した回数繰り返す。

・推論結果の確認

・最後のパラメータの確認

W_nk,H_kmの最後の更新値から、元のデータ$\bf{X}$をどれだけ表現できているのかを見てみる。

# W と H の積を計算
NMF_WH = np.dot(W_nk,H_km)

元のデータと同様に可視化する。

元のデータと比較する。

おおよその特徴は掴めていることがわかる。

・潜在変数のパラメータの確認

$\bf{W,H}$がどれほど生成したものを推定できたか確認する。
最後の更新値W_nk,H_kmと、データ$\bf{X}$を生成するために用意したtruth_W,truth_Hを比較する。

まずW_nkを可視化する。

次にtruth_Wを可視化する。

同様にH_km,truth_Hの比較も行う。

適宜推論の結果は更新していきます。

・さいごに

　本記事は初学者による『ベイズ推論による機械学習入門』の学習時ノートです。前回の記事でギブスサンプリングによるアルゴリズム導出を行ったため、その実装例となります。コーディングについても初学者のため、解釈の違い等ございましたら、是非コメント頂ければと思います。

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

　これからも勉強した内容をあげていければと思いますので、よろしくお願いいたします。

・参考文献

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