\newcommand{\vf}{\mathbf{f}}
\newcommand{\vg}{\mathbf{g}}
\newcommand{\vj}{\mathbf{j}}
\newcommand{\vu}{\mathbf{u}}
\newcommand{\vw}{\mathbf{w}}
\newcommand{\vx}{\mathbf{x}}
\newcommand{\vy}{\mathbf{y}}
\newcommand{\vz}{\mathbf{z}}
\newcommand{\vmu}{\boldsymbol{\mu}}
\newcommand{\vzero}{\mathbf{0}}
\newcommand{\vone}{\mathbf{1}}
\newcommand{\vwi}[1]{\mathbf{w}_{(#1)}}
\newcommand{\Ma}{\mathbf{A}}
\newcommand{\Mb}{\mathbf{B}}
\newcommand{\Mc}{\mathbf{C}}
\newcommand{\Ml}{\mathbf{L}}
\newcommand{\Mm}{\mathbf{M}}
\newcommand{\Mo}{\mathbf{O}}
\newcommand{\Mr}{\mathbf{R}}
\newcommand{\Ms}{\mathbf{S}}
\newcommand{\Mu}{\mathbf{U}}
\newcommand{\Mk}{\mathbf{K}}
\newcommand{\Mv}{\mathbf{V}}
\newcommand{\Mw}{\mathbf{W}}
\newcommand{\Mx}{\mathbf{X}}
\newcommand{\Mz}{\mathbf{Z}}
\newcommand{\Mgamma}{\boldsymbol{\Gamma}}
\newcommand{\Msigma}{\boldsymbol{\Sigma}}
\newcommand{\Momega}{\boldsymbol{\Omega}}
\newcommand{\Mi}{\mathbf{I}}
\newcommand{\diag}{\mathrm{diag}}
\newcommand{\norm}{\mathcal{N}}
\newcommand{\tr}[1]{#1 ^ \mathrm{T}}
\newcommand{\trace}[1]{\mathrm{Tr} \left[ #1 \right]}
\newcommand{\pderiv}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\lik}{\mathcal{L}}
\newcommand{\vecsize}[1]{\left\| #1 \right\|}
この記事は、統計学・機械学習の教科書である C. M. Bishop 著「パターン認識と機械学習」（略称：PRML）において、式変形の行間を独自に埋めたものです。
前提：確率的主成分分析の概説
確率的主成分分析では、長さ $D$ の観測変数ベクトル $\vx$ を、長さ $M$ の潜在変数ベクトル $\vz$ で近似する。まず、潜在変数の事前分布にはガウス分布
p(\vz) = \norm (\vz | \vzero, \Mi)
\tag{12.31}
を与える。そして、潜在変数に $D \times M$ 行列 $\Mw$ で表される因子負荷量をかけることで観測変数が近似できるとする。実際の観測変数は、近似値にノイズを加えたものとして説明される。
p(\vx | \vz) = \norm (\vx | \Mw \vz + \vmu, \sigma^2 \Mi)
\tag{12.32}
この記事で解説すること
以上のモデルにおいて、尤度を最大化するパラメータ $\vmu, \Mw, \sigma^2$ を求める。まず、$\vmu$ の最尤解はデータ平均 $\overline{\vx}$ であることがすぐにわかる（本文参照）。これを代入すると、対数尤度関数は
\ln p(\vx | \Mw, \vmu, \sigma^2)
=
- \frac{N}{2}
\left\{
D \ln(2 \pi) + \ln |\Mc| + \trace{\Mc^{-1} \Ms}
\right\}
\tag{12.44}
と表される。なお、$N$ はデータ数、$\mathbf{C} = \Mw \Mw^{\mathrm{T}} + \sigma^2 \Mi$、$\Ms$ はデータの共分散行列である。これをもとに $\Mw$ の最尤解を求めたいのだが、その過程で対数尤度の微分が $0$ となるような $\Mw$ は
\Mw
=
\Mu_M (\Ml_M - \sigma^2 \Mi) ^ {1/2} \Mr
\tag{12.45}
を必ず満たすことが紹介されている。なお、$\Mu_M$ は $D \times M$ 行列で、その列ベクトルは $\Ms$ の固有ベクトルの任意の部分集合で与えられる。$\Ml_M$は、固有ベクトルに対応する固有値 $\lambda_m \ (m=1,...,M)$ 対角要素に持つ $M \times M$ 対角行列である。$\Mr$は任意の $M \times M$ 直交行列である。
この記事の目的は、式 (12.44) と (12.45) の間の計算過程を埋めることである。
参考文献
本文では、PRMLの著者自身が確率的主成分分析を提唱した論文 (Tipping and Bishop, 1999) が引用されており、この記事の対象とする部分についても Appendix A.1 に解説されている。元論文の記述で理解できる方なら、そちらを参照したほうが正確だろう。
ただ私は行間が理解できなかったので、こちらのブログ記事（https://masamunetogetoge.com/ppca）を参考にした。$\Mw$ の最尤解の導出についても、元論文より丁寧に解説されているが、行列計算の知識の乏しい私には計算過程を追うことができなかった。
そこで、「The Matrix Cookbook」という行列計算の公式集を参考に、自力で計算してみたのがこの記事である。
解説
対数尤度を $\Mw$ で偏微分し、それを $\Mo$ とおくと、式 (12.44) より
\pderiv{}{\Mw} \ln |\mathbf{C}|
+ \pderiv{}{\Mw} \trace{\mathbf{C}^{-1} \mathbf{S}}
= \Mo
\tag{1}
となる。これ以降、
- 式 (1) の第二項を計算する
- 式 (1) の第一項を計算する
- 以上の結果を合わせる
- $\Mw$ を固有値分解する
という順で進んでいく。
式 (1) の第二項を計算する
まずは第二項にとりかかろう。とはいえそのままだと難しいので、$\Mw$ の $(d,m)$ 成分 $W_{d,m}$ に注目して式を展開する。一般に成り立つ行列の公式
\begin{align}
&\pderiv{}{X} \trace{\Ma} = \trace{\pderiv{\Ma}{X}}
\\
&\pderiv{\Ma^{-1}}{X} = -\Ma^{-1} \pderiv{\Ma}{X} \Ma^{-1}
\end{align}
を使うと以下のようになる。
\begin{align}
\pderiv{}{W_{d,m}} \trace{\Mc^{-1} \Ms}
&=
\trace{\pderiv{}{W_{d,m}} (\Mc^{-1} \Ms)}
\\
&=
\trace{ \left( \pderiv{\Mc^{-1}}{W_{d,m}} \right) \Ms}
\\
&=
- \trace{ \Mc^{-1} \pderiv{\Mc}{W_{d,m}} \Mc^{-1} \Ms}
\tag{2}
\end{align}
ここで $\partial \Mc / \partial W_{d,m}$ に注目する。$\Mc = \Mw \tr{\Mw} + \sigma^2 \mathrm{I}$ なので、$\Mw \tr{\Mw}$ を $W_{d,m}$ で微分すればよい。$\Mw \tr{\Mw}$ は $D \times D$ 行列で、その $(d_1, d_2)$ 成分は
\sum_{m^{'}=1}^M W_{d_1, m^{'}} W_{d_2, m^{'}}
で表されるから、$\partial \Mc / \partial W_{d,m}$ の $(d_1, d_2)$ 成分は
- $d_1 = d, d_2 = d$ のとき、$2W_{d,m}$
- $d_1 \ne d, d_2 = d$ のとき、$W_{d,m}$
- $d_1 = d, d_2 \ne d$ のとき、$W_{d,m}$
- $d_1 \ne d, d_2 \ne d$ のとき、$0$
つまり、$d$ 行目と $d$ 列目にのみ値が入った行列となる。表現を変えると、$d$ 列目の要素が $W_{1,m},...,W_{D,m}$、それ以外の要素が $0$ の $D \times D$ 行列 $\Momega$ を用いて
\pderiv{\Mc}{W_{d,m}} = \Momega + \tr{\Momega}
\tag{3}
と表せる。従って、
\begin{align}
\pderiv{}{W_{d,m}} \trace{\Mc^{-1} \Ms}
&=
- \trace{ \Mc^{-1} (\Momega + \tr{\Momega}) \Mc^{-1} \Ms}
\\
&=
- \trace{ \Mc^{-1} \Momega \Mc^{-1} \Ms}
- \trace{ \Mc^{-1} \tr{\Momega} \Mc^{-1} \Ms}
\\
&=
- 2 \trace{ \Mc^{-1} \Momega \Mc^{-1} \Ms}
\\
&=
- 2 \trace{ \Mc^{-1} \Ms \Mc^{-1} \Momega}
\tag{4}
\end{align}
が得られる。
式 (1) の第一項を計算する
次に、式 (1) の第一項 $\partial \ln{|\Mc|} / \partial \Mw$ の計算に進む。微分の対象は $\Mc$ の関数で、$\Mc$ が $\Mw$ の関数という階層構造なので、連鎖律 (chain rule) を適用したい。ただ、行列の微分における連鎖律は簡単に書けず、要素ごとに計算することになる。一般に、関数 $f(\Ma)$ を行列 $\Mx$ で微分した結果は、 $\Mx$ の $(i, j)$ 成分 $X_{ij}$ ごとに書ける。
\pderiv{f(\Ma)}{X_{ij}}
=
\trace{
\tr{\left(\pderiv{f(\Ma)}{\Ma}\right)}
\pderiv{\Ma}{X_{ij}}
}
従って、$\Mw$ の $(d,m)$ 成分 $W_{d,m}$ において
\pderiv{}{W_{d,m}} \ln |\mathbf{C}|
=
\trace{
\tr{\left(\pderiv{}{\Mc} \ln |\mathbf{C}| \right)}
\pderiv{\Mc}{W_{d,m}}
}
\tag{5}
となる。$\trace{\ }$ の中の前半は公式
\pderiv{\ln|\Ma|}{\Ma} = {\left( \Ma^{-1} \right)}^{\mathrm{T}}
が適用でき、後半には式 (3) が代入できるので、式 (5) は以下のように計算できる。
\begin{align}
\pderiv{}{W_{d,m}} \ln |\mathbf{C}|
&=
\trace{\Mc^{-1} (\Momega + \tr{\Momega})}
\\
&=
\trace{\Mc^{-1} \Momega} + \trace{\Mc^{-1} \tr{\Momega}}
\\
&=
2 \trace{\Mc^{-1} \Momega}
\tag{6}
\end{align}
以上の結果を合わせる
従って、式 (1) の $(d,m)$ 成分において
\begin{align}
- 2 \trace{ \Mc^{-1} \Ms \Mc^{-1} \Momega} + 2 \trace{\Mc^{-1} \Momega} &= 0
\\
\therefore \
\trace{ \Mc^{-1} \Ms \Mc^{-1} \Momega} &= \trace{\Mc^{-1} \Momega}
\tag{7}
\end{align}
が成り立つことが分かった。よく見ると、両辺とも $\trace{\Ma \Momega}$ の形をしている。$\Momega$ は $d$ 列目以外の要素が $0$ なので、$\Ma \Momega$ の対角成分は $d$ 個目以外 $0$ である。従ってこれは
\trace{\Ma \Momega}
= \sum_{d^{'}=1}^D A_{d,d^{'}} W_{d^{'},m}
= \{ \Ma \Mw \}_{d,m}
すなわち、行列 $\Ma \Mw$ の $(d,m)$ 成分であることが分かる。よって、式 (7) の $d=1,...,D, \ \ m=1,...,M$ の場合を行列形式で並べることができ、
\begin{align}
\Mc^{-1} \Ms \Mc^{-1} \Mw &= \Mc^{-1} \Mw
\\
\therefore \ \Ms \Mc^{-1} \Mw &= \Mw
\tag{8}
\end{align}
という単純な形で書くことができる。この式を満たす $\Mw$ を探せばよい。
Wを固有値分解する
式 (8) の自明な解として $\Mw = \Mo$ や $\Ms = \Mc$ があるが、それ以外の場合を考える。そこで、$\Mw$ を固有値分解し、
\Mw = \Mu \Ml \tr{\Mv} \tag{9}
と書く。$\Mu$ は $D \times D$ 直交行列、$\Mv$ は $M \times D$ 直交行列である。$\Ml$ は $D \times D$ 対角行列であり、対角成分の $M$ 番目までの要素は $\Mw$ の特異値 ($l_1, ..., l_M$)、それ以降の要素は $0$ である 。以下では、式 (8) を満たす $\Mw$ の代わりに、式 (8) を満たす $\Mu, \Ml, \Mv$ について考える。
式 (9) を $\mathbf{C} = \Mw \Mw^{\mathrm{T}} + \sigma^2 \Mi$ に代入すると、直交行列の性質 $\Mu \tr{\Mu} = \Mi, \Mv \tr{\Mv} = \Mi$ から
\begin{align}
\Mc
&= \Mu \Ml \tr{\Mv} \Mv \Ml \tr{\Mu} + \sigma^2 \Mi
\\
&= \Mu \Ml^2 \tr{\Mu} + \Mu (\sigma^2 \Mi) \tr{\Mu}
\\
&= \Mu (\Ml^2 + \sigma^2 \Mi) \tr{\Mu}
\tag{10}
\end{align}
となる。一般に成り立つ公式 $(\Ma \Mb)^{-1} = \Mb^{-1} \Ma^{-1}$ 、および直交行列の性質 $\tr{\Mu} = \Mu^{-1}$ を用いると、
\begin{align}
\Mc^{-1}
&= (\tr{\Mu})^{-1} (\Ml^2 + \sigma^2 \Mi)^{-1} \Mu^{-1}
\\
&= \Mu (\Ml^2 + \sigma^2 \Mi)^{-1} \tr{\Mu}
\tag{11}
\end{align}
となり、これらを式 (8) に代入すると
\begin{align}
\Ms \Mu (\Ml^2 + \sigma^2 \Mi)^{-1} \tr{\Mu} \Mu \Ml \tr{\Mv}
&= \Mu \Ml \tr{\Mv}
\\
\Ms \Mu (\Ml^2 + \sigma^2 \Mi)^{-1} \Ml \tr{\Mv}
&= \Mu \Ml \tr{\Mv}
\\
\Ms \Mu (\Ml^2 + \sigma^2 \Mi)^{-1}
&= \Mu
\\
\Ms \Mu &= \Mu (\Ml^2 + \sigma^2 \Mi)
\tag{12}
\end{align}
となる。これを $\Mu$ の列ベクトル $\vu_d \ (d=1,...,D)$ それぞれについて見ると、
- $d \leq M$ のとき $\Ms \vu_d = (l_d^2 + \sigma^2) \vu_d$
- $d > M$ のとき $\Ms \vu_d = \sigma^2 \vu_d$
となっている。これは$\vu_d$ が $\Ms$ の固有ベクトルであり、かつ、$l_m \ (m = 1, ..., M)$ が対応する $\Ms$ の固有値 $\lambda_m$ を用いて
l_m = (\lambda_m - \sigma^2)^{1/2}
と書けることを意味する。以上が式 (8) を満たす $\Mu, \Ml, \Mv$ の性質である。なお、$\Mv$ については特に条件が得られないので、任意の直交行列を与えればよい。
結果をまとめると、対数尤度を最大化する $\Mw$ は
\Mw
=
\Mu_M (\Ml_M - \sigma^2 \Mi) ^ {1/2} \Mr
\tag{12.45}
という形で書ける。ここで、$\Mu_M$ は $D \times M$ 行列で、その列ベクトルは $\Ms$ の固有ベクトルの任意の部分集合で与えられる。$\Ml_M$は、固有ベクトルに対応する固有値 $\lambda_m \ (m=1,...,M)$ 対角要素に持つ $M \times M$ 対角行列である。$\Mr$は任意の $M \times M$ 直交行列である。