確率的主成分分析とEMアルゴリズム(その1)

確率的主成分分析で、各パラメータを期待値最大化法(EMアルゴリズム)により求める式の導出の仕方を備忘録として書いてみる。確率的主成分分析では、観測値$x$が、潜在変数$z$によって、確率的に生成されたと考えるらしい。$\mu$は$x$の平均値を表している($z=0$のときに$\mu$に移る)。観測値$x$を$M$次元、潜在変数$z$を$m$次元とする。$z$の事前分布$p\left(z\right)$は、取り敢えず正規分布としておく。$p\left(x|z\right)$のパラメータ$W,\mu ,{\sigma }^2$を、観測データ$D=\left\{x_1,\cdots ,x_N\right\}$から求めるのが目標となる。

$$p\left(x|z\right)=N(x|Wz+\mu ,{\sigma }^2{I}_M)$$ $$p\left(z\right)=N(z|0,I_m)$$

平均ベクトル$\mu$は簡単に求められる。多変数正規分布のベイズの定理を用いると、$p\left(x\right)$は次のようになる。

$$p\left(x\right)=N(x|\mu ,C)(C=W{W}^T+{\sigma }^2{I}_M)$$

最尤推定を行うために対数尤度関数$L(\mu ,W,{\sigma }^2|D)$を求める。

$$\begin{array}{cc}L(\mu ,W,{\sigma }^2|D)& =\ln \prod _{i=1}^{N}p\left(x_i\right)\\ & =\ln \prod _{i=1}^{N}N(x_i|\mu ,C)\\ & =\sum _{i=1}^N\ln N(x_i|\mu ,C)\\ & =\sum _{i=1}^N\ln \left((2\pi {)}^{-\frac{M}{2}}|C{|}^{-\frac{1}{2}}\exp \left(-\frac{1}{2}(x_i-\mu {)}^T{C}^{-1}(x_i-\mu )\right)\right)\\ & =-\frac{MN}{2}\ln \left(2\pi \right)-\frac{N}{2}\ln |C|-\frac{1}{2}\sum _{i=1}^N(x_i-\mu {)}^T{C}^{-1}(x_i-\mu )\end{array}$$

$L(\mu ,W,{\sigma }^2|D)$を$\mu$で偏微分して$0$とおくことによって、$\mu$の最尤推定値$\hat{\mu }$が得られる。

$$\begin{array}{cc}\frac{\partial L}{\partial \mu }& =-\frac{1}{2}\sum _{i=1}^N\frac{\partial }{\partial \mu }(x_i-\mu {)}^T{C}^{-1}(x_i-\mu )\\ & =-\frac{1}{2}\sum _{i=1}^N\left((-1)(C^{-1}+(C^{-1}{)}^T\left)\right(x_i-\mu )\right)\\ & =-\frac{1}{2}\sum _{i=1}^N\left(-2C^{-1}(x_i-\mu )\right)\\ & =C^{-1}\sum _{i=1}^N(x_i-\mu )\\ & =0\end{array}$$

$$\sum _{i=1}^N{x}_i=\sum _{i=1}^N\mu =N\mu$$ $$\mu =\frac{1}{N}\sum _{i=1}^N{x}_i$$

$\hat{\mu }$は単なる標本平均となっている。これ以降$\mu$は既知として扱うことにして、$\mu =\overline{x}$と書く。パラメータ$W,{\sigma }^2$を求めるために、対数尤度関数を改めて書き直す。

$$\begin{array}{cc}L(W,{\sigma }^2|D)& =\ln \prod _{i=1}^{N}p\left(x_i\right)\\ & =\sum _{i=1}^N\ln p\left(x_i\right)\\ & =\sum _{i=1}^N\ln \int p\left(x_i|z_i\right)p\left(z_i\right)d{z}_i\end{array}$$

ここでイェンセンの不等式を用いると次のようになる($q\left(z_i\right)$は何らかの確率分布)。

$$\begin{array}{cc}L(W,{\sigma }^2|x_i)& =\ln \int p(x_i|z_i,W,{\sigma }^2)p\left(z_i\right)d{z}_i\\ & \ge \int q\left(z_i\right)\ln \frac{p(x_i|z_i,W,{\sigma }^2)p\left(z_i\right)}{q\left(z_i\right)}d{z}_i\\ & =\int q\left(z_i\right)\ln \frac{p(x_i,z_i|W,{\sigma }^2)}{q\left(z_i\right)}d{z}_i\end{array}$$

