確率的主成分分析で、各パラメータを期待値最大化法(EMアルゴリズム)により求める式の導出の仕方を備忘録として書いてみる。確率的主成分分析では、観測値\(\boldsymbol{x}\)が、潜在変数\(\boldsymbol{z}\)によって、確率的に生成されたと考えるらしい。\(\boldsymbol{\mu}\)は\(\boldsymbol{x}\)の平均値を表している(\(\boldsymbol{z} = \boldsymbol{0}\)のときに\(\boldsymbol{\mu}\)に移る)。観測値\(\boldsymbol{x}\)を\(M\)次元、潜在変数\(\boldsymbol{z}\)を\(m\)次元とする。\(\boldsymbol{z}\)の事前分布\(p(\boldsymbol{z})\)は、取り敢えず正規分布としておく。\(p(\boldsymbol{x} | \boldsymbol{z})\)のパラメータ\(\boldsymbol{W}, \boldsymbol{\mu}, \sigma^2\)を、観測データ\(\mathcal{D} = \left\{ \boldsymbol{x}_1, \cdots, \boldsymbol{x}_N \right\}\)から求めるのが目標となる。
\[p(\boldsymbol{x} | \boldsymbol{z}) = \mathcal{N}(\boldsymbol{x} | \boldsymbol{W} \boldsymbol{z} + \boldsymbol{\mu}, \sigma^2 \boldsymbol{I}_M)\] \[p(\boldsymbol{z}) = \mathcal{N}(\boldsymbol{z} | \boldsymbol{0}, \boldsymbol{I}_m)\]
平均ベクトル\(\boldsymbol{\mu}\)は簡単に求められる。多変数正規分布のベイズの定理を用いると、\(p(\boldsymbol{x})\)は次のようになる。
\[p(\boldsymbol{x}) = \mathcal{N}(\boldsymbol{x} | \boldsymbol{\mu}, \boldsymbol{C}) \qquad (\boldsymbol{C} = \boldsymbol{W} \boldsymbol{W}^T + \sigma^2 \boldsymbol{I}_M)\]
最尤推定を行うために対数尤度関数\(L(\boldsymbol{\mu}, \boldsymbol{W}, \sigma^2 | \mathcal{D})\)を求める。
\[\begin{align} L(\boldsymbol{\mu}, \boldsymbol{W}, \sigma^2 | \mathcal{D}) &= \ln \prod_{i = 1}^N p(\boldsymbol{x}_i) \\ &= \ln \prod_{i = 1}^N \mathcal{N}(\boldsymbol{x}_i | \boldsymbol{\mu}, \boldsymbol{C}) \\ &= \sum_{i = 1}^N \ln \mathcal{N}(\boldsymbol{x}_i | \boldsymbol{\mu}, \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 - \boldsymbol{\mu})^T \boldsymbol{C}^{-1} (\boldsymbol{x}_i - \boldsymbol{\mu}) \right) \right) \\ &= -\frac{MN}{2} \ln(2\pi) - \frac{N}{2} \ln|\boldsymbol{C}| - \frac{1}{2} \sum_{i = 1}^N (\boldsymbol{x}_i - \boldsymbol{\mu})^T \boldsymbol{C}^{-1} (\boldsymbol{x}_i - \boldsymbol{\mu}) \\ \end{align}\]
\(L(\boldsymbol{\mu}, \boldsymbol{W}, \sigma^2 | \mathcal{D})\)を\(\boldsymbol{\mu}\)で偏微分して\(\boldsymbol{0}\)とおくことによって、\(\boldsymbol{\mu}\)の最尤推定値\(\hat{\boldsymbol{\mu}}\)が得られる。
\[\begin{align} \frac{\partial L}{\partial \boldsymbol{\mu}} &= -\frac{1}{2} \sum_{i = 1}^N \frac{\partial}{\partial \boldsymbol{\mu}} (\boldsymbol{x}_i - \boldsymbol{\mu})^T C^{-1} (\boldsymbol{x}_i - \boldsymbol{\mu}) \\ &= -\frac{1}{2} \sum_{i = 1}^N \left( (-1) (C^{-1} + (C^{-1})^T)(\boldsymbol{x}_i - \boldsymbol{\mu}) \right) \\ &= -\frac{1}{2} \sum_{i = 1}^N \left( -2 C^{-1} (\boldsymbol{x}_i - \boldsymbol{\mu}) \right) \\ &= C^{-1} \sum_{i = 1}^N (\boldsymbol{x}_i - \boldsymbol{\mu}) \\ &= \boldsymbol{0} \end{align}\]
\[\sum_{i = 1}^N \boldsymbol{x}_i = \sum_{i = 1}^N \boldsymbol{\mu} = N \boldsymbol{\mu}\] \[\boldsymbol{\mu} = \frac{1}{N} \sum_{i = 1}^N \boldsymbol{x}_i\]
\(\hat{\boldsymbol{\mu}}\)は単なる標本平均となっている。これ以降\(\boldsymbol{\mu}\)は既知として扱うことにして、\(\boldsymbol{\mu} = \bar{\boldsymbol{x}}\)と書く。パラメータ\(\boldsymbol{W}, \sigma^2\)を求めるために、対数尤度関数を改めて書き直す。
\[\begin{align} L(\boldsymbol{W}, \sigma^2 | \mathcal{D}) &= \ln \prod_{i = 1}^N p(\boldsymbol{x}_i) \\ &= \sum_{i = 1}^N \ln p(\boldsymbol{x}_i) \\ &= \sum_{i = 1}^N \ln \int p(\boldsymbol{x}_i | \boldsymbol{z}_i) p(\boldsymbol{z}_i) d\boldsymbol{z}_i \end{align}\]
ここでイェンセンの不等式を用いると次のようになる(\(q(\boldsymbol{z}_i)\)は何らかの確率分布)。
\[\begin{align} L(\boldsymbol{W}, \sigma^2 | \boldsymbol{x}_i) &= \ln \int p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) p(\boldsymbol{z}_i) d\boldsymbol{z}_i \\ &\ge \int q(\boldsymbol{z}_i) \ln \frac{p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) p(\boldsymbol{z}_i)}{q(\boldsymbol{z}_i)} d\boldsymbol{z}_i \\ &= \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 \\ \end{align}\]
\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2)\)とすると上手くいく。
\[\begin{align} \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 p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2) \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 \\ &= \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{\frac{p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) p(\boldsymbol{z}_i)}{p(\boldsymbol{x}_i)}} d\boldsymbol{z}_i \\ &= \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2) \ln \frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{\frac{p(\boldsymbol{x}_i, \boldsymbol{z}_i | \boldsymbol{W}, \sigma^2)}{p(\boldsymbol{x}_i)}} d\boldsymbol{z}_i \\ &= \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2) \ln \left( p(\boldsymbol{x}_i) \right) d\boldsymbol{z}_i \\ &= \ln \left( p(\boldsymbol{x}_i) \right) \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i \\ &\qquad (\because \lnの中身が\boldsymbol{z}_iに依存しないため, 係数として積分の外側に括り出す) \\ &= \ln \left( \int p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) p(\boldsymbol{z}_i) d \boldsymbol{z}_i \right) \int p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2) d\boldsymbol{z}_i \\ &\qquad (\because p(\boldsymbol{x}_i)を積分の形で表現する) \\ &= \ln \int p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) p(\boldsymbol{z}_i) d \boldsymbol{z}_i \\ &\qquad (\because 確率分布p(\boldsymbol{z}_i | \boldsymbol{x}_i)に関する規格化条件) \\ &= L(\boldsymbol{W}, \sigma^2 | \boldsymbol{x}_i) \end{align}\]
上式のように、対数関数の中身が積分変数(上の場合は\(\boldsymbol{z}_i\))に関係なくなるときに、イェンセンの不等式による近似の精度が最も良くなる(元々の対数尤度と等しくなっている)。\(q(\boldsymbol{z}_i)\)は、多変数正規分布のベイズの定理から次のようになる。
\[q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}, \sigma^2) = \mathcal{N}(\boldsymbol{z}_i | \boldsymbol{M}^{-1} \boldsymbol{W}^T (\boldsymbol{x} - \bar{\boldsymbol{x}}), \sigma^2 \boldsymbol{M}^{-1})\] \[\boldsymbol{M} = \boldsymbol{W}^T \boldsymbol{W} + \sigma^2 \boldsymbol{I}_m\]
\(q(\boldsymbol{z}_i)\)は上式に、パラメータ\(\boldsymbol{W}, \sigma^2\)を代入すれば求められる。代入されるパラメータは、\(q(\boldsymbol{z}_i)\)を求める前に既に用意してあった仮のものであるため、これらを\(\boldsymbol{W}^{\prime}, \sigma^{2^{\prime}}\)と表記する\((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{Z} = \left\{ \boldsymbol{z}_1, \cdots, \boldsymbol{z}_N \right\})\)。なぜこの式の期待値を最大化すればよいかについてはまた後で考える。
\[\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) p(\boldsymbol{z}_i) \\ &= \sum_{i = 1}^N (\ln p(\boldsymbol{x}_i | \boldsymbol{z}_i) + \ln p(\boldsymbol{z}_i)) \\ \end{align}\]
\(q(\boldsymbol{z}_i) = p(\boldsymbol{z}_i | \boldsymbol{x}_i, \boldsymbol{W}^{\prime}, \sigma^{2^{\prime}})\)(潜在変数の事後分布)についての期待値を取ると、上式は次のようになる。
\[\begin{align} & \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \left\{ \ln p(\boldsymbol{x}_i | \boldsymbol{z}_i, \boldsymbol{W}, \sigma^2) + \ln p(\boldsymbol{z}_i) \right\} d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) \left\{ \ln \mathcal{N}(\boldsymbol{x}_i | \boldsymbol{W} \boldsymbol{z}_i + \bar{\boldsymbol{x}}, \sigma^2 \boldsymbol{I}_M) + \ln \mathcal{N}(\boldsymbol{z}_i | \boldsymbol{0}, \boldsymbol{I}_m) \right\} d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) (-\frac{M}{2} \ln(2\pi) - \frac{1}{2} \ln|\sigma^2 \boldsymbol{I}_M| - \frac{1}{2} (\boldsymbol{x}_i - (\boldsymbol{W} \boldsymbol{z}_i + \bar{\boldsymbol{x}}))^T (\sigma^2 \boldsymbol{I}_m)^{-1} (\boldsymbol{x}_i - (\boldsymbol{W} \boldsymbol{z}_i + \bar{\boldsymbol{x}})) \\ & \qquad - \frac{m}{2} \ln(2\pi) - \frac{1}{2} \boldsymbol{z}_i^T \boldsymbol{z}_i) d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) (-\frac{M}{2} \ln(2\pi) - \frac{1}{2} \ln((\sigma^2)^M) - \frac{1}{2 \sigma^2} (\boldsymbol{x}_i - (\boldsymbol{W} \boldsymbol{z}_i + \bar{\boldsymbol{x}}))^T (\boldsymbol{x}_i - (\boldsymbol{W} \boldsymbol{z}_i + \bar{\boldsymbol{x}})) \\ & \qquad - \frac{m}{2} \ln(2\pi) - \frac{1}{2} \boldsymbol{z}_i^T \boldsymbol{z}_i) d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) (-\frac{M}{2} \ln(2\pi) - \frac{M}{2} \ln(\sigma^2) \\ & \qquad - \frac{1}{2 \sigma^2} (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) + \frac{1}{\sigma^2} (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \boldsymbol{z}_i - \frac{1}{2 \sigma^2} (\boldsymbol{W} \boldsymbol{z}_i)^T (\boldsymbol{W} \boldsymbol{z}_i) \\ & \qquad - \frac{m}{2} \ln(2\pi) - \frac{1}{2} \boldsymbol{z}_i^T \boldsymbol{z}_i) d\boldsymbol{z}_i \\ &= \sum_{i = 1}^N \int q(\boldsymbol{z}_i) (-\frac{M}{2} \ln(2 \pi \sigma^2) - \frac{m}{2} \ln(2\pi) \\ & \qquad - \frac{1}{2 \sigma^2} ||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 + \frac{1}{\sigma^2} (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \boldsymbol{z}_i - \frac{1}{2 \sigma^2} \mathrm{Tr}((\boldsymbol{W} \boldsymbol{z}_i)^T (\boldsymbol{W} \boldsymbol{z}_i)) - \frac{1}{2} \mathrm{Tr}(\boldsymbol{z}_i^T \boldsymbol{z}_i)) d\boldsymbol{z}_i \\ &= -\frac{MN}{2} \ln(2 \pi \sigma^2) - \frac{mN}{2} \ln(2\pi) - \frac{1}{2 \sigma^2} \sum_{i = 1}^N \int q(\boldsymbol{z}_i) (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \boldsymbol{z}_i + \mathrm{Tr}(\boldsymbol{z}_i^T \boldsymbol{W}^T \boldsymbol{W} \boldsymbol{z}_i) + \sigma^2 \mathrm{Tr}(\boldsymbol{z}_i \boldsymbol{z}_i^T)) d\boldsymbol{z}_i \\ &= -\frac{MN}{2} \ln(2 \pi \sigma^2) - \frac{mN}{2} \ln(2\pi) - \frac{1}{2 \sigma^2} \sum_{i = 1}^N \int q(\boldsymbol{z}_i) (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \boldsymbol{z}_i + \mathrm{Tr}(\boldsymbol{z}_i \boldsymbol{z}_i^T \boldsymbol{W}^T \boldsymbol{W}) + \sigma^2 \mathrm{Tr}(\boldsymbol{z}_i \boldsymbol{z}_i^T)) d\boldsymbol{z}_i \end{align}\]
\(q(\boldsymbol{z}_i)\)による期待値を\(\langle \cdot \rangle\)で表すことにすると、上式は次のようになる。
\[\begin{align} & -\frac{MN}{2} \ln(2 \pi \sigma^2) - \frac{mN}{2} \ln(2\pi) - \frac{1}{2 \sigma^2} \sum_{i = 1}^N (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}) + \sigma^2 \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle)) \end{align}\]
ここで\(\langle \boldsymbol{z}_i \rangle\)と、\(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle\)は、更新前の仮のパラメータ\(\boldsymbol{W}^{\prime}, \sigma^{2^{\prime}}\)を使って次のように計算できることが分かる。以下の2つの式が期待値ステップで利用される。
\[\langle \boldsymbol{z}_i \rangle = \boldsymbol{M}^{\prime^{-1}} W^{\prime^{T}} (\boldsymbol{x}_i - \bar{\boldsymbol{x}})\] \[\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle = \sigma^{2^\prime} \boldsymbol{M}^{\prime^{-1}} + \langle \boldsymbol{z}_i \rangle \langle \boldsymbol{z}_i \rangle^T\] \[\boldsymbol{M}^{\prime} = \boldsymbol{W}^{\prime^T} \boldsymbol{W}^{\prime} + \sigma^{2^\prime} \boldsymbol{I}_m\]
最大化ステップで使用する式(パラメータ\(\boldsymbol{W}, \sigma^2\)の更新式)は、上述の対数尤度の近似式を、\(\boldsymbol{W}, \sigma^2\)について微分して\(0\)とおくことによって得られる。
\[\begin{align} & \frac{\partial}{\partial \boldsymbol{W}} \{ -\frac{MN}{2} \ln(2 \pi \sigma^2) - \frac{mN}{2} \ln(2\pi) - \frac{1}{2 \sigma^2} \sum_{i = 1}^N (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}) + \sigma^2 \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle))\} \\ &= -\frac{1}{2 \sigma^2} \frac{\partial}{\partial \boldsymbol{W}} \sum_{i = 1}^N (-2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W})) \\ &= -\frac{1}{2 \sigma^2} \frac{\partial}{\partial \boldsymbol{W}} \sum_{i = 1}^N (-2 \mathrm{Tr}((\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle) + \mathrm{Tr}(\boldsymbol{W} \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T)) \\ &= -\frac{1}{2 \sigma^2} \frac{\partial}{\partial \boldsymbol{W}} \sum_{i = 1}^N (-2 \mathrm{Tr}(\boldsymbol{W} \langle \boldsymbol{z}_i \rangle (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T) + \mathrm{Tr}(\boldsymbol{W} \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T)) \\ &= -\frac{1}{2 \sigma^2} \sum_{i = 1}^N (-2 (\langle \boldsymbol{z}_i \rangle (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T)^T + \boldsymbol{W} (\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle + \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle^T)) \\ &= -\frac{1}{2 \sigma^2} \sum_{i = 1}^N (-2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) \langle \boldsymbol{z}_i \rangle^T + 2 \boldsymbol{W} \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle) \\ &= \boldsymbol{0} \end{align}\]
これより\(\boldsymbol{W}\)の更新式は次のようになる。
\[\begin{align} \therefore \sum_{i = 1}^N \boldsymbol{W} \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle &= \sum_{i = 1}^N (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) \langle \boldsymbol{z}_i \rangle^T \\ \boldsymbol{W} \left( \sum_{i = 1}^N \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \right) &= \sum_{i = 1}^N (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) \langle \boldsymbol{z}_i \rangle^T \\ \boldsymbol{W} &= \left\{ \sum_{i = 1}^N (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) \langle \boldsymbol{z}_i \rangle^T \right\} \left( \sum_{i = 1}^N \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \right)^{-1} \end{align}\]
\(\sigma^2\)についても同様に行う。
\[\begin{align} & \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 (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}) + \sigma^2 \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle))\} \\ &= \frac{\partial}{\partial \sigma^2} \{ -\frac{MN}{2} \ln(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_{i = 1}^N (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}))\} \\ &= -\frac{MN}{2} \frac{2\pi}{2 \pi \sigma^2} + \frac{1}{2 \sigma^4} \sum_{i = 1}^N (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}))\} \\ &= -\frac{MN}{2} \frac{1}{\sigma^2} + \frac{1}{2 \sigma^4} \sum_{i = 1}^N (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}))\} \\ &= \boldsymbol{0} \end{align}\]
これより\(\sigma^2\)の更新式は次のようになる。
\[\begin{align} \therefore \frac{MN}{2} \frac{1}{\sigma^2} &= \frac{1}{2 \sigma^4} \sum_{i = 1}^N (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}))\} \\ \sigma^2 &= \frac{1}{MN} \sum_{i = 1}^N (||\boldsymbol{x}_i - \bar{\boldsymbol{x}}||^2 \\ & \qquad - 2 (\boldsymbol{x}_i - \bar{\boldsymbol{x}})^T \boldsymbol{W} \langle \boldsymbol{z}_i \rangle + \mathrm{Tr}(\langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle \boldsymbol{W}^T \boldsymbol{W}))\} \end{align}\]
従って、\(\boldsymbol{W}, \sigma^2\)の初期値を適当に決めてから、両者のパラメータが収束するまで、\(\langle \boldsymbol{z}_i \rangle, \langle \boldsymbol{z}_i \boldsymbol{z}_i^T \rangle\)の計算と、\(\boldsymbol{W}, \sigma^2\)の更新を交互に繰り返していけばよいことが分かる。
あまり理解できていないな…