データ解析の手法の一つである主成分分析(principal component analysis ; PCA)について,それなりに利用頻度が高いものの,そのたびに勉強しなおしていて効率が悪かったので,その基礎をまとめておくことにしました.

=================================================================================
[ 1. 準備 ]
データ行列を準備します.{displaystyle P } 個の変数についてそれぞれ {displaystyle N } 個のサンプルがある場合を考え,測定値を {displaystyle { x_{np}^* }_{n=1,ldots,N ;  p=1,ldots,P} } とします.議論を簡単にするために,各変数についてその平均

(1.1){displaystyle ;;;;;; ar{x}_p = frac{1}{N} sum_{n=1}^N x_{np} , ;;;;;; p=1,ldots,P }

からの偏差

(1.2){displaystyle ;;;;;; x_{np} = x_{np}^* - ar{x}_p, ;;;;;;  n=1,ldots,N  ;  p=1,ldots,P }

を導入します.すると測定データ全体は以下の {displaystyle N 	imes P } 行列であるデータ行列 {displaystyle X } で表現できます.各サンプルを意味する {displaystyle N } 個の測定データベクトル {displaystyle { oldsymbol{x}_{n} }_{n=1,ldots,N} } を定義しておきます.

(1.3){displaystyle ;;;;;; X=  egin{bmatrix} x_{11} & x_{12} & cdots & x_{1P}  x_{21} & x_{22} & cdots &  x_{2P}  vdots & vdots & vdots & vdots   vdots & vdots  & vdots & vdots   x_{N1} & x_{N2}  & cdots & x_{NP} end{bmatrix} = egin{bmatrix} oldsymbol{x}_{1}^T   oldsymbol{x}_{2}^T  vdots  vdots  oldsymbol{x}_{N}^T end{bmatrix} , ;;;;;; oldsymbol{x}_{n} = egin{bmatrix} x_{n1}  x_{n2}   vdots  x_{nP} end{bmatrix}, ;;; n=1,ldots,N }


分散共分散行列を準備します.分散共分散行列は以下の {displaystyle P 	imes P } 行列 {displaystyle V } で表現できます,その {displaystyle (i,j) } 要素 {displaystyle v_{ij} } も示します.

(1.4){displaystyle ;;;;;; V = frac{1}{N-1} X^T  X  }

(1.5){displaystyle ;;;;;; v_{ij} = frac{1}{N-1} sum_{n=1}^N x_{ni} x_{nj} = frac{1}{N-1} sum_{n=1}^N ( x_{ni}^* - ar{x}_i) (x_{nj}^* - ar{x}_j) = v_{ji} }

(1.5)より分散共分散行列 {displaystyle V= V^T } なので実対称行列です.


分散共分散行列の固有値,固有ベクトルとその性質を準備します.{displaystyle V }{displaystyle P } 個の固有値を {displaystyle { lambda_p }_{p=1,ldots,P} },長さ1の固有ベクトルを {displaystyle { oldsymbol{v}_p^* }_{p=1,2,ldots,P} } とします.つまり以下が成り立つとします. {displaystyle langle cdot , cdot 
angle  } は通常の内積, {displaystyle || cdot||_2 = ( langle cdot , cdot 
angle )^{1/2} } はユークリッドノルムです.

(1.6){displaystyle ;;;;;; V oldsymbol{v}_p^* = lambda_p oldsymbol{v}_p^*, ;;; ||oldsymbol{v}_p^*||_2=1, ;;;;;; p=1,2,ldots,P }

測定データから計算されていることから,現実のデータ解析において固有方程式が重解を持つことはあまりないと思われます.よって,以下を仮定します.

仮定.
{displaystyle V } のすべての固有値は相異なる.すなわち {displaystyle lambda_i 
eq lambda_j,  i 
eq j, ;; i,j =1,2,ldots,P } である.

