はじめに

　『トピックモデルによる統計的潜在意味解析』の学習時のメモです。本と併せて読んでください。

　この記事では、3.4.3項のLDAの確率的変分ベイズ法について書いています。図3.6の疑似コードを基にR言語で実装していきます。

　プログラムからアルゴリズムの理解を目指します。できるだけ基本的な(直感的に理解できる)関数を用いて組んでいきます。

【数理編】

www.anarchive-beta.com

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

はじめに
3.4.3 LDAの確率的変分ベイズ法
参考文献
おわりに

3.4.3 LDAの確率的変分ベイズ法

　図3.6の疑似コードを参考にLDAに対する確率的変分ベイズ法を実装していきます。

# 利用パッケージ
library(tidyverse)

・簡易文書データ(クリックで展開)

## 簡易文書データ
# 文書数
M <- 10

# 語彙数
V <- 20

# 文書ごとの各語彙の出現回数
n_dv <- matrix(sample(0:10, M * V, replace = TRUE), M, V)

# 文書ごとの単語数
n_d <- apply(n_dv, 1, sum)

・パラメータの設定

　始めにパラメータ類の初期値を設定していきます。

# サンプリング回数を指定
S <- 1000

# イタレーション数を指定
InIter <- 25

# ステップサイズ(学習率)を指定
nu <- 0.01

# トピック数を指定
K <- 5

　サンプリング回数(試行回数)をS、内側ループの回数をInIterとして指定します。

　また学習率をnu、トピック数をKとします。

# 事前分布のパラメータを指定
alpha_k <- rep(2, K)
beta_v  <- rep(2, V)

　トピック分布のパラメータ$\boldsymbol{\alpha} = (\alpha_1, \cdots, \alpha_K)$、単語分布のパラメータ$\boldsymbol{\beta} = (\beta_1, \cdots, \beta_V)$の値をそれぞれ指定します。どちらもディリクレ分布のパラメータなので、$\alpha_k > 0,\ \beta_v > 0$である必要があります。

　この例では値を一様に設定します。

# 事後分布のパラメータを指定
xi_theta_dk <- rep(n_d / K, K) + rep(alpha_k, each = M) %>% # E_n_dk + alpha_dk
  matrix(nrow = M, ncol = K)
xi_phi_kv <- seq(1, 5) %>% # とり得る値の範囲を指定
  sample(size = K * V, replace = TRUE) %>% 
  matrix(nrow = K, ncol = V)

　パラメータ$\boldsymbol{\theta},\ \boldsymbol{\phi}$(の分布)の初期値を設定するイメージで(?)、事後分布のパラメータ$\xi_{d,k}^{\theta},\ \xi_{k,v}^{\phi}$の初期値を設定します。

　$\xi_{d,k}^{\theta}$は、定義式(3.89)

$$ \xi_{d,k}^{\theta} = \mathbb{E}_{q(\boldsymbol{z}_d)}[n_{d,k}] + \alpha_k \tag{3.89} $$

に従い設定します。ただし$\mathbb{E}_{q(\boldsymbol{z}_d)}[n_{d,k}]$の初期値については

$$ \mathbb{E}_{q(\boldsymbol{z}_d)}[n_{d,k}] = \frac{n_d}{K} $$

各文書でどのトピックに対しても等しくなるようにします。

　$\xi_{d,k}^{\theta}$の初期値は、値の範囲を決めてランダムに設定します。

# 潜在トピック集合の分布
z_dv_k <- array(0, dim = c(M, V, K))
for(d in 1:M) {
  for(v in 1:V) {
    if(n_dv[d, v] > 0) {
      # ランダムに値を生成
      tmp_q_z <- sample(seq(0, 1, by = 0.01), size = K, replace = TRUE)
      # 正規化
      z_dv_k[d, v, ] <- tmp_q_z / sum(tmp_q_z)
    }
  }
}

　各単語に潜在トピックを割り当てる確率$q(\boldsymbol{z})$の初期値をランダムに設定します。

　$q(z_{d,i} = k)$の計算式(3.99)より同一単語の確率は同じになることから、ここではz_dv_kとして、文書・(単語$i$ではなく)語彙・トピックの3次元配列として扱うことにします。

# 文書dの語彙vにおいてトピックkが割り当てられた単語数の期待値の受け皿
E_n_dkv <- matrix(0, nrow = K, ncol = V) ## (K次元ベクトルで扱うことに注意)

