はじめに
機械学習プロフェッショナルシリーズの『トピックモデル』の勉強時に自分の理解の助けになったことや勉強会資料のまとめです。トピックモデルの各種アルゴリズムを「数式」と「プログラム」から理解することを目指します。
この記事は、3.4節「変分ベイズ推定」の内容です。
変分ベイズ推定による混合ユニグラムモデルにおけるパラメータ推定を実装します。
【数式編】
【前節の内容】
【他の節一覧】
【この節の内容】
3.4 変分ベイズ推定：ループなし版
混合ユニグラムモデル(混合カテゴリ分布)に対する変分ベイズ推定を実装する。
「【R】3.4：混合ユニグラムモデルの変分ベイズ推定の実装【『トピックモデル』の勉強ノート】 - からっぽのしょこ」では、アルゴリズム3.2を参考にして、3重のループを使って実装した。この記事では、行列計算を使って、ループを使わずに実装する。
生成モデルの設定と文書データの生成、分布の可視化については「ループあり版」の記事を参照のこと。
利用するパッケージを読み込む。
# 利用パッケージ library(tidyverse) library(gganimate) library(gtools)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はない。ただし、作図コードに関してはパッケージ名を省略するので、ggplot2を読み込む必要がある。
また、magrittrパッケージのパイプ演算子%>%ではなく、ネイティブパイプ演算子|>を使う。%>%に置き換えても処理できるが、その場合はmagrittrを読み込む必要がある。
文書データの簡易生成
人工的に生成した次の文書データ(bag-of-words)を利用する。
各文書の語彙ごとの出現度数N_dvの作成については「ループあり版」の記事を参照のこと。
推論処理
生成した文書データを用いて、トピック分布と単語分布を推定する。
パラメータの初期化
トピック数を指定して、文書数・語彙数を設定する。
# トピック数を指定 K <- 4 # 文書数と語彙数を取得 D <- nrow(N_dv) V <- ncol(N_dv) D; V
## [1] 20 ## [1] 26
トピック数$K$を指定する。
文書数$D$はN_dvの行数、語彙数$V$は列数に対応する。
各文書の単語数を作成する。
# 文書ごとの単語数を計算 N_d <- rowSums(N_dv) N_d[1:5]
## [1] 150 109 116 105 148
各文書の語彙ごとの出現回数$N_{dv}$の語彙$v$について(行ごと)の和を計算する。
事前分布のパラメータ(ハイパーパラメータの初期値)を指定する。
# トピック分布θの事前分布のパラメータαを指定 alpha <- 1 # 単語分布Φの事前分布のパラメータβを指定 beta <- 1
トピック分布$\boldsymbol{\theta}$の事前分布(ディリクレ分布)$p(\boldsymbol{\theta} | \alpha)$のパラメータ$\alpha > 0$と、単語分布$\boldsymbol{\Phi}$の事前分布(ディリクレ分布)$p(\boldsymbol{\phi}_k | \beta)$のパラメータ$\beta > 0$を指定する。それぞれ、次ステップの変分事後分布のパラメータの初期値として使用する。
トピック分布の変分事後分布のパラメータの初期値を作成する。
# トピック分布θの変分事後分布のパラメータαの初期化 alpha_k <- runif(n = K, min = 0.01, max = 2) # 範囲を指定 alpha_k
## [1] 1.4205099 0.1397011 0.3297111 0.5594082
K個の一様乱数をrunif()で生成し、$\boldsymbol{\theta}$の変分事後分布(ディリクレ分布)$q(\boldsymbol{\theta} | \boldsymbol{\alpha})$のパラメータ$\boldsymbol{\alpha} = (\alpha_1, \cdots, \alpha_K)$、$\alpha_k > 0$とする。
初期値によるトピック分布の期待値(変分事後分布の期待値)を計算する。
# トピック分布の期待値を計算:式(1.15) E_theta_k <- alpha_k / sum(alpha_k) E_theta_k
## [1] 0.57995848 0.05703645 0.13461276 0.22839231
ディリクレ分布の期待値は、次の式で計算できる。
$\alpha_k > 0$を総和で割っているので、$0 \leq \mathbb{E}[\theta_k] \leq 1$、$\sum_{k=1}^K \mathbb{E}[\theta_k] = 1$となり、カテゴリ分布のパラメータの条件を満たす。
作図用に、生成した値をデータフレームに格納する。
# 作図用のデータフレームを作成 E_topic_df <- tibble::tibble( topic = factor(1:K), # トピック番号 probability = E_theta_k # 割り当て確率 ) E_topic_df
## # A tibble: 4 × 2 ## topic probability ## <fct> <dbl> ## 1 1 0.580 ## 2 2 0.0570 ## 3 3 0.135 ## 4 4 0.228
因子型の1からKの整数をトピック番号列として、対応する確率とあわせて格納する。
初期値によるトピック分布の期待値をグラフで確認する。
# 推定トピック分布を作図 ggplot() + geom_bar(data = E_topic_df, mapping = aes(x = topic, y = probability, fill = topic), stat = "identity", show.legend = FALSE) + # トピック分布の期待値 labs(title = "Topic Distribution", subtitle = "initial", x = "k", y = expression(E(theta[k])))
x軸をトピック番号、y軸を各トピックとなる確率として棒グラフを描画する。
同様に、単語分布の変分事後分布のパラメータの初期値を作成する。
# 単語分布Φの変分事後分布のパラメータβの初期化 beta_kv <- runif(n = K*V, min = 0.01, max = 2) |> # 範囲を指定 matrix(nrow = K, ncol = V) beta_kv[, 1:5]
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.3131834 1.8139104 0.3900829 1.1481236 0.5579850 ## [2,] 1.8351323 1.6963435 1.4327709 0.1698565 0.2275119 ## [3,] 0.1153063 0.1751875 1.4162313 1.4793283 1.3684136 ## [4,] 0.2835802 0.6060316 1.9112660 1.9239297 0.3576422
K×V個の一様乱数を生成してK行V列のマトリクスに整形し、K個の$\boldsymbol{\phi}_k$の変分事後分布(ディリクレ分布)$q(\boldsymbol{\phi}_k | \boldsymbol{\beta}_k)$のパラメータ$\boldsymbol{\beta} = \{\boldsymbol{\beta}_1, \cdots, \boldsymbol{\beta}_K\}$、$\boldsymbol{\beta}_k = (\beta_{k1}, \cdots, \beta_{kV})$、$\beta_{kv} > 0$とする。
初期値による単語分布の期待値(変分事後分布の期待値)を計算する。
# 単語分布の期待値を計算:式(1.15) E_phi_kv <- beta_kv / rowSums(beta_kv) E_phi_kv[, 1:5]
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.012235897 0.070868457 0.01524032 0.044856542 0.021800158 ## [2,] 0.060771948 0.056175841 0.04744741 0.005624941 0.007534247 ## [3,] 0.004485943 0.006815595 0.05509789 0.057552654 0.053237562 ## [4,] 0.010936031 0.023371098 0.07370636 0.074194726 0.013792171
トピック分布のときと同様にして、トピック$k$(行)ごとに$\mathbb{E}[\phi_{kv}] = \frac{\beta_{kv}}{\sum_{v=1}^V \beta_{kv'}}$を計算する。
生成した値をデータフレームに格納する。
# 作図用のデータフレームを作成 E_word_df <- E_phi_kv |> tibble::as_tibble(.name_repair = ~as.character(1:V)) |> tibble::add_column( topic = factor(1:K) # トピック番号 ) |> tidyr::pivot_longer( cols = !topic, names_to = "vocabulary", # 列名を語彙番号に変換 names_ptypes = list(vocabulary = factor()), values_to = "probability" # 出現確率列をまとめる ) E_word_df
## # A tibble: 104 × 3 ## topic vocabulary probability ## <fct> <fct> <dbl> ## 1 1 1 0.0122 ## 2 1 2 0.0709 ## 3 1 3 0.0152 ## 4 1 4 0.0449 ## 5 1 5 0.0218 ## 6 1 6 0.0557 ## 7 1 7 0.00382 ## 8 1 8 0.0579 ## 9 1 9 0.0110 ## 10 1 10 0.0159 ## # … with 94 more rows
phi_true_kvはK行V列のマトリクスなので、as_tibble()でデータフレームに変換して、トピック番号列を追加する。列名を1からVの整数にしておく。
語彙ごと(V列)の出現確率列をpivot_longer()で1列にまとめる。names_ptypes引数で列名を因子型に変換して、語彙番号列とする。
初期値による単語分布の期待値をグラフで確認する。
# 推定単語分布を作図 ggplot() + geom_bar(data = E_word_df, mapping = aes(x = vocabulary, y = probability, fill = vocabulary), stat = "identity", show.legend = FALSE) + # 単語分布の期待値 facet_wrap(topic ~ ., labeller = label_bquote(k==.(topic)), scales = "free_x") + # 分割 labs(title = "Word Distribution", subtitle = "initial", x = "v", y = expression(E(phi[kv])))
x軸を語彙番号、y軸を各語彙となる確率として棒グラフを描画する。
facet_wrap()にトピック番号列を指定して、トピックごとにグラフを分割して描画する。
最尤推定では、トピック分布$\boldsymbol{\theta}$と単語分布$\boldsymbol{\Phi}$をそれぞれ点推定した。変分ベイズ推定では、ハイパーパラメータ$\boldsymbol{\alpha}, \boldsymbol{\beta}$を求めることで、$\boldsymbol{\theta}, \boldsymbol{\Phi}$を分布推定する。上の2つのグラフは、それぞれ推定した分布の期待値を可視化したものである。
ここまでで、推論処理の準備ができた。続いて、推論処理を行う。
コード全体
試行回数を指定して、トピック分布と単語分布の更新を繰り返す。
# 試行回数を指定 MaxIter <- 20 # 推移の確認用の受け皿を作成 trace_q_dki <- array(NA, dim = c(D, K, MaxIter)) trace_alpha_ki <- matrix(NA, nrow = K, ncol = MaxIter+1) trace_beta_kvi <- array(NA, dim = c(K, V, MaxIter+1)) # 初期値を保存 trace_alpha_ki[, 1] <- alpha_k trace_beta_kvi[, , 1] <- beta_kv # 変分ベイズ推定 for(i in 1:MaxIter) { # 負担率qの計算:式(3.22) term_alpha_k <- digamma(alpha_k) - digamma(sum(alpha_k)) term_beta_kd <- digamma(beta_kv) %*% t(N_dv) - digamma(rowSums(beta_kv)) %*% t(N_d) log_q_dk <- t(term_alpha_k + term_beta_kd) log_q_dk <- log_q_dk - apply(log_q_dk, MARGIN = 1, FUN = min) # アンダーフロー対策 log_q_dk <- log_q_dk - apply(log_q_dk, MARGIN = 1, FUN = max) # オーバーフロー対策 q_dk <- exp(log_q_dk) #q_dk[rowSums(q_dk) == 0, ] <- 1 / K # 全ての要素が0の場合は等確率(任意の同じ値)に設定 q_dk <- q_dk / rowSums(q_dk) # 正規化 # トピック分布θの変分事後分布のパラメータαの計算:式(3.19') alpha_k <- alpha + colSums(q_dk) # 単語分布Φの変分事後分布のパラメータβの計算:式(3.20') beta_kv <- beta + t(q_dk) %*% N_dv # i回目の更新値を保存 trace_q_dki[, , i] <- q_dk trace_alpha_ki[, i+1] <- alpha_k trace_beta_kvi[, , i+1] <- beta_kv # 動作確認 message("\r", i, "/", MaxIter, appendLF = FALSE) }
推移の確認用に、負担率とトピック分布、単語分布の初期値と更新値をそれぞれtrace_***に格納していく。
コードの解説
繰り返し行う更新処理をそれぞれ確認する。
負担率の計算
先に、前ステップで更新した(または初期値の)トピック分布$\boldsymbol{\theta}$のハイパーパラメータ$\boldsymbol{\alpha}$と単語分布$\boldsymbol{\Phi}$のハイパーパラメータ$\boldsymbol{\beta}$を用いて、負担率$\mathbf{q}$を更新する。各文書の負担率$\mathbf{q}_d = (q_{d1}, \cdots, q_{dK})$をまとめて、$\mathbf{q} = \{\mathbf{q}_1, \cdots, \mathbf{q}_D\}$とする。
次の処理によって、負担率を更新(計算)する。
# 負担率qの計算:式(3.22') term_alpha_k <- digamma(alpha_k) - digamma(sum(alpha_k)) term_beta_kd <- digamma(beta_kv) %*% t(N_dv) - digamma(rowSums(beta_kv)) %*% t(N_d) log_q_dk <- t(term_alpha_k + term_beta_kd) log_q_dk[1:5, ]
## [,1] [,2] [,3] [,4] ## [1,] -379.9799 -624.4040 -628.3267 -552.6445 ## [2,] -297.5354 -446.5964 -421.2790 -412.6085 ## [3,] -294.5431 -464.4651 -480.4457 -436.4662 ## [4,] -268.7081 -412.2376 -431.8243 -391.6616 ## [5,] -380.6398 -604.2793 -570.5513 -572.5941
負担率の各要素の更新式は、次の式である。
ここで、$\Psi(x)$はディガンマ関数で、digamma()で計算できる。
正規化前の値を$q_{dk} \propto \tilde{q}_{dk}$とおき、対数をとった(指数をとっていない)$\log \tilde{q}_{dk}$を計算する。
3つ目の項$\sum_{v=1}^V N_{dv}\ (d = 1, \dots, D,\ k = 1, \dots, K)$について、$(D \times K)$の行列とすると、2つの行列の積に分解できる。
文書と語彙ごとの度数$N_{dv}\ (d = 1, \dots, D,\ v = 1, \dots, V)$を並べた$(D \times V)$の行列と、$\Psi(\beta_{kv})\ (k = 1, \dots, K,\ v = 1, \dots, V)$を並べた$(V \times K)$の行列の積で計算できる。後の項は、$\boldsymbol{\beta}$を$(K \times V)$の行列として、転置して各要素を$\Psi(x)$で計算した行列と言える。
また、4つ目の項$N_d \Psi(\sum_{v=1}^V \beta_{kv})$についても、$(D \times K)$の行列とすると、2つのベクトルの積に分解できる。
文書ごとの度数$N_d\ (d = 1, \dots, D)$を並べた$D$次元の縦ベクトルと、$\Psi(\sum_{v=1}^V \beta_{kv})\ (k = 1, \dots, K)$を並べた$K$次元の横ベクトルの積で計算できる。
よって、$\log \tilde{q}_{dk}\ (d = 1, \dots, D,\ k = 1, \dots, K)$について、$(D \times K)$の行列とすると、次の式で計算できる。
1つ目の項は、$\Psi(\alpha_k)\ (k = 1, \dots, K)$をD行に複製した行列で、2つ目の項は、$\Psi(\sum_{k=1}^K \alpha_k)$をD×K個に複製した行列である。
各要素の指数を計算する。
# 負担率qの計算:式(3.22') log_q_dk <- log_q_dk - apply(log_q_dk, MARGIN = 1, FUN = min) # アンダーフロー対策 log_q_dk <- log_q_dk - apply(log_q_dk, MARGIN = 1, FUN = max) # オーバーフロー対策 q_dk <- exp(log_q_dk) q_dk[1:5, ]
## [,1] [,2] [,3] [,4] ## [1,] 1 7.046356e-107 1.394302e-108 1.029719e-75 ## [2,] 1 1.834969e-65 1.814912e-54 1.057802e-50 ## [3,] 1 1.598773e-74 1.834506e-81 2.309821e-62 ## [4,] 1 4.633825e-63 1.443895e-71 3.999483e-54 ## [5,] 1 7.492395e-98 3.330435e-83 4.318486e-84
$x = \exp(\log x)$より、各要素の$\log$を外す。
その際に、行ごとに最小値を引くことでアンダーフロー、最大値を引くことでオーバーフローが起きにくくする。アンダーフロー対策については「ソフトマックス関数のオーバーフロー対策【ゼロつく1のノート(数学)】 - からっぽのしょこ」を参照のこと。
各行を和で割って正規化する。
# 負担率qの計算:式(3.22) #q_dk[rowSums(q_dk) == 0, ] <- 1 / K # 全ての要素が0の場合は等確率(任意の同じ値)に設定 q_dk <- q_dk / rowSums(q_dk) # 正規化 q_dk[1:5, ]
## [,1] [,2] [,3] [,4] ## [1,] 1 7.046356e-107 1.394302e-108 1.029719e-75 ## [2,] 1 1.834969e-65 1.814912e-54 1.057802e-50 ## [3,] 1 1.598773e-74 1.834506e-81 2.309821e-62 ## [4,] 1 4.633825e-63 1.443895e-71 3.999483e-54 ## [5,] 1 7.492395e-98 3.330435e-83 4.318486e-84
各要素をトピックごとの和で割って正規化する。
q_dkの行ごとの和が$\sum_k$の計算に対応する。
アンダーフローによって、$\tilde{q}_{dk}\ (k = 1, \dots, K)$の全ての値が0になると、正規化の計算において0を0で割ることになり、NaNになる。そこで、q_dk[d, ]の要素が全て0(行の和が0)のときは、$q_{dk} = \frac{1}{K}$とする(全ての値を0や、正規化したランダムな値の方がいいのかもしれない)。この例では、明示的に1 / Kを代入しているが、K個の要素を同じ値にすれば正規化時に1 / Kになる。
トピック分布と単語分布の計算
次に、更新した負担率$\mathbf{q}$を用いて、トピック分布$\boldsymbol{\theta}$のハイパーパラメータ$\boldsymbol{\alpha}$と単語分布$\boldsymbol{\Phi}$のハイパーパラメータ$\boldsymbol{\beta}$を更新する。
次の処理によって、トピック分布のハイパーパラメータを更新する。
# トピック分布θの変分事後分布のパラメータαの計算:式(3.19') res_alpha_k <- alpha + colSums(q_dk) res_alpha_k
## [1] 12 2 4 6
トピック分布の変分事後分布のパラメータの各要素の更新式は、次の式である。
$\alpha_k\ (1, \dots, K)$について、$K$次元のベクトルとすると、2つのベクトルの和に分解できる。
q_dkの列ごとの和が$\sum_d$の計算に対応する。
(次にやる作図処理に影響しないように変数名をres_***にしている。)
続いて、次の処理によって、単語分布のハイパーパラメータを更新する。
# 単語分布Φの変分事後分布のパラメータβの計算:式(3.20') res_beta_kv <- beta + t(q_dk) %*% N_dv res_beta_kv[, 1:5]
## [,1] [,2] [,3] [,4] [,5] ## [1,] 8 5 33 21 23 ## [2,] 1 3 9 7 10 ## [3,] 8 4 31 15 5 ## [4,] 9 2 6 23 70
単語分布の変分事後分布のパラメータの各要素の更新式は、次の式である。
$\beta_{kv}\ (k = 1, \dots, K,\ v = 1, \dots, V)$について、$(K \times V)$の行列とすると、3つの行列の積と和に分解できる。
負担率$q_{dk}\ (d = 1, \dots, D,\ k = 1, \dots, K)$を並べた$(K \times D)$の行列と、度数$N_{dv}\ (d = 1, \dots, D,\ v = 1, \dots, V)$を並べた$(D \times V)$の行列の積に、事前分布のパラメータ$\beta$を加えて計算できる。2つ目の項は、$\mathbf{q}$を$(D \times K)$の行列として、転置した行列と言える。
以上が、1ステップでの更新処理である。この処理を指定した回数繰り返し行う。
推論結果の可視化
推定したトピック分布と単語分布を可視化する。詳しくは「ループあり版」の記事を参照のこと。
トピック分布の期待値(変分事後分布の期待値)を計算して、推定したトピック分布を真の分布と重ねてグラフを作成する。
# トピック分布の期待値を計算:式(1.15) E_theta_k <- alpha_k / sum(alpha_k) E_theta_k
## [1] 0.50000000 0.08333333 0.16666667 0.25000000
推定した分布を塗りつぶし、真の分布を破線で示す。
単語分布の期待値(変分事後分布の期待値)を計算して、推定した単語分布を真の分布と重ねてグラフを作成する。
# 単語分布の期待値を計算:式(1.15) E_phi_kv <- beta_kv / rowSums(beta_kv) E_phi_kv[, 1:5]
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.005372733 0.003357958 0.02216253 0.01410343 0.01544661 ## [2,] 0.004484305 0.013452915 0.04035874 0.03139013 0.04484305 ## [3,] 0.021505376 0.010752688 0.08333333 0.04032258 0.01344086 ## [4,] 0.013975155 0.003105590 0.00931677 0.03571429 0.10869565
トピック分布と単語分布のアニメーションを作成する。
この記事では、変分ベイズ推定を実装した。次の記事では、ギブスサンプリングを実装する。
参考書籍
- 岩田具治(2015)『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社
おわりに
初めて本を読んで実装してた時にアレコレ悩んでたことが、数行でできてしまい感動です。ついでに、なんでNaNの対策もできました。
行列恐怖症だった自分が、行列計算の形にした方が分かりやすくなるとは、辛かったけど特訓してよかった辛かったけど。記事中の計算式はゴチャゴチャしていますが、行列表記用の記号をそれぞれ用意したらスッキリすると思います。4章まで修正できたら共通して使える表記を考えてみます。
【次節の内容】
つづく
