PCA と固有値問題

主成分分析

主成分分析(PCA)は、多次元のデータをより少次元に置き換える、次元縮約の方法の一つである。多次元のデータに対して、分散が最も大きくなるような次元を探し当てれば、その次元だけで、データが持つ情報の多くを説明できるようになる。PCA はそのような分散が最も大きくなるような次元(主成分)を探す方法である。

PCA と固有値問題

n 個のサンプルのそれぞれに対して、p 項目の調査データ(特徴)があると仮定したとき、データ全体は n × p の行列で表せる。ここで、i (1 ≤ i ≤ p) 番目の特徴を \(\mathbf{x}_{i} = (x_{1i}, x_{2i}, \cdots, x_{ni})^{T} \) で表し、全サンプルのデータからなる n × p の行列をベクトル表示で \(\mathbf{X} = (\mathbf{x}_{1}, \mathbf{x}_{2}, \cdots, \mathbf{x}_{n})\) と表す。

分散共分散行列

次に、p 個の特徴(調査項目)のうち、k 番目と l 番目の特徴の共分散は、定義式に従って、次のように計算することができる。

\[ \begin{eqnarray} Cov(\mathbf{x}_{k}, \mathbf{x}_{l}) &=& E\left[\left(\mathbf{x}_{k}-E[\mathbf{x}_{k}]\right)\left(\mathbf{x}_{l}-E[\mathbf{x}_{l}]\right)\right] \\ &=& \frac{1}{n}\sum_{j=1}^{n}\left( x_{jk} - \mu_{k} \right) \left( x_{jl} - \mu_{l} \right) \\ &=& \frac{1}{n}\mathbf{x'}_{k}^{T}\mathbf{x'}_{l} \end{eqnarray} \]