定義より {displaystyle V } は半正定値であり,またエルミート行列でもあるのですべての固有値は実数です.さらにすべての固有値は相異なると仮定したので,{displaystyle V } の固有ベクトルはそれぞれ互いに直交します.それらの証明は以下の過去記事にあります.

分散共分散行列(と相関行列)は半正定値であることを証明する - エンジニアを目指す浪人のブログ

エルミート行列のすべての固有値は実数であることの証明をメモする - エンジニアを目指す浪人のブログ

エルミート行列の相異なる固有値に対する固有ベクトルは直交することの証明をメモする - エンジニアを目指す浪人のブログ

以下の2つの事実を証明なしで用います.
・半正定値行列のすべての固有値は非負である
{displaystyle n 	imes n } エルミート行列 {displaystyle A }{displaystyle n } 個の固有ベクトルからなるユニタリ空間 {displaystyle mathbb{C}^n } の正規直交基底をとることができる.

よって {displaystyle { oldsymbol{v}_1^*,oldsymbol{v}_2^*, ldots ,oldsymbol{v}_P^* } }{displaystyle mathbb{R}^P } の正規直交基底です.また固有値は以下のように順序付けできます.

(1.7){displaystyle ;;;;;; 0 le lambda_mathrm{min} = lambda_1 lt lambda_2 lt cdots lt lambda_{P-1} lt lambda_P =lambda_mathrm{max}  }


結合係数を準備します.ある正の整数 {displaystyle M  (le P) } について結合係数 {displaystyle { w_{pm} }_{ p=1,ldots,P  ;  m=1,ldots,M}  } を導入します.これを表現する {displaystyle M } 個の結合係数ベクトルを定義しておきます.

(1.8){displaystyle ;;;;;; oldsymbol{w}_m = egin{bmatrix} w_{1m}  w_{2m}  vdots  w_{Pm} end{bmatrix}, ;;;  m=1,ldots,M  }



[ 2. 主成分分析の概要 ]
主成分分析とは,{displaystyle X } の持つ情報を,情報の損失を最小限に抑えながら,{displaystyle { x_{cdot p} }_{p=1,ldots,P} } の一次結合として与えられる互いに独立な {displaystyle M  (le P)  } 個の主成分 {displaystyle { z_m }_{m=1,ldots,M} },

(2.1){displaystyle ;;;;;; z_m = sum_{p=1}^P w_{pm} x_{cdot p}, ;;;;;; m=1,ldots,M }

を用いて表現する手法です.行列で表現すると以下となります.

(2.2){displaystyle egin{bmatrix} z_1  z_2  vdots  z_M end{bmatrix} = egin{bmatrix} w_{11} & w_{21} & cdots & cdots & w_{P1}  w_{12} & w_{22} & cdots & cdots & w_{P2}  cdots & cdots & cdots & cdots & cdots  w_{1M} & w_{2M} & cdots & cdots & w_{PM} end{bmatrix} egin{bmatrix} x_{cdot 1}  x_{cdot 2}  vdots  vdots  x_{cdot P} end{bmatrix} = egin{bmatrix}  oldsymbol{w}_{1}^T    oldsymbol{w}_{2}^T    vdots   oldsymbol{w}_{M}^T   end{bmatrix} oldsymbol{x}_{cdot} = egin{bmatrix} langle oldsymbol{w}_{1} , oldsymbol{x}_{cdot} 
angle   langle oldsymbol{w}_{2} ,  oldsymbol{x}_{cdot} 
angle  vdots    langle oldsymbol{w}_{M} , oldsymbol{x}_{cdot} 
angle end{bmatrix} = egin{bmatrix} langle oldsymbol{x}_{cdot},oldsymbol{w}_{1} 
angle   langle oldsymbol{x}_{cdot},oldsymbol{w}_{2} 
angle  vdots    langle oldsymbol{x}_{cdot},oldsymbol{w}_{M} 
angle end{bmatrix}   }

