はじめに
PRML解答例まとめを参照
演習 12.1
この演習問題では，帰納法を使って，射影されたデータの分散を最大化するようなM次元部分空間の上への線形写像がデータ共分散行列\mathbf{S},
\mathbf{S}=\frac{1}{N} \sum_{n=1}^{N}\left(\mathbf{x}_{n}-\overline{\mathbf{x}}\right)\left(\mathbf{x}_{n}-\overline{\mathbf{x}}\right)^{\mathrm{T}} \tag{12.3}
の上位M個の固有値に属するM本の固有ベクトルにより定義されることを証明する．12.1節では，M=1に対してこの結果を証明した．今度は，ある一般的な値Mに対してこの結果が成り立つと仮定して，その下でM+1次元に対しても成り立つことを示す．これを行うため，最初に，射影されたデータの分散の，ベクトル\mathbf{u}_{M+1}に対する微分を0とおく．\mathbf{u}_{M+1}はデータ空間における新しい方向を定義する．このとき，次の2つの制約を同時に満足しなければならない．ひとつは，\mathbf{u}_{M+1}がすでに求めたベクトル\mathbf{u}_{1},\ldots,\mathbf{u}_{M}と直交するという制約であり，もうひとつは単位長さに規格化しておかなければならないという制約である．この制約を取り込むためにラグランジュ乗数(付録E)を使ってみよ．そうして，新しいベクトル\mathbf{u}_{M+1}が\mathbf{S}の固有ベクトルであることを示すために，ベクトル\mathbf{u}_{1},\ldots,\mathbf{u}_{M}の正規直交性を利用せよ．最後に，固有値が大きい順に並べられているときに，その固有ベクトル\mathbf{u}_{M+1}を\lambda_{M+1}に対応する固有ベクトルに選べば，分散が最大化されることを示せ．
※P.278下部の「一般の場合としてM次元の射影空間を考えればデータ分散行列\mathbf{S}の, 大きい順にM個の固有値\lambda_{1},\ldots,\lambda_{M}に対応するM個の固有ベクトル\mathbf{u}_{1},\ldots,\mathbf{u}_{M}により，射影されたデータの分散を最大にする最適な線形射影が得られる」ことを帰納法で示す。
(i) M=1のとき、P.278の(12.4)-(12.6)の手続きによって\mathbf{Su}_{1} = \lambda_{1}\mathbf{u}_{1}となり、このとき\mathbf{u}_{1}^{\mathrm T}\mathbf{Su}_{1}は\lambda_{1}で最大となることが示される（本文参照）。
(ii) M次元についてもP.278の議論が成立していると仮定する。すなわち、
\mathbf{S}\mathbf{u}_{M} = \lambda_{M}\mathbf{u}_{M}
である。この条件下で、\mathbf{u}_{M+1}\mathbf{S}\mathbf{u}_{M+1}を\mathbf{u}_{M+1}に対して最大化したとき、\mathbf{S}\mathbf{u}_{M+1} = \lambda_{M+1}\mathbf{u}_{M+1}が成立することを示す。（この式から\lambda_{M+1}は\lambda_{1},\ldots,\lambda_{M}に次ぐ最大の\mathbf{S}の固有値となり、\mathbf{u}_{M+1}\mathbf{S}\mathbf{u}_{M+1} = \lambda_{M+1}となることは示される。）
2つの制約をラグランジュ未定乗数法を用いて導入する。1つは\mathbf{u}_{M+1}が正規直交基底であることの
\mathbf{u}_{M+1}^{\mathrm T}\mathbf{u}_{M+1} = 1
であり、もう1つは\mathbf{u}_{M+1}がすでに求めたベクトル\mathbf{u}_{1},\ldots,\mathbf{u}_{M}と直交することの
\mathbf{u}_{i}^{\mathrm T}\mathbf{u}_{M+1} = 0 \quad \textrm{for} \quad i=1,\ldots, M
である。これらの制約をそれぞれ未定乗数\lambda, \eta_{i}を用いてラグランジュ関数にすると
\mathbf{u}_{M+1}\mathbf{S}\mathbf{u}_{M+1} + \lambda_{M+1}\left( 1-\mathbf{u}_{M+1}^{\mathrm T}\mathbf{u}_{M+1} \right) + \sum_{i=1}^{M}\eta_{i}\mathbf{u}_{i}^{\mathrm T}\mathbf{u}_{M+1}
この関数の停留点を求める。\mathbf{u}_{M+1}についてこのラグランジュ関数を微分すると
2\mathbf{S}\mathbf{u}_{M+1} - 2\lambda_{M+1}\mathbf{u}_{M+1} + \sum_{i=1}^{M}\eta_{i}\mathbf{u}_{i} = 0
この式に左から\mathbf{u}_{j}^{\mathrm T}\ (j=1,\ldots, M)をかけると、第2項は正規直交基底の性質から0となる。第1項については
\begin{aligned}
\mathbf{u}_{j}^{\mathrm T}\mathbf{S}\mathbf{u}_{M+1} &= \mathbf{u}_{M+1}^{\mathrm T}\mathbf{S}\mathbf{u}_{j} \quad (\because \textrm{scholar}) \\
&=\mathbf{u}_{M+1}^{\mathrm T}\lambda_{j}\mathbf{u}_{j} \quad (\because \textrm{仮定}) \\
&=0
\end{aligned}
より0となる。よって第3項について、i \neq jでは0, i=jでは\eta_{j}となるので、
\eta_{i} = 0 \quad \textrm{for} \quad i=1,\ldots, M
となる。すなわち停留点は
\mathbf{S}\mathbf{u}_{M+1} = \lambda_{M+1}\mathbf{u}_{M+1}
のときに得られる。以上(i), (ii)から任意のM次元について\mathbf{S}\mathbf{u}_{M} = \lambda_{M}\mathbf{u}_{M}が成立するときに射影された分散\mathbf{u}_{M}^{\mathrm T}\mathbf{S}\mathbf{u}_{M}が最大化されることが帰納的に示された。
最後に、\mathbf{u}_{M+1}^{\mathrm T}を左からかければ
\mathbf{u}_{M+1}^{\mathrm T}\mathbf{S}\mathbf{u}_{M+1} = \lambda_{M+1}
がただちに得られ、\mathbf{u}_{M+1}ベクトルに対して射影されたデータの分散\mathbf{u}_{M+1}^{\mathrm T}\mathbf{S}\mathbf{u}_{M+1}は\lambda_{M+1}で最大値をとることが得られる。
演習 12.2
J=\frac{1}{N} \sum_{n=1}^{N} \sum_{i=M+1}^{D}\left(\mathbf{x}_{n}^{\mathrm{T}} \mathbf{u}_{i}-\overline{\mathbf{x}}^{\mathrm{T}} \mathbf{u}_{i}\right)^{2}=\sum_{i=M+1}^{D} \mathbf{u}_{i}^{\mathrm{T}} \mathbf{S} \mathbf{u}_{i} \tag{12.15}
で与えられる主成分分析の歪み尺度Jの，正規直交条件
\mathbf{u}_{i}^{\mathrm T}\mathbf{u}_{j} = \delta_{ij} \tag{12.7}
の下での\mathbf{u}_{i}に対する最小値は，\mathbf{u}_{i}がデータ共分散行列\mathbf{S}の固有ベクトルであるときに得られることを示せ．これを行うために，ラグランジュ乗数の行列\mathbf{H}を導入し制約条件のそれぞれを取り込む．その結果，歪み尺度の式は修正を受け，行列形式で表すと
\tilde{J}=\operatorname{Tr}\left\{\widehat{\mathbf{U}}^{\mathrm{T}} \mathbf{S} \widehat{\mathbf{U}}\right\}+\operatorname{Tr}\left\{\mathbf{H}\left(\mathbf{I}-\widehat{\mathbf{U}}^{\mathrm{T}} \widehat{\mathbf{U}}\right)\right\} \tag{12.93}
のようになる．ここで，\widehat{\mathbf{U}}はD \times (D - M)行列で，その列ベクトルは\mathbf{u}_{i}で与えられる．\widehat{\mathbf{U}}についてこの\tilde{J}を最小化しその解が\mathbf{S\widehat{U}}=\widehat{\mathbf{U}} \mathbf{H}を満たすことを示せ．明らかに，可能なひとつの解は，\widehat{\mathbf{U}}の列ベクトルが\mathbf{S}の固有ベクトルとなっている場合である．その場合， \mathbf{H}は対応する固有値を持った対角行列となる．一般的な解を得るために，\mathbf{H}が対称行列と仮定できることを示し，その固有ベクトル展開を用いることで，\mathbf{S\widehat{U}}=\widehat{\mathbf{U}} \mathbf{H}の一般解が，\mathbf{U}の列ベクトルを\mathbf{S}の固有ベクトルに選ぶという特解と同じ\tilde{J}の値を与えることを示せ．これらの解は等価なので固有ベクトルの方の解を選んだ方が楽である．
※ P.280下部の議論はM=1, D=2の特別な場合であり、これを任意のDとD<M条件下に拡張する。
（前半）
演習12.1のように、ラグランジュ未定乗数法を用いてラグランジュ関数\tilde{J}を定義する。
\tilde{J} = \underbrace{\sum_{i=M+1}^D \mathbf{u}_{i} \mathbf{S} \mathbf{u}_{i}}_{歪み尺度J} +
\underbrace{\sum_{i=M+1}^D \lambda_i (1-\mathbf{u}_{i}^{\mathrm T} \mathbf{u}_{i})}_{\mathbf{u}_{i} が規格化されている条件} +
\underbrace{\sum_{i=M+1}^{D-1}\sum_{j=i+1}^{D}\mu_{ij}(-\mathbf{u}_{j}\mathbf{u}_{i})}_{\mathbf{u}_{j} と \mathbf{u}_{i} が直交している条件}
これを展開していくと(12.93)が得られることを示す。問題文の設定から
\widehat{\mathbf{U}} = \begin{pmatrix} \mathbf{u}_{M+1} & \mathbf{u}_{M+2} & \cdots & \mathbf{u}_{D} \end{pmatrix}
ラグランジュ乗数の行列\mathbf{H}を以下のように設定する。
\mathbf{H}=\begin{pmatrix}
\lambda_{M+1} & \frac{1}{2} \mu_{M+1, M+2} & \ldots & \frac{1}{2} \mu_{M+1, D} \\
\frac{1}{2} \mu_{M+1, M+2} & \lambda_{M+2} & \ldots & \vdots \\
\vdots & \vdots & \ddots & \vdots \\
\frac{1}{2} \mu_{M+1, D} & \ldots & \ldots & \lambda_{D}
\end{pmatrix} \quad (\mathbf{H}は対称行列)
まず(12.93)の第1項の歪み尺度\displaystyle \sum_{i=M+1}^D \mathbf{u}_{i} \mathbf{S} \mathbf{u}_{i}が\operatorname{Tr}\left\{\widehat{\mathbf{U}}^{\mathrm{T}} \mathbf{S} \widehat{\mathbf{U}}\right\}と一致することを示す。これは
\begin{aligned}
\hat{\mathbf{U}}^{\mathrm T} \mathbf{S} \hat{\mathbf{U}} &=\begin{pmatrix}
\mathbf{u}_{M+1}^{\mathrm T} \\
\mathbf{u}_{M+2}^{\mathrm T} \\
\vdots \\
\mathbf{u}_{D}^{\mathrm T}
\end{pmatrix} \mathbf{S} \begin{pmatrix} \mathbf{u}_{M+1} & \mathbf{u}_{M+2} & \ldots &\mathbf{u}_{D}
\end{pmatrix}\\
&=\begin{pmatrix}
\mathbf{u}_{M+1}^{\mathrm T} \\
\mathbf{u}_{M+2}^{\mathrm T} \\
\vdots \\
\mathbf{u}_{D}^{\mathrm T}
\end{pmatrix}
\begin{pmatrix}
\mathbf{S} \mathbf{u}_{M+1} & \mathbf{S} \mathbf{u}_{M+2} & \ldots & \mathbf{S} \mathbf{u}_{D}
\end{pmatrix}
\\
&=\begin{pmatrix}
\mathbf{u}_{M+1}^{\mathrm T} \mathbf{S u}_{M+1} & \mathbf{u}_{M+1}^{\mathrm T} \mathbf{S u}_{M+2} & \ldots & \mathbf{u}_{M+1}^{\mathrm T} \mathbf{S u}_{D} \\
\mathbf{u}_{M+2}^{\mathrm T} \mathbf{S u}_{M+1} & \mathbf{u}_{M+2}^{\mathrm T} \mathbf{S u}_{M+2} & \ldots & \vdots \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{u}_{D}^{\mathrm T} \mathbf{S u}_{M+1} & \ldots & \ldots & \mathbf{u}_{D}^{\mathrm T} \mathbf{S u}_{D}
\end{pmatrix}
\end{aligned}
これより\displaystyle \sum_{i=M+1}^D \mathbf{u}_{i} \mathbf{S} \mathbf{u}_{i} = \operatorname{Tr}\left\{\widehat{\mathbf{U}}^{\mathrm{T}} \mathbf{S} \widehat{\mathbf{U}}\right\}が示された。
続いて残りを計算する
\begin{aligned}
&\mathbf{H}\left(\mathbf{I}-\widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}}\right) \\
=\ &\mathbf{H}\left(\mathbf{I}-\begin{pmatrix}
\mathbf{u}_{M+1}^{\mathrm T} \\
\mathbf{u}_{M+2}^{\mathrm T} \\
\vdots \\
\mathbf{u}_{D}^{\mathrm T}
\end{pmatrix}\begin{pmatrix}
\mathbf{u}_{M+1} & \mathbf{u}_{M+2} & \ldots & \mathbf{u}_{D}
\end{pmatrix}\right)\\
=\ & \mathbf{H}\left(\begin{pmatrix}
1 & 0 & \ldots & 0 \\
0 & 1 & \ldots & \vdots \\
\vdots & \vdots & \ddots & \vdots \\
0 & \ldots & \ldots & 1
\end{pmatrix}-\begin{pmatrix}
\mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{M+1} & \mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{M+2} & \ldots & \mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{D} \\
\mathbf{u}_{M+2}^{\mathrm T} \mathbf{u}_{M+1} & \mathbf{u}_{M+2}^{\mathrm T} \mathbf{u}_{M+2} & \ldots & \vdots \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{u}_{D}^{\mathrm T} \mathbf{u}_{M+1} & \ldots & \ldots & \mathbf{u}_{D}^{\mathrm T} \mathbf{u}_{D}
\end{pmatrix}\right) \\
=\ &\begin{pmatrix}
\lambda_{M+1} & \frac{1}{2} \mu_{M+1, M+2} & \cdots & \frac{1}{2} \mu_{M+1, D} \\
\frac{1}{2} \mu_{M+1, M+2} & \lambda_{M+2} & \cdots & \vdots \\
\vdots & \vdots & \ddots & \vdots \\
\frac{1}{2} \mu_{M+1, D} & \cdots & \cdots & \lambda_{D}
\end{pmatrix}\begin{pmatrix}
1-\mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{M+1} & -\mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{M+2} & \cdots & -\mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{D} \\
-\mathbf{u}_{M+2}^{\mathrm T} \mathbf{u}_{M+1} & 1-\mathbf{u}_{M+2}^{\mathrm T} \mathbf{u}_{M+2} & \cdots & \vdots \\
\vdots & \vdots & \ddots & \vdots \\
-\mathbf{u}_{D}^{\mathrm T} \mathbf{u}_{M+1} & \cdots & \cdots & 1-\mathbf{u}_{D}^{\mathrm T} \mathbf{u}_{D}
\end{pmatrix}
\end{aligned}
これより
\begin{aligned}
& \operatorname{Tr}\left(\mathbf{H}\left(\mathbf{I}-\widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}}\right)\right) \\
=\ & \lambda_{M+1}\left(1-\mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{M+1}\right)+\frac{1}{2} \mu_{M+1, M+2}\left(-\mathbf{u}_{M+2}^{\mathrm T} \mathbf{u}_{M+1}\right)+\ldots+\frac{1}{2} \mu_{M+1, D}\left(-\mathbf{u}_{D}^{\mathrm T} \mathbf{u}_{M+1}\right) \\
&+\frac{1}{2} \mu_{M+1, M+2}\left(-\mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{M+2}\right)+\lambda_{M+2}\left(1-\mathbf{u}_{M+2}^{\mathrm T} \mathbf{u}_{M+2}\right)+\ldots \\
& \quad\quad\quad \vdots \\
&+\frac{1}{2} \mu_{M+1, D}\left(-\mathbf{u}_{M+1}^{\mathrm T} \mathbf{u}_{D}\right)+\ldots+\lambda_{D}\left(1-\mathbf{u}_{D}^{\mathrm T} \mathbf{u}_{D}\right) \\
=\ & \sum_{i=M+1}^{D} \lambda_{i}\left(1-\mathbf{u}_{i}^{\mathrm T} \mathbf{u}_{i}\right)+\sum_{i=M+1}^{D-1} \sum_{j=i+1}^{D} \mu_{i j}\left(-\mathbf{u}_{j}^{\mathrm T} \mathbf{u}_{i}\right)
\end{aligned}
以上から(12.93)式が示された。
（2）「\widehat{\mathbf{U}}についてこの\tilde{J}を最小化しその解が\mathbf{S\widehat{U}}=\widehat{\mathbf{U}} \mathbf{H}を満たすことを示せ」
(12.93)式を\widehat{\mathbf{U}}で微分すると、Matrix Cookbookの公式(108),(112)を使って
\begin{aligned}
\frac{\partial\tilde{J}}{\partial \widehat{\mathbf{U}}}&=
\frac{\partial}{\partial \widehat{\mathbf{U}}}\operatorname{Tr}\left\{\widehat{\mathbf{U}}^{\mathrm{T}} \mathbf{S} \widehat{\mathbf{U}}\right\}+\frac{\partial}{\partial \widehat{\mathbf{U}}}\operatorname{Tr}\left\{\mathbf{H}\left(\mathbf{I}-\widehat{\mathbf{U}}^{\mathrm{T}} \widehat{\mathbf{U}}\right)\right\} \\
&=(\mathbf{S} \widehat{\mathbf{U}}+\mathbf{S}^{\mathrm T} \widehat{\mathbf{U}}) - (\widehat{\mathbf{U}}\mathbf{H}^{\mathrm T}+\widehat{\mathbf{U}}\mathbf{H}) \\
&=2\mathbf{S} \widehat{\mathbf{U}} - 2\widehat{\mathbf{U}}\mathbf{H}
\end{aligned}
これを0とすれば、\mathbf{S} \widehat{\mathbf{U}} = \widehat{\mathbf{U}}\mathbf{H}を満たす\widehat{\mathbf{U}}が解となり、このとき\tilde{J}は最小値
\begin{aligned} \tilde{J} &=\operatorname{Tr}\left(\widehat{\mathbf{U}}^{\mathrm T} \mathbf{S} \widehat{\mathbf{U}}\right)+\operatorname{Tr}\left(\mathbf{H}\left(\mathbf{I}-\widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}}\right)\right) \\ &=\operatorname{Tr}\left(\widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}} \mathbf{H}\right)+\operatorname{Tr}\left(\mathbf{H}\left(\mathbf{I}-\widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}}\right)\right) \\ &=\operatorname{Tr}(\mathbf{H})=\sum_{i=M+1}^{D} \lambda_{i} \end{aligned}
を得る。これはすなわち\mathbf{S}の固有値のうち大きい順からM+1\sim D番目の固有値となるときに最小となる。
(3) 一般的な解を得るために，\mathbf{H}が対称行列と仮定できることを示し，その固有ベクトル展開を用いることで，\mathbf{S\widehat{U}}=\widehat{\mathbf{U}} \mathbf{H}の一般解が，\mathbf{U}の列ベクトルを\mathbf{S}の固有ベクトルに選ぶという特解と同じ\tilde{J}の値を与えることを示せ．
\mathbf{H}の一般解を議論するために、すでに上で設定した対称行列\mathbf{H}を用いる。\mathbf{\Phi}を(D-M)\times(D-M)の対角行列、すなわち\mathbf{\Phi}^{\mathrm T}\mathbf{\Phi} = \mathbf{\Phi}\mathbf{\Phi}^{\mathrm T}=\mathbf{I}を満たす行列とすれば、任意の対称行列の異なる固有値に対応する固有ベクトルは直交するので、
\mathbf{H\Psi} = \mathbf{\Psi L}
と書ける。ここで\mathbf{L}は対角成分が\mathbf{H}の固有値となっている対角行列である（※ここでは固有値は縮退している可能性がある）。
\mathbf{S} \widehat{\mathbf{U}} = \widehat{\mathbf{U}}\mathbf{H}に右から\mathbf{\Psi}をかけ、新たに\widetilde{\mathbf{U}} = \widehat{\mathbf{U}}\mathbf{\Psi}を定義すると
\mathbf{SU\Psi} =\widehat{\mathbf{U}} \mathbf{H\Psi} =\widehat{\mathbf{U}} \mathbf{\Psi L} \\
\mathbf{S}\widetilde{\mathbf{U}} = \widetilde{\mathbf{U}}\mathbf{L}
となる。よって\widetilde{\mathbf{U}}の列ベクトルは\mathbf{S}の固有ベクトルになり、\mathbf{L}の対角成分が固有値となる。
最後に\mathbf{S}\widetilde{\mathbf{U}} = \widetilde{\mathbf{U}}\mathbf{L}を用いたときに\tilde{J}の最小値が\mathbf{S} \widehat{\mathbf{U}} = \widehat{\mathbf{U}}\mathbf{H}を用いたときと等価になることを示す。
\begin{aligned}
\tilde{J} &=\operatorname{Tr}\left(\widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}}\right)+\operatorname{Tr}\left(\mathbf{H}\left(\mathbf{I}-\widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}}\right)\right) \\
&=\operatorname{Tr}\left(\widehat{\mathbf{U}}^{\mathrm T} \mathbf{S} \widehat{\mathbf{U}} \mathbf{\Psi \Psi}^{\mathrm T}\right)+\operatorname{Tr}\left(\mathbf{H \Psi \Psi}^{\mathrm T}\right)-\operatorname{Tr}\left(\mathbf{H} \widehat{\mathbf{U}}^{\mathrm T} \widehat{\mathbf{U}} \mathbf{\Psi \Psi}^{\mathrm T}\right) \\
&=\operatorname{Tr}\left(\widetilde{\mathbf{U}}^{\mathrm T} \mathbf{S} \widetilde{\mathbf{U}}\right)+\operatorname{Tr}\left(\mathbf{\Psi S \Psi}^{\mathrm T}\right)-\operatorname{Tr}\left(\mathbf{L \Psi}^{\mathrm T} \widehat{\mathbf{U}}^{\mathrm T} \widetilde{\mathbf{U}}\right) \\
&=\operatorname{Tr}(\mathbf{L})+\operatorname{Tr}(\mathbf{L})-\operatorname{Tr}(\mathbf{L}) \\
&=\operatorname{Tr}(\mathbf{L}) \\
&=\sum_{i=1}^{D-m} \lambda_{i}
\end{aligned}
これより、\tilde{J}を最小化するためには小さい順に\mathbf{S}のD-m個の固有値を選べば良い。この結果は一般的な対称行列\mathbf{H}を用いた場合の値と等価である。
演習 12.3
\mathbf{u}_{i}=\frac{1}{\left(N \lambda_{i}\right)^{1 / 2}} \mathbf{X}^{\mathrm{T}} \mathbf{v}_{i} \tag{12.30}
で定義される固有ベクトルが単位長さに規格化されていることを示せ．ただし固有ベクトル\mathbf{v}_{i}は単位長さを持っていると仮定する．
\| \mathbf{u}_{i}\|^{2}=1であることを示せば良い。P.286の\mathbf{v}_{i} = \mathbf{Xu}_{i}の定義と
\frac{1}{N} \mathbf{X X}^{\mathrm{T}} \mathbf{v}_{i}=\lambda_{i} \mathbf{v}_{i} \tag{12.28}
を利用する。
\begin{aligned}
\mathbf{u}_{i}^{\mathrm T} \mathbf{u}_{i} &=\left(\frac{1}{\left(N \lambda_{i}\right)^{1/2}}\right)^{2} \mathbf{v}_{i}^{\mathrm T} \mathbf{X} \mathbf{X}^{\mathrm T} \mathbf{v}_{i} \\
&=\frac{1}{N \lambda_{i}} \mathbf{v}_{i}^{\mathrm T}\left(N \lambda_{i}\right) \mathbf{v}_{i} \quad (\because(12.28)) \\
&=\mathbf{v}_{i}^{\mathrm T} \mathbf{v}_{i}=1 \quad \left(\because\left\|\mathbf{v}_{i}\right\|^{2}=1\right)
\end{aligned}
以上で示された。
演習 12.4
確率的主成分分析モデルの潜在変数の空間分布
p(\mathbf{z})=\mathcal{N}(\mathbf{z} \mid \mathbf{0}, \mathbf{I}) \tag{12.31}
では，平均0で単位共分散のものが仮定されている．それを，より一般的なガウス分布\mathcal{N}(\mathbf{z}\mid \mathbf{m},\mathbf{\Sigma})で置き換えるものとしよう．モデルのパラメータの再定義により，\mathbf{m}と\mathbf{\Sigma}の任意の選択が，観測変数についての同ーの周辺分布p(\mathbf{x})を導くことを示せ．
(2.113)-(2.115)式より、
\begin{aligned}
p(\mathbf{x}) &= \mathcal{N}(\mathbf{Wm} + \boldsymbol{\mu}, \sigma^2\mathbf{I} + \mathbf{W\Sigma W}^T)
\end{aligned}
ここで、
\begin{aligned}
\tilde{\boldsymbol{\mu}} &= \mathbf{Wm} + \boldsymbol{\mu} \\
\tilde{\mathbf{W}} &= \mathbf{W\Sigma}^{1/2}
\end{aligned}
とおくと、p(\mathbf{x})は
\begin{aligned}
p(\mathbf{x}) &= \mathcal{N}(\tilde{\boldsymbol{\mu}}, \sigma^2\mathbf{I} + \tilde{\mathbf{W}}\tilde{\mathbf{W}}^{T})
\end{aligned}
と表すことができ、これは
p(\mathbf{x})=\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{C}) \tag{12.35}
と同じ形である。
演習 12.5
\mathbf{x}をD次元のガウス分布\mathcal{N}(\mathbf{x}\mid \boldsymbol{\mu},\mathbf{\Sigma})に従う確率変数とし，\mathbf{y}をM \times D行列\mathbf{A}によって\mathbf{y} = \mathbf{Ax} + \mathbf{b}で定義されるM次元確率変数とする．\mathbf{y}もまたガウス分布を持つことを示し，その平均と分散に対する表現を見出せ．M \lt D，M = DおよびM \gt Dの各場合に対し，そのガウス分布の形がどうなるか考察せよ．
\mathbf{A}の「第i行、第j列」成分をa_{ij}、\mathbf{\Sigma}の「第i行、第j列」成分を\sigma_{ij}とおく。以降は、\mathbf{A}が正方行列ではないケースも包含した議論（のつもり）。
ガウス分布の線型性から、線型結合で得られる\mathbf{y}がガウス分布に従うことは明らか。
また、期待値に関する線型性から、
\begin{aligned}
\mathbb{E}[\mathbf{A}\mathbf{x}+\mathbf{b}]
= \mathbf{A} \mathbb{E}[\mathbf{x} ] +\mathbf{b}
= \mathbf{A} \boldsymbol{\mu} +\mathbf{b}
\end{aligned}
である。
共分散行列\mathbb{V}[\mathbf{Ax}+\mathbf{b}]=\mathbb{V}[\mathbf{Ax}]の「第i行、第j列」成分に関しては、
\begin{aligned}
\big( \mathbb{V}[\mathbf{Ax}]\big)_{ij}
&= {\rm E} \left[ \sum_{k=1}^D a_{ik}x_k \sum_{l=1}^D a_{jl}x_l \right] - {\rm E} \left[ \sum_{k=1}^D a_{ik}x_k \right] E\left[ \sum_{l=1}^D a_{jl}x_l \right] \\
&= \sum_{k,l=1}^D a_{ik}a_{jl}{\rm E} \left[x_k x_l \right] - \sum_{k,l=1}^D a_{ik} a_{jl}{\rm E} \left[ x_k \right]E\left[ x_l \right] \\
&= \sum_{k,l=1}^D a_{ik}a_{jl} \left( {\rm E}\left[x_k x_l \right] - {\rm E} \left[ x_k \right]E\left[ x_l \right]
\right)\\
&= \sum_{k,l=1}^D a_{ik}a_{jl} \sigma_{kl}\\
&= \left( \mathbf{A}\boldsymbol{\Sigma} \mathbf{A}^{\rm T} \right)_{ij}
\end{aligned}
であり、この式変形は\mathbf{A}が正方行列でなくても成立する・・・ような気がするのですが。
演習 12.6
確率的主成分分析モデルに対する有向確率グラフを描け．有向確率グラフについては12.2節に説明がある．観測変数\mathbf{x}の各要素は明示的に分離したノードとして示されるはずであり，したがって確率的主成分分析モデルでは，8.2.2節で議論したナイーブベイズモデルと同様の独立構造を持っていることを確かめよ．
確率的主成分分析モデルの有向確率グラフは図12.10のようになり，この図中のパラメータを省略し確率変数z, xのみに着目すると，8.2.2節の図8.24のナイーブベイズの有向グラフと同一の構造になることがわかる．このことから同様な独立構造を持っていることがわかる．
演習 12.7
一般的な分布の平均と分散に対する
\begin{aligned}
\mathbb{E}[x] &=\mathbb{E}_{y}\left[\mathbb{E}_{x}[x \mid y]\right] & (2.270)\\ \operatorname{var}[x] &=\mathbb{E}_{y}\left[\operatorname{var}_{x}[x \mid y]\right]+\operatorname{var}_{y}\left[\mathbb{E}_{x}[x \mid y]\right] & (2.271)
\end{aligned}
の結果を利用して。確率的主成分分析モデルの周辺分布p(\mathbf{x})に対する結果
p(\mathbf{x})=\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{C}) \tag{12.35}
を導け．
\begin{aligned}
\mathbb{E}[\mathbf{x}] &= \mathbb{E}_z[\mathbb{E}_x[\mathbf{x}|\mathbf{z}]] \\
&= \mathbb{E}_z[\mathbf{Wz} + \mathbf{\mu}] \\
&= \mathbf{\mu}.
\end{aligned}
\begin{aligned}
cov[\mathbf{x}] &= \mathbb{E}_z[cov_x[\mathbf{x}|\mathbf{z}]] + cov_z[\mathbb{E}_x[\mathbf{x}|\mathbf{z}]] \\
&= \mathbb{E}_z[\sigma^2\mathbf{I}] + cov_z[\mathbf{Wz} + \mathbf{\mu}] \\
&= \sigma^2\mathbf{I} + \mathbb{E}_z[\mathbf{Wz}\mathbf{z}^T\mathbf{W}^T] \\
&= \sigma^2\mathbf{I} + \mathbf{W}\mathbf{W}^T
\end{aligned}
演習 12.8
p(\mathbf{x} \mid \mathbf{y})=\mathcal{N}\left(\mathbf{x} \mid \mathbf{\Sigma}\left\{\mathbf{A}^{\mathrm{T}} \mathbf{L}(\mathbf{y}-\mathbf{b})+\mathbf{\Lambda} \boldsymbol{\mu}\right\}, \mathbf{\Sigma}\right) \tag{2.116}
の結果を利用して確率的主成分分析モデルで出てくる事後分布p(\mathbf{z}\mid \mathbf{x})が
p(\mathbf{z} \mid \mathbf{x})=\mathcal{N}\left(\mathbf{z} \mid \mathbf{M}^{-1} \mathbf{~W}^{\mathbf{T}}(\mathbf{x}-\boldsymbol{\mu}), \sigma^{2} \mathbf{M}^{-1}\right) \tag{12.42}
で与えられることを示せ．
(2.113), (2.114), (2.116), (2.117), (12.31), (12.32)式より、
\begin{aligned}
p(\mathbf{z}|\mathbf{x}) &= \mathcal{N}(\mathbf{z}|(\mathbf{I} + \sigma^{-2}\mathbf{W}^T\mathbf{W})^{-1}\mathbf{W}^T\sigma^{-2}\mathbf{I}(\mathbf{x} - \mathbf{\mu}), (\mathbf{I} + \sigma^{-2}\mathbf{W}^T\mathbf{W})^{-1})
\end{aligned}
ここで、
\begin{aligned}
(\mathbf{I} + \sigma^{-2}\mathbf{W}^T\mathbf{W})^{-1} &= (\sigma^{-2}\mathbf{M})^{-1}\ \ (\because \mathrm{Eq}(12.41)) \\
&= \sigma^2\mathbf{M}^{-1},
\end{aligned}
であることから、
\begin{aligned}
p(\mathbf{z}|\mathbf{x}) &= \mathcal{N}(\mathbf{z}|\mathbf{M}^{-1}\mathbf{W}^T(\mathbf{x} - \mathbf{\mu}), \sigma^2\mathbf{M}^{-1})
\end{aligned}
演習 12.9
確率的主成分分析モデルの対数尤度
\begin{aligned}
\ln p\left(\mathbf{X} \mid \boldsymbol{\mu}, \mathbf{W}, \sigma^{2}\right) &= \sum_{n=1}^{N} \ln p\left(\mathbf{x}_{n} \mid \mathbf{W}, \boldsymbol{\mu}, \sigma^{2}\right) \\
&=-\frac{N D}{2} \ln (2 \pi)-\frac{N}{2} \ln |\mathbf{C}|-\frac{1}{2} \sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)^{\mathrm{T}} \mathbf{C}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)
\end{aligned}
\tag{12.43}
をパラメータ\boldsymbol{\mu}に対して最大化すると，\boldsymbol{\mu}_{\mathrm{ML}}=\overline{\mathbf{x}}の結果になることを確かめよ．ただし\overline{\mathbf{x}}はデータベクトルの平均である．
対数尤度を\ln L = \ln p(\mathbf{X}\mid \boldsymbol{\mu}, \mathbf{W}, \sigma^2)とし、\frac{\partial \ln L}{\partial \boldsymbol{\mu}} = 0を求める。(12.36)の定義より\mathbf{C}=\mathbf{W} \mathbf{W}^{\mathrm{T}}+\sigma^{2} \mathbf{I}である。
\begin{aligned}
\frac{\partial \ln L}{\partial \boldsymbol{\mu}} &=-\frac{1}{2} \sum_{n=1}^{N} \frac{\partial}{\partial \boldsymbol{\mu}}\left(\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)^{\mathrm T} \mathbf{C}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)\right) \\
&=-\frac{1}{2} \sum_{n=1}^{N}\left(-2 \mathbf{C}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)\right) \\
&=-\sum_{n=1}^{N} \mathbf{C}^{-1} \boldsymbol{\mu}+\sum_{n=1}^{N} \mathbf{C}^{-1} \mathbf{x}_{n} \\
&= -N\mathbf{C}^{-1} \boldsymbol{\mu}+\sum_{n=1}^{N} \mathbf{C}^{-1} \mathbf{x}_{n}=0
\end{aligned}
これを解くと
\boldsymbol{\mu} = \frac{1}{N}\sum_{n=1}^{N}\mathbf{x}_{n} = \overline{\mathbf{x}}
すなわちデータベクトルの平均が得られた。
演習 12.10
確率的主成分分析モデルの対数尤度関数
\begin{aligned}
\ln p\left(\mathbf{X} \mid \boldsymbol{\mu}, \mathbf{W}, \sigma^{2}\right) &= \sum_{n=1}^{N} \ln p\left(\mathbf{x}_{n} \mid \mathbf{W}, \boldsymbol{\mu}, \sigma^{2}\right) \\
&=-\frac{N D}{2} \ln (2 \pi)-\frac{N}{2} \ln |\mathbf{C}|-\frac{1}{2} \sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)^{\mathrm{T}} \mathbf{C}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)
\end{aligned}
\tag{12.43}
のパラメータ\boldsymbol{\mu}に対する2次微分を求めることにより，停留点\boldsymbol{\mu}_{\mathrm{ML}}=\overline{\mathbf{x}}が唯一の最大値を与える点となることを示せ．
演習問題12.9の続きで、もう一度\boldsymbol{\mu}で微分すると
\frac{\partial}{\partial \boldsymbol{\mu}}\left(\frac{\partial \ln L}{\partial \boldsymbol{\mu}}\right) = -N\mathbf{C}^{-1}
となる。ここで\mathbf{C}は正定値対称行列なので-N\mathbf{C}^{-1}は負定値対称行列になる。これはすなわち対数尤度関数が\boldsymbol{\mu}の値によらず上に凸であることを表すため、12.9で求めた点は極大かつ最大となる。
演習 12.11
\sigma^{2}\to 0の極限において，確率的主成分分析モデルの事後平均が，通常の主成分分析と同様に主部分空間の上への直交射影となることを示せ．
\sigma^2 \rightarrow 0のとき、 (12.41)式を(12.48)式に代入すると、
\begin{aligned}
\mathbb{E}[\mathbf{z}|\mathbf{x}] = (\mathbf{W}_{ML}^T\mathbf{W}_{ML})^{-1}\mathbf{W}_{ML}^T (\mathbf{x} - \bar{\mathbf{x}})
\end{aligned}
を得る。\mathbf{R} = \mathbf{I}として、(12.45)式を代入する。
\begin{aligned}
(\mathbf{W}_{ML}^T\mathbf{W}_{ML})^{-1}\mathbf{W}_{ML}^T (\mathbf{x} - \bar{\mathbf{x}}) &= \left[\left\{(\mathbf{L}_M - \sigma^2\mathbf{I})^{1/2}\right\}^T\mathbf{U}_M^T\mathbf{U}_M(\mathbf{L}_M - \sigma^2\mathbf{I})^{1/2}\right]\mathbf{W}^T (\mathbf{x} - \bar{\mathbf{x}}) \\
&= (\mathbf{L}_M - \sigma^2\mathbf{I})^{-1}(\mathbf{L}_M - \sigma^2\mathbf{I})^{1/2}\mathbf{U}_M^T (\mathbf{x} - \bar{\mathbf{x}}) \\
&= (\mathbf{L}_M - \sigma^2\mathbf{I})^{-1/2}\mathbf{U}_M^T (\mathbf{x} - \bar{\mathbf{x}}).
\end{aligned}
ここで1行目から2行目は(\mathbf{L}_M - \sigma^2\mathbf{I})が対称行列であること及び定義より\mathbf{U}_M^T\mathbf{U}_M = \mathbf{I}であることを用いた。\sigma^2 \rightarrow 0のとき上式は
\begin{aligned}
\mathbf{L}_M^{-1/2}\mathbf{U}_M^T (\mathbf{x} - \bar{\mathbf{x}})
\end{aligned}
となるが、これは非確率的な主成分分析について得られた結果 ((12.24)式) に等しい。
演習 12.12
\sigma^{2}\gt 0のとき，確率的主成分分析モデルの事後平均が，直交射影の結果と比べると，原点に向かってシフトしていることを示せ．
演習12.11の結果を用い、\sigma^2 \rightarrow 0のとき及び\sigma^2 > 0のときの\mathbf{z}の期待値 (\mathbf{z}_0及び\mathbf{z}_1と表記する) のノルム (の2乗) を比較する。
\begin{aligned}
\|\mathbf{z}_0\|^2 - \|\mathbf{z}_1\|^2 &= (\mathbf{x} - \bar{\mathbf{x}})^T(\mathbf{U}_M\mathbf{L}_M^{-1}\mathbf{U}_M^T - \mathbf{W}_{ML}\mathbf{M}^{-1}\mathbf{M}^{-1}\mathbf{W}_{ML})(\mathbf{x} - \bar{\mathbf{x}}).
\end{aligned}
ここで、
\begin{aligned}
\mathbf{M}^{-1} &= (\mathbf{W}_{ML}^T\mathbf{W}_{ML} + \sigma^2\mathbf{I})^{-1} \\
&= ((\mathbf{L}_M - \sigma^2\mathbf{I}) + \sigma^2\mathbf{I})^{-1} \\
&= \mathbf{L}_M^{-1}.
\end{aligned}
よって、
\begin{aligned}
\|\mathbf{z}_0\|^2 - \|\mathbf{z}_1\|^2 &= (\mathbf{x} - \bar{\mathbf{x}})^T(\mathbf{U}_M\mathbf{L}_M^{-1}\mathbf{U}_M^T - \mathbf{W}_{ML}\mathbf{L}_M^{-2}\mathbf{W}_{ML})(\mathbf{x} - \bar{\mathbf{x}}) \\
&= (\mathbf{x} - \bar{\mathbf{x}})^T\mathbf{U}_M(\mathbf{L}_M^{-1} - (\mathbf{L}_M - \sigma^2\mathbf{I})^{1/2}\mathbf{L}_M^{-2}(\mathbf{L}_M - \sigma^2\mathbf{I})^{1/2})\mathbf{U}_M^T(\mathbf{x} - \bar{\mathbf{x}}) \tag{*}
\end{aligned}
となる。ここで、行列(\mathbf{L}_M^{-1} - (\mathbf{L}_M - \sigma^2\mathbf{I})^{1/2}\mathbf{L}_M^{-2}(\mathbf{L}_M - \sigma^2\mathbf{I})^{1/2})は、定義よりi番目の対角要素が
\begin{aligned}
\frac{1}{\lambda_i} - \frac{\lambda_i - \sigma^2}{\lambda_i^2} = \frac{\sigma^2}{\lambda_i^2} > 0
\end{aligned}
に等しい対角行列である。よって二次形式(*)は正。つまり、\|\mathbf{z}_0\| > \|\mathbf{z}_1\|となり、\sigma^{2}\gt 0のときに事後平均が原点に向かってシフトしていることが示された。
演習 12.13
確率的主成分分析の下でもともとのデータ点を再現する際，通常の主成分分析の最小2乗射影のコスト関数を使うと，最適な再現点は以下で与えられることを示せ．
\widetilde{\mathbf{x}}=\mathbf{W}_{\mathrm{ML}}\left(\mathbf{W}_{\mathrm{ML}}^{\mathrm{T}}\mathbf{W}_{\mathrm{ML}}\right)^{-1} \mathbf{M} \mathbb{E}[\mathbf{z} \mid \mathbf{x}] \tag{12.94}
\mathbb{E}[\mathbf{z} \mid \mathbf{x}]=\mathbf{M}^{-1} \mathbf{W}_{\mathrm{ML}}^{\mathrm{T}}(\mathbf{x}-\overline{\mathbf{x}}) \tag{12.48}
の右辺を(12.94)に代入すると、次のようになる。
\widetilde{\mathbf{x}}_{n}=\mathbf{W}_{\mathrm{ML}}\left(\mathbf{W}_{\mathrm{ML}}^{\mathrm{T}} \mathbf{W}_{\mathrm{ML}}\right)^{-1} \mathbf{W}_{\mathrm{ML}}^{\mathrm{T}}\left(\mathbf{x}_{n}-\overline{\mathbf{x}}\right)
\mathbf{W}_{\mathrm{ML}}=\mathbf{U}_{M}\left(\mathbf{L}_{M}-\sigma^{2} \mathbf{I}\right)^{1 / 2} \mathbf{R} \tag{12.45}
より
\left(\mathbf{W}_{\mathrm{ML}}^{\mathrm{T}} \mathbf{W}_{\mathrm{ML}}\right)^{-1}=\left(\mathbf{L}_{M}-\sigma^{2} \mathbf{I}\right)^{-1}
そして
\widetilde{\mathbf{x}}_{n}=\mathbf{U}_{M} \mathbf{U}_{M}^{\mathrm{T}}\left(\mathbf{x}_{n}-\overline{\mathbf{x}}\right)
これは、M個の最大固有値に対応するM個の固有ベクトルを使ってx_n - xを再構成したもので、12.1.2節から最小二乗射影コスト(12.11)を最小にすることがわかる。
演習 12.14
確率的主成分分析モデルの共分散行列の中にある独立なパラメータの数は，
D M+1-M(M-1) / 2 \tag{12.51}
で与えられる．ただし，潜在変数空間がM次元データ空間がD次元である．M = D-1の場合に，独立なパラメータの数が一般的な共分散行列を持つガウス分布と一致することを確かめよ．一方，M=0のときにそれが等方的な共分散を持つガウス分布と同じであることも確かめよ．
M = D - 1のとき、パラメータの数は
\begin{aligned}
DM + 1 - \frac{M(M - 1)}{2} &= D(D - 1) + 1 - \frac{(D - 1)(D - 2)}{2} \\
&= \frac{1}{2}(D^2 + D) \\
&= \frac{1}{2}(D^2 - D) + D.
\end{aligned}
となる。また、式(12.51)にM = 0を代入すると、このときパラメータの数は1となる (i.e., パラメータは\sigma^2のみ) ことがわかる。
演習 12.15
\begin{aligned}\mathbb{E}\left[\ln p\left(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \mathbf{W}, \sigma^{2}\right)\right]=&-\sum_{n=1}^{N}\left\{\frac{D}{2} \ln \left(2 \pi \sigma^{2}\right)+\frac{1}{2} \operatorname{Tr}\left(\mathbb{E}\left[\mathbf{z}_{n} \mathbf{z}_{n}^{\mathrm{T}}\right]\right)\right. \\
&+\frac{1}{2 \sigma^{2}}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}\right\|^{2}-\frac{1}{\sigma^{2}} \mathbb{E}\left[\mathbf{z}_{n}\right]^{\mathrm{T}} \mathbf{W}^{\mathrm{T}}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right) \\
&+\frac{1}{2 \sigma^{2}} \operatorname{Tr}\left(\mathbb{E}\left[\mathbf{z}_{n} \mathbf{z}_{n}^{\mathrm{T}}\right] \mathbf{W}^{\mathrm{T}} \mathbf{W}\right)+\left.\frac{M}{2} \ln (2 \pi)\right\}
\end{aligned} \tag{12.53}
で与えられる完全データの対数尤度関数の期待値を最大化することにより、確率的主成分分析モデルのMステップの式
\mathbf{W}_{\text {new }}=\left[\sum_{n=1}^{N}\left(\mathbf{x}_{n}-\overline{\mathbf{x}}\right) \mathbb{E}\left[\mathbf{z}_{n}\right]^{\mathrm{T}}\right]\left[\sum_{n=1}^{N} \mathbb{E}\left[\mathbf{z}_{n} \mathbf{z}_{n}^{\mathrm{T}}\right]\right]^{-1} \tag{12.56}
と
\begin{aligned} \sigma_{\text {new }}^{2}=& \frac{1}{N D} \sum_{n=1}^{N}\left\{\left\|\mathbf{x}_{n}-\overline{\mathbf{x}}\right\|^{2}-2 \mathbb{E}\left[\mathbf{z}_{n}\right]^{\mathrm{T}} \mathbf{W}_{\text {new }}^{\mathrm{T}}\left(\mathbf{x}_{n}-\overline{\mathbf{x}}\right)\right.\\ &\left.+\operatorname{Tr}\left(\mathbb{E}\left[\mathbf{z}_{n} \mathbf{z}_{n}^{\mathrm{T}}\right] \mathbf{W}_{\text {new }}^{\mathrm{T}} \mathbf{W}_{\text {new }}\right)\right\} \end{aligned} \tag{12.57}
を導け．
Mステップの更新式を求めるためには、(12.53)式を\mathbf{W}及び\sigma^2について微分する必要がある。
まず、
\begin{aligned}
\frac{\partial}{\partial \mathbf{W}}\mathbb{E}[\ln p (\mathbf{X}, \mathbf{Z}|\mathbf{\mu}, \mathbf{W}, \sigma^2)] &= \sum_{n = 1}^N \left\{\frac{1}{\sigma^2}(\mathbf{x}_n - \mathbf{\mu})\mathbb{E}[\mathbf{z}_n]^T - \frac{1}{\sigma^2}\mathbf{W}\mathbb{E}[\mathbf{z}_n\mathbf{z}_n^T]\right\}.
\end{aligned}
ここではThe Matrix Cookbookの(71)式及びPRML(C.24)式を用いた。
また、
\begin{aligned}
\frac{\partial}{\partial \mathbf{\sigma^2}}\mathbb{E}[\ln p (\mathbf{X}, \mathbf{Z}|\mathbf{\mu}, \mathbf{W}, \sigma^2)] &= \sum_{n = 1}^N \left\{\frac{D}{2\sigma^2} - \frac{1}{2\sigma^4}\|\mathbf{x}_n - \mathbf{\mu}\|^2 + \frac{1}{\sigma^4}\mathbb{E}[\mathbf{z}_n]\mathbf{W}^T(\mathbf{x}_n - \mathbf{\mu}) - \frac{1}{2\sigma^4}\mathrm{Tr}(\mathbb{E}[\mathbf{z}_n\mathbf{z}_n^T]\mathbf{W}^T\mathbf{W}) \right\}.
\end{aligned}
これらの式を0とおき、\mathbf{\mu} = \bar{\mathbf{x}}を代入したうえで方程式を解くと(12.56)式及び(12.57)式を得る。
Discussion