　　$\mathbb{E}_{q(\boldsymbol{z})}[n_{d,v,k}]$は、ループ処理の中で添字を使って代入していくため、予め変数を用意しておく必要があります。ただしこの変数についてもサンプリングされた文書以外の値を保存しておく必要がないため、(命名規則には反していますが)初期値が0の$M$行$K$列のマトリクスをE_n_dkvとします。

・確率的変分ベイズ法

　続いて推論部分を実装していきます。

・コード全体(クリックで展開)

# 推移の確認用
trace_xi_theta <- array(0, dim = c(M, K, S + 1))
trace_xi_phi   <- array(0, dim = c(K, V, S + 1))
# 初期値を代入
trace_xi_theta[, , 1] <- xi_theta_dk
trace_xi_phi[, , 1]   <- xi_phi_kv

for(s in 1:S) { ## (サンプリング)
  
  # 動作確認用に開始時間を記録
  start_time <- Sys.time()
  
  # 文書(番号)をサンプリング
  d <- sample(1:M, size = 1)
  
  for(II in 1:InIter) { ## (内側ループ)
    
    for(v in 1:V) { ## (各語彙)
      if(n_dv[d, v] > 0) {
        
        # 潜在トピック集合の近似事後分布を計算:式(3.99)
        term1 <- digamma(xi_phi_kv[, v]) - digamma(apply(xi_phi_kv, 1, sum))
        term2 <- digamma(xi_theta_dk[d, ]) - digamma(sum(xi_theta_dk[d, ]))
        tmp_q_z <- exp(term1 + term2)
        z_dv_k[d, v, ] <- tmp_q_z / sum(tmp_q_z) # 正規化
        
      }
    } ## (/各語彙)
    
    # 文書dにおいてトピックkが割り当てられた単語数の期待値を計算
    tmp_n_vk <- z_dv_k[d, , ] * n_dv[d, ]
    E_n_dk <- apply(tmp_n_vk, 2, sum) ## (K次元ベクトルで扱うことに注意)
    
    # トピック分布の近似事後分布のパラメータを更新:式(3.89)
    xi_theta_dk[d, ] <- E_n_dk + alpha_k
    
  } ## (/内側ループ)
  
  for(k in 1:K) { ## (各トピック)
    
    # 文書dの語彙vにおいてトピックkが割り当てられた単語数の期待値を計算
    E_n_dkv[k, ] <- tmp_n_vk[, k]
    
    # 単語分布の近似事後分布のパラメータを更新:式(3.159)
    xi_phi_kv[k, ] <- xi_phi_kv[k, ] + nu * (M * E_n_dkv[k, ] + beta_v - xi_phi_kv[k, ])
    
  } ## (/各トピック)
  
  # 推移の確認用に値を保存
  trace_xi_theta[, , s + 1] <- xi_theta_dk
  trace_xi_phi[, , s + 1]   <- xi_phi_kv
  
  # 動作確認
  print(paste0("s=", s, "(", round(s / S * 100, 1), "%)", "...", round(Sys.time() - start_time, 3)))
}

・$q(\boldsymbol{z})$の更新

# 潜在トピック集合の近似事後分布を計算:式(3.99)
term1 <- digamma(xi_phi_kv[, v]) - digamma(apply(xi_phi_kv, 1, sum))
term2 <- digamma(xi_theta_dk[d, ]) - digamma(sum(xi_theta_dk[d, ]))
tmp_q_z <- exp(term1 + term2)
z_dv_k[d, v, ] <- tmp_q_z / sum(tmp_q_z) # 正規化

　語彙ごとに、式(3.99)の計算を行い値を更新します。

$$ q(z_{d,i} = k) \propto \frac{ \exp \left[ \Psi( \xi_{k,w_{d,i}}^{\phi} ) \right] }{ \exp \left[ \Psi \left( \sum_{v'=1}^V \xi_{k,v'}^{\phi} \right) \right] } \frac{ \exp \left[ \Psi( \xi_{d,k}^{\theta} ) \right] }{ \exp \left[ \Psi \left( \sum_{k'=1}^K \xi_{d,k'}^{\theta} \right) \right] } \tag{3.99} $$

　またこの計算結果を総和が1になるように正規化します。

　項ごとにexp()の計算を行う必要もないので、最後にまとめてすることにします。それによって「*」が「+」に、「/」が「-」に変わっています。

　ちなみに式通りだとこうですね。

# 潜在トピック集合の近似事後分布を計算:式(3.99)
term1 <- exp(digamma(xi_phi_kv[, v])) / exp(digamma(apply(xi_phi_kv, 1, sum)))
term2 <- exp(digamma(xi_theta_dk[d, ])) / digamma(sum(xi_theta_dk[d, ]))
tmp_q_z <- term1 * term2
z_dv_k[d, v, ] <- tmp_q_z / sum(tmp_q_z) # 正規化