{displaystyle { z_m } } は第 {displaystyle m } 主成分と呼ばれ,結合係数 {displaystyle { w_{pm} }_{ p=1,ldots,  p ;  m=1,ldots,M} } は以下の条件を満たすように決定します.

条件.
第1主成分 {displaystyle z_1 }{displaystyle { x_{cdot p} }_{p=1,ldots,P} } のあらゆる1次式のなかで分散が最大となるものであり,第 {displaystyle m } 主成分 {displaystyle { z_m }_{m=2,ldots,M} } の分散は {displaystyle { z_l }_{l=1,ldots,m-1} } のすべてと無相関な {displaystyle { x_{cdot p} }_{p=1,ldots,P} } の1次式のなかで分散が最大となるものである.ただし({displaystyle oldsymbol{w}_m  } を大きくすると {displaystyle z_m  } の分散はいくらでも大きくなってしまうので),

(2.3) {displaystyle ;;;;;; sum_{p=1}^P (w_{pm})^2 =  || oldsymbol{w}_m ||_2^2 =1,  ;;;;;; m=1,ldots,M }

とする.

条件.について補足します.分散を最大にすることは以下の意味があります(文献[3]にわかりやすい図があります).
・射影したデータのバラツキが大きいほど,もとのデータの情報を多く含んでいると考えられる.
・もとのデータの情報の損失ができるだけ小さくなるような軸を探したい.



[ 3. 第 {displaystyle m } 主成分得点の分散の導出 ]
(1.8)(2.1)(2.2)より,{displaystyle n } 番目のサンプルに対応する第 {displaystyle m } 主成分の値,いいかえると第 {displaystyle m } 主成分得点 {displaystyle t_{nm}  }

(3.1){displaystyle ;;;;;; t_{nm} = sum_{p=1}^P w_{pm} x_{np} = egin{bmatrix} w_{1m} & w_{2m} & cdots & w_{Pm} end{bmatrix} egin{bmatrix} x_{n1}  x_{n2}  vdots  x_{nP} end{bmatrix} = oldsymbol{w}_m^T oldsymbol{x}_n = langle oldsymbol{w}_m , oldsymbol{x}_n 
angle }

とし,それらを {displaystyle N } 個のサンプルについてまとめた第 {displaystyle m } 主成分得点ベクトルを(1.3)も用いて

(3.2){displaystyle ;;;;;; oldsymbol{t}_m = egin{bmatrix} t_{1m}  t_{2m}  vdots  vdots  t_{Nm} end{bmatrix} = egin{bmatrix} oldsymbol{w}_m^T oldsymbol{x}_1  oldsymbol{w}_m^T oldsymbol{x}_2  vdots  vdots  oldsymbol{w}_m^T oldsymbol{x}_N end{bmatrix} = egin{bmatrix} oldsymbol{x}_1^T oldsymbol{w}_m   oldsymbol{x}_2^T oldsymbol{w}_m  vdots  vdots  oldsymbol{x}_N^T oldsymbol{w}_m end{bmatrix} = egin{bmatrix} oldsymbol{x}_1^T   oldsymbol{x}_2^T  vdots  vdots  oldsymbol{x}_N^T end{bmatrix} oldsymbol{w}_m = X oldsymbol{w}_{m}  }

とします.第 {displaystyle m } 主成分得点 {displaystyle t_{nm} } の平均 {displaystyle ar{t}_m } は(1.1)(1.2)(3.1)を用いて

(3.3){displaystyle ;;;;;; ar{t}_m = frac{1}{N} sum_{n=1}^N t_{nm} = frac{1}{N} sum_{n=1}^N  left( sum_{p=1}^P w_{pm} x_{np} 
ight) =  sum_{p=1}^P  w_{pm} left( frac{1}{N} sum_{n=1}^N x_{np} 
ight) = 0 }

であるので,第 {displaystyle m } 主成分得点 {displaystyle t_{nm} } の分散 {displaystyle sigma_{z_m}^2 } は(1.4)(3.2)も用いて以下となります.