$q\left(z_i\right)=p(z_i|x_i,W,{\sigma }^2)$とすると上手くいく。

$$\begin{array}{cc}\int q\left(z_i\right)\ln \frac{p(x_i,z_i,W,{\sigma }^2)}{q\left(z_i\right)}d{z}_i& =\int p(z_i|x_i,W,{\sigma }^2)\ln \frac{p(x_i,z_i|W,{\sigma }^2)}{p(z_i|x_i,W,{\sigma }^2)}d{z}_i\\ & =\int p(z_i|x_i,W,{\sigma }^2)\ln \frac{p(x_i,z_i|W,{\sigma }^2)}{\frac{p(x_i|z_i,W,{\sigma }^2)p\left(z_i\right)}{p\left(x_i\right)}}d{z}_i\\ & =\int p(z_i|x_i,W,{\sigma }^2)\ln \frac{p(x_i,z_i|W,{\sigma }^2)}{\frac{p(x_i,z_i|W,{\sigma }^2)}{p\left(x_i\right)}}d{z}_i\\ & =\int p(z_i|x_i,W,{\sigma }^2)\ln \left(p\left(x_i\right)\right)d{z}_i\\ & =\ln \left(p\left(x_i\right)\right)\int p(z_i|x_i,W,{\sigma }^2)d{z}_i\\ & (\because \ln の中身がz_iに依存しないため,係数として積分の外側に括り出す)\\ & =\ln \left(\int p(x_i|z_i,W,{\sigma }^2)p\left(z_i\right)d{z}_i\right)\int p(z_i|x_i,W,{\sigma }^2)d{z}_i\\ & (\because p(x_i\left)を積分の形で表現する\right)\\ & =\ln \int p(x_i|z_i,W,{\sigma }^2)p\left(z_i\right)d{z}_i\\ & (\because 確率分布p(z_i|x_i\left)に関する規格化条件\right)\\ & =L(W,{\sigma }^2|x_i)\end{array}$$

上式のように、対数関数の中身が積分変数(上の場合は$z_i$)に関係なくなるときに、イェンセンの不等式による近似の精度が最も良くなる(元々の対数尤度と等しくなっている)。$q\left(z_i\right)$は、多変数正規分布のベイズの定理から次のようになる。

$$q\left(z_i\right)=p(z_i|x_i,W,{\sigma }^2)=N\left(z_i|M^{-1}{W}^T\right(x-\overline{x}),{\sigma }^2{M}^{-1})$$ $$M=W^{T}W+{\sigma }^2{I}_m$$

$q\left(z_i\right)$は上式に、パラメータ$W,{\sigma }^2$を代入すれば求められる。代入されるパラメータは、$q\left(z_i\right)$を求める前に既に用意してあった仮のものであるため、これらを$W^{\prime },{\sigma }^{2^{\prime }}$と表記する$\left(q\right(z_i)=p(z_i|x_i,W^{\prime },{\sigma }^{2^{\prime }}\left)\right)$。さて、最大化しようとしている対数尤度$L(W,{\sigma }^2|D,Z)$は次のようになる$(Z=\left\{z_1,\cdots ,z_N\right\})$。なぜこの式の期待値を最大化すればよいかについてはまた後で考える。

$$\begin{array}{cc}L(W,{\sigma }^2|D,Z)& =\ln p(D,Z|W,{\sigma }^2)\\ & =\ln \prod _{i=1}^{N}p\left(x_i|z_i\right)p\left(z_i\right)\\ & =\sum _{i=1}^N\left(\ln p\right(x_i|z_i)+\ln p(z_i\left)\right)\end{array}$$

$q\left(z_i\right)=p(z_i|x_i,W^{\prime },{\sigma }^{2^{\prime }})$(潜在変数の事後分布)についての期待値を取ると、上式は次のようになる。

