機械学習・自然言語処理の勉強メモ

学んだことのメモやまとめ

Stan:NMF

はじめに



今回はNMF(Non-negative Matrix Factorization)について勉強する。
日本語では「非負値行列因子分解」と言い、線形次元手法の1つ。
その名前の通り「観測データの非負性」を仮定したモデル。
「非負性」の制約を設けることで、これに適用可能な観測データはより良い表現ができるらしい。
推薦アルゴリズムなどでも使われている。

詳しいことは省略して実装するうえで大事なところを押さえていく。MCMCではなく、最適化により求める方法は「非負値行列因子分解 - NTTコミュニケーション科学基礎研究所」がとても丁寧だった。

尚、今回の記事のいくつかは「ベイズ推論による機械学習入門 」を参考とした。

NMF



まず、観測データ(各要素が非負値の行列)行列\rm{X}\in \mathbb{R}^{D\times N}を以下のように表す。
\rm{X}\approx\rm{WH}
尚、\rm{W}\in \mathbb{R}^{D\times M}\rm{H}\in \mathbb{M}^{M\times N}
Mは観測データに出現しない潜在変数である。
\rm{X}を要素ごとに書くと次のようになる。
X_{d,n}\approx\sum_{m=1}^M W_{d,m}H_{m,n}
例えば、行列\rm{X}がユーザー×アイテムのレーティング値の場合、\rm{W}はユーザー×潜在変数(ユーザーの嗜好)を表し、行列\rm{H}は潜在変数×アイテム(アイテムの特徴)のようなものを表している。
この場合、ユーザーdのアイテムnへの評価は、ユーザーdの嗜好ベクトルとアイテムnの特徴ベクトルの内積で決定される。2つのベクトルの方向が近い(ユーザーdの嗜好とアイテムnの特徴が近い)ほど、レーティング値は大きくなるようになる。

モデル



今回は観測データの要素が非負の整数を仮定したモデルを考える。
非負の整数なので、データ生成にはポアソン分布が使える。
各要素は以下のように表せたので、
X_{d,n}\approx\sum_{m=1}^M W_{d,m}H_{m,n}
これをそのままポアソン分布にあてはめる。
p(X_{d,n})\sim{\rm Poi} (X_{d,n}|\sum_{m=1}^M W_{d,m}H_{m,n})

今回はこれが尤度となる。(計算では対数を取る。)
NMFで目指すのは、\rm{X}\rm{WH}で近似させること。
各要素X_{d,n}sum_{m=1}^M W_{d,m}H_{m,n}が近似できているかどうかは「尤度」で測ることができる。
ポアソン分布の場合、観測値が平均値と一致する値が尤度最大となる値である。この場合、X_{d,n}=sum_{m=1}^M W_{d,m}H_{m,n}である。観測値が平均から離れるほど尤度は小さくなる
簡単であるが、ポアソン分布での尤度を比べてみた。

dpois(1,1,log=TRUE)
[1] -1
dpois(1,5,log=TRUE)
[1] -3.390562

観測値と平均値が近いほど、尤度が大きいことが分かる。
NMFではポアソン分布の尤度を大きくすることが、\rm{X}\rm{WH}させることと同じ意味なのである。

行列全体の尤度は以下のように各要素の積で書ける。
\begin{eqnarray}
L &=& \prod_{d=1}^D\prod_{n=1}^N {\rm Poi} (X_{d,n}|\sum_{m=1}^M W_{d,m}H_{m,n}) \\
\log L  &=& \sum_{d=1}^D\sum_{n=1}^N \log{\rm Poi} (X_{d,n}|\sum_{m=1}^M W_{d,m}H_{m,n})
\end{eqnarray}

W,Hの各要素は、ガンマ分布を取っておくのが良いと思われる。
p(W_{d,m})\sim{\rm Gam}(W_{d,m}|a,b)
p(H_{m,n})\sim{\rm Gam}(H_{d,m}|a,b)

コード



コードは以下の通り。
(以前、ネットに落ちていたのを少し改良した。)

data{
  int<lower=1> D; // row dimension of input
  int<lower=1> N; // column dimension of input
  int<lower=1> M; // hidden dimension
  
  int X[D,N]; // original data
  
  real<lower=0> alpha; //hyperparameter
  real<lower=0> beta; //hyperparameter
}

parameters{
  matrix<lower=0>[D, M] W;
  matrix<lower=0>[M, N] H;
}


model{
  //estimate W matrix dist
  for(d in 1:D){
    for(m in 1:M){
      W[d,m] ~ gamma(alpha,beta);
    }
  }
  //estimate H matrix dist
  for(m in 1:M){
    for(n in 1:N){
      H[m,n] ~ gamma(alpha,beta);
    }
  }

  for(d in 1:D){
      for(n in 1:N){
        target += poisson_lpmf(X[d,n] | sum(W[d,:])*sum(H[:,n]));
      }
  }
}

以下のようにしてキックする。

stan.dat <- list(D=nrow(data),D=ncol(data),M=2,X=data,alpha=0.1,beta=0.1)
stan.fit <- stan("NMF.stan", data=stan.dat)

ちなみにデータはこんな感じ。

data
  item1 item2 item3 item4 item5 item6 item7 item8 item9 item10
1     0     2     4     1     0     0     0     1     2      5
2     0     2     3     0     0     0     1     2     2      3
3     1     0     1     4     3     0     1     1     1      0
4     0     0     1     4     5     1     0     1     2      1
5     1     1     5     4     1     1     1     0     1      5

収束はするがデータが小さいので面白い発見は出来ない・・

今回は観測データが非負の整数を取ったが、取る範囲が広いほど推論が難しい気がする。観測データが二値でベルヌーイ分布でやってみるのも面白いかもしれない。