# 機械学習初心者が確率的主成分分析を頑張って理解しようとしてみた
この記事は[古川研究室 Workout_calendar](https://qiita.com/flab5420/private/aee38a24be700d146817) 16日目の記事です。
本記事は古川研究室の学生が学習の一環として書いたものです。内容が曖昧であったり表現が多少異なったりする場合があります。

# はじめに
私は記事の一本目に線形代数の重要な要素である固有値、二本目にディープラーニングの応用例であるVAEを用いたマッチョVAEを執筆したわけですが、VAEの論文を読んでいるときに確率への理解が乏しいことに気づきました。そこで確率力(造語)を高めるべく、今まで執筆した記事の内容とシナジーのありそうな確率的主成分分析について記事を書こうと思いいたりました。
この記事の目的は自分で確率的主成分分析を噛み砕きながら理解することで、確率への理解を深めることです(純度100%のエゴ)。

# 確率的主成分分析とは
確率的主成分分析は、主成分分析を確率的に解釈しなおしたものです。
主成分分析で用いる値(共分散行列とか潜在空間とか、入力データとか)を確率的に表現することによって、最尤推定とか、EMアルゴリズムとか、ベイズ推定によってパラメータを決めれるようになったことが大きな利点っぽいです。

確率的主成分分析では、観測データ$\mathbf{x}\in \mathbb{R}^D$が、条件付分布$p(\mathbf{x}|\mathbf{z})$で生成されると仮定します。ここで、$\mathbf{x}$の条件である$\mathbf{z}$を潜在変数と呼びます。($\mathbf{z}\in\mathbb{R}^M$)。
潜在変数の次元は観測データの次元より低い$M<D$とします。 潜在変数$\mathbf{z}$の事前分布は、平均が0で共分散行列が単位行列のガウス分布

$$
p(z)=N(z|0,I)
$$

で与えられると仮定します。ここで共分散行列が単位行列ということは、お互いの次元の共分散が0であることを表し、潜在変数$z$の各次元がお互いに独立であることを意味します。
求めたい条件付分布$p(\mathbf{x}|\mathbf{z})$は$\mathbf{z}$によって条件つけられた、ガウス分布だと仮定します。

$$
p(\mathbf{x}|\mathbf{z})=N(\mathbf{x}|W\mathbf{z}+\mu_b,\sigma ^2 I)
$$
と表します。ここで行列$W$は$D×M$次元の行列で、$M$次元の$\mathbf{z}$と行列積を取ることにより$D$次元となります。そこに$\mathbf{x}$の各次元の平均$\mu_b\in \mathbb{R}^D$を加えることで$D$次元の平均としています。
また、$\mathbf{x}$の各次元の分散は、パラメータ$\sigma ^2$によってのみ決まっています。

この式を生成モデル的に見ていくと、$D$次元の観測データ$\mathbf{x}$は、$M$次元の潜在変数$\mathbf{z}$を$W$と$\mu_b$で線形変換したものに、ガウス分布のノイズ（大きさは$\sigma ^2$によって決まる)が加えられたものと考えることができます。
$$
\mathbf{x}=W\mathbf{z}+\mu_b+\epsilon
$$

ここで$\epsilon$は$D$次元の平均0で共分散$\sigma ^2 I$のガウスノイズを表しています。つまり我々は、潜在変数と観測データ$\mathbf{x}$を繋いでいる$W$,$\mu_b$,$\sigma ^2$を求めればいいわけです。

# 予測分布の導出
最尤法などでパラメータを決定するためには、事後分布$p(\mathbf{z|x})$を導出しなければなりません。事後分布を導出するためには、観測変数の周辺分布$p(\mathbf{x})$を表す必要があります。周辺分布$p(\mathbf{x})$は確率の加法定理、乗法定理から
$$
p(\mathbf{x}) = \int p(\mathbf{x|z})p(\mathbf{z})d\mathbf{z}
$$
となります。
$p(\mathbf{x|z}),p(\mathbf{z})$共にガウス分布ですので、積分の結果もガウス分布となります。
よって$p(\mathbf{x})$は次のように表されます。
$$
p(\mathbf{x})=N(\mathbf{x}|\mathbf{\mu} ,C) , C=WW^T+\sigma ^2 I
$$
$p(\mathbf{x})$のガウス分布を定めるパラメータである$C,\mu$は以下に示す周辺分布の公式より導くことができます。
周辺分布の公式
>xの周辺ガウス分布とxが与えられた時のyの条件付ガウス分布が以下のように表されるとき
>\begin{equation}
>p(x)=N(x|\mu,\Lambda ^{-1})
>\end{equation}
>\begin{equation}
p(y|x)=N(y|Ax+b,L^{-1})
\end{equation}
yの周辺分布は
\begin{equation}
p(y)=N(y|A\mu+b,L^{-1}+A\Lambda ^{-1}A^T)
\end{equation}
で与えられる

