確率論的主成分分析(PPCA)の話

マサムネのイメージ図 ベイズ統計学

主成分分析(PCA)を、確率論的に行う記事です。普通のPCAはデータの次元を削減したり、クラスタリングしたりに使用します。一方、確率論的主成分分析(PPCA)は、次元を削減するだけでなく、データを生成する事も出来ます。その意味で、PPCAは生成モデルになっています。
PCAの復習から初めて、PPCAの問題設定を解説し、最尤法でパラメーターを推定します。

主成分分析

普通のPCAを復習しておきましょう。詳しくは、解説記事を読んでください。

主成分分析からカーネル主成分分析へ

データの次元を削減する事を目標として計算をしていきます。
\(D \)次元のデータが\(N\)個から、\(D \times N \)行列Xが出来ているとします。各データの平均で正規化した行列\( \tilde{X} \)を考えます。
$$\begin{eqnarray}
\mu &=& \frac{1}{n} \sum \vec{x_i } \\
\tilde{X} &=& (\vec{x_1 } -\mu , \cdots , \vec{x_n } -\mu _m )
\end{eqnarray} $$
\(\tilde{X} \)を正方行列にするために、転置したものと掛け合わせ、固有値と固有ベクトルを求めます。
$$\begin{eqnarray}
\tilde{X} \tilde{X} ^{T} u_i &=& \lambda _i u_i \\
\tilde{X}^{T} \tilde{X} v_i &=& \mu _i v_i
\end{eqnarray} $$
\(\tilde{X} ^T \tilde{X} \)の対称性から、\( \lambda , \mu , u , v\)には関係があり、結局以下の等式が成り立つ事が分かります。
$$\begin{eqnarray}
\tilde{X} ^{T} u_i&=& \sqrt{\lambda _i } v_i \\
\tilde{X} v_i &=&\sqrt{ \lambda _i } u_i
\end{eqnarray} $$
ここで、\( \tilde{X} ^T \tilde{X} \)の相異なる固有値の数を\(r \)とすると、\(u_1 , \cdots , u_r \)という正規直交基底が得られます。正規直交基底を適当に延長する事で、正規直交基底\(u_1 , \cdots , u_r , u_{r+1} , \cdots , u_{D} \)が得られます。
\(V_D= (v_1, \cdots , v_{D} ) \)も作っておくことで、\(X \)を分解できます。
$$\begin{eqnarray}
\tilde{X} = U_D \Lambda ^{1/2} V_D ^{T}
\end{eqnarray}$$
ただし、\( \Lambda ^{1/2} \)は\( \tilde{X} ^T \tilde{X} \)の固有値の1/2乗が対角成分に並んだ対角行列です。この分解を特異値分解と呼びます。
\(U_D \)に属する正規直交基底を使ってデータの次元を圧縮する事が出来ます。
具体的には、圧縮したい次元\( k\)決めて、\(U_k = ( u_1 , \cdots , u_k ) \)を作り、データ\( x \)に対して
$$\begin{eqnarray}
X_k = U_k^{T} (x – \mu )
\end{eqnarray} $$
とします。これが、普通のPCAです。D次元のデータを、任意の次元のデータに圧縮する事が出来ます。データの情報は固有値が全て持っていると思い、
$$\begin{eqnarray}
L_ k = \frac{\sum_{i=1} ^k \lambda _k}{\sum_{j=1}^D \lambda _j }
\end{eqnarray}$$
が80%を超えるくらいのkを選ぶ事が良くあります。

確率論的主成分分析

