廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

変分ベイズによるトピックモデル(GaP; Gamma-Poisson Model)のパラメータ推定の高速化

今日の川柳

行列の分解がトピックモデルの一種として解釈できることは以下に書いた:
行列の知識ゼロからはじめてトピックモデル(GaP)の結果だけ理解する - 廿TT

ここでは
[math/0604410] Discrete Component Analysis
の Gamma-Poisson モデル(GaP)の変分推論がちょっと工夫するとより簡単になることを示す。

モデルと変分推論

カウントデータの行列  Y = (y_{n,k}) があるとする。これを2つの正の実数の行列  W=(w_{n,l}) H=(h_{l,k}) で、
 Y \approx WH
と近似することを考える。

そのために、以下のような生成モデルを考え、 w_{n,l}h_{l,k} を推定する。

 y_{n,k} \sim \mathrm{Poisson}(\sum_{l=1}^{L}w_{n,l}h_{l,k})
 w_{n,l} \sim \mathrm{Gamma}(a,b)
 h_{l,1:K} \sim \mathrm{Dirichlet}(\alpha_{1:K})

このままだと対数尤度にしたとき対数の中に足し算が入る形になり、計算しづらい。なので、ポアソン分布の部分は以下のように書き換える。

 s_{n,l,k} \sim \mathrm{Poisson}(w_{n,l}h_{l,k})
 y_{n,k}= \sum_{l=1}^{L}s_{n,l,k}

こうやって潜在変数 s_{n,l,k}かますことで平均場近似による近似事後分布は以下のように導ける。

\log q(w_{n,l}) =( \sum_{k=1}^{K}s_{n,l,k} +a-1) \log w_{n,l} -(b+1)w_{n,l} + \mathrm{const.}
\log q(h_{l,1:K}) =\sum_{k=1}^{K}( \sum_{n=1}^{N}s_{n,l,k} +\alpha_k-1) \log h_{l,k} + \mathrm{const.}

これらはそれぞれガンマ分布、ディリクレ分布に対応する。

また、ポアソン分布を和で条件付けてやると多項分布になるので、潜在変数については以下の更新式を用いる。

E[s_{n,l,k}] = y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n
 Z_n = \sum_{l=1}^{L}\exp(E[\log w_{n,l}]E[\log h_{l,k}])

以上をそのまま実装してもアルゴリズムは動くが、更新式をもう少し簡単にすることができる。

潜在変数 s_{n,l,k} についての結果を \sum_{k=1}^{K}s_{n,l,k} に代入すると、

\sum_{k=1}^{K}s_{n,l,k} \\= \sum_{k=1}^{K} y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/\{ \sum_{l=1}^{L}\exp(E[\log w_{n,l}]E[\log h_{l,k}]) \}

\sum_{n=1}^{N}s_{n,l,k} に代入すると、
\sum_{n=1}^{N}s_{n,l,k} \\= \sum_{n=1}^{N} y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/\{ \sum_{l=1}^{L}\exp(E[\log w_{n,l}]E[\log h_{l,k}]) \}

これにより潜在変数 s_{n,l,k} (3次元配列)を保存することなく、パラメータ推定を行うことができる。

添字に注意して整理すると、\sum_{k=1}^{K}s_{n,l,k} の行列は、

E[\exp(\log W)] \circ \{(Y/(\exp(E[\log W])\exp(E[\log H])))(E[\log H]^\top)\}

\sum_{n=1}^{N}s_{n,l,k} の行列は、

E[\exp(\log H)] \circ \{(E[\log W]^\top)(Y/(\exp(E[\log W])\exp(E[\log H])))\}

となる。ここで \circ は要素ごとのかけ算、/ は要素ごとのわり算を表す。

R による実装例

NMFVB <-function(Y,L=2,a=1,b=1,alpha=rep(1,ncol(Y)),maxit=1000){
  N <- nrow(Y)
  K <- ncol(Y)
  EelW <- matrix(rgamma(N*L,shape=a,rate=1),N,L)
  unnorm <- matrix(alpha+colSums(Y),L,K,byrow = TRUE)
  EelH <- unnorm/rowSums(unnorm)
  log1pb <- log1p(b)
  for(iter in 1:maxit){
    Sw <- EelW * (((Y)/(EelW %*% EelH)) %*% t(EelH))
    Sh <- EelH * (t(EelW) %*% (Y/(EelW %*% EelH)))
    alpha_W <- a + Sw
    alpha_H <- alpha + Sh
    den <- rowSums(alpha_H)
    EelW <-exp(digamma(alpha_W)-log1pb)
    EelH <-exp(digamma(alpha_H)-digamma(den))
  }
  EW <- alpha_W/(b+1)
  EH <- alpha_H/den
  return(list(W=EW,H=EH,alpha_W=alpha_W,beta_W=1+b))
}

簡単なシミュレーションを行ってみる。

N <- 100
K <- 10
L <- 3
W <- matrix(rgamma(N*L,shape=1,scale=1000),N,L)
W <- W[,order(W[1,])]
unnormH <- matrix(rgamma(L*K,shape=1,rate=1),L,K)
H <- unnormH/rowSums(unnormH)

Y <- matrix(rpois(N*K,W%*%H),N,K)

out <- NMFVB(Y,L=L,a=0.5,b=0)


plot(out$W[,order(out$W[1,])],W, main="W", xlab = "estimates", ylab = "true value")
abline(0,1,lty=2)
plot(out$H[order(out$W[1,]),],H, main="H", xlab = "estimates", ylab = "true value")
abline(0,1,lty=2)

f:id:abrahamcow:20190323162911p:plain

f:id:abrahamcow:20190323162922p:plain

Wのスケールが大きい場合、シミュレーションで設定した真値と推定値がそこそこ一致することがわかる。

ただしラベルスイッチングの問題があるので、 w_{11} < w_{12} < w_{13} となるようソートし直している。

まったく同じ議論が Dirichlet-multinomial モデルのパラメータ推定にも適応できる:
変分ベイズによるトピックモデル(Dirichlet-multinomial Model)のパラメータ推定の高速化 - 廿TT