クラスタリングで省略した、混合ガウス分布のEMアルゴリズムの計算過程を記述します。
EMアルゴリズム
観測データ \(\mathbf{X}=\{\mathbf{x}_1,\ldots,\mathbf{x}_N\}\) に対する対数尤度関数は、\((9.7)\) から以下のように書けます。（※明示的にパラメータを記載しています。）
$$\ln p\left(\mathbf{X}|\boldsymbol{\pi,\mu,\Sigma}\right) = \sum_{n=1}^N\ln\left\{\sum_{k=1}^K\pi_k\mathcal{N}\left(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k\right)\right\}\qquad (9.14)$$
E step
$$\begin{aligned}
p(z_{k}=1|\mathbf{x})
&\equiv \frac{p(z_{k}=1)p(\mathbf{x}|z_{k}=1)}{p(\mathbf{x})}\qquad (\because\text{Bayes' theorem})\\
&\propto p(z_{k}=1)p(\mathbf{x}|z_{k}=1) \\
&=\pi_k\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)
\end{aligned}$$
なので、総和が \(1\) になるように正規化を行えば負担率が以下のように求まります。
$$\gamma(z_{k})\equiv p(z_{k}=1|\mathbf{x})
=\frac{\pi_{k}\mathcal{N}(\mathbf{x}|\mathbf{\mu_{\rm{k}}}, \mathbf{\Sigma_{\rm{k}}})}{\sum_{j=1}^{K}\pi_{j}\mathcal{N}(\mathbf{x}|\mathbf{\mu_{\rm{j}}}, \mathbf{\Sigma_{\rm{j}}})}\qquad (9.13)$$
※ なお、この時分母である \(\sum_{j=1}^{K}\pi_{j}\mathcal{N}(\mathbf{x}|\mathbf{\mu_{\rm{j}}}, \mathbf{\Sigma_{\rm{j}}})\) が \(p(\mathbf{x})\) であることは有用です。（生物データマイニング論 第1回では、これを用いて可視化していました。）
Maximization step
※ 対数尤度関数 \((9.14)\) を \(\ln L\) と記述します。
\(\boldsymbol{\mu}_k\)
対数尤度 \((9.14)\) を \(\boldsymbol{\mu}_k\) で微分します。この時正規分布は \(e^{\text{hoge}}\) という形をしており、微分しても形が変わらないので、
$$\begin{aligned}
\frac{\partial\ln L}{\partial\boldsymbol{\mu}_k}
&= \sum_{i=n}^N \frac{\pi_k\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j)}\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)\\
&= \boldsymbol{\Sigma}_k^{-1}\sum_{n=1}^N \gamma_{nk}(\mathbf{x}_n-\boldsymbol{\mu}_k)\\
&= 0\\
\therefore\boldsymbol{\mu}_k^{\star}
&= \frac{\sum_{n=1}^N\gamma(z_{nk})\mathbf{x}_n}{\sum_{n=1}^N\gamma(z_{nk})} \qquad (9.17)\\
\end{aligned}$$
と最適解が求まりました。
\(\pi_k\)
パラメータ \(\pi_k\) については、\(\sum_{k=1}^K \pi_k = 1\) という制約に注意する必要があります。ですが、隠れマルコフモデルの最尤推定で行なった能登同様に、ラグランジュの未定乗数 \(\lambda\) を導入すればこの問題は解けます。
$$L_{\lambda} = \log L + \lambda\left(\sum_{c=1}^K \pi_c - 1\right)$$
の導関数が \(0\) となる条件を求める事になるので、
$$\begin{aligned}
\frac{\partial L_{\lambda}}{\partial\pi_k}
&= \sum_{n=1}^N \frac{\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\sum_k \pi_k \mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} + \lambda\\
&= \sum_{n=1}^N \frac{\gamma_{nk}}{\pi_k} + \lambda = 0 &\therefore \pi_k = -\frac{1}{\lambda}\sum_{n=1}^N \gamma_{ic}\\
\frac{\partial L_{\lambda}}{\partial\lambda}
&= \sum_{k=1}^K \pi_k - 1\\
&= \left(-\frac{1}{\lambda}\sum_{k=1}^K\sum_{n=1}^N \gamma_{nk}\right) - 1 = 0 &\therefore \lambda = -n
\end{aligned}$$
より、
$$\pi_k^{\star} = \frac{\sum_{n=1}^N\gamma(z_{nk})}{N} \qquad(9.22)$$
と最適解が求まりました。
\(\boldsymbol{\Sigma}_k\)
\(\boldsymbol{\mu}_k\) と同様に、対数尤度 \((9.14)\) を \(\boldsymbol{\Sigma}_k\) で微分します。すると、
$$\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\log L = \sum_{n=1}^N\frac{\pi_k\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j)}$$
となります。ここで、\(\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\) を計算すると、
$$\begin{aligned}
&\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\\
=&\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\frac{1}{(2\pi)^{D/2}}\frac{1}{|\boldsymbol{\Sigma}_k|^{1/2}}\exp\left\{-\frac{1}{2}(\mathbf{x}_n-\boldsymbol{\mu}_k)^T\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)\right\}\\
=&\frac{1}{(2\pi)^{D/2}}\biggl(\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\frac{1}{|\boldsymbol{\Sigma}_k|^{1/2}}\biggr)\exp\left\{-\frac{1}{2}(\mathbf{x}_n-\boldsymbol{\mu}_k)^T\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)\right\}\\
&+\frac{1}{(2\pi)^{D/2}}\frac{1}{|\boldsymbol{\Sigma}_k|^{1/2}}\left(\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\exp\left\{-\frac{1}{2}(\mathbf{x}_n-\boldsymbol{\mu}_k)^T\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)\right\}\right)\\
=&\frac{1}{(2\pi)^{D/2}}\biggl(-\frac{1}{2}\biggr)|\boldsymbol{\Sigma}_k|^{-1/2}\boldsymbol{\Sigma}_k^{-1}\exp\left\{-\frac{1}{2}(\mathbf{x}_n-\boldsymbol{\mu}_k)^T\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)\right\}\quad (\ast1)\\
&+\frac{1}{(2\pi)^{D/2}}\frac{1}{|\boldsymbol{\Sigma}_k|^{1/2}}\frac{1}{2}\exp\left\{-\frac{1}{2}(\mathbf{x}_n-\boldsymbol{\mu}_k)^T\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)\right\}\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)(\mathbf{x}_n-\boldsymbol{\mu}_k)^{T}\boldsymbol{\Sigma}_k^{-1}\quad (\ast2)\\
=&-\frac{1}{2}\boldsymbol{\Sigma}_k^{-1}\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)+\frac{1}{2}\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)(\mathbf{x}_n-\boldsymbol{\mu}_k)^{T}\boldsymbol{\Sigma}_k^{-1}\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\\
=&\frac{1}{2}\left\{\boldsymbol{\Sigma}_k^{-1}-\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)(\mathbf{x}_n-\boldsymbol{\mu}_k)^{T}\boldsymbol{\Sigma}_k^{-1}\right\}\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)
\end{aligned}$$
したがって、
$$
\boldsymbol{\Sigma}_k^{\star}= \frac{\sum_{n=1}^N\gamma(z_{nk})\left(\mathbf{x}_n-\boldsymbol{\mu}_k\right)\left(\mathbf{x}_n-\boldsymbol{\mu}_k\right)^{\mathrm{T}}}{\sum_{n=1}^N\gamma(z_{nk})} \qquad (9.19)
$$
と最適解が求まりました。
※ \(\ast1,\ast2\) の式変形は以下でより詳細に説明します。が、その前にいくつか基本事項の定義を振り返ります。
行列式
$$\det A=\displaystyle\sum_{\sigma\in S_n}\mathrm{sgn}(\sigma)\prod_{i=1}^na_{i\sigma(i)}=\displaystyle\sum_{\sigma\in S_n}\mathrm{sgn}(\sigma)a_{1\sigma(1)}a_{2\sigma(2)}\cdots a_{n\sigma(n)}$$
- \(\sigma\) は \(1\) から \(n\) の置換（順列）を表す。
- \(\mathrm{sgn}(\sigma)\) は置換の符号を表す。なお、置換の符号は互換の数に注目しており、奇置換（互換の数が奇数個）なら \(-1\)、偶置換なら \(+1\)
余韻子
\(n\) 次正方行列 \(A := (a_{ij})\) に対し、\(i\) 行と \(i\) 列を\(1\)つずつ取り去って作られる小行列を \(M_{ij}\) とする。つまり、
$$
M_{ij} = \small \left[
\begin{array}{ccc}
a_{11} & \cdots &a_{1\hspace{1mm}j-1} & a_{1\hspace{1mm}j+1} & \cdots & a_{1n}
\\
\vdots & &\vdots & \vdots & &\vdots
\\
a_{i-1 \hspace{1mm} 1} & \cdots &a_{i-1\hspace{1mm}j-1} & a_{i-1\hspace{1mm}j+1} & \cdots & a_{i-1 \hspace{1mm}n}
\\
a_{i+1 \hspace{1mm} 1} & \cdots &a_{i+1\hspace{1mm}j-1} & a_{i+1\hspace{1mm}j+1} & \cdots & a_{i+1 \hspace{1mm}n}
\\
\vdots & &\vdots & \vdots & &\vdots
\\
a_{n1} &\cdots &a_{n\hspace{1mm}j-1} & a_{n\hspace{1mm}j+1} &\cdots &a_{nn}
\end{array}
\right]
$$
です。ここで、\(\Delta_{ij}=(-1)^{i+j}|M_{ij}|\) とすると、以下の余因子展開ができます。
- \(A\) の行列式は \(j\) 列に関して、以下のように展開されます。
$$\det A=\Delta_{1j}a_{1j}+\Delta_{2j}a_{2j}+\cdots+\Delta_{nj}a_{nj}$$
- \(A\) の行列式は \(i\) 行に関して、以下のように展開されます。
$$\det A=\Delta_{i1}a_{i1}+\Delta_{i2}a_{i2}+\cdots+\Delta_{in}a_{in}$$
余韻子行列
\(n\) 次正方行列 \(A := (a_{ij})\) に対し、\((i, j)\) 余因子を \((j, i)\) 成分に持つ行列
$$\tilde{A} := \begin{pmatrix}&\Delta_{11}&\Delta_{21}&\cdots&\Delta_{n1}\\&\Delta_{21}&\Delta_{22}&\ldots&\Delta_{n2}\\&\vdots&\vdots&\ddots&\vdots\\&\Delta_{n1}&\Delta_{n2}&\cdots&\Delta_{nn}\end{pmatrix}$$
を余韻子行列と呼びます。ここで、余韻子展開を考えれば、余韻子行列に関して、
$$\tilde{A}A=A\tilde{A}=\det(A)I_n$$
が成り立つことがわかります。
\(\ast1\) の変形
以上を踏まえると、\(\frac{\partial\det(A)}{\partial a_{ij}} = \Delta_{ij}\) なので、
$$\frac{\partial|\Sigma|}{\partial\Sigma}
=\tilde\Sigma = |\Sigma|\Sigma^{-1}$$
が成り立ち、
$$\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\frac{1}{|\boldsymbol{\Sigma}_k|^{1/2}}
=-\frac{1}{2}|\boldsymbol{\Sigma}_k|^{-\frac{3}{2}}\frac{\partial}{\partial\boldsymbol{\Sigma}_k}\boldsymbol{\Sigma_k}=-\frac{1}{2}|\boldsymbol{\Sigma}_k|^{-1/2}\boldsymbol{\Sigma_k}^{-1}$$
となることがわかります。
トレースと固有値の関係性
【命題】
トレースには、\(Tr(A) = \sum_{k=1}^{n}\lambda_{k}\) という関係がある。
- トレース： \(n\times n\) の正方行列 \(A\) に対して、対角成分の和 \(\sum_{k=1}^{n}a_{kk}\) を \(A\) のトレースと呼び、\(\mathrm{Tr}(A),\mathrm{tr}A\) と表す。
【証明】
まず、固有方程式は、
$$\phi(t) = |A-tI| = \left|\begin{array}{cccc}
a_{11}-t & a_{12} & \ldots & a_{1n} \\
a_{21} & a_{22}-t & \ldots & a_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{n1} & a_{n2} & \ldots & a_{nn}-t
\end{array}
\right|$$
である。ここで、\(t\) の係数に着目する。
- \(t^n\) の係数
これは、対角成分を全て掛け合わせた次の多項式$$(a_{11}-t)(a_{22}-t)\ldots(a_{nn}-t)$$における \(t^n\) の係数と等しくなる。よって、係数は「\((-1)^n\)」
- \(t^{n-1}\) の係数 これも、先ほどの多項式の \(t^{n-1}\) の係数と等しくなる。（\(\because\) 余韻子展開を考える。\(\sigma\) が全単射のため、１行（列）対角成分ではないものがあると、少なくともあと１つは体格成分でないものがあることがわかる。）
- 定数項 これは、余韻子展開を考えれば \(|A|\) そのもの。
よって、係数は「\((-1)^{n-1}(a_{11}+a_{22}+\ldots+a_{nn})\)」であり、これは対角和を用いて「$(-1)^{n-1}\mathrm{Tr}(A) $」とも書ける。
一方先ほどの固有方程式は固有値を解に持つ。
つまり、\(A\)の固有値を\(\lambda_1\)〜\(\lambda_n\)として、
$$\phi(t)=(-1)^n(t-\lambda_1)(t-\lambda_2)\cdots(t-\lambda_n)$$
とかける。（ \((-1)^n\) によって、\(t^n\) の係数を合わせている。）
この式を展開すると、\(t^{n-1}\) の係数が「\((-1)^{n-1}(\lambda_1+\lambda_2+\ldots+\lambda_n)\)」であることから
$$\mathrm{Tr}(A) = \lambda_1+\lambda_2+\ldots+\lambda_n$$
※ちなみに、\(\phi(t)=(-1)^n(t-\lambda_1)(t-\lambda_2)\cdots(t-\lambda_n)\) の定数項を考えることで、
$$|A| = \lambda_1\lambda_2\ldots\lambda_n$$
であることもわかる。
トレースの循環性
\(A\) を \(m \times n\)、\(B\) を \(n \times m\) の行列とすると、\(AB\) は \(m \times m\) の行列であり、
$$
\begin{aligned}
\mathrm{Tr}[AB]&= \sum_{i=1}^{m} (AB)_{ii}
\\&= \sum_{i=1}^{m}\sum_{j=1}^{n} A_{ij}B_{ji}
\\&= \sum_{j=1}^{n} \sum_{i=1}^{m} B_{ji} A_{ij}
\\&= \sum_{j=1}^{n} (BA)_{jj}
\\&= \mathrm{Tr}[BA]
\end{aligned}
$$
となる。これを応用すれば、以下の循環性が証明できる。
\(A,B,C\) をそれぞれ \(m\times n,n \times l,l \times m\) の行列とするとき、
$$\begin{aligned}
\mathrm{Tr}[ABC] =& \mathrm{Tr}[BCA]
\\=& \mathrm{Tr}[CAB]
\end{aligned}
$$
が成り立つ。（ \(\because\) ２つの行列積をセットで考えれば明らか）
逆行列の微分
正則行列 \(A\) に対して、\(A^{-1}A=I\) が成立するので、この等式の両辺を \(A\) で微分して、
$$\begin{aligned}
\biggr(\frac{\partial }{\partial A}A^{-1}\biggl)A + A^{-1}\biggr(\frac{\partial }{\partial A}A\biggl) = 0\\
\biggr(\frac{\partial A^{-1}}{\partial A}\biggl) = -A^{-1}\biggr(\frac{\partial }{\partial A}A\biggl)A^{-1}
\end{aligned}$$
シングルエントリ行列
\((i,j)\) 成分のみが\(1\) で、残りが全て \(0\) の行列 \(\mathbf{J}^{ij}\) をシングルエントリ行列と呼ぶ。なお、以下の式が成り立つ。
$$\mathrm{Tr}\bigr(\mathbf{AJ}^{ij}\bigl)=\mathbf{A}_{ji}$$
\(2 \times 2\) の行列でこれを示す。
$$\begin{aligned}
\mathrm{Tr}\bigr(\mathbf{AJ}^{12}\bigl)&=\mathrm{Tr}\biggr(\begin{pmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{pmatrix}\begin{pmatrix}0&1\\0&0\end{pmatrix}\biggl)\\
&=\mathrm{Tr}\begin{pmatrix}0&a_{11}\\0&a_{21}\end{pmatrix}\\
&=a_{21}\\
&=\mathbf{A}_{21}
\end{aligned}$$
\(\ast2\) の変形
以上を踏まえれば、
$$\begin{aligned}
\frac{\partial}{\partial\boldsymbol{\Sigma}_{ij}}(\mathbf{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})
&=\frac{\partial}{\partial\boldsymbol{\Sigma}_{ij}}\mathrm{Tr}\biggr((\mathbf{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\biggl)\\
&=\frac{\partial}{\partial\boldsymbol{\Sigma}_{ij}}\mathrm{Tr}\biggr(\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\biggl)\\
&=\mathrm{Tr}\biggr(\frac{\partial}{\partial\boldsymbol{\Sigma}_{ij}}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\biggl)\\
&=\mathrm{Tr}\biggr(-\boldsymbol{\Sigma}^{-1}\frac{\partial\boldsymbol{\Sigma}}{\partial\boldsymbol{\Sigma}_{ij}}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\biggl)\\
&=\mathrm{Tr}\biggr(-\boldsymbol{\Sigma}^{-1}\mathbf{J}^{ij}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\biggl)\\
&=-\mathrm{Tr}\biggr(\boldsymbol{\Sigma}^{-1}\mathbf{J}^{ij}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\biggl)\\
&=-\mathrm{Tr}\biggr(\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}\mathbf{J}^{ij}\biggl)\\
&=-\left\{\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}\right\}_{ji}\\
&=-\left\{\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}\right\}_{ij}\\
\end{aligned}$$
よって、
$$\frac{\partial}{\partial\boldsymbol{\Sigma}}(\mathbf{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})=-\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}$$