廿TT

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

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

今日の川柳

変分ベイズによるトピックモデル(GaP; Gamma-Poisson Model)のパラメータ推定の高速化 - 廿TT とまったく同じ議論により、
[math/0604410] Discrete Component Analysis
の Dirichlet-multinomial モデルの変分推論もより簡単にすることができる。

Dirichlet-multinomial モデルにより推定されるパラメータ W は Gamma-Poisson モデルの W が行ごとに足して1になるようノーマライズされたもの、と捉えればよい。

モデルと変分推論

以下の生成モデルを考え、 w_{n,l}h_{l,k} を推定する。

 w_{n,1:L} \sim \mathrm{Dirichlet}(\beta_{1:L})
 h_{l,1:K} \sim \mathrm{Dirichlet}(\alpha_{1:K})
 v_{n,1:L} \sim \mathrm{Multinomial}(T_{n}, w_{n,1:L})
 s_{n,l,1:K} \sim \mathrm{Multinomial}(v_{n,l}, h_{l,1:K})
 y_{n,k} = \sum_{l=1}^{L}s_{n,l,k}

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

\log q(w_{n,1:L}) =\sum_{l=1}^{L}(\sum_{k=1}^{K}s_{n,l,k} +\beta_l-1)\log 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 による実装例

library(gtools)
NMFVB <-function(Y,L=2,alpha=1,beta=1,maxit=1000){
  N <- nrow(Y)
  K <- ncol(Y)
  EelW <- gtools::rdirichlet(N,rep(beta,L))
  EelH <- gtools::rdirichlet(L,rep(alpha,K))
  for(iter in 1:maxit){
    Sw <- EelW * (((Y)/(EelW %*% EelH)) %*% t(EelH))
    Sh <- EelH * (t(EelW) %*% (Y/(EelW %*% EelH)))
    beta_W <- beta + Sw
    alpha_H <- alpha + Sh
    EelW <-exp(sweep(digamma(beta_W),1,digamma(rowSums(beta_W))))
    EelH <-exp(sweep(digamma(alpha_H),1,digamma(rowSums(alpha_H))))
  }
  EW <- sweep(beta_W,1,rowSums(beta_W),"/")
  EH <- sweep(alpha_H,1,rowSums(alpha_H),"/")
  return(list(W=EW,H=EH,alpha_H=alpha_H,beta_W=beta_W))
}

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

W <- rdirichlet(100,rep(1,3))
W <- W[,order(W[1,])]
H <- rdirichlet(3,rep(1,100))
Y <- t(apply(W%*%H,1,function(p){rmultinom(1,10000,p)}))
L=3
out <- NMFVB(Y,L=L)
plot(out$W[,order(out$W[1,])],W, main="W", xlab = "estimates", ylab = "true value")
abline(0,1,col="red")
plot(out$H[order(out$W[1,]),],H, main="H", xlab = "estimates", ylab = "true value")
abline(0,1,col="red")

f:id:abrahamcow:20190323203924p:plain

f:id:abrahamcow:20190323203905p:plain

シミュレーションで設定した真値と推定値がそこそこ一致することがわかる。

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