・$q(\boldsymbol{\theta})$の更新

# 文書dにおいてトピックkが割り当てられた単語数の期待値を計算
tmp_n_vk <- z_dv_k[d, , ] * n_dv[d, ]
E_n_dk <- apply(tmp_n_vk, 2, sum) ## (K次元ベクトルで扱うことに注意)

　$\mathbb{E}_{q(\boldsymbol{z}_d)} [n_{d,k}]$の計算式

$$ \mathbb{E}_{q(\boldsymbol{z}_d)} [n_{d,k}] = \sum_{i=1}^{n_d} q(z_{d,i} = k) $$

に従い、E_n_dkを更新します。ただしここでは語彙ごとに求めているので、各語彙に対して出現回数(n_dv[d, ])を掛けることで$\sum_{i=1}^{n_d}$の計算となります。

　初期値の設定時にE_n_dkをM行K列のマトリクスとして作成しましたが、このアルゴリズムでは今サンプリング(更新)している文書以外の$\mathbb{E}_{q(\boldsymbol{z}_d)} [n_{d,k}]$の値を保存する必要がないことが良さの1つです。なので更新値を(文書3の推定中であればE_n_3kのようなイメージで)要素がK個のベクトルとして再代入することにします。)

　次にこの統計量E_n_dkを用いて、ハイパーパラメータxi_theta_dk[d, ]を更新します。また更新にはtmp_n_vkも使います。

# トピック分布の近似事後分布のパラメータを更新:式(3.89)
xi_theta_dk[d, ] <- E_n_dk + alpha_k

　$q(\boldsymbol{\theta}_d)$のパラメータの更新式(3.89)の計算を行い、ハイパーパラメータを更新します。

　これを文書1からMまで繰り返し、$q(\boldsymbol{\theta})$を全て更新します。

・$q(\boldsymbol{\phi})$の更新

# 文書dの語彙vにおいてトピックkが割り当てられた単語数の期待値を計算
E_n_dkv[k, ] <- tmp_n_vk[, k]

　$\mathbb{E}_{q(\boldsymbol{z})}[n_{d,k,v}]$の定義式

$$ \mathbb{E}_{q(\boldsymbol{z})}[n_{d,k,v}] = \sum_{i=1}^{n_d} q(z_{d,i} = k) \delta(w_{d,i} = v) $$

に従い、E_n_dkv[k, ]を更新します。$\delta(w_{d,i} = v)$は、単語を語彙に変換するための項のようなものですが、この例ではそもそも語彙として扱っているので、特に計算の必要はありません。

　次にこの統計量E_n_dkv[k, ]を用いて、ハイパーパラメータxi_phi_kv[k, ]を更新します。

# 単語分布の近似事後分布のパラメータを更新:式(3.159)
xi_phi_kv[k, ] <- xi_phi_kv[k, ] + nu * (M * E_n_dkv[k, ] + beta_v - xi_phi_kv[k, ])

　$q(\boldsymbol{\phi}_k)$のパラメータの更新式(3.95)

$$ \xi_{k,v}^{\phi (s+1)} = \xi_{k,v}^{\phi (s)} + \nu_s ( M \mathbb{E}_{q(\boldsymbol{z})}[ n_{d,k,v} ] + \beta_{v} - \xi_{k,v}^{\phi (s)} ) \tag{3.159} $$

の計算を行い、ハイパーパラメータを更新します。

　これをトピック1からKまで繰り返し、$q(\boldsymbol{\phi})$を全て更新します。

　分かりやすさのためにも極力疑似コードに忠実に組むのがコンセプトですが、この部分はfor(k in 1:K)によってトピック1つずつ更新するよりも次のようにひとまとめに処理した方がむしろ分かりやすいかもしれません。

# 文書dの語彙vにおいてトピックkが割り当てられた単語数の期待値を計算
E_n_dkv <- t(tmp_n_vk)

# 単語分布の近似事後分布のパラメータを更新:式(3.159)
xi_phi_kv <- xi_phi_kv + nu * (M * E_n_dkv + beta_v - xi_phi_kv)

・推定結果の確認

　ggplot2パッケージを利用して、推定結果を可視化しましょう。

　ディリクレ分布の期待値(2.10)の計算を行いパラメータの期待値を作図します。

・トピック分布(クリックで展開)

## トピック分布(期待値)
# thetaの期待値を計算:式(2.10)
theta_dk <- xi_theta_dk / apply(xi_theta_dk, 1, sum)