今回は$p(\mathbf{z})=N(\mathbf{z}|0,I)$,$p(\mathbf{x|z})=N(\mathbf{x|Wz+\mu_b},\sigma ^2 I)$なので、各パラメータを上の公式に当てはめると$p(x)$を導けます。公式があるなんてびっくり。

さて、新しいデータ点$\mathbf{x}^*$に対しての予測分布$p(\mathbf{x}^*|\mathbf{X})$を導出します。(PRMLではなぜか周辺分布と同じ$p(\mathbf{x})$と表記されている。)
予測分布を求めるためには周辺分布を示す共分散行列$C$の逆数$C^{-1}$が必要になりますが、これは$D×D$次元の計算を必要とするので、$D$が大きい場合、非常に計算コストがかかります。そこで以下に示すウッドべリーの公式を用いて変換します。
>Woodburyの公式
>$$
>(A+BD^{-1}C)^{-1} = A^{-1}-A^{-1}B(D+CA^{-1}B)^{-1}CA^{-1}
>$$

$$
C^{-1}=\sigma ^{-2}I-\sigma ^{-2}WM^{-1}W^T
$$

ここで$M×M$行列$M$は
$$
M=W^TW+\sigma ^{2}I
$$
で表されます。

※Woodburyによる導出と事後分布を導出する意味を理解し次第、以下に事後分布$p(\mathbf{z|x})$の導出を記述する予定
# 最尤法による確率的主成分分析
前章により、$p(\mathbf{x})$の周辺分布が分かったため、それを尤度関数として、最尤法を適用します。周辺分布$p(\mathbf{x})$は

$$
p(\mathbf{x})=N(\mathbf{x}|\mu ,C) , C=WW^T+\sigma ^2 I
$$

$\mathbf{x}$は$D$次元の多変量ですから、多変量ガウス分布で上記を書き下します。

多変量ガウス分布は

$$
f(\mathbf{x})=\frac{1}{(\sqrt{2\pi})^D \sqrt{|\Sigma}|} \exp{(-\frac{1}{2}(\mathbf{x}-\mu)^T \Sigma ^{-1}(\mathbf{x}-\mu))}
$$

で表されます。Σはデータの分散共分散行列です。
なので$p(\mathbf{x})$の各パラメータ$\mu,C$を代入すると

$$
f(\mathbf{x})=\frac{1}{(\sqrt{2\pi})^D \sqrt{|C}|} \exp{(\sum^N_{n=1}-\frac{1}{2}(\mathbf{x}_n-\mu)^T C^{-1}(\mathbf{x}_n-\mu))}
$$

となります。最尤法を適用する際に計算を簡単にするために対数尤度を取ります。

$$
\ln p(X|\mu,C)=\sum^N_{n=1}\ln p(\mathbf{x}_n|\mu,C)=-\frac{ND}{2}\ln-\frac{N}{2}\ln (|C|)-\sum^N_{n=1}\frac{1}{2}(\mathbf{x}_n-\mu)^TC^{-1}(\mathbf{x}_n-\mu)
$$

これを$\mu$で偏微分することで平均を導けます。
第1、第2項に関しては偏微分すると0になることは自明なので第3項について考えます。以下に示す公式を用いて第3項を偏微分します。

>$$
\frac{\partial}{\partial x}(x^TAx)=(A+A^T)x
$$

$$
\frac{\partial}{\partial \mu}\ln p(X|\mu,C)=\frac{1}{2}(C^{-1}+C^{-T})\sum^N_{n=1}(\mathbf{x}_n-\mu)=0
$$

$$
\sum^N_{n=1}\mathbf{x}_n=N\mu
$$