(3.4){displaystyle ;;;;;; sigma_{z_m}^2 = frac{1}{N-1}sum_{n=1}^N (t_{nm})^2 = frac{1}{N-1} oldsymbol{t}_m^T oldsymbol{t}_m }
{displaystyle = frac{1}{N-1} (X oldsymbol{w}_m)^T  X oldsymbol{w}_m =  oldsymbol{w}_m^T left( frac{1}{N-1} X^T  X 
ight) oldsymbol{w}_m =  oldsymbol{w}_m^T V oldsymbol{w}_m (ge 0), ;;;;;; m=1,2,ldots,M }

最後の {displaystyle ge } は分散共分散行列 {displaystyle V } が半正定値であることによります.



[ 4. 第1主成分得点の分散を最大化する結合係数の導出 ]
いよいよ結合係数 {displaystyle oldsymbol{w}_1 } を決めていきます.条件.より,制約(2.3)を満たし,(3.4)で {displaystyle  m = 1 } とした {displaystyle sigma_{z_1}^2 = oldsymbol{w}_1^T V oldsymbol{w}_1 } が最大となるような {displaystyle oldsymbol{w}_1^* } を決めればよいことになります.そのためにレイリー商の性質を用います.以下の過去記事の定理4.2.2.(c)をそのまま用いることができます.(1.6)(1.7)も用います.

レイリー商についての定理を証明する - エンジニアを目指す浪人のブログ

具体的に式で書くと以下となります.

(4.1){displaystyle ;;;;;; S = S_{P} = mathrm{span} { oldsymbol{v}_1^*,oldsymbol{v}_2^*, ldots ,oldsymbol{v}_{P}^* } = mathbb{R}^P }

(4.2){displaystyle  ;;;;;; max_{ { oldsymbol{w}_1  ;  oldsymbol{w}_1 in S_{P},  ||oldsymbol{w}_1||_2 = 1 } } sigma_{z_1}^2 = max_{ { oldsymbol{w}_1  ;  oldsymbol{w}_1 in S_{P},  ||oldsymbol{w}_1||_2 = 1 } } oldsymbol{w}_1^T V oldsymbol{w}_1 = lambda_P = lambda_mathrm{max} }

(4.3){displaystyle ;;;;;; oldsymbol{w}_1^* = { oldsymbol{w}_1 in mathbb{R}^P  ;   V oldsymbol{w}_1  =  lambda_mathrm{max} oldsymbol{w}_1 } = oldsymbol{v}_P^* }

{displaystyle oldsymbol{w}_1^* }{displaystyle  lambda_mathrm{max} = lambda_P } に対する長さ1の固有ベクトルであることがわかります.(3.2)で {displaystyle  oldsymbol{w}_m = oldsymbol{w}_1^* } としたとき {displaystyle oldsymbol{t}_1^* } と書くことにします.



[ 5. 第 {displaystyle m } 主成分得点の分散を最大化する結合係数の導出 ]
次に,第2主成分の結合係数 {displaystyle oldsymbol{w}_2 } を決めます.条件.を満たすように {displaystyle z_2 }{displaystyle z_1 } と無相関となるためには,(3.3)(4.3)を用いて

(5.1){displaystyle ;;;;;; 0 = frac{1}{N-1} oldsymbol{t}_2^T oldsymbol{t}_1^* = frac{1}{N-1} (X oldsymbol{w}_2)^T X oldsymbol{w}_1^* = oldsymbol{w}_2^T left( frac{1}{N-1} X^T X 
ight) oldsymbol{w}_1^* = oldsymbol{w}_2^T V oldsymbol{w}_1^* }
{displaystyle ;;;;;;;;;;;;;; = oldsymbol{w}_2^T lambda_mathrm{max} oldsymbol{w}_1^* = lambda_mathrm{max} oldsymbol{w}_2^T oldsymbol{w}_1^* = lambda_mathrm{max} langle oldsymbol{w}_2 , oldsymbol{w}_1^*  
angle }