# 作図用のデータフレームを作成
theta_df_wide <- cbind(
  as.data.frame(theta_dk), 
  doc = as.factor(1:M) # 文書番号
)

# データフレームをlong型に変換
theta_df_long <- pivot_longer(
  theta_df_wide, 
  cols = -doc, # 変換せずにそのまま残す現列名
  names_to = "topic", # 現列名を格納する新しい列の名前
  names_prefix = "V", # 現列名から取り除く文字列
  names_ptypes = list(topic = factor()), # 現列名を要素とする際の型
  values_to = "prob" # 現要素を格納する新しい列の名前
)

# 作図
ggplot(theta_df_long, aes(x = topic, y = prob, fill = topic)) + 
  geom_bar(stat = "identity", position = "dodge") + # 棒グラフ
  facet_wrap(~ doc, labeller = label_both) + # グラフの分割
  labs(title = "Stochastic Variational Bayes for LDA", 
       subtitle = expression(xi^theta)) # ラベル

・単語分布(クリックで展開)

## 単語分布(期待値)
# phiの期待値を計算:式(2.10)
phi_kv <- xi_phi_kv / apply(xi_phi_kv, 1, sum)

# 作図用のデータフレームを作成
phi_df_wide <- cbind(
  as.data.frame(phi_kv), 
  topic = as.factor(1:K)
)

# データフレームをlong型に変換
phi_df_long <- pivot_longer(
  phi_df_wide, 
  cols = -topic, # 変換せずにそのまま残す現列名
  names_to = "word", # 現列名を格納する新しい列の名前
  names_prefix = "V", # 現列名から取り除く文字列
  names_ptypes = list(word = factor()), # 現列名を要素とする際の型
  values_to = "prob" # 現要素を格納する新しい列の名前
)

# 作図
ggplot(phi_df_long, aes(x = word, y = prob, fill = word, color = word)) + 
  geom_bar(stat = "identity", position = "dodge") + # 棒グラフ
  facet_wrap(~ topic, labeller = label_both) + # グラフの分割
  theme(legend.position = "none") + # 凡例
  scale_x_discrete(breaks = seq(0, V, by = 10)) + # x軸目盛(連続値)
  labs(title = "Stochastic Variational Bayes for LDA", 
       subtitle = expression(xi^phi)) # ラベル

・更新値の推移の確認

　続いて、各分布のパラメータの更新値の推移を可視化します。

・トピック分布のパラメータ(クリックで展開)

## トピック分布のパラメータ
# 作図用データフレームを作成
trace_xi_theta_df_wide <- data.frame()
for(s in 1:(S + 1)) {
  # データフレームに変換
  tmp_theta_df <- cbind(
    as_tibble(trace_xi_theta[, , s]), 
    doc = as.factor(1:M), 
    sampling = s - 1
  )
  # 結合
  trace_xi_theta_df_wide <- rbind(trace_xi_theta_df_wide, tmp_theta_df)
}

# データフレームをlong型に変換
trace_xi_theta_df_long <- pivot_longer(
  trace_xi_theta_df_wide, 
  cols = -c(doc, sampling), # 変換しない列
  names_to = "topic", # 現列名を格納する列の名前
  names_prefix = "V", # 現列名から取り除く文字列
  names_ptypes = list(topic = factor()), # 現列名を格納する際のデータ型
  values_to = "value" # 現セルを格納する列の名前
)

# 文書番号を指定
doc_num <- 16

# 作図
trace_xi_theta_df_long %>% 
  filter(doc == doc_num) %>% 
  ggplot(aes(x = sampling, y = value, color = topic)) + 
    geom_line() + # 折れ線グラフ
    labs(title = "Stochastic Variational Bayes for LDA", 
         subtitle = paste0("d=", doc_num)) # ラベル

・単語分布のパラメータ(クリックで展開)

## 単語分布のパラメータ
# 作図用データフレームを作成
trace_xi_phi_df_wide <- data.frame()
for(s in 1:(S + 1)) {
  # データフレームに変換
  tmp_trace_xi_phi_df <- cbind(
    as_tibble(trace_xi_phi[, , s]), 
    topic = as.factor(1:K), 
    sampling = s - 1
  )
  # 結合
  trace_xi_phi_df_wide <- rbind(trace_xi_phi_df_wide, tmp_trace_xi_phi_df)
}

# データフレームをlong型に変換
trace_xi_phi_df_long <- pivot_longer(
  trace_xi_phi_df_wide, 
  cols = -c(topic, sampling), # 変換しない列
  names_to = "word", # 現列名を格納する列の名前
  names_prefix = "V", # 現列名から取り除く文字列
  names_ptypes = list(word = factor()), # 現列名を格納する際のデータ型
  values_to = "value" # 現セルを格納する列の名前
)

