[data analysis] 직교행렬과 주성분

직교행렬과 주성분(Principal Component)

식 1의 행렬 X를 평균-편차 형태로 가정합니다.

$$\tag{식 1}X=\left[\begin{matrix} x_1 & x_2 & \cdots &x_n \end{matrix}\right]$$

주성분분석의 목적은 식 2와 같이 X=PY 형태로 변환가능한 p×p차원의 직교 행렬 P를 발견하는 것입니다.

\begin{align}\tag{식 2} \begin{bmatrix}x_1 \\ x_2 \\ \vdots \\ x_p \end{bmatrix} & = \begin{bmatrix}u_1 & u_2 & \cdots & u_p \end{bmatrix}\begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_p \end{bmatrix}\\ P&=\begin{bmatrix} u_1 & u_2 & \cdots & u_p \end{bmatrix} \end{align}

식 2에서 새로운 변수 $y_1, y_2, \cdots, y_p$는 상관성이 없고 내림차순으로 정렬된 형태입니다. Y는 P에 대해 X의 좌표벡터가 됩니다. P가 직교행렬이므로 식 3이 성립합니다(정규직교의 특징 참조).

$$\tag{식 3}P^{-1}PY=P^TPY=P^TX \rightarrow Y=P^TX$$

식 3으로 부터 새로운 변수인 Y 역시 평균-편차 형태이므로 Y의 공분산 행렬을 유도 할 수 있습니다(식 4).

$$\tag{식 4}YY^T=(P^TX)(P^TX)^T=P^TXX^TP=P^TSP$$

S는 변수 X에 대한 공분산 행렬로서 대칭행렬입니다. 식 4에서 $YY^T=A$라 하고 S에 대해 정리하면 식 5와 같습니다.

$$\tag{식 5}S=PAP^T$$

S가 대칭행렬이고 P가 직교행렬이므로 위 식은 스펙트럴 분해와 동일한 형태가 됩니다. 그러므로 행렬 $A=YY^T$는 행렬 S의 고유값들을 내림차순으로 정렬한 값들을 대각요소로 하는 대각행렬이 됩니다. 결과적으로 새로운 변수의 공분산 행렬은 원래의 변수의 공분산 행렬의 고유값에 의한 대각행렬이 됩니다.

공분산 행렬 S의 고유값에 대응하는 단위 고유벡터들 $u_1, u_2, \cdots, u_p$은 자료의 주성분(principal component)라고 합니다. 고유값들을 내림차순으로 정렬하기 때문에 제1 주성분은 S의 가장 큰 고유값에 대응하는 고유벡터, 두번째 주성분은 두번째로 큰 고유값에 대응하는 고유벡터가 됩니다. 즉, P는 S의 단위고유벡터들을 성분으로 하는 고유행렬이 되므로 제1 주성분에 의한 새로운 변수는 식 6과 같이 나타낼 수 있습니다.

$$\tag{식 6}y_1 = u_1^TX$$

import numpy as np
import numpy.linalg as la
import pandas as pd
from sympy import *
np.set_printoptions(precision=4, suppress=True)

예 1)

어떤 데이터의 공분산 행렬 S가 다음과 같다면 그 자료의 주성분?

$$S=\begin{bmatrix}2382.78 & 2611.84 & 2136.2\\2611.84 & 3160.47 & 2553.9\\2136.2 & 2553.9 & 2650.71\end{bmatrix}$$

S=np.array([[2382.78, 2611.84, 2136.20],[2611.84, 3160.47, 2553.90],[2136.20, 2553.90, 2650.71]])
print(S)

[[2382.78 2611.84 2136.2 ]
 [2611.84 3160.47 2553.9 ]
 [2136.2  2553.9  2650.71]]

S를 식 5에 대입하기 위해 고유행렬과 고유벡터를 계산하기 위해 np.eig() 함수를 적용합니다.

d, p = la.eig(S)
print(d.round(2))

[7635.72  125.24  433.  ]

print(p.round(2))

[[-0.54 -0.7  -0.46]
 [-0.63  0.7  -0.33]
 [-0.56 -0.11  0.82]]

식 5의 A 즉, 새로운 변수들의 공분산 행렬(YY^T)은 다음의 결과와 같이 S의 고유값을 대각요소로 하는 대각행렬입니다.

D=np.diag(d)  #새로운 변수의 공분산 행렬
print(D.round(2))

[[7635.72    0.      0.  ]
 [   0.    125.24    0.  ]
 [   0.      0.    433.  ]]