でなければなりません.(1.7)より {displaystyle lambda_mathrm{max} 
eq 0 } であるので

(5.2){displaystyle ;;;;;; langle oldsymbol{w}_2 , oldsymbol{w}_1^*  
angle = 0 }

を得ます.したがって制約(2.3)(5.2)を満たし,(3.4)で {displaystyle  m = 2 } とした {displaystyle sigma_{z_2}^2 = oldsymbol{w}_2^T V oldsymbol{w}_2 } が最大となるような {displaystyle oldsymbol{w}_2^* } を決めればよいことになります.そのために第1主成分のときと同様にレイリー商の性質を用います.先ほどの過去記事の定理4.2.2.(a)(b)を用いることができます.具体的に式で書くと以下となります.

(5.3){displaystyle ;;;;;; S = S_{P-1} = mathrm{span} { oldsymbol{v}_1^*,oldsymbol{v}_2^*, ldots ,oldsymbol{v}_{P-1}^*,oldsymbol{0} } }

(5.4){displaystyle ;;;;;; max_{ { oldsymbol{w}_2  ;  oldsymbol{w}_2 in S_{P-1},  ||oldsymbol{w}_2||_2 = 1 } } sigma_{z_2}^2 = max_{ { oldsymbol{w}_2  ;  oldsymbol{w}_2 in S_{P-1},  ||oldsymbol{w}_2||_2 = 1 } } oldsymbol{w}_2^T V oldsymbol{w}_2 = lambda_{P-1} }

(5.5){displaystyle ;;;;;; oldsymbol{w}_2^* = { oldsymbol{w}_2 in mathbb{R}^P  ;   V oldsymbol{w}_2  =  lambda_{P-1} oldsymbol{w}_2 } = oldsymbol{v}_{P-1}^* }

{displaystyle oldsymbol{w}_2^* }{displaystyle lambda_{P-1} } に対する長さ1の固有ベクトルであることがわかります.


以下同様に,第 {displaystyle m } 主成分の結合係数 {displaystyle oldsymbol{w}_{m=3,cdots,M} } を決めることができます.条件.を満たすように {displaystyle z_m }{displaystyle { z_l }_{l=1,2,cdots,m-1}   } と無相関となるためには(5.1)と同様の計算により

(5.6) {displaystyle ;;;;;; oldsymbol{w}_i^T oldsymbol{w}_j =0, ;;; i 
eq j, ;;;;;; i,j=1,ldots,M }

を得ます.したがって制約(2.3)(5.6)を満たし,(3.4)の {displaystyle sigma_{z_m}^2 = oldsymbol{w}_m^T V oldsymbol{w}_m } が最大となるような {displaystyle oldsymbol{w}_m^* } を決めればよいことになります.そのために第2主成分のときと同様に先ほどの過去記事の定理4.2.2.(a)(b)を用いることができます.

(5.7){displaystyle ;;;;;; S = S_{P-m+1} = mathrm{span} { oldsymbol{v}_1^*,oldsymbol{v}_2^*, ldots ,oldsymbol{v}_{P-m+1}^*,oldsymbol{0},ldots,oldsymbol{0} } }

(5.8){displaystyle ;;;;;; max_{ { oldsymbol{w}_m  ;  oldsymbol{w}_m in S_{P-m+1},  ||oldsymbol{w}_m||_2 = 1 } } sigma_{z_m}^2 = max_{ { oldsymbol{w}_m  ;  oldsymbol{w}_m in S_{P-m+1},  ||oldsymbol{w}_m||_2 = 1 } } oldsymbol{w}_2^T V oldsymbol{w}_2 = lambda_{P-m+1} }