# トピック番号を指定
topic_num <- 1

# 作図
trace_xi_phi_df_long %>% 
  filter(topic == topic_num) %>% 
  ggplot(aes(x = sampling, y = value, color = word)) + 
    geom_line(alpha = 0.5) + # 折れ線グラフ
    theme(legend.position = "none") + # 凡例
    labs(title = "Stochastic Variational Bayes for LDA", 
         subtitle = paste0("k=", topic_num)) # ラベル

・gif画像を作成

・コード(クリックで展開)

　更新ごとの分布の変化をgif画像でみるためのコードです。gganimateパッケージを利用してgif画像を作成します。作図には上で作成したデータフレームを使います。

# 利用パッケージ
library(gganimate)

## トピック分布のパラメータ
# 作図
graph_xi_theta <- ggplot(trace_xi_theta_df_long, aes(x = topic, y = value, fill = topic)) + 
  geom_bar(stat = "identity", position = "dodge") + # 棒グラフ
  facet_wrap(~ doc, labeller = label_both) + # グラフの分割
  transition_manual(sampling) + # フレーム
  labs(title = "Stochastic Variational Bayes for LDA", 
       subtitle = "s={current_frame}") # ラベル

# gif画像を作成
animate(graph_xi_theta, nframes = S + 1, fps = 20)

## 単語分布のパラメータ
# 作図
graph_xi_phi <- ggplot(trace_xi_phi_df_long, aes(x = word, y = value, fill = word, color = word)) + 
  geom_bar(stat = "identity", position = "dodge") + # 棒グラフ
  facet_wrap(~ topic, labeller = label_both) + # グラフの分割
  theme(legend.position = "none") + # 凡例
  scale_x_discrete(breaks = seq(0, V, by = 10)) + # x軸目盛(連続値)
  transition_manual(sampling) + # フレーム
  labs(title = "Stochastic Variational Bayes for LDA", 
       subtitle = "s={current_frame}") # ラベル

# gif画像を作成
animate(graph_xi_phi, nframes = S + 1, fps = 20)

　transition_manual()に時系列値の列を指定します。

(やっぱりファイルサイズが大きいとダメなのかアップできないのでTwitterに上げときましたのでそれで…)

LDAの確率的変分ベイズ法：単語分布のパラメータ pic.twitter.com/2GBnpAHdKM
— しょこ📚 (@anemptyarchive) 2020年5月22日

参考文献

佐藤一誠『トピックモデルによる統計的潜在意味解析』(自然言語処理シリーズ 8)奥村学監修,コロナ社,2015年.

おわりに

　なんとなくですがこの手法が一番好きです。組むのが楽しかったとか、gif画像で見たときの分布の変化の仕方がなんかかわいいとかそんな理由です。
　パラメータの推移を折れ線グラフでみると、ギブスサンプリングはギザギザで、変分ベイズは滑らかで、これはカクカクしてるというのだけでも、(アルゴリズム的に当然だと思えたのも含めて)自分で組んで実感して面白いものです。

　楽しかったとは言うもののいくつかの点で詰まりました。

　1つ目は逐次学習の意味を捉え違っていたことです。取得した文書データはその場限りで推定するのだと理解していました。いくつかのアルゴリズムを並行して組んでたのですが、次の粒子フィルタでいうところの逐次学習は理解の通りだったので気付くのに時間がかかりました。

　2つ目は学習率の値ですね。学習率が1のとき更新式(3.159)は変分ベイズの更新式と等しくなるとのことだったので(合わせた方が比較しやすかったりするのかなぁと)$\nu = 1$に設定して組んでいたのですが、中々収束せずというか分布が変な動きしているので、組み方が間違っているのだと思いアレコレ弄ってました。が、何かの拍子で学習率を変えてみたらとっくに組めていたというオチ。あぁこれが値が小さすぎると学習に時間がかかり、逆に大きすぎると正しい値を飛び越えてしまうというヤツか～

　理解してしまえばどちらも本にちゃんと書いてあるんですよねぇ。でも解らないと分からないんですよね。

　という訳で、LDA基本編は次でラストです！ちなみに実装できませんでした！！が、前回と同様に一旦諦めてアップします！！！

【次節の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【R】3.4.3：LDAの確率的変分ベイズ法【白トピックモデルのノート】

はじめに

3.4.3 LDAの確率的変分ベイズ法

・パラメータの設定

・確率的変分ベイズ法

・推定結果の確認

・更新値の推移の確認

参考文献

おわりに