PCA
개요
주성분 분석(PCA)는 데이터를 대표하는 성분을 찾고 설명되는 분산이 큰 상위 N개의 성분을 추려내 데이터의 차원을 축소하는 기법입니다. 다차원의 자료를 2차원으로 시각화하거나 차원을 축소하여 소위 ‘차원의 저주’를 피하려는 목적으로 사용합니다.
제1 주성분 벡터로 각 샘플을 projection시키면, 각 projection vectors 크기의 분산은 최대화합니다. 분산이 최대화한다는 것은 projection vectors의 크기가 제각각이라는 뜻이며 그만큼 각 샘플이 제1 성분으로 projection되었을 때 서로 잘 구분된다는 의미입니다. 샘플 벡터를 주성분으로 사영한 벡터 크기의 분산이 최대화함을 이 링크에서 시각적으로 확인할 수 있습니다.
따라서 주성분 벡터를 찾고 사용 목적에 적합한 수만큼 주성분을 고르는 것이 PCA의 핵심입니다. 일반적으로 선형변환을 통해 주성분 벡터를 얻습니다. 이번 포스트에서는 고유값 분해와 SVD로 주성분 벡터를 구하고 데이터 샘플을 차원-축소해보겠습니다. 그 결과를 라이브러리를 썼을 때와 비교하여 검증하려 합니다.
필요한 라이브러리와 데이터 불러오기
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import time
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler
from packages.GramSchmidt import gram_schmidt
from packages.SVD import svd_thin
검증에 사용할 샘플 데이터를 생성합니다. 샘플 개수는 1,000,000개 차원은 20-차원, 레이블은 0과 1입니다.
data = make_classification(n_redundant=0, n_samples=10**6, weights=[0.9], random_state= 42, )
X = data[0]
y = data[1]
print('샘플 개수:', X.shape[0])
print('차원 수:', X.shape[1])
샘플 개수: 1000000
차원 수: 20
먼저 샘플을 scaling해줍니다. 특히 centering 즉 샘플 행렬의 열(features)의 평균을 0으로 만들어주지 않으면 링크와 같이 주성분이 잘못 산출됩니다. 샘플의 개수를 n, 차원 수를 p라 하면, 샘플 행렬은 아래와 같이 표현됩니다.
$\large \mathbf{X} = \begin{pmatrix} x_{1,1} & \cdots & x_{1,p} \\ \vdots & \ddots & \vdots \\ x_{n,1} & \cdots & x_{n,p} \end{pmatrix} = \begin{pmatrix} X_1 & \cdots&X_p \end{pmatrix} _{n\times p}$
Standard scaling을 통해 샘플 행렬 각 열(features)의 평균을 0, 표준편차를 1로 만들어줍니다.
$\large \mathbf{X}_{\normalsize std} = \begin{pmatrix} \frac{X_1 - \mu_1}{\sigma_1} & \frac{X_2 - \mu_2}{\sigma_2} & \cdots &\frac{X_p - \mu_p}{\sigma_p} \end{pmatrix} $, where $\mu_i$ and $\sigma_i$ are mean and standard deviation of a vector $X_i$ for $i \in {1,2,\cdots,p}$.
scaler = StandardScaler()
X_std = scaler.fit_transform(X)
print('means of each columns are all 0?:', np.allclose(X_std.mean(axis=0), np.zeros(X_std.shape)))
print('standard deviations of each columns are all 1?:', np.allclose(X_std.std(axis=0), np.ones(X_std.shape)))
means of each columns are all 0?: True
standard deviations of each columns are all 1?: True
본 포스팅에서는 샘플 데이터의 차원을 20에서 2차원으로 줄여보겠습니다.
dim_reduced = 2
Scikit-learn 라이브러리 사용
먼저 라이브러리를 사용해 PCA를 실행해봅시다.
start = time.time()
pca = PCA(n_components= dim_reduced, svd_solver='full')
X_pca = pca.fit_transform(X_std)
end = time.time()
print('runtime:', end - start)
runtime: 0.9079928398132324
백만 개 샘플을 PCA하는 데 1초 가량 걸립니다.
2차원으로 축소하였으므로 2개의 주성분이 주어집니다. components_ attribute로 확인할 수 있습니다.
print('shape:',pca.components_.T.shape)
print(pca.components_.T)
shape: (20, 2)
[[-0.70599712 -0.02179824]
[-0.03820457 -0.29716737]
[-0.00511981 0.24191311]
[ 0.00316976 -0.1552799 ]
[-0.00607306 0.19315766]
[ 0.01378161 -0.31379142]
[ 0.01395261 -0.19963659]
[-0.01303674 -0.10695125]
[ 0.00663818 0.12910551]
[-0.01650511 -0.02872533]
[ 0.00849686 0.44358762]
[-0.01663547 -0.19563718]
[-0.00496126 0.17365063]
[-0.01661978 0.05243992]
[ 0.02078452 -0.28187502]
[-0.01667663 0.14918825]
[-0.00635741 0.06491868]
[ 0.7051932 -0.00858833]
[-0.00136236 -0.16749529]
[-0.02155668 0.47118897]]
$\large Comp = \begin{bmatrix} V_{comp1} & V_{comp2} \end{bmatrix}$
scaling을 해주었기 때문에 각 성분의 크기는 1이고 서로 다른 성분끼리는 직교합니다.
$\large \lVert V_{comp1} \rVert = \lVert V_{comp2} \rVert = 1$
$\large V_{comp1} \perp V_{comp2}$
print('perpendicular to each component?:',np.allclose((pca.components_).dot(pca.components_.T), np.eye(2)))
perpendicular to each component?: True
scaling된 샘플 행렬에 주성분을 곱하면 차원-축소된 샘플을 얻습니다.
$\large X_{pca} = X_{std} \times Comp = X_{std} \times \begin{bmatrix} V_{comp1} & V_{comp2} \end{bmatrix}$
print('X_pca == X_std.dot(pca.components_.T)?:', np.allclose(X_pca, X_std.dot(pca.components_.T)))
X_pca == X_std.dot(pca.components_.T)?: True
앞서 설명한 바와 같이 주성분으로 사영된 벡터 크기의 분산 즉 ‘설명되는 분산’이 큰 순으로 주성분이 정렬됩니다.
pca = PCA()
X_pca = pca.fit_transform(X_std)
pca.explained_variance_ == X_pca.var(ddof=1, axis=0)
print('pca.explained_variance_ == X_pca.var(ddof=1, axis=0)?:', np.allclose(pca.explained_variance_, X_pca.var(ddof=1, axis=0)))
pca.explained_variance_ == X_pca.var(ddof=1, axis=0)?: True
print('explained variance:', pca.explained_variance_)
explained variance: [1.07743561 1.00654863]
분산 설명량을 비율로 따지면 아래와 같습니다.
print('explained variance ratio:', pca.explained_variance_ratio_)
explained variance ratio: [0.05387173 0.05032738]
주성분 두 개의 분산 설명량은 전체 분산의 약 10%입니다. 차원이 축소되면서 부득이 90% 가량의 정보를 상실하게 됩니다. 대신 2차원으로 축소되었으므로 평면에 시각화하여 나타낼 수 있습니다. 1,000개 샘플만 뽑아서 평면 위에 표시해봅시다. 두 주성분 벡터 역시 2차원으로 축소할 수 있으며 아래와 같이 서로 직교합니다.
fig, ax = plt.subplots(1,1, figsize = (16,12))
ax.set_aspect(1)
ax.axline((.0, .0), (pca.transform(pca.components_)[0,0], pca.transform(pca.components_)[0,1]), c='b', label='1st component')
ax.axline((.0, .0), (pca.transform(pca.components_)[1,0], pca.transform(pca.components_)[1,1]), c='r', label='2nd component')
show_n_samp = 1000
sns.scatterplot(x = X_pca[:show_n_samp, 0], y=X_pca[:show_n_samp, 1], hue=y[:show_n_samp], alpha=.7 , s=20)
plt.title('PCA of the standardized dataset', fontsize = 16)
plt.show()
numpy로 구현
with EVD of covariance matrix
샘플의 공분산 행렬을 고윳값 분해(eigen decomposition)하여 주성분 벡터를 구할 수 있습니다. 공식은 아래와 같습니다.
$\displaystyle \mathrm{Cov}(X) = \begin{pmatrix} \mathrm{cov}[X_1, X_1] & \cdots & \mathrm{cov}[X_1, X_p] \\ \vdots & \ddots & \vdots \\ \mathrm{cov}[X_p, X_1] & \cdots & \mathrm{cov}[X_p, X_p] \end{pmatrix} _{p\times p} = \begin{pmatrix} \mathrm{var}[X_1] & \cdots & \mathrm{cov}[X_1, X_p] \\ \vdots & \ddots & \vdots \\ \mathrm{cov}[X_p, X_1] & \cdots & \mathrm{var}[X_p] \end{pmatrix} _{p\times p}$
cov = np.cov(X_std, rowvar=False)
scaling된 행렬의 공분산 행렬은 본래 행렬과 전치 행렬의 곱으로 계산할 수 있습니다.
$\large \mathrm{Cov}(X_{std}) = \frac{1}{n-1} X_{std}^T X_{std}$
np.allclose(cov, (X_std.T).dot(X_std)/(X_std.shape[0]-1))
True
$\large \mathrm{cov}[X_i, X_j] = \mathrm{cov}[X_j, X_i] \ \ \normalsize \forall i, j \in {1,2,\dots,p}$ 이고 공분산 값은 실수이므로,
위에서 계산한 공분산 행렬 $\large \mathrm{Cov}(X_{std})$은 대칭행렬(symmetric)입니다.
즉 $\large \mathrm{Cov}(X_{std}) = \left [ \mathrm{Cov}(X_{std}) \right ]^T$ 입니다.
print('cov is symmetric:',np.alltrue(cov == cov.T))
cov is symmetric: True
실수 대칭행렬의 고윳값은 실수이며 서로 다른 고윳값에 대응하는 고유벡터는 직교합니다.
위 공분산 행렬의 고유 벡터를 열 벡터로 하는 행렬을 $\large Q$, 내림차순으로 정렬된 고윳값을 대각 성분으로 하는 대각행렬을 $\large \Lambda$라 하면,
$\large Q Q^T = I_{p\times p}$
$\large \mathrm{Cov}(X_{std}) = Q \Lambda Q^{-1} = Q \Lambda Q^T$
evalues, evectors = np.linalg.eig(cov)
evalues_mat = np.sort(evalues)[::-1] * np.eye(evalues.size)
evectors = evectors[:,evalues.argsort()[::-1]]
앞서 구한 공분산 행렬은 그 고유벡터 $\large Q$와 고윳값 $\large \Lambda$로 대각화(diagonalizing)할 수 있습니다.
np.allclose(cov, evectors.dot(evalues_mat.dot(evectors.T)))
True
행렬 $\large Q$의 앞 N개 열이 바로 제1 ~ 제N 주성분 벡터입니다. 즉
$\large X_{pca} = X_{std} \times \begin{bmatrix} \mathbf{q_{1}} & \mathbf{q_{2}} \end{bmatrix}$
X_pca_eig = X_std.dot(evectors[:,:dim_reduced])
comp_reduced = (evectors[:,:dim_reduced].T).dot(evectors[:,:dim_reduced])
고윳값 분해로 구한 주성분으로 PCA한 결과를 scikit-learn을 쓴 결과와 비교하면, 상하 반전된 점(y-성분의 부호가 서로 반대)을 제외하고 수치가 동일함을 확인할 수 있습니다.
print('x-components:',np.allclose(X_pca_eig[:,0],X_pca[:,0]))
print('y-components:',np.allclose(X_pca_eig[:,1],-X_pca[:,1]))
x-components: True
y-components: True
- Visualization: the same as the x-axis reflection of above scatter plot
fig, ax = plt.subplots(1,1, figsize = (16,12))
ax.set_aspect(1)
ax.axline((.0, .0), (comp_reduced[0,0], comp_reduced[0,1]), c='b', label='1st component')
ax.axline((.0, .0), (comp_reduced[1,0], comp_reduced[1,1]), c='r', label='2nd component')
show_n_samp = 1000
sns.scatterplot(x = X_pca_eig[:show_n_samp, 0], y=X_pca_eig[:show_n_samp, 1], hue=y[:show_n_samp], alpha=.7 , s=20)
plt.title('PCA of the standardized dataset: with numpy(EVD)', fontsize = 16)
plt.show()
with SVD of dataset
지난 포스트에서 다루었던 SVD를 이용하여 주성분 벡터를 찾아봅시다. 고윳값 분해 즉 EVD를 쓸 때는 샘플 행렬의 공분산 행렬을 계산해야 했지만, SVD는 샘플 행렬 자체에 적용됩니다. 지난 포스트에서 구현한 SVD 알고리즘으로 scaling된 샘플 행렬을 (reduced) SVD합니다.
$\large X_{standardized} = U S V^*$
# with numpy SVD
u, s, v = svd_thin(X_std, )
계산된 좌측 특이벡터 $U$와 그 켤레 전치 행렬은 직교합니다. 즉
$\large U^* U = I$ 입니다.
np.allclose(np.conjugate(u).T.dot(u), np.eye(u.shape[1]))
True
이를 공분산 행렬 계산식에 대입하여 보면 우측 특이벡터와 앞서 구한 고유 벡터 $Q$가 같음을 알 수 있습니다.
$\mathrm{Cov}(X_{std}) =\frac{1}{n-1} X_{std}^T X_{std} =\frac{1}{n-1} V S^* U^* U S V^* =\frac{1}{n-1} V S^* I S V^* =\frac{1}{n-1} V S^2 V^* = Q \Lambda Q^T $
$\large V = Q$
print(np.allclose(v, evectors))
True
그렇다면 EVD를 활용할 때와 마찬가지로 우측 특이벡터 $V$의 앞 N개 열 벡터가 제1 ~ 제N 주성분입니다.
$\large X_{pca} = X_{std} \times \begin{bmatrix} \mathbf{v_{1}} & \mathbf{v_{2}} \end{bmatrix}$
X_pca_svd = X_std.dot(v[:,:dim_reduced])
comp_reduced_svd = (v[:,:dim_reduced].T).dot(v[:,:dim_reduced])
SVD로 구한 주성분으로 PCA한 결과를 scikit-learn으로 PCA한 결과와 비교하면, 상하 반전된 점(y-성분의 부호가 서로 반대)을 제외하고 수치가 동일함을 확인할 수 있습니다.
print('x-components:',np.allclose(X_pca_svd[:,0],X_pca[:,0]))
print('y-components:',np.allclose(X_pca_svd[:,1],-X_pca[:,1]))
x-components: True
y-components: True
fig, ax = plt.subplots(1,1, figsize = (16,12))
ax.set_aspect(1)
ax.axline((.0, .0), (comp_reduced_svd[0,0], comp_reduced_svd[0,1]), c='b', label='1st component')
ax.axline((.0, .0), (comp_reduced_svd[1,0], comp_reduced_svd[1,1]), c='r', label='2nd component')
show_n_samp = 1000
sns.scatterplot(x = X_pca_svd[:show_n_samp, 0], y=X_pca_svd[:show_n_samp, 1], hue=y[:show_n_samp], alpha=.7 , s=20)
plt.title('PCA of the standardized dataset: with numpy(SVD)', fontsize = 16)
plt.show()
Leave a comment