ただし、\( \mu_{k} = \sum_{j=1}^{n}x_{jk} \) とし、また、\( \mathbf{x'}_{k} = (x_{1k}-\mu_{k}, x_{2k}-\mu_{k}, \cdots, x_{nk}-\mu_{k})^{T} \) とした。

p 個の特徴に対して、特徴のペア同士の共分散を、すべての組み合わせに対して計算すると、次のような分散共分散行列 Σ が得られる。

\[ \mathbf{\Sigma}_{\mathbf{X}} = \begin{bmatrix} Cov(\mathbf{x}_{1}, \mathbf{x}_{1}) & Cov(\mathbf{x}_{1}, \mathbf{x}_{2}) & \cdots & Cov(\mathbf{x}_{1}, \mathbf{x}_{n}) \\ Cov(\mathbf{x}_{2}, \mathbf{x}_{1}) & Cov(\mathbf{x}_{2}, \mathbf{x}_{2}) & \cdots & Cov(\mathbf{x}_{2}, \mathbf{x}_{n}) \\ \vdots & \vdots & \ddots & \vdots \\ Cov(\mathbf{x}_{n}, \mathbf{x}_{1}) & Cov(\mathbf{x}_{n}, \mathbf{x}_{2}) & \cdots & Cov(\mathbf{x}_{n}, \mathbf{x}_{n}) \end{bmatrix} = \frac{1}{n}\mathbf{X'}^{T}\mathbf{X'} \]

ただし、\(\mathbf{X'} = (\mathbf{x'}_{1}, \mathbf{x'}_{2}, \cdots, \mathbf{x'}_{n} )\) とした。

固有値問題への帰着

PCA の目標は、特徴の分散が最大となるような次元を見つけることである。そのため、既存のデータ \(\mathbf{X}\) に対して座標変換を行い、変換後のデータの分散が最大となるような変換を見つければよい。ここで、\(\mathbf{X}\) に対して、座標変換を行うために、ある単位ベクトル \( \mathbf{w} \) をかけて、新しい座標を \(\mathbf{Z}\) を計算する。

\[ \mathbf{Z} = \mathbf{X}\mathbf{w} \]

こうして得られた変換後の新しい座標 \(\mathbf{Z} = (\mathbf{z}_{1}, \mathbf{z}_{2}, \cdots, \mathbf{z}_{n})\) の分散を最大にする単位ベクトル \(\mathbf{w}\) を求めればよい。そのために、\(\mathbf{Z}\) の分散を計算する必要がある。上と同様にして、\(\mathbf{Z}\) の分散共分散行列は次のように計算できる。

\[ \mathbf{\Sigma}_{\mathbf{Z}} = \frac{1}{n}\mathbf{Z'}^{T}\mathbf{Z'} \]

\( \mathbf{Z} = \mathbf{X}\mathbf{w} \) だから、\(\mathbf{Z}\) の分散共分散行列は、\(\mathbf{X}\) の分散共分散行列を使って次のように置き換えることができる。

\[ \begin{eqnarray} \mathbf{\Sigma}_{\mathbf{Z}} = \frac{1}{n}\mathbf{Z'}^{T}\mathbf{Z'} &=& \frac{1}{n}(\mathbf{X'}\mathbf{w})^{T}(\mathbf{X'}\mathbf{w}) \\ &=& \frac{1}{n}\mathbf{w}^{T}\mathbf{X'}^{T}\mathbf{X'}\mathbf{w} \\ &=& \mathbf{w}^{T}\mathbf{\Sigma}_{\mathbf{X}}\mathbf{w} \end{eqnarray} \]

変換後のデータの分散を最大にしたいので、\(\mathbf{\Sigma}_{\mathbf{Z}}\) となるようなベクトル \(\mathbf{w}\) を求めればよい。ベクトル \(\mathbf{w}\) は単位ベクトルなので、 \(|\mathbf{w}| = 1\) という制約を利用すれば、ラグランジェの未定乗数法で \(\mathbf{w}\) を解くことができる。

\[ L(\mathbf{w}) = \mathbf{w}^{T}\mathbf{\Sigma}\mathbf{w} - \lambda (\mathbf{w}^{T}\mathbf{w} - 1) \]

この式を最大化する \(\mathbf{w}\) は、関数 L を \(\mathbf{w}\) で偏微分した一次導関数を 0 と置くことで計算される。

\[ \frac{\partial L(\mathbf{w})}{\partial \mathbf{w}} = 2 \mathbf{\Sigma}\mathbf{w} - 2 \lambda \mathbf{w} = 0 \] \[ \Sigma_{\mathbf{X'}}\mathbf{w} = \lambda\mathbf{w} \]

このように PCA は固有値問題へと導かれる。このとき、 \(\mathbf{w}\) を固有ベクトル、 \(\lambda\) を固有値という。

固有ベクトル \(\mathbf{w}\) がもとまると、変換後のデータの座標 \(\mathbf{Z}\)も得られ、これが主成分スコアにあたる。

R による PCA

R には 2 つの PCA を行う関数 prcompprincomp がある。2 つの関数の入力パラメーターと出力フォーマットが異なるが、実際の解析ではともに PCA を行なう。

各特徴の単位が揃っているときは、データをそのまま PCA 行うことができる。しかし、各特徴の単位が揃っていないとき(cm だったり、g だったり)、各特徴を正規化(標準化)しておく必要がある。標準化を行う場合、prcomp 関数では scale = TRUE を指定し、princomp 関数では cor = TRUE を指定する。各関数の実行例を次に示す。(このデータはフィッシャーのアヤメのデータであり、測定単位はすべて cm となっているため、本来は正規化する必要はない。)

data(iris)

iris.data <- iris[, 1:4]
iris.name <- iris[, 5]

p1 <- prcomp(iris.data, scale = TRUE)
summary(p1)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000


p2 <- prcomp(iris.data, scale = FALSE)
summary(p2)
## Importance of components:
##                           PC1     PC2    PC3     PC4
## Standard deviation     2.0563 0.49262 0.2797 0.15439
## Proportion of Variance 0.9246 0.05307 0.0171 0.00521
## Cumulative Proportion  0.9246 0.97769 0.9948 1.00000


p3 <- princomp(iris.data, cor = TRUE)
summary(p3)
## Importance of components:
##                           Comp.1    Comp.2     Comp.3      Comp.4
## Standard deviation     1.7083611 0.9560494 0.38308860 0.143926497
## Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
## Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000


p4 <- princomp(iris.data, cor = FALSE)
summary(p4)
## Importance of components:
##                           Comp.1     Comp.2     Comp.3      Comp.4
## Standard deviation     2.0494032 0.49097143 0.27872586 0.153870700
## Proportion of Variance 0.9246187 0.05306648 0.01710261 0.005212184
## Cumulative Proportion  0.9246187 0.97768521 0.99478782 1.000000000

prcomp

prcomp 関数の出力結果には標準偏差 sdev、固有ベクトル rotation、特徴量平均値 center、正規化パラメーター scale、主成分スコア x が保存されている。

data(iris)
iris.data <- iris[, 1:4]
iris.name <- iris[, 5]

p <- prcomp(iris.data, scale = TRUE)

names(p)
## [1] "sdev"     "rotation" "center"   "scale"    "x"

head(p$x)
##            PC1        PC2         PC3          PC4
## [1,] -2.684126 -0.3193972  0.02791483  0.002262437
## [2,] -2.714142  0.1770012  0.21046427  0.099026550
## [3,] -2.888991  0.1449494 -0.01790026  0.019968390
## [4,] -2.745343  0.3182990 -0.03155937 -0.075575817
## [5,] -2.728717 -0.3267545 -0.09007924 -0.061258593
## [6,] -2.280860 -0.7413304 -0.16867766 -0.024200858

princomp

princomp 関数の出力結果には標準偏差 sdev、固有ベクトル loadings、特徴量平均値 center、正規化パラメーター scale、サンプル数 n.obs、主成分スコア scores、関数実行コマンドログ call が保存されている。

data(iris)
iris.data <- iris[, 1:4]
iris.name <- iris[, 5]

p <- prcomp(iris.data, scale = TRUE)

names(p)
## [1] "sdev"     "loadings" "center"   "scale"    "n.obs"    "scores"   "call"

head(p$scores)
##         Comp.1     Comp.2      Comp.3       Comp.4
## [1,] -2.684126  0.3193972  0.02791483  0.002262437
## [2,] -2.714142 -0.1770012  0.21046427  0.099026550
## [3,] -2.888991 -0.1449494 -0.01790026  0.019968390
## [4,] -2.745343 -0.3182990 -0.03155937 -0.075575817
## [5,] -2.728717  0.3267545 -0.09007924 -0.061258593
## [6,] -2.280860  0.7413304 -0.16867766 -0.024200858