確率的主成分分析とEMアルゴリズム(その1)の続きで、しっくりこない部分を再度考えてみる。平均ベクトル\(\boldsymbol{\mu}\)は、最尤推定で簡単に求めることができた。対数尤度関数\(L(\boldsymbol{\mu}, \boldsymbol{W}, \sigma^2 | \mathcal{D})\)を\(\boldsymbol{\mu}\)について微分し、\(0\)と等値するのみであった。しかし、未知のパラメータ\(\boldsymbol{W}, \sigma^2\)を最尤推定によって求めるのは難しい。\(L(\boldsymbol{W}, \sigma^2 | \mathcal{D})\)は以下のように、\(\boldsymbol{W}, \sigma^2\)についての複雑な式になるためである。\(\boldsymbol{W}, \sigma^2\)が行列\(\boldsymbol{C} = \boldsymbol{W} \boldsymbol{W}^T + \sigma^2 \boldsymbol{I}_M\)の中に埋め込まれているので、\(L(\boldsymbol{W}, \sigma^2 | \mathcal{D})\)を\(\boldsymbol{W}, \sigma^2\)で微分して\(0\)とおき、解を直接求めるのは大変そうである。
\[\begin{align} L(\boldsymbol{W}, \sigma^2 | \mathcal{D}) &= \ln p(\mathcal{D} | \boldsymbol{W}, \sigma^2) \\ &= \ln \prod_{i = 1}^N p(\boldsymbol{x}_i | \boldsymbol{W}, \sigma^2) \\ &= \sum_{i = 1}^N \ln p(\boldsymbol{x}_i | \boldsymbol{W}, \sigma^2) \\ &= \sum_{i = 1}^N \ln \mathcal{N}(\boldsymbol{x}_i | \boldsymbol{C}) \\ &= \sum_{i = 1}^N \ln \left( (2\pi)^{-\frac{M}{2}} |\boldsymbol{C}|^{-\frac{1}{2}} \exp \left( -\frac{1}{2} (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{C}^{-1} (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) \right) \right) \\ &= -\frac{MN}{2} \ln(2\pi) - \frac{N}{2} \ln|\boldsymbol{C}| - \frac{1}{2} \sum_{i = 1}^N (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{C}^{-1} (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) \\ \end{align}\]
平均ベクトル\(\boldsymbol{\mu}\)の最尤推定量は、観測値\(\boldsymbol{x}\)についての単なる標本平均であるため、\(\boldsymbol{\mu}\)は既知として扱い、\(\boldsymbol{\mu} = \bar{\boldsymbol{x}}\)と書いている。
\(L(\boldsymbol{W}, \sigma^2 | \mathcal{D})\)を直接最大化するのは難しい。そこで次のように、潜在変数\(\boldsymbol{z}\)についての分布\(q(\boldsymbol{z})\)を導入し、2つの項に分解してみることにする。
\[\begin{align} L(\boldsymbol{W}, \sigma^2 | \mathcal{D}) &= \sum_{i = 1}^N \ln p(\boldsymbol{x}_i | \boldsymbol{W}, \sigma^2) \\ &= \sum_{i = 1}^N \underbrace{\left( \int q(\boldsymbol{z}_i) d\boldsymbol{z}_i \right)}_{= 1} \ln p(\boldsymbol{x}_i | \boldsymbol{W}, \sigma^2) \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \ln p(\boldsymbol{x}_i | \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)} d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} \frac{q(\boldsymbol{z}_i)}{p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)} d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \left( \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} - \ln \frac{p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} \right) d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \left\{ \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} d\boldsymbol{z}_i - \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} d\boldsymbol{z}_i \right\} \\ &= \sum_{i = 1}^N \left\{ \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} d\boldsymbol{z}_i + \mathrm{KL}(q_i || p_i) \right\} \\ &= \sum_{i = 1}^N \left\{ \mathcal{L}(q_i, \boldsymbol{W}, \sigma^2) + \mathrm{KL}(q_i || p_i) \right\} \end{align}\]
上式では、\(q_i = q(\boldsymbol{z}_i), p_i = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)\)と略した。\(\mathrm{KL}(q_i || p_i)\)は、\(q(\boldsymbol{z}_i)\)と事後分布\(p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)\)間の、カルバック-ライブラーダイバージェンス(KLダイバージェンス)と呼ばれる量であり、\(q_i\)と\(p_i\)がどれだけ分布として違っているのかを表している。\(\mathrm{KL}(q_i || p_i) \ge 0\)であり、\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)\)のときに限り等号が成立する。従って、以下の不等式が得られる(前回はイェンセンの不等式から導出した)。
\[L(\boldsymbol{W}, \sigma^2 | \mathcal{D}) \ge \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} d\boldsymbol{z}_i = \sum_{i = 1}^N \mathcal{L}(q_i, \boldsymbol{W}, \sigma^2)\]
さて、現在のパラメータが\(\boldsymbol{W}^{\prime}, \sigma^{2^\prime}\)であるとする。このとき、\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime})\)とすれば、上式において等号が成立することは前回確認した。対数の中身は\(p(\boldsymbol{x}_i | \boldsymbol{W}^{\prime}, \sigma^{2^\prime})\)のようになり、積分変数\(\boldsymbol{z}_i\)に依存しなくなるためだった。\(\mathrm{KL}(q_i || p_i) = 0\)となることからも、等号の成立が明らかである。これは、\(\boldsymbol{W}^{\prime}, \sigma^{2^\prime}\)の2つのパラメータを固定したまま、\(q_i = q(\boldsymbol{z}_i)\)をうまく調整することによって、\(\mathcal{L}(q_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime})\)を最大化したことに相当する。
現在のパラメータ\(\boldsymbol{W}^{\prime}, \sigma^{2^\prime}\)を元に、\(q(\boldsymbol{z}_i)\)を更新する方法が分かった。あとは、現在の\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime})\)を元に、パラメータ\(\boldsymbol{W}, \sigma^2\)を更新すればよい。これは、\(q_i = q(\boldsymbol{z}_i)\)を固定しつつ、\(\boldsymbol{W}, \sigma^2\)について\(\mathcal{L}(q_i, \boldsymbol{W}, \sigma^2)\)を最大化するだけである。対数尤度\(L(\boldsymbol{W}, \sigma^2 | \mathcal{D})\)を直接最大化するのが困難であっても、このような2段階のステップによって\(\mathcal{L}(q_i, \boldsymbol{W}, \sigma^2)\)を徐々に大きくしていけば、\(\mathcal{L}(q_i, \boldsymbol{W}, \sigma^2)\)の合計(上の不等式の右辺)が\(L(\boldsymbol{W}, \sigma^2 | \mathcal{D})\)の最大値(上の不等式の左辺の最大値)に近づいていくので、上手く行きそうである。\(\mathcal{L}(q_i, \boldsymbol{W}, \sigma^2)\)は次のようになる。
\[\begin{align} \mathcal{L}(q_i, \boldsymbol{W}, \sigma^2) &= \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{q(\boldsymbol{z}_i)} d\boldsymbol{z}_i \\ &= \int q(\boldsymbol{z}_i) \ln p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i - \int q(\boldsymbol{z}_i) \ln q(\boldsymbol{z}_i) d\boldsymbol{z}_i \end{align}\]
\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime})\)を代入すると以下のようになる。
\[\begin{align} & \int q(\boldsymbol{z}_i) \ln p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i - \int q(\boldsymbol{z}_i) \ln q(\boldsymbol{z}_i) d\boldsymbol{z}_i \\ &= \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime}) \ln p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i - \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime}) \ln p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime}) d\boldsymbol{z}_i \\ &= \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime}) \ln p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i + \mathrm{Const.} \\ &= \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime}) \left\{ \ln p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) + \ln p(\boldsymbol{z}_i) \right\} d\boldsymbol{z}_i + \mathrm{Const.} \end{align}\]
上式の第2項は、\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime})\)のエントロピーであり、\(\boldsymbol{W}, \sigma^2\)には関係しない定数として扱うことができる。上式を\(\boldsymbol{W}, \sigma^2\)について最大化するために\(\boldsymbol{W}, \sigma^2\)で微分しても、第2項には\(\boldsymbol{W}, \sigma^2\)が含まれないため単に無視される。結局、\(\boldsymbol{W}, \sigma^2\)によって最大化されるのは上式の第1項のみだと分かる。従って、第1項を\(\boldsymbol{W}, \sigma^2\)で微分して\(0\)と等値すれば、更新後の新たなパラメータ\(\boldsymbol{W}, \sigma^2\)が求まる。因みに前回の例では、\(p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2)\)は\(\mathcal{N}(\boldsymbol{x}_i | \boldsymbol{W} \boldsymbol{z}_i + \bar{\boldsymbol{x}}, \sigma^2 \boldsymbol{I}_M)\)、また\(p(\boldsymbol{z}_i)\)は\(\mathcal{N}(\boldsymbol{z}_i | \boldsymbol{0}, \boldsymbol{I}_m)\)であった。
上式は\(i\)番目の\(\mathcal{L}(q_i, \boldsymbol{W}, \sigma^2)\)ついてのみ考えたものだったので、\(i = 1, \cdots, N\)について足し合わせて、以下のように変形してみる。定数項を除いた部分が、\(\boldsymbol{W}, \sigma^2\)について最大化する対象となる。
\[\begin{align} \sum_{i = 1}^N \mathcal{L}(q_i, \boldsymbol{W}, \sigma^2) &= \sum_{i = 1}^N \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^\prime}) \ln p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i + \mathrm{Const.} \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \left\{ \ln p(\boldsymbol{z}_i | \boldsymbol{x}_i , \boldsymbol{W}, \sigma^2) + \ln p(\boldsymbol{z}_i) \right\} d\boldsymbol{z}_i + \mathrm{Const.} \end{align}\]
ここで、
\[\begin{align} L(\boldsymbol{W}, \sigma^2 | \mathcal{D}, \boldsymbol{Z}) &= \ln p(\mathcal{D}, \boldsymbol{Z} | \boldsymbol{W}, \sigma^2) \\ &= \ln \prod_{i = 1}^N p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) p(\boldsymbol{z}_i) \\ &= \sum_{i = 1}^N (\ln p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) + \ln p(\boldsymbol{z}_i)) \\ \end{align}\]
であるから、結局\(\boldsymbol{W}, \sigma^2\)によって最大化していたのは、対数尤度\(L(\boldsymbol{W}, \sigma^2 | \mathcal{D}, \boldsymbol{Z})\)の、潜在変数の事後分布\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^{\prime}})\)についての期待値だということが分かった。
これより、対数尤度\(L(\boldsymbol{W}, \sigma^2 | \mathcal{D}, \boldsymbol{Z})\)を最初に用意し、現在のパラメータ\(\boldsymbol{W}^{\prime}, \sigma^{2^\prime}\)により計算された潜在変数の事後分布\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^{\prime}})\)における期待値を求めた後に、この期待対数尤度を最大化することによって、新たなパラメータ\(\boldsymbol{W}, \sigma^2\)が得られることが分かった。
これで前回は曖昧になってしまった部分がやっと分かったので一安心した。