$$
\frac{1}{N}\sum^N_{n=1}\mathbf{x}_n=\mu
$$

$\mu$は結局与えられた全データの平均となりました
あとはこれを各パラメータ$W,\sigma ^2$で同様に行うのですが、偉い人たち(Tipping&Bishop)が既に答えを出してくれています(実際に最尤法で解こうとするとえらい複雑らしい)。

$$
W_{ML}=U_M(L_M-\sigma ^2I)^{\frac{1}{2}}R
$$

ここで$U_M$は$D×M$行列で、その各列ベクトルは、$X$から得られるデータ共分散行列$S$の固有ベクトルの任意の部分集合からなります。$L_M$は固有ベクトルに対応する固有値$\lambda_i$で構成される$M×M$の対角行列です。$R$は任意の$M×M$の直交行列です。
また、$W_{ML}$のパラメータ$\sigma$の最尤解は以下のように求められます。

$$
\sigma^2_{ML}=\frac{1}{D-M}\sum^D_{i=M+1}\lambda_i
$$

さらにTipping&Bishopらは固有値の上位$M$個を選ぶことで最も尤度が高い$W_{ML}$と$\sigma^2_{ML}$が得られることを示しています。
次節ではこの値を使って実際に主成分分析を実装してみます。
# 実装

import numpy as np#numpyをインストール
import matplotlib.pyplot as plt#グラフを書くやつ
import scipy.linalg
from mpl_toolkits.mplot3d import Axes3D#3次元plotを行うためのツール
from sklearn.preprocessing import Normalizer
from sklearn.preprocessing import StandardScaler
from sklearn import decomposition#PCAのライブラリ
from sklearn import datasets#datasetsを取ってくる
from numpy import linalg as LA
UM=[]
n_component=2
np.set_printoptions(suppress=True)
fig = plt.figure(2, figsize=(14, 6))
iris = datasets.load_iris()
X = iris.data
y = iris.target
N=150
D=4
M=2
step=100
W =np.random.randn(4, 2)
sigma = 0
for i in range(4):#各データの平均を元のデータから引いている
    mean = np.mean(X[:,i])
    X[:,i]=(X[:,i]-mean)
    
    
X_cov=np.dot(X.T,X)#共分散行列を生成
w,v=np.linalg.eig(X_cov)#共分散行列の固有値、固有ベクトルを算出(固有値は既に大きさ順に並べられている)
for i in range(n_component):#固有値の大きい順に固有ベクトルを取り出す
    UM.append(v[:,i])
UM=np.array(UM)
a = np.array([w[0], w[1]])
Lm = np.diag(a)
m = float(D-M)

for i in range(D-M+1):
  sigma += w[i]/m
  
Xpc = scipy.linalg.sqrtm(Lm-sigma*np.eye(M))
Wml=np.dot(UM.T,Xpc)#WMLの計算
Xafter=np.dot(X,W)#取り出した固有ベクトルでデータを線形写像する


pca = decomposition.PCA(n_components=2)
pca.fit(X)
Xlib = pca.transform(X)    


ax1 = fig.add_subplot(121)
for label in np.unique(y):
    ax1.scatter(Xlib[y == label, 0],
                Xlib[y == label, 1])
ax1.set_title('library pca')
ax1.set_xlabel("X_axis")
ax1.set_ylabel("Y_axis")
ax2 = fig.add_subplot(122)
for label in np.unique(y):
    ax2.scatter(Xafter[y == label, 0],
                Xafter[y == label, 1])
ax2.set_title('original pca')
ax2.set_xlabel("X_axis")
ax2.set_ylabel("Y_axis")
plt.show()


左の図は主成分分析の結果を、右の図は確率的主成分分析の結果を示しています。$W$の初期値を乱数としているため、実行するごとに結果は変化しますが、おおむね主成分を残して次元削減出来ていると思います。

# まとめ
今回の記事では最尤推定での実装までしか行えませんでした。最尤推定では主成分分析を確率的に解釈しなおすことのメリットはあまり感じることはできなかったです。今後はEMアルゴリズムやベイズ推定を用いた確率的主成分分析を実装することで、確率的主成分分析の真髄を体感してみようと思います。

## 参考文献
パターン認識と機械学習(下巻)