$$\begin{array}{cc}& \sum _{i=1}^N\int q\left(z_i\right)\left\{\ln p(x_i|z_i,W,{\sigma }^2)+\ln p\left(z_i\right)\right\}d{z}_i\\ & =\sum _{i=1}^N\int q\left(z_i\right)\left\{\ln N(x_i|W{z}_i+\overline{x},{\sigma }^2{I}_M)+\ln N(z_i|0,I_m)\right\}d{z}_i\\ & =\sum _{i=1}^N\int q\left(z_i\right)(-\frac{M}{2}\ln (2\pi )-\frac{1}{2}\ln |{\sigma }^2{I}_M|-\frac{1}{2}(x_i-(W{z}_i+\overline{x}){)}^T\left({\sigma }^2{I}_m{)}^{-1}\right(x_i-(W{z}_i+\overline{x}))\\ & -\frac{m}{2}\ln \left(2\pi \right)-\frac{1}{2}{z}_i^T{z}_i)d{z}_i\\ & =\sum _{i=1}^N\int q\left(z_i\right)(-\frac{M}{2}\ln (2\pi )-\frac{1}{2}\ln (\left({\sigma }^2{)}^M\right)-\frac{1}{2{\sigma }^2}(x_i-(W{z}_i+\overline{x}\left){)}^T\right(x_i-(W{z}_i+\overline{x}))\\ & -\frac{m}{2}\ln \left(2\pi \right)-\frac{1}{2}{z}_i^T{z}_i)d{z}_i\\ & =\sum _{i=1}^N\int q\left(z_i\right)(-\frac{M}{2}\ln (2\pi )-\frac{M}{2}\ln ({\sigma }^2)\\ & -\frac{1}{2{\sigma }^2}(x_i-\overline{x}{)}^T(x_i-\overline{x})+\frac{1}{{\sigma }^2}(x_i-\overline{x}{)}^{T}W{z}_i-\frac{1}{2{\sigma }^2}\left(W{z}_i{)}^T\right(W{z}_i)\\ & -\frac{m}{2}\ln \left(2\pi \right)-\frac{1}{2}{z}_i^T{z}_i)d{z}_i\\ & =\sum _{i=1}^N\int q\left(z_i\right)(-\frac{M}{2}\ln (2\pi {\sigma }^2)-\frac{m}{2}\ln (2\pi )\\ & -\frac{1}{2{\sigma }^2}||x_i-\overline{x}|{|}^2+\frac{1}{{\sigma }^2}(x_i-\overline{x}{)}^{T}W{z}_i-\frac{1}{2{\sigma }^2}Tr(\left(W{z}_i{)}^T\right(W{z}_i\left)\right)-\frac{1}{2}Tr\left(z_i^T{z}_i\right))d{z}_i\\ & =-\frac{MN}{2}\ln \left(2\pi {\sigma }^2\right)-\frac{mN}{2}\ln \left(2\pi \right)-\frac{1}{2{\sigma }^2}\sum _{i=1}^N\int q\left(z_i\right)(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W{z}_i+Tr(z_i^T{W}^{T}W{z}_i)+{\sigma }^{2}Tr(z_i{z}_i^T\left)\right)d{z}_i\\ & =-\frac{MN}{2}\ln \left(2\pi {\sigma }^2\right)-\frac{mN}{2}\ln \left(2\pi \right)-\frac{1}{2{\sigma }^2}\sum _{i=1}^N\int q\left(z_i\right)(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W{z}_i+Tr(z_i{z}_i^T{W}^{T}W)+{\sigma }^{2}Tr(z_i{z}_i^T\left)\right)d{z}_i\end{array}$$

$q\left(z_i\right)$による期待値を$⟨\cdot ⟩$で表すことにすると、上式は次のようになる。

$$\begin{array}{cc}& -\frac{MN}{2}\ln \left(2\pi {\sigma }^2\right)-\frac{mN}{2}\ln \left(2\pi \right)-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W)+{\sigma }^{2}Tr(⟨z_i{z}_i^T⟩\left)\right)\end{array}$$

ここで$⟨z_i⟩$と、$⟨z_i{z}_i^T⟩$は、更新前の仮のパラメータ$W^{\prime },{\sigma }^{2^{\prime }}$を使って次のように計算できることが分かる。以下の2つの式が期待値ステップで利用される。

