ジョンとヨーコのイマジン日記

想像してください。「あなたはぼくをプラグマティストだと言うかもしれない」と歌う、逆イマジンです。

トピックモデル(GaP; Gamma-Poisson Model)の ELBO の導出

ここでは
[math/0604410] Discrete Component Analysis
の Gamma-Poisson モデル(GaP)の ELBO (evidence lower bound) を導出する。

まず、行列の分解がトピックモデルの一種として解釈できること説明する。

次に、モデルのパラメータ推定方法について述べる。

そしてELBOを導出し、最後にRによるシミュレーションを行う。

行列の復習からはじめてトピックモデルの結果だけ理解する

トピックモデルと総称されるモデルの中には様々なものがありますが、ここでは
[math/0604410] Discrete Component Analysis
の Gamma-Poisson モデル(GaP)を紹介します。

トピックモデルはカウントデータ(なにかの数を数えたデータ)の行列が与えられたとき、それをより小さい行列で近似して理解しやすくするためのモデルです。

行列というのはこんな感じの表のことです。

f:id:abrahamcow:20190211093407p:plain

多くのトピックモデルでは、単語の並び順は無視して、単語の出現頻度にのみ注目し、文書ごとに集計した表を利用します。
ここでは簡単のため3単語のみを集計したことにしていますが、もっと数が増えても大丈夫です。
その場合行列が横に長くなっていきます。

トピックモデルではこの行列がいくつかのトピックが混ざった結果生成されたと考えます。

トピックは単語の出現確率の違いで特徴づけられます。

例えば「刀」とか「印籠」という単語が出てくる確率は時代劇だったら高そうですが、現代劇にはめったに出てこなさそうですよね。

このような単語の出現確率の違いが文書のトピックを特徴づけていると考えるのがトピックモデルです。

トピックも行列で表しておくと便利です。

f:id:abrahamcow:20190211095029p:plain

今回の行列はトピックごとの単語の出現確率を表しているので、横合計が1になります。

トピック1は単語3が出る確率が高め、トピック2は単語2が出る確率が高めです。

ここでは簡単のため2種類のトピックのみを考えていますが、もっと数が増えても大丈夫です。

その場合行列が縦に長くなっていきます。

さて、各文書はトピックが混ざった結果出現すると書きました。

たとえば文書1はトピック1が40、トピック2が10出現したとしましょう。

表と見比べながら慎重に計算していくと、この文書の単語の出現頻度は次のようになります。

f:id:abrahamcow:20190211100755p:plain

トピック1に40をかけて、トピック2に10をかけています。単語はどちらのトピックにおいても出現確率を持っているので、トピック1から出現した数とトピック2から出現した数を足しています。

このような計算を行列では以下のように省略して書きます。

\begin{pmatrix}
40 & 10\\
\end{pmatrix}
\begin{pmatrix}
0.1 & 0.2 & 0.7\\
0.3 & 0.5 & 0.2
\end{pmatrix}=\begin{pmatrix}
7 & 13 & 30\\
\end{pmatrix}

意味は

40 \times 0.1 + 10 \times 0.3 =7\\
40 \times 0.2 + 10 \times 0.5 = 13\\
40 \times 0.7 + 10 \times 0.2 = 30

とまったく一緒です。

次に文書2はトピック1が20、トピック2が90出現したとしましょう。

行列を使って省略して書くと

\begin{pmatrix}
20 & 90\\
\end{pmatrix}
\begin{pmatrix}
0.1 & 0.2 & 0.7\\
0.3 & 0.5 & 0.2
\end{pmatrix}=\begin{pmatrix}
29 & 49 & 32\\
\end{pmatrix}

となります。

行列ではこの2回の計算をさらに省略して、

\begin{pmatrix}
40 & 10\\
20 & 90\\
\end{pmatrix}
\begin{pmatrix}
0.1 & 0.2 & 0.7\\
0.3 & 0.5 & 0.2
\end{pmatrix}
= \begin{pmatrix}
7 & 13 & 30 \\
29 & 49 & 32
\end{pmatrix}

と書きます。右辺は1行目に文書1についての計算結果、2行目に文書2について計算結果が入ります。

文書の数が増えても左辺の一個目の行列が長くなるだけで、まったく同じように書くことができます。


\begin{pmatrix}
40 & 10\\
20 & 90\\
\vdots & \vdots\\
40 & 70
\end{pmatrix}
\begin{pmatrix}
0.1 & 0.2 & 0.7\\
0.3 & 0.5 & 0.2
\end{pmatrix}
= \begin{pmatrix}
7 & 13 & 30 \\
29 & 49 & 32 \\
\vdots & \vdots & \vdots\\
25 & 43 & 42 \\
\end{pmatrix}

以上がトピックモデルのデータ生成過程です。トピックモデルではこれの逆問題をやります。つまり文書と単語の行列がつぎのように与えられたとき、


\begin{pmatrix}
7 & 13 & 30 \\
29 & 49 & 32 \\
\vdots & \vdots & \vdots\\
25 & 43 & 42 \\
\end{pmatrix}

この行列の分解


\begin{pmatrix}
40 & 10\\
20 & 90\\
\vdots & \vdots\\
40 & 70
\end{pmatrix}
\begin{pmatrix}
0.1 & 0.2 & 0.7\\
0.3 & 0.5 & 0.2
\end{pmatrix}

を探すのです。

トピックモデルは確率モデルですので本当は


\begin{pmatrix}
7 & 13 & 30 \\
29 & 49 & 32 \\
\vdots & \vdots & \vdots\\
25 & 43 & 42 \\
\end{pmatrix} \approx
\begin{pmatrix}
40 & 10\\
20 & 90\\
\vdots & \vdots\\
40 & 70
\end{pmatrix}
\begin{pmatrix}
0.1 & 0.2 & 0.7\\
0.3 & 0.5 & 0.2
\end{pmatrix}