S의 고유값 중 가장 큰 7635.72에 대응하는 고유벡터 p[:,0]이 제1 주성분이 됩니다. 즉, 행렬 p의 각 열이 주성분이 됩니다.

예 2)

다음은 5명의 남아의 몸무게와 키에 대한 자료입니다.

	0	1	2	3	4
Weight	120	125	125	135	145
Height	61	60	64	68	72

자료 내에서 대부분의 변화를 설명하는 새로운 단일한 변수를 발견하기 위한 주성분 분석?

data=np.array([[120, 125, 125, 135, 145], [61, 60, 64, 68, 72]])
print(data)

[[120 125 125 135 145]
 [ 61  60  64  68  72]]

S=np.cov(data, rowvar=True) # 위 data의 공분산 행렬
print(S)

[[100.   47.5]
 [ 47.5  25. ]]

d,P=la.eig(S)
print(d.round(2))

[123.02   1.98]

print(P.round(2)) # data 공분산행렬에 대한 주성분

[[ 0.9  -0.44]
 [ 0.44  0.9 ]]

제1 주성분에 의한 새로운 단일한 변수 y₁은 다음과 같이 결정됩니다.

\begin{align} Y& =u_1^TX\\ u_1 & = \text{P[:,0]}\end{align}

weight, height=symbols("w h")
x=Matrix(2,1,[weight, height]); x

$\left[\begin{matrix}w\\h\end{matrix}\right]$

pc1=Matrix(P[:,0].round(2))
pc1

$\left[\begin{matrix}0.9\\0.44\end{matrix}\right]$

y1=pc1.T*x
y1

$\left[\begin{matrix}0.44 h + 0.9 w\end{matrix}\right]$

주성분에 의해 새로운 변수를 생성할 수 있습니다. 즉, weight, height 두 변수는 하나의 새로운 변수 y로 대체 될 수 있습니다. 결과적으로 자료의 차원이 축소됩니다. 위의 분석은 python의 모듈 sklearn.decomposition.PCA() 클래스로 실행될 수 있습니다.

• sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False, svd_solver=’auto’, tol=0.0, iterated_power=’auto’, random_state=None)
  • n_components : 산출하고자 하는 주성분의 수로 정수 입니다. 이 값을 정하지 않는다면 모든 주성분은 다음과 같이 지정됩니다. n_components=min(n_sample, n_features), 결국 원 자료의 공분산행렬의 모든 고유값에 대응되는 모든 고유값을 사용합니다.

PCA() 클래스를 예 2에 적용합니다.

from sklearn.decomposition import PCA

data1=data.T
print(data1)

[[120  61]
 [125  60]
 [125  64]
 [135  68]
 [145  72]]

pca=PCA() #PCA 클래스(pca) 생성
pca.fit(data1) #pca 모형생성
pc=pca.components_#고유벡터, 주성분 벡터
print(pc)

[[ 0.89990119  0.43609386]
 [ 0.43609386 -0.89990119]]

n_components를 지정하지 않았으므로 모든 주성분을 고려합니다.

eigvalY=pca.explained_variance_#고유값
print(eigvalY)

[123.01859218   1.98140782]

위 결과에 의한 대각행렬은 weight와 height의 공분산 행렬(XX^T)을 기반으로 하는 새로운 변수(Y)의 공분산행렬이 됩니다.

\begin{align} X&=PY \Leftrightarrow P^{-1}X=P^TY=X\\ \tag{식 7}YY^T& =P^TX(P^TX)^T\\ & =P^TXX^TP\\& =P^TSP\\ S & =P(YY^T)P^T \end{align}

D=np.diag(eigvalY)
print(D)

[[123.01859218   0.        ]
 [  0.           1.98140782]]

위 결과 D는 새로운 변수 Y의 공분산행렬로 Y를 구성하는 여러 벡터들의 분산만 존재하며 그들사이의 공분산은 0입니다. 이 모델의 속성 .explained_variance_ratio_을 사용하여 각각의 주성분이 전체 분산에 영향을 주는 비율을 확인할 수 있습니다.

위 결과는 $S=PDP^T=P(YY^T)P^T$에서 새로운 변수 Y 의 공분산 행렬로 각 주성분에 대한
분산을 나타냅니다. (D는 대각 행렬외에 모두 0이므로 공분산은 존재하지 않습니다. 즉, 원래의 자료를 직교행렬을 사용하여 변환하는 것으로서 변수들 간의 상관성은 0이 됩니다.)

pca.explained_variance_ratio_

array([0.98414874, 0.01585126])

위의 경우 제 1 주성분에 의한 분산은 전체의 약 98%를 차지합니다.