$$⟨z_i⟩=M^{{\prime }^{-1}}{W}^{{\prime }^T}(x_i-\overline{x})$$ $$⟨z_i{z}_i^T⟩={\sigma }^{2^{\prime }}{M}^{{\prime }^{-1}}+⟨z_i⟩⟨z_i{⟩}^T$$ $$M^{\prime }=W^{{\prime }^T}{W}^{\prime }+{\sigma }^{2^{\prime }}{I}_m$$

最大化ステップで使用する式(パラメータ$W,{\sigma }^2$の更新式)は、上述の対数尤度の近似式を、$W,{\sigma }^2$について微分して$0$とおくことによって得られる。

$$\begin{array}{cc}& \frac{\partial }{\partial W}\{-\frac{MN}{2}\ln (2\pi {\sigma }^2)-\frac{mN}{2}\ln (2\pi )-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W)+{\sigma }^{2}Tr(⟨z_i{z}_i^T⟩\left)\right)\}\\ & =-\frac{1}{2{\sigma }^2}\frac{\partial }{\partial W}\sum _{i=1}^N(-2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W))\\ & =-\frac{1}{2{\sigma }^2}\frac{\partial }{\partial W}\sum _{i=1}^N(-2Tr((x_i-\overline{x}{)}^{T}W⟨z_i⟩)+Tr(W⟨z_i{z}_i^T⟩W^T))\\ & =-\frac{1}{2{\sigma }^2}\frac{\partial }{\partial W}\sum _{i=1}^N(-2Tr(W⟨z_i⟩(x_i-\overline{x}{)}^T)+Tr(W⟨z_i{z}_i^T⟩W^T))\\ & =-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(-2(⟨z_i⟩(x_i-\overline{x}{)}^T{)}^T+W(⟨z_i{z}_i^T⟩+⟨z_i{z}_i^T{⟩}^T\left)\right)\\ & =-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(-2(x_i-\overline{x})⟨z_i{⟩}^T+2W⟨z_i{z}_i^T⟩)\\ & =0\end{array}$$

これより$W$の更新式は次のようになる。

$$\begin{array}{cc}\therefore \sum _{i=1}^{N}W⟨z_i{z}_i^T⟩& =\sum _{i=1}^N(x_i-\overline{x})⟨z_i{⟩}^T\\ W\left(\sum _{i=1}^N⟨z_i{z}_i^T⟩\right)& =\sum _{i=1}^N(x_i-\overline{x})⟨z_i{⟩}^T\\ W& =\left\{\sum _{i=1}^N(x_i-\overline{x})⟨z_i{⟩}^T\right\}{\left(\sum _{i=1}^N⟨z_i{z}_i^T⟩\right)}^{-1}\end{array}$$

${\sigma }^2$についても同様に行う。

$$\begin{array}{cc}& \frac{\partial }{\partial {\sigma }^2}\{-\frac{MN}{2}\ln (2\pi {\sigma }^2)-\frac{mN}{2}\ln (2\pi )-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W)+{\sigma }^{2}Tr(⟨z_i{z}_i^T⟩\left)\right)\}\\ & =\frac{\partial }{\partial {\sigma }^2}\{-\frac{MN}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W\left)\right)\}\\ & =-\frac{MN}{2}\frac{2\pi }{2\pi {\sigma }^2}+\frac{1}{2{\sigma }^4}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W\left)\right)\}\\ & =-\frac{MN}{2}\frac{1}{{\sigma }^2}+\frac{1}{2{\sigma }^4}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W\left)\right)\}\\ & =0\end{array}$$

これより${\sigma }^2$の更新式は次のようになる。

$$\begin{array}{cc}\therefore \frac{MN}{2}\frac{1}{{\sigma }^2}& =\frac{1}{2{\sigma }^4}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W\left)\right)\}\\ {\sigma }^2& =\frac{1}{MN}\sum _{i=1}^N(||x_i-\overline{x}|{|}^2\\ & -2(x_i-\overline{x}{)}^{T}W⟨z_i⟩+Tr(⟨z_i{z}_i^T⟩W^{T}W\left)\right)\}\end{array}$$

従って、$W,{\sigma }^2$の初期値を適当に決めてから、両者のパラメータが収束するまで、$⟨z_i⟩,⟨z_i{z}_i^T⟩$の計算と、$W,{\sigma }^2$の更新を交互に繰り返していけばよいことが分かる。

あまり理解できていないな…