- 講師：津田宏治
- 資料：生物データマイニング論(2019)
主成分分析(PCA)
- 次元削減・非可逆データ圧縮・特徴抽出・データの可視化などに使われる。
- 主部分空間(principal subspace)と呼ばれる低次元の線形空間の上への、データ点の直交射影として、以下の2つの定義方法がある。
- 「射影されたデータの分散が最大化されるような直交射影」として定義。→分散最大化による定式化
- 「データ点と射影した点の間の2乗距離の平均値で定義される射影コスト関数の期待値を最小化するような線形射影」として定義。（要はデータ点からの距離を最小化）→誤差最小化による定式化
- Pythonによる実装
分散最大化による定式化
次元 $D$ のユークリッド空間内のデータ集合 $\{\mathbf{x}_n\}(n=1,\ldots,N)$ を考える。ここで、「射影されたデータ点の分散を最大化しながら、データを次元 $M
1次元空間の上への射影（$M=1$）¶
この空間の方向を $D$ 次元ベクトル $\mathbf{u}_1$ として表す。ここで、興味があるのは $\mathbf{u}_1$ の方向だけであるので、$\mathbf{u}_1^T\mathbf{u}_1 = 1$ としても一般性を失わない。
したがって、各データ点 $\mathbf{x}_n$ はスカラー値 $\mathbf{u}_1^T\mathbf{x}_n$ の上に射影される。ここで、射影されたデータの平均・分散は以下のようになる。
| 統計量 | |
|---|---|
| 平均値 | $$\mathbf{u}_1^T\bar{\mathbf{x}}\qquad\bar{\mathbf{x}} = \frac{1}{N}\sum_{n=1}^N\mathbf{x}_n\qquad (12.1)$$ |
| 分散 | $$\mathbf{u}_1^T\mathbf{Su}_1\qquad\mathbf{S} = \frac{1}{N}\sum_{n=1}^N\left(\mathbf{x}_n-\bar{\mathbf{x}}\right)\left(\mathbf{x}_n-\bar{\mathbf{x}}\right)^T\qquad (12.3)$$ |
ここで、射影されたデータの分散を、$\mathbf{u}_1$ に関して最大化することを考える。
制約条件 $\mathbf{u}_1^T\mathbf{u}_1 = 1$ をラグランジュ乗数 $\lambda_1$ を利用して導入して
$$\mathbf{u}_1^T\mathbf{Su}_1 + \lambda_1\left(1-\mathbf{u}^T_1\mathbf{u}_1\right)\qquad (12.4)$$を最大化すれば良い。$\mathbf{u}_1$ に関する微分を $0$ とおくことにより、停留条件は
$$\mathbf{Su}_1 = \lambda_1\mathbf{u}_1\qquad (12.5)$$となることがわかる。これは、$\mathbf{u}_1$ が $\mathbf{S}$ の固有ベクトルであることを示している。左から $\mathbf{u}^T_1$ をかけると、
$$\text{Variance} = \mathbf{u}^T_1\mathbf{Su}_1 = \lambda_1 \qquad (12.6)$$となる。よって、分散は $\mathbf{u}_1$ を、最大固有値 $\lambda_1$ に属する固有ベクトルに選んだ時に最大になることがわかる。
$M$ 次元空間の上への射影（一般解）¶
演習12.1
帰納法を用いて、射影されたデータの分散を最大化するようなM次元部分空間の上への線形写像が、データ共分散行列Sの上位M個の固有値に属するM本の固有ベクトルによって定義されることを示す。仮定¶
- $M$ に対して上記が成り立つ。
- データ空間に対する新しい方向ベクトル $\mathbf{u}_{M+1}$ は、次を満たす。
- $\mathbf{u}_{M+1}$ は、$\mathbf{u}_1,\ldots,\mathbf{u}_M$ と直交する。（$\mathbf{u}_{M+1}^T\mathbf{u}_{j}=0$）
- 単位長に規格化されている。（$\mathbf{u}_{M+1}^T\mathbf{u}_{M+1}=1$）
制約条件を踏まえると、
$$L = \sum_{i=1}^{M+1}\mathbf{u}^T_i\mathbf{Su}_i + \lambda_{M+1}\left(1-\mathbf{u}^T_{M+1}\mathbf{u}_{M+1}\right) + \sum_{j=1}^M\lambda^{\prime}_j\mathbf{u}_{M+1}^T\mathbf{u}_j$$両辺を $\mathbf{u}_{M+1}$ で微分すると、
$$0 = \frac{\partial L}{\partial \mathbf{u}_{M+1}} = 2\left(\mathbf{Su}_{M+1} - \lambda_{M+1}\mathbf{u}_{M+1}\right) + \sum_{j=1}^M\lambda^{\prime}_j\mathbf{u}_j\cdots(\ast)$$これに左から $\mathbf{u}^T_{j}$ をかけると、
$$0 = \lambda^{\prime}_j\quad(\because \mathbf{u}_{M+1}^T\mathbf{u}_{j}=0)$$となるので、これを $(\ast)$ に代入して、
$$\mathbf{Su}_{M+1} = \lambda_{M+1}\mathbf{u}_{M+1}$$また、この時分散は
$$\begin{aligned} \sum_{i=1}^{M+1}\mathbf{u}^T_i\mathbf{Su}_i &= \sum_{i=1}^{M}\mathbf{u}^T_i\mathbf{Su}_i + \mathbf{u}^T_{M+1}\mathbf{Su}_{M+1}\\ &= \sum_{i=1}^{M}\lambda_i + \lambda_{M+1} \end{aligned}$$ここで、仮定より $\lambda_i$ は、データ共分散行列 $\mathbf{S}$ の上位 $M$ 個の固有値なので、この分散を最大する $\lambda_{M+1}$ は、$M+1$ 番目に大きい固有値となる。
誤差最小化による定式化
続いて、射影誤差の最小化に基づいた主成分分析の定式化を考える。まず、$D$ 次元の基底ベクトル $\{\mathbf{u}_i\}$ から成る完全正規直交系を導入する。（$\mathbf{u}_i^T\mathbf{u}_i = \delta_{ij}$）
すると、各データ点はこれらの基底ベクトルの線形結合なので、以下の形で表すことができる。
$$\mathbf{x}_n = \sum_{i=1}^D\alpha_{ni}\mathbf{u}_i\qquad (12.8)$$なお、ここで式 $(12.8)$ は、各データ点 $\{\mathbf{x}_n\}$ を、もともとの座標系から $\{\mathbf{u}_i\}$ で表される新しい座標系への回転を表していると解釈できる（もともと $D$ 個の成分 $\{x_{n1},\ldots,x_{nD}\}$ が、等価な集合 $\{a_{n1},\ldots,a_{nD}\}$ に置き換えられている）。
上式と $\mathbf{u}_j$ との内積を取り、正規直交性を使うと、$\alpha_{nj} = \mathbf{x}_n^T\mathbf{u}_j$ を得るので、
$$\mathbf{x}_n = \sum_{i=1}^D\left(\mathbf{x}_n^T\mathbf{u}_i\right)\mathbf{u}_i\qquad (12.9)$$$D$ 次元から $M$ 次元への近似¶
ここで、$M<D$ 次元の線形部分空間のデータ点でこれらを近似することが目的となる。（$M$ 次元空間で完全に $\mathbf{x}_n$ を表現するのは一般に不可能なので、近似する。）
したがって、各データ点を
$$\tilde{\mathbf{x}}_n = \sum_{i=1}^Mz_{ni}\mathbf{u}_i+\sum_{i=M+1}^Db_i\mathbf{u}_i\qquad (12.10)$$- $\{z_{ni}\}$ は特定のデータ点に依存。
- $\{b_i\}$ は全てのデータ点に共通な定数。
と近似して表す。つまり、$\{\mathbf{u}_i\}$ と $\{z_{ni}\}$ と $\{b_i\}$ を調整することで、次元が減ることによってもたらされる歪みを最小化することがここでの目的となる。
近似による歪み尺度 $J$¶
近似による歪みの尺度として、元々のデータ点 $\mathbf{x}_n$ とその近似 $\tilde{\mathbf{x}}_n$ の間の2乗距離をデータ集合に渡って平均したものを採用する。
$$\begin{aligned} J &= \frac{1}{N}\sum_{n=1}^N\|\mathbf{x}_n-\tilde{\mathbf{x}}_n\|^2\qquad (12.11)\\ &= \frac{1}{N}\sum_{n=1}^N\left\|\mathbf{x}_n-\left(\sum_{i=1}^Mz_{ni}\mathbf{u}_i+\sum_{i=M+1}^Db_i\mathbf{u}_i\right)\right\|^2 \end{aligned}$$そこで、各変数についての微分を $0$ としてこれを最小化する変数を求める。
$\{z_{ni}\}$ に対する最小化¶
$$ \begin{aligned} \frac{\partial J}{\partial z_{ni}} &= \frac{1}{N}\sum_{n=1}^N\frac{\partial}{\partial z_{ni}}\left\{\left\|\sum_{i=1}^Mz_{ni}\mathbf{u}_i\right\|^2 -2\left(\mathbf{x}_n-\sum_{i=M+1}^Db_i\mathbf{u}_i\right)^T\left(\sum_{i=1}^Mz_{ni}\mathbf{u}_i\right) + C \right\}\\ &= \frac{1}{N}\sum_{n=1}^N\frac{\partial}{\partial z_{ni}}\left\{\sum_{i=1}^Mz_{ni}^2-2\sum_{i=1}^Mz_{ni}\mathbf{x}_n^T\mathbf{u}_i\right\}\\ &= \frac{1}{N}\sum_{n=1}^N\sum_{i=1}^M2\left\{z_{ni}-\mathbf{x}_n^T\mathbf{u}_i\right\} = 0\\ \therefore z_{ni}^{\star} &= \mathbf{x}_n^T\mathbf{u}_i \qquad (12.12) \end{aligned} $$$b_i$ に対する最小化¶
$$ \begin{aligned} \frac{\partial J}{\partial b_{i}} &= \frac{1}{N}\sum_{n=1}^N\frac{\partial}{\partial b_{i}}\left\{\left\|\sum_{i=M+1}^Db_{i}\mathbf{u}_i\right\|^2 -2\left(\mathbf{x}_n-\sum_{i=1}^Mz_{ni}\mathbf{u}_i\right)^T\left(\sum_{i=M+1}^Db_{i}\mathbf{u}_i\right) + C \right\}\\ &= \frac{1}{N}\sum_{n=1}^N\frac{\partial}{\partial b_{i}}\left\{\sum_{i=M+1}^Db_{i}^2-2\sum_{i=M+1}^Db_{i}\mathbf{x}_n^T\mathbf{u}_i\right\}\\ &= \frac{1}{N}\sum_{n=1}^N\sum_{i=M+1}^D2\left\{b_{i}-\mathbf{x}_n^T\mathbf{u}_i\right\} = 0\\ \therefore b_{i}^{\star} &= \bar{\mathbf{x}}_n^T\mathbf{u}_i \qquad (12.13) \end{aligned} $$$J$ の整頓¶
ここで、$(12.10)$ の $z_{ni},b_i$ を置き換えると、
$$ \begin{aligned} \mathbf{x}_n-\tilde{\mathbf{x}}_n &= \mathbf{x}_n-\left(\sum_{i=1}^M\mathbf{x}_n^T\mathbf{u}_i\mathbf{u}_i + \sum_{i=M+1}^D\bar{\mathbf{x}}^T\mathbf{u}_i\mathbf{u}_i\right)\\ &= \sum_{i=1}^D\mathbf{x}_n^T\mathbf{u}_i\mathbf{u}_i - \left(\sum_{i=1}^M\mathbf{x}_n^T\mathbf{u}_i\mathbf{u}_i + \sum_{i=M+1}^D\bar{\mathbf{x}}^T\mathbf{u}_i\mathbf{u}_i\right)\\ &= \sum_{i=M+1}^D\left\{\left(\mathbf{x}_n-\bar{\mathbf{x}}\right)^T\mathbf{u}_i\right\}\mathbf{u}_i\qquad (12.14) \end{aligned} $$と表すことができる。この式から、$\tilde{\mathbf{x}}_n$ から $\mathbf{x}_n$ への変位を表すベクトルは、主部分空間に直交する空間にあることがわかる。
以上より、歪み尺度 $J$ に対する表現を、純粋に $\{\mathbf{u}_i\}$ の関数として、以下のような形で得られることがわかる。
$$J = \frac{1}{N}\sum_{n=1}^N\sum_{i=M+1}^D\left(\mathbf{x}_n^T\mathbf{u}_i-\bar{\mathbf{x}}^T\mathbf{u}_i\right)^2 = \sum_{i=M+1}^D\mathbf{u}^T_i\mathbf{Su}_i\qquad (12.15)$$したがって、あとは$J$ を $\{\mathbf{u}_i\}$ に対して最小化すれば良い。
$M=1,D=2$ の場合¶
$M=1,D=2$ の場合、「規格化条件 $\mathbf{u}^T_2\mathbf{u}_2 = 1$ の下で、方向 $\mathbf{u}_2$ を、$J=\mathbf{u}^T_2\mathbf{Su}_2$ が最小化されるように選ぶ」問題となる。
つまり、ラグランジュ乗数を用いて、
$$\tilde{J} =\mathbf{u}^T_2\mathbf{Su}_2 + \lambda_2\left(1-\mathbf{u}^T_2\mathbf{u}_2\right)\qquad (12.16)$$を最小化する問題となる。
$\mathbf{u}_i$ についての微分を $0$ とおくと、$\mathbf{Su}_2 = \lambda_2\mathbf{u}_2$ を得る。つまり、$\mathbf{u}_2$ は固有値 $\lambda_2$ に属する $\mathbf{S}$ の固有ベクトルであるので、任意の固有ベクトルは、この歪み尺度の停留点を定義することとなる。
$(12.16)$ に $\mathbf{Su}_2 = \lambda_2\mathbf{u}_2$ を代入すると、$J=\lambda_2$ となる。ゆえに、「$\mathbf{u}_2$ を、2つある固有値のうち小さい方に対応した固有ベクトルを選ぶことにより、$J$ の最小値を得ることができる」（↔︎主部分空間では大きい方）
任意の $D$ と任意の $M<D$ の場合¶
任意の $D$ と任意の $M<D$ に対する $J$ の最小化問題の一般解は、$\{\mathbf{u}_i\}$ を
$$\mathbf{Su}_i = \lambda_i\mathbf{u}_i\qquad (12.17)$$で与えられる共分散行列の固有ベクトルに選ぶことによって得られる。ここで $i=1,\ldots,D$ であり、固有ベクトル $\{\mathbf{u}_i\}$ は正規直交するように選ばれる。
また、対応する歪み尺度の値は、
$$J = \sum_{i=M+1}^D\lambda_i\qquad (12.18)$$のように、主部分空間に直交する固有ベクトルの固有値についての単なる和で与えられる。したがって、$J$ の最小値は、小さい順から $D-M$ 個の固有値に対応する固有ベクトルを選ぶことによって得られ、つまり主部分空間を定義する固有ベクトルは大きい順から $M$ 個の固有値に対応する固有ベクトルになる。
実装
Wine Data Set (UCI Machine Learning Repository) のデータを利用し、13次元、178データ、3種類のワインの属性データから、2つの主成分を取って可視化を行います。
また、PCAは自作のモジュール（kerasy.ML.decomposition.PCA）を用います。
# データ取得。
! wget https://raw.githubusercontent.com/maskot1977/ipython_notebook/master/toydata/wine.txt
from kerasy.ML.decomposition import PCA
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
colors = [(1, 0, 0), (0, 0, 1), (0, 1, 0)]
n_bin = len(colors)
cmap_name = 'my_list'
cm = LinearSegmentedColormap.from_list(cmap_name, colors, N=n_bin)
df = pd.read_csv("wine.txt", sep="\t", index_col=0)
df.head(3)
class カラムがワインの種類を表しているので、この列は主成分分析から取り除きます。
cls = df.iloc[:, 0].values
df_use = df.iloc[:, 1:]
また、データの標準化を行います。
df_std = df_use.apply(lambda x: (x-x.mean())/x.std(), axis=0)
data = np.array(df_std)
データの準備ができたので、主成分分析を行います。
pca = PCA()
pca.fit(data)
features = pca.transform(data)
plt.figure(figsize=(6, 6))
plt.scatter(features[:, 0], features[:, 1], c=cls, cmap=cm, edgecolors='black')
plt.title("Same fig with slide :)")
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.grid()
plt.show()