PCAを確率論的に行います。隠れ変数を導入し、データはその隠れ変数から生成されていると考えます。
つまり、D次元のデータ\(x \)はM次元の隠れ変数\( z \sim \mathcal{N} (z| 0, I_{M} ) \)と1 、ノイズ\( \epsilon \sim \mathcal{N} (\epsilon | 0, \sigma ^2 I_{D} ) \)によって生成されていると考えます。2
$$\begin{eqnarray}
x = Wz + \mu + \epsilon
\end{eqnarray} $$
\(W \)は\(D\times M \)行列で、\( \mu \)は\(D \)次元の定数ベクトルです。これから、
$$\begin{eqnarray}
p(x ) &=& \mathcal{N}(x | \mu , \Sigma ) \\
\Sigma &=& W W^T + \sigma ^2 I_D
\end{eqnarray}$$
となります。
ここで、\(x \)を生成するには、\( \Sigma ^{-1} \)を計算しなくてはいけないのですが、ウッドベリーの公式から、
$$\begin{eqnarray}
\Sigma ^{-1} &=& \sigma ^{-2} I_{D} -\sigma ^{-2} WK^{-1} W^T \\
K&=& W^T W +\sigma ^2 I_{M}
\end{eqnarray}$$
と書けることに注意しましょう。この計算は、\( M\times M \)行列の逆行列の計算が一番大変なわけですが、これは\( D\times D \)行列の逆行列の計算よりは軽くなっています。
隠れ変数が従う確率分布を出すために、\(p(x|z ) \)が必要です。\(x \)と\(z \)の関係から、\(z \)を定数だと思う事で
$$\begin{eqnarray}
p(x|z ) &=& \mathcal{N}(x | Wz+\mu , \sigma ^2 I_{D} )
\end{eqnarray}$$
となります。これから、ベイズの定理を使う事で、隠れ変数が従う確率分布は、
$$\begin{eqnarray}
K&=& W^T W +\sigma ^2 I_{M} \\
p(z|x ) &=& \mathcal{N}(z |K^{-1} W^T (x-\mu ) , \sigma^{-2} K^{-1} )
\end{eqnarray}$$
となります。計算は、正規分布\( \mathcal{N} (x|\mu , \Sigma ) \) の\(\exp \)の肩の部分が
$$\begin{eqnarray}
-\frac{ (x-\mu)^T \Sigma ^{-1}(x- \mu ) }{2}
\end{eqnarray}$$
となる事を利用して、\(p(x|z) p(z) \)の\( \exp \)の肩の部分をまとめると出て来ます。
確率論的に考える事のメリットは、隠れ変数\(z \)をサンプルすることで、\(p(x|z ) \)からデータを生成出来る事です。この意味で、PPCAは生成モデルになっています。

最尤法によるパラメーターの決定

Mは適当な数字を与えているとして、\( \mu , \sigma ^2 , W \)を最尤推定で決めます。
初めに、全データの対数尤度は以下のようになります。
$$\begin{eqnarray}
\log p(X|\mu, W, \sigma ^2 ) &=& \sum \log p(x_n |\mu, W, \sigma ^2 ) \\
&=& -\frac{ND}{2} \log \pi -\frac{N}{2}\log |\Sigma | -\frac{1}{2} \sum (x_n – \mu ) \Sigma ^{-1} (x_n – \mu )
\end{eqnarray}$$
\( \mu \)で偏微分して、\( \mu \)の最尤推定値を求めるのは簡単です。
$$\begin{eqnarray}
\hat{ \mu } = \frac{x_n}{N}
\end{eqnarray}$$
\( S= \sum (x_n – \hat{ \mu } ) (x_n – \hat{ \mu } ) ^T \)として、\(\hat{\mu} \)を代入すると、対数尤度は以下のように書き直せます。
$$\begin{eqnarray}
L= \log p(X|\mu, W, \sigma ^2 ) =
-\frac{ND}{2} \log \pi -\frac{N}{2}\log |\Sigma | -\frac{1}{2} Tr (\Sigma ^{-1} S)
\end{eqnarray}$$
次に、\(W , \sigma ^2 \)の最尤推定値を求めましょう。3

Wの最尤推定