(5.9){displaystyle ;;;;;; oldsymbol{w}_m^* = { oldsymbol{w}_m in mathbb{R}^P  ;   V oldsymbol{w}_m  =  lambda_{P-m+1} oldsymbol{w}_m } = oldsymbol{v}_{P-m+1}^* }

{displaystyle oldsymbol{w}_m^* }{displaystyle lambda_{P-m+1} } に対する長さ1の固有ベクトルであることがわかります.


まとめると,条件.をみたす結合係数 {displaystyle { w_{pm} }_{ p=1,ldots,  p ;  m=1,ldots,M} } は以下となります.

(5.10){displaystyle ;;;;;; oldsymbol{w}_m = oldsymbol{w}_m^* = oldsymbol{v}_{P-m+1}^*, ;;; m=1,2,ldots,M }



[ 6.考察 ]
(2.2)において,{displaystyle M=P } とすると,

(6.1){displaystyle ;;;;;; oldsymbol{x}_{cdot} = z_1 oldsymbol{w}_{1} + z_2 oldsymbol{w}_{2} + cdots + z_P oldsymbol{w}_{P} }
{displaystyle ;;;;;;;;;;;;;;;; = langle oldsymbol{x}_{cdot},oldsymbol{w}_{1} 
angle oldsymbol{w}_{1} + langle oldsymbol{x}_{cdot},oldsymbol{w}_{2} 
angle oldsymbol{w}_{2} + cdots + langle oldsymbol{x}_{cdot},oldsymbol{w}_{P} 
angle oldsymbol{w}_{P} }

となります.これは主成分 {displaystyle z_1,z_2,cdots,z_P } は各サンプル {displaystyle  { oldsymbol{x}_{n} }_{n=1,cdots,N} }{displaystyle oldsymbol{w}_{1},oldsymbol{w}_{2},cdots,oldsymbol{w}_{P} } に射影したベクトルの長さであることを意味しています.{displaystyle M lt P } であれば,(6.1)で {displaystyle M+1 } 以降の項を打ち切り {displaystyle oldsymbol{x}_{cdot} } を以下のように近似することになります.

(6.2){displaystyle ;;;;;; oldsymbol{x}_{cdot} approx z_1 oldsymbol{w}_{1} + z_2 oldsymbol{w}_{2} + cdots + z_M oldsymbol{w}_{M} }
{displaystyle ;;;;;;;;;;;;;;;; = langle oldsymbol{x}_{cdot},oldsymbol{w}_{1} 
angle oldsymbol{w}_{1} + langle oldsymbol{x}_{cdot},oldsymbol{w}_{2} 
angle oldsymbol{w}_{2} + cdots + langle oldsymbol{x}_{cdot},oldsymbol{w}_{M} 
angle oldsymbol{w}_{M} }

主成分分析とは,「{displaystyle mathbb{R}^P } において測定データから計算される分散共分散行列 {displaystyle V } の固有ベクトル {displaystyle { oldsymbol{w}_1^* , oldsymbol{w}_2^*, cdots ,oldsymbol{w}_P^*  } } (固有値の大きい順)による直交座標系に線形変換し,そのうちの {displaystyle M (le P)} 次元のみを用いて(すなわち次元削減をおこない)測定データを近似すること」であることがわかります.
=================================================================================

以上,主成分分析の基礎をまとめました.



参考文献
[1] 京都大学 加納 学 先生のノート http://manabukano.brilliant-future.net/document/text-PCA.pdf
[2] 東京工業大学 渡辺澄夫先生のノート 第4回 主成分分析 http://www.ocw.titech.ac.jp/index.php?module=General&action=T0300&GakubuCD=4&GakkaCD=342200&KeiCD=&KougiCD=201602395&Nendo=2016&lang=JA&vid=05
[3] 統計科学研究所のページ http://www.statistics.co.jp/reference/software_R/statR_9_principal.pdf
[4] Cross Validated Stack Exchange http://stats.stackexchange.com/questions/153928/why-are-principal-component-scores-uncorrelated