はじめに
機械学習プロフェッショナルシリーズの『トピックモデル』の勉強時に自分の理解の助けになったことや勉強会資料のまとめです。トピックモデルの各種アルゴリズムを「数式」と「プログラム」から理解することを目指します。
この記事は、3.3節「EMアルゴリズム」の内容です。
EMアルゴリズを用いた最尤推定による混合ユニグラムモデルにおけるパラメータ推定を導出します。
この記事で用いる記号類や混合ユニグラムモデルの定義については前節の記事をご確認ください。
【プログラム編】
www.anarchive-beta.com
www.anarchive-beta.com
【前節の内容】
www.anarchive-beta.com
【他の節一覧】
www.anarchive-beta.com
【この節の内容】
3.3 EMアルゴリズム
ユニグラムモデルでは、対数尤度が最大となるパラメータを解析的に求めることができた。しかし混合モデルでは解析的に得ることができない。そこでEMアルゴリズムを用いて最尤推定を行う。
EMアルゴリズムとは、E(期待値)ステップとM(最大化)ステップを収束するまで繰り返すことでパラメータを推定する手法のことである。この節では対数尤度の下限を最大化することで、パラメータの局所最適値を求める。
3.1節の混合ユニグラムモデルの文書集合の生成過程(生成モデル)$p(\mathbf{W} | \boldsymbol{\theta}, \boldsymbol{\Phi})$より、混合ユニグラムモデルの対数尤度は次のように変形できる。
$$
\begin{align}
\log p(\mathbf{W} | \boldsymbol{\theta}, \boldsymbol{\Phi})
&= \sum_{d=1}^D
\log p(\mathbf{w}_d | \boldsymbol{\theta}, \boldsymbol{\Phi})
\\
&= \sum_{d=1}^D
\log \sum_{k=1}^K
p(\mathbf{w}_d, z_d = k | \boldsymbol{\theta}, \boldsymbol{\phi}_k)
\\
&= \sum_{d=1}^D
\log \sum_{k=1}^K
p(z_d = k | \boldsymbol{\theta})
p(\mathbf{w}_d | \boldsymbol{\phi}_k)
\tag{3.1}
\end{align}
$$
この対数尤度$\log p(\mathbf{W} | \boldsymbol{\theta}, \boldsymbol{\Phi})$に、負担率を導入する。負担率とは各文書が各トピックに属する確率のことである。文書$d$がトピック$k$となる確率を$q_{dk}$とし、次の定義に従う。
$$
0 \leq q_{dk} \leq 1,\
\sum_{k=1}^K q_{dk} = 1
$$
対数尤度(3.1)の式は対数$\log$の計算の中に総和$\sum$の計算を含むことから、解析的に計算できない。そこでイェンゼンの不等式を用いて、対数尤度の下限を求める。
$$
\begin{align}
\log p(\mathbf{W} | \boldsymbol{\theta}, \boldsymbol{\Phi})
&= \sum_{d=1}^D
\log \sum_{k=1}^K
p(z_d = k | \boldsymbol{\theta})
p(\mathbf{w}_d | \boldsymbol{\phi}_k)
\tag{3.1}\\
&= \sum_{d=1}^D \log \sum_{k=1}^K
q_{dk}
\frac{
p(z_d = k | \boldsymbol{\theta})
p(\mathbf{w}_d | \boldsymbol{\phi}_k)
}{
q_{dk}
}
\\
&\geq
\sum_{d=1}^D \sum_{k=1}^K
q_{dk}
\log \frac{
p(z_d = k | \boldsymbol{\theta})
p(\mathbf{w}_d | \boldsymbol{\phi}_k)
}{
q_{dk}
}
\\
&= \sum_{d=1}^D \sum_{k=1}^K
q_{dk}
\log \frac{
\theta_k
\prod_{v=1}^V \phi_{kv}^{N_{dv}}
}{
q_{dk}
}
\\
&= \sum_{d=1}^D \sum_{k=1}^K
q_{dk} \left(
\log \theta_k
+ \sum_{v=1}^V N_{dv} \log \phi_{kv}
- \log q_{dk}
\right)
\equiv F
\tag{3.2}
\end{align}
$$
【途中式の途中式】
- 負担率$q_{dk}$を$\frac{q_{dk}}{q_{dk}} = 1$として、分母分子を分割して式(3.1)に掛ける。
- 対数関数$y = \log x$は上に凸な関数であることから、イェンゼンの不等式(1.7)を用いて下限みる。
- 文書集合の生成過程(3.1節)より、それぞれ具体的な式に置き換える
- 各文書に対してトピック$k$が割り当てられる確率は、$p(z_d = k | \boldsymbol{\theta}) = \theta_k$である。
- トピック$k$の単語分布に従う単語集合の生成確率は、$p(\mathbf{w}_d | \boldsymbol{\phi}_k) = \prod_{v=1}^V \phi_{kv}^{N_{dv}}$である。
- $\log \frac{A B}{C} = \log A + \log B - \log C$であり、また$\log x^a = a \log x$であることから式を変形する。
対数尤度の下限を$F$と置く。これにより$\sum$を$\log$の外に出すことができた(解析的に計算できる)。
次からは、下限$F$を最大化するパラメータをEMアルゴリズムにより推定していく。
・Eステップ
Eステップでは、下限$F$が最大となる負担率$(q_{11}, \cdots, q_{DK})$を求める。
・負担率の更新式の導出
$\sum_{k=1}^K q_{dk} = 1$の制約の下で下限$F$が最大となる$q_{dk}$を、ラグランジュの未定乗数法を用いて求める。ラグランジュの未定乗数$\lambda$を用いて式を立て、(ここでは負担率$q_{dk}$に関心があるためこの式を$q_{dk}$の関数として)$F(q_{dk})$と置く。
$$
F(q_{dk})
= q_{dk} \left(
\log \theta_k
+ \sum_{v=1}^V
N_{dv} \log \phi_{kv}
- \log q_{dk}
\right)
+ \lambda \left(
\sum_{k=1}^K q_{dk} - 1
\right)
\tag{3.3}
$$
この式を$q_{dk}$に関して微分して、$\frac{\partial F(q_{dk})}{\partial q_{dk}} = 0$となる$q_{dk}$を求める。
$$
\begin{align}
\frac{\partial F(q_{dk})}{\partial q_{dk}}
= \log \theta_k
&+ \sum_{v=1}^V
N_{dv} \log \phi_{kv}
- \log q_{dk}
- 1 + \lambda
= 0
\\
\log q_{dk}
&= \log \theta_k
+ \sum_{v=1}^V
N_{dv} \log \phi_{kv}
+ \lambda - 1
\\
q_{dk}
&= \exp(\lambda - 1)
\theta_k
\prod_{v=1}^V \phi_{kv}^{N_{dv}}
\tag{3.5}
\end{align}
$$
【途中式の途中式】
- $F(q_{dk})$を$q_{dk}$に関して微分する。
各項の微分について順番に確認する。
1つ目の因子は、括弧の中を1つの係数とみて
$$
\left\{
q_{dk} \left(
\log \theta_k
+ \sum_{v=1}^V
N_{dv} \log \phi_{kv}
\right)
\right\}'
= \log \theta_k
+ \sum_{v=1}^V
N_{dv} \log \phi_{kv}
$$
になる。
2つ目の項は、積の微分$\{f(x)g(x)\}' = f'(x)g(x) + f(x)g'(x)$と対数関数の微分$(\log x)' = \frac{1}{x}$より
$$
\begin{aligned}
(-q_{dk} \log q_{dk})'
&= (-q_{dk})' \log q_{dk}
+ (-q_{dk}) (\log q_{dk})'
\\
&= - 1 \log q_{dk}
- q_{dk} \frac{1}{q_{dk}}
\\
&= - \log q_{dk} - 1
\end{aligned}
$$
になる。
3つ目の因子は、$q_{dk}$以外の$q_{d1}$から$q_{dK}$は定数として扱うため$q_{dk}$で微分すると0になる。
$$
\begin{aligned}
\left\{
\lambda \left(
\sum_{k=1}^K q_{dk} - 1
\right)
\right\}'
&= \Bigl\{
(\lambda q_{d1} - \lambda) + \cdots + (\lambda q_{dk} - \lambda) + \cdots + (\lambda q_{dK} - \lambda)
\Bigr\}'
\\
&= 0 + \cdots + \lambda + \cdots + 0
\\
&= \lambda
\end{aligned}
$$
このように、ある変数に関してのみ行う微分を偏微分と呼ぶ。
- $\log q_{dk}$に関して式を整理する。
- 両辺で$\log$を外す。
この式の両辺で$k = 1$から$K$まで和をとると、制約条件より
$$
\begin{align}
\sum_{k=1}^K q_{dk}
&= \sum_{k=1}^K
\exp(\lambda - 1) \theta_k
\prod_{v=1}^V
\phi_{kv}^{N_{dv}}
\\
1 &= \exp(\lambda - 1)
\sum_{k=1}^K
\theta_k
\prod_{v=1}^V
\phi_{kv}^{N_{dv}}
\\
\exp(\lambda - 1)
&= \frac{
1
}{
\sum_{k=1}^K
\theta_k
\prod_{v=1}^V
\phi_{kv}^{N_{dv}}
}
\tag{3.6}
\end{align}
$$
となる。これを式(3.5)に代入すると
$$
q_{dk}
= \frac{
\theta_k
\prod_{v=1}^V \phi_{kv}^{N_{dv}}
}{
\sum_{k'=1}^K \theta_{k'}
\prod_{v=1}^V \phi_{k'v}^{N_{dv}}
}
\tag{3.3}
$$
が得られる。
この式が負担率$q_{dk}$の更新式となる。またこの式は、文書$d$とトピック$k$の同時尤度$p(\mathbf{w}_d, z_d = k | \theta_k, \boldsymbol{\phi}_k) = \theta_k \prod_{v=1}^V \phi_{kv}^{N_dv}$を$k$についての和が1となるように正規化したものになっている。
次は、更新した負担率を用いて(固定して)他のパラメータを更新する。
・Mステップ
Mステップでは、負担率が与えられた下での下限$F$が最大化するようにパラメータ$\boldsymbol{\theta},\ \boldsymbol{\Phi}$を更新する。
・トピック分布の更新式の導出
Eステップと同様に、$\sum_{k=1}^K \theta_k = 1$の制約の下で式(3.2)が最大となるパラメータ$\boldsymbol{\theta}$を、ラグランジュの未定乗数法を用いて求める。ここでは$\theta_k$に注目するため、式を$F(\theta_k)$とする。
$$
\begin{aligned}
F(\theta_k)
= \sum_{d=1}^D
q_{dk} \left(
\log \theta_k
+ \sum_{v=1}^V N_{dv} \log \phi_{kv} - \log q_{dk}
\right)
+ \lambda \left(
\sum_{k=1}^K \theta_k - 1
\right)
\end{aligned}
$$
この式を$\theta_k$に関して微分して、$\frac{\partial F(\theta_k)}{\partial \theta_k} = 0$となる$\theta_k$を求める。
$$
\begin{align}
\frac{\partial F(\theta_k)}{\partial \theta_k}
= \frac{
\sum_{d=1}^D q_{dk}
}{
\theta_k
}
+ \lambda
&= 0
\\
\theta_k
&= - \frac{
\sum_{d=1}^D q_{dk}
}{
\lambda
}
\tag{1}
\end{align}
$$
【途中式の途中式】
$\log \theta_k$を$g(x)$とすると、合成関数の微分$\{f(g(x))\}' = f'(g(x)) g'(x)$と対数関数の微分$(\log x)' = \frac{1}{x}$より
$$
\left(
\sum_{d=1}^D
q_{dk} \log \theta_k
\right)'
= \sum_{d=1}^D
q_{dk}
\frac{1}{\theta_k}
$$
になる。
この式の両辺で$k = 1$から$K$まで和をとると、$\sum_{k=1}^K \theta_k = 1$より
$$
\begin{aligned}
\sum_{k=1}^K \theta_k
&= - \sum_{k=1}^K \sum_{d=1}^D
\frac{q_{dk}}{\lambda}
\\
1 &= - \sum_{k=1}^K \sum_{d=1}^D
\frac{q_{dk}}{\lambda}
\\
\lambda
&= - \sum_{k=1}^K \sum_{d=1}^D
q_{dk}
\end{aligned}
$$
となる。これを(1)に代入すると
$$
\begin{align}
\theta_k
&= - \frac{
\sum_{d=1}^D q_{dk}
}{
- \sum_{k=1}^K \sum_{d=1}^D
q_{dk}
}
\\
&= \frac{
\sum_{d=1}^D q_{dk}
}{
\sum_{k'=1}^K \sum_{d=1}^D q_{dk'}
}
\tag{3.7}
\end{align}
$$
が得られる。
この式がパラメータ$\theta_k$の更新式となる。この式の分子は、全ての文書に関して和をとった負担率である。また更に全てのトピックに関して和をとったもので割ることで、$\sum_{k=1}^K \theta_k = 1$となるように正規化している。
・単語分布の更新式の導出
$\boldsymbol{\phi}_k$についても同様に、$\sum_{v=1}^V \phi_{kv} = 1$の制約の下で式(3.2)が最大となるパラメータ$\boldsymbol{\phi}_k$を、ラグランジュの未定乗数法を用いて求める。
$$
\begin{aligned}
F(\phi_{kv})
= \sum_{d=1}^D
q_{dk} \left(
\log \theta_k
+ \sum_{v=1}^V N_{dv} \log \phi_{kv} - \log q_{dk}
\right)
+ \lambda \left(
\sum_{k=1}^K \phi_{kv} - 1
\right)
\end{aligned}
$$
$\phi_{kv}$に関して微分して、$\frac{\partial F(\phi_{kv})}{\partial \theta_k} = 0$となる$\theta_k$を求める。
$$
\begin{align}
\frac{\partial F(\phi_{kv})}{\partial \phi_{kv}}
= \frac{
\sum_{d=1}^D q_{dk} N_{dv}
}{
\phi_{kv}
}
+ \lambda
&= 0
\\
\phi_{kv}
&= - \frac{
\sum_{d=1}^D q_{dk} N_{dv}
}{
\lambda
}
\tag{2}
\end{align}
$$
両辺で$v = 1$から$V$まで和をとると、$\sum_{v=1}^V \phi_{kv} = 1$より
$$
\begin{aligned}
\sum_{v=1}^V \phi_{kv}
&= - \sum_{v=1}^V
\frac{
\sum_{d=1}^D q_{dk} N_{dv}
}{
\lambda
}
\\
1 &= - \frac{
\sum_{v=1}^V \sum_{d=1}^D
q_{dk} N_{dv}
}{
\lambda
}
\\
\lambda
&= - \sum_{v=1}^V \sum_{d=1}^D
q_{dk} N_{dv}
\end{aligned}
$$
となる。これを式(2)に代入すると
$$
\begin{align}
\phi_{kv}
&= - \frac{
\sum_{d=1}^D q_{dk} N_{dv}
}{
- \sum_{v=1}^V \sum_{d=1}^D q_{dk} N_{dv}
}
\\
&= \frac{
\sum_{d=1}^D
q_{dk} N_{dv}
}{
\sum_{v=1}^V \sum_{d=1}^D
q_{dk} N_{dv}
}
\tag{3.8}
\end{align}
$$
が得られる。
この式がパラメータ$\phi_{kv}$の更新式となる。この式の分子は、語彙の出現回数によって重み付けした負担率を全ての文書に関して和をとったものである。また更にそれを全ての語彙に関して和をとったもので割ることで、$\sum_{v=1}^V \phi_{kv} = 1$となるように正規化している。
トピック分布も単語分布も(重み付けした)負担率に比例していることが分かる。
次は更新したパラメータ$\boldsymbol{\Phi},\ \boldsymbol{\theta}$を用いて(固定して)、負担率を更新する(Eステップを行う)。EMアルゴリズムは、この操作を繰り返すことでパラメータを推定する。
・対数尤度と下限$F$の関係性
対数尤度$\log p(\mathbf{W} | \boldsymbol{\theta}, \boldsymbol{\Phi})$と下限$F$との差を求めることで、2つの関係をみる。
$$
\begin{align}
\log p(\mathbf{W} | \boldsymbol{\theta}, \boldsymbol{\Phi})
- F
&= \sum_{d=1}^D \sum_{k=1}^K
q_{dk}
\log p(\mathbf{w}_d | \boldsymbol{\theta}, \boldsymbol{\Phi})
- \sum_{d=1}^D \sum_{k=1}^K
q_{dk}
\log \frac{
p(z_d = k | \boldsymbol{\theta})
p(\mathbf{w}_d | z_d = k, \boldsymbol{\phi}_k)
}{
q_{dk}
}
\\
&= \sum_{d=1}^D \sum_{k=1}^K
q_{dk} \left(
\log p(\mathbf{w}_d | \boldsymbol{\theta}, \boldsymbol{\Phi})
- \log \frac{
p(z_d = k | \boldsymbol{\theta})
p(\mathbf{w}_d | z_d = k, \boldsymbol{\phi}_k)
}{
q_{dk}
}
\right)
\\
&= \sum_{d=1}^D \sum_{k=1}^K
q_{dk}
\log \left(
p(\mathbf{w}_d | \boldsymbol{\theta}, \boldsymbol{\Phi})
\frac{
q_{dk}
}{
p(\mathbf{w}_d, z_d = k | \boldsymbol{\theta}, \boldsymbol{\phi}_k)
}
\right)
\\
&= \sum_{d=1}^D \sum_{k=1}^K
q_{dk}
\log \frac{
q_{dk}
}{
p(z_d = k |\mathbf{w}_d, \boldsymbol{\theta}, \boldsymbol{\phi}_k)
}
\\
&= \sum_{d=1}^D
\mathrm{KL}(
\boldsymbol{q}_d, p(z_d |\mathbf{w}_d, \boldsymbol{\theta}, \boldsymbol{\Phi})
)
\tag{3.9}
\end{align}
$$
【途中式の途中式】
- 両式の形を揃えるために、対数尤度$\log p(\mathbf{W} | \boldsymbol{\theta}, \boldsymbol{\Phi}) = \sum_{d=1}^D \log p(\mathbf{w}_d | \boldsymbol{\theta}, \boldsymbol{\Phi})$に$\sum_{k=1}^K q_{dk} = 1$を掛ける。
- $\sum_{d=1}^D \sum_{k=1}^K q_{dk} = 1$で括る。
- $- \log \frac{A}{B} = - \log A + \log B = \log \frac{B}{A}$と$p(A | B) p(B) = P(A, B)$の変形を行い、$\log$の項をまとめる。
- $\frac{p(A, B)}{p(A)} = p(B | A)$の変形を行う。
- 式の形状が(離散確率分布の)KLダイバージェンス(1.6)となっている。
ここで、文書$d$に関する負担率を$\boldsymbol{q}_d = (q_{d1}, q_{d2}, \cdots, q_{dK})$と置く。
対数尤度は負担率の影響を受けないことから、負担率に対して対数尤度は定数となる。つまり下限$F$の最大化は、対数尤度と下限$F$の差の最小化を意味する。
またこの差の式は、文書$d$の負担率$\boldsymbol{q}_d$とトピックの事後分布$p(z_d |\mathbf{w}_d, \boldsymbol{\theta}, \boldsymbol{\Phi})$のKLダイバージェンスであることが分かった。つまり下限$F$の最大化は、このKLタイバージェンスの最小化でもある。
これは、下限$F$を最大化することで、解析的に求められないトピックの事後分布を負担率で近似できるということである。
参考書籍
- 岩田具治(2015)『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社
おわりに
2019/07/26:加筆修正しました。
2020/07/24:加筆修正しました。
【次節の内容】
www.anarchive-beta.com