対数尤度をWで偏微分します。
行列の微分には以下の公式を使います。
$$\begin{eqnarray}
\frac{\partial Tr(ABA^T ) }{\partial A} &=& A(B+ B^T ) \\
\frac{\partial Tr(AB) }{\partial A} &=&B^T \\
\frac{\partial \log |A| }{\partial A} &=& (A^{-1} )^{T}
\end{eqnarray}$$
これらを用いて、
$$\begin{eqnarray}
\frac{\partial L }{\partial W } =N(\Sigma ^{-1} S\Sigma ^{-1} W -\Sigma ^{-1} W )
\end{eqnarray}$$
です。上の式を0と置く事で、
$$\begin{eqnarray}
S\Sigma ^{-1} W =W
\end{eqnarray}$$
です。W=0や、\( S =\Sigma \)の時以外の\(W \)を求めます。そのために、普通のPCAのように、Wを特異値分解します。
$$\begin{eqnarray}
W= U \Lambda ^{1/2} V^T
\end{eqnarray}$$
ここで、\(U \)と\(\Lambda ^{1/2} \)は\(D \times D \)行列、\(V \)は\(M\times D \)行列です。注意なのは、\( \Lambda \)は対角行列で、\( m \)列目までは0以外の値が入り得ますが、\( m+1 \)列目以降は必ず0になっているという点です。\( \Sigma = WW^T + \sigma ^2 I_M \)なので、特異値分解したWを代入すると
$$\begin{eqnarray}
\Sigma = U (\Lambda +\sigma ^2 )U^T
\end{eqnarray}$$
となります。これを \( S\Sigma ^{-1} W =W \)に代入して整理すると4
$$\begin{eqnarray}
SU = (\Lambda + \sigma ^2)U
\end{eqnarray}$$
となります。これは行列としての等式ですが、Uの列成分毎に見ると、
$$\begin{eqnarray}
Su_j = (\Lambda + \sigma ^2)u_j
\end{eqnarray}$$
となり、初めの等式から、Sに関する固有値方程式が出てくることが分かります。\(u_j \)たちは、Wに付随する固有ベクトルでしたが、式\( S\Sigma ^{-1} W =W \) によって、Sの固有ベクトルにもなりました。5Sの固有値を大きい順に\(s_1 , \cdots , s_D \)と並べると、固有値についての等式
$$\begin{eqnarray}
s_j = \lambda _j + \sigma ^2
\end{eqnarray}$$
が成り立てば、\( S\Sigma ^{-1} W =W \)が成り立ちます。つまり、固有値が
$$\begin{eqnarray}
\lambda _j =
\begin{cases}
s_j – \sigma ^2 & j\leq m \\
0 & j>m
\end{cases}
\end{eqnarray}$$
となるような行列\(WW^T \)ならなんでも良いという事です。ただし、mは、非ゼロな\(\lambda _j \)の個数です。
この行列は、一般に対角行列Kと、\(D\times M \)行列Rで、
$$\begin{eqnarray}
W_{ML} &=& U(K- \sigma ^2 I_{D} )R
\end{eqnarray}$$
と書けます。ただし、
$$\begin{eqnarray}
{\rm diag} K&=&
\begin{cases}
s_j & j\leq m \\
\sigma ^2 & j>m
\end{cases} \\
R^T R &=& I_{M} \\
R^T R&=&I_{D}
\end{eqnarray}$$
です。Rは任意とは言え、どうやって取るんだという疑問がありますが、特異値分解で出てくる\(V \)を使えば良いです。6使いやすい形としては、
$$\begin{eqnarray}
W_{ML} = U(K- \sigma ^2 I_{D} )V^T
\end{eqnarray}$$
となります。また、mはWの特異値のうち、0出ないものの数でしたが、特異値はSの固有値と共有されているので、m=Mとしても良いことになります。つまり、まとめると
$$\begin{eqnarray}
W_{ML} &=& U(K- \sigma ^2 I_{D} )V^T \\
{\rm diag} K&=&
\begin{cases}
s_j & j\leq M \\
\sigma ^2 & j>M
\end{cases} \\
U &=&(u_1 , \cdots , u_D ) \\
Su_j &=& s_j u_j
\end{eqnarray}$$
となります。

\( \sigma ^2\)の最尤推定

\(W \)が最尤推定で決まっているとすると、\(\Sigma \)の固有値が分かっているので、\( \log |\Sigma | \)の値も書き下す事が出来ます。
$$\begin{eqnarray}
\log | \Sigma | = \sum _{j=1} ^{M} \log (\lambda _j + \sigma ^2 ) + \sum_{j=M+1}^D \log \sigma ^2
\end{eqnarray}$$
対数尤度の式は、
$$\begin{eqnarray}
L = -\frac{N}{2} ( D\log (2 \pi ) +\sum_{i=1} ^M \log s_j +\frac{1}{\sigma ^2 } \sum_{j=M+1} ^D s_j +(D-M)\log ( \sigma ^2 ) +M )
\end{eqnarray}$$
となります。
これから、\( \sigma ^2 \)の最尤推定値は
$$\begin{eqnarray}
\sigma ^2 _{ML} = \frac{1}{D-M} \sum_{j=M+1} ^{D}s_j
\end{eqnarray}$$
です。

このようにパラメーターを指定する事で\(x\)や\(z \)をサンプルすることが出来ます。また、M<D を適当に決める事で、次元の削減が出来ます。別の記事で、実際にPPCAでデータを生成してみます。

まとめ

  • PCAの説明をした
  • PPCAの説明をした
  • 最尤法の計算を説明した
  1. M<D としておけば、データの次元を削減できます。
  2. \(I_k \)はk次元の単位行列を表しています。
  3. 元論文のリンクを貼っておきます。
  4. \(U^T =U^{-1} \)に注意です。
  5. j>mでは、Sの固有値=\(\sigma ^2 \)となることに注意しましょう。
  6. \(D\times M \)行列を適当に生成して、SVD で出てくるVを使ったりしても良いです。
タイトルとURLをコピーしました