となる行列を探します。\approx はニアリーイコールです。トピックモデルで得られる行列のかけ算の結果は、もとの行列と厳密には一致しませんが、一致しない分は確率的な誤差とみなすのです。

モデルと変分推論

カウントデータの行列  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})

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

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

ここで  u_{n,l} = \sum_{k=1}^{K} s_{n,l,k} である。

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

\log q(w_{n,l}) =( \sum_{k=1}^{K}\hat a-1) \log w_{n,l} -\hat bw_{n,l} + \mathrm{const.}
\log q(h_{l,1:K}) =\sum_{k=1}^{K}( \hat \alpha_{l,k}-1) \log h_{l,k} + \mathrm{const.}

これらはそれぞれガンマ分布、ディリクレ分布である。

ここで
 \hat a =s_{n,l,k} +a
 \hat b =b+1
 \hat \alpha_{l,k} = \sum_{n=1}^{N} s_{n,l,k} +\alpha_k
とおいた。

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

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 は要素ごとのかけ算、/ は要素ごとのわり算を表す。

ELBO

潜在変数を Z、観測されるデータを Y としたとき、ELBO は次の式で与えられる.

 \operatorname{ELBO} = E[ \log p (Y|Z) + \log p (Z) -\log q(Z) ]

ここでの期待値は変分事後分布 q(Z) による期待値。

GaP に出てくる確率(密度)関数の対数をすべて書き下す。

 \log p(u_{n,l}) = u_{n,l} \log w_{n,l} - w_{n,l} - \log \Gamma(u_{n,l}+1)

 \log p(s_{n,l,1:K}) = \sum_k s_{n,l,k} \log h_{l,k} + \log \Gamma(u_{n,l}+1)- \sum_k  \log \Gamma(s_{n,l,k}+1)

 \log q(s_{n,l,1:K}) = s_{n,l,1:K} \log (\exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n) + \log \Gamma (y_{n,k}+1) \\- \sum_l \log \Gamma (s_{n,l,k}+1)

 \log p(w_{n,l}) = (a-1) \log w_{n,l} -b w_{n,l} + a\log b -\log \Gamma(a)

 \log q(w_{n,l}) = (\hat a-1) \log w_{n,l} -\hat b w_{n,l} + a\log b -\log \Gamma(a)

 \log p (h_{l,1:K}) = \sum_{k} (\alpha_k-1) \log h_{l,k} + \log \Gamma (\sum_{k} \alpha_k) - \sum_{k} \log \Gamma(\alpha_k)

 \log q (h_{l,1:K}) = \sum_{k} (\hat \alpha_{l,k}-1) \log h_{l,k} + \log \Gamma (\sum_{k} \hat \alpha_{l,k}) \\- \sum_{k} \log \Gamma(\hat \alpha_{l,k})

期待値の計算の難しい  \log \Gamma(s_{n,l,k}+1) \log \Gamma(u_{n,l}+1) の項はキャンセルされて消える。

また、
 E[u_{n,l}] + a = \hat a
  b+1 = \hat b
 E[\sum_{n=1}^{N} s_{n,l,k} +\alpha_k] = \hat \alpha_{l,k}
であるから、
 (\hat a-1) \log w_{n,l} と、
 \hat b w_{n,l} と、
 \sum_{k} (\hat \alpha_{l,k}-1) \log h_{l,k}
の項もキャンセルされて消える。

また、
 E[s_{n,l,k} ]\log (\exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n)\\
= y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n )  \log (\exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n)\\
= y_{n,k} [ (\{ \exp( E[\log w_{n,l}] ) E[\log w_{n,l}] \} \exp(E[\log h_{l,k}]) +\\  \exp( E[\log w_{n,l}] ) \{\exp(E[\log h_{n,l}]) E[\log h_{l,k}]\})/ Z_n -\log Z_n ]

R による実装例

GaPVB <-function(Y,L=2,a=0.5,b=1e-8,alpha=rep(1,ncol(Y)),maxit=500){
  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
  ELBO <- -sum(lfactorial(Y))+
    L*sum(lgamma(sum(alpha))-sum(lgamma(alpha))) +
    sum(-lgamma(rowSums(alpha_H))+rowSums(lgamma(alpha_H))) + 
    sum(-alpha_W*log1pb+a*log(b)+lgamma(alpha_W)-lgamma(a)) + 
    sum(-Y*(((EelW*log(EelW))%*%EelH + EelW%*%(EelH*log(EelH)))/(EelW%*%EelH)-log(EelW%*%EelH))) 
  return(list(W=EW,H=EH,ELBO=ELBO))
}

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

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} となるようソートし直している。

シミュレーション

トピック数を5にしてデータを生成した。

library(parallel)
simfunc <- function(i){
  set.seed(i)
  W <- matrix(rgamma(100*5,1,0.001),100,5)
  H <- matrix(rgamma(5*20,1,1),5,20)
  H <- H/rowSums(H)
  Y <- matrix(rpois(100*20,W%*%H),100,20)
  sapply(1:10, function(L)GaPVB(Y,L)$ELBO)
}
out <- mclapply(1:100,simfunc,mc.cores = detectCores())
print(table(apply(simplify2array(out),2,which.max)))
boxplot(t(simplify2array(out)))

f:id:abrahamcow:20190830221703p:plain

100回のシミュレーション中100回とも正しいトピック数を選んでいる。

> print(table(apply(simplify2array(out),2,which.max)))

  5 
100