[재업]차원 축소와 주성분분석(Principal Component Anaylse, PCA)

실제 현상은 여러 변수가 영향을 미치지만, 그래프로 표현할 수 있는 것은 3차원 까지이다. 그렇다면 적당히 데이터의 분포를 잘 보여주도록 차원을 축소하여 그래프로 한눈에 볼 수는 없을까? 이에 대한 답변중 하나는 주성분 분석 (Pricipal Component Anaylse, PCA)이다. 이번 포스트에서는 PCA의 수학적인 원리를 알아보고, 이를 파이썬을 이용해 구현해볼 것이다. (scikit-learn에도 주성분 분석할 수 있도록 추상화된 함수를 제공하지만, 원리를 공부하기 위해 사용하지 않았다.)

 

0.Prequestion

아래 간단하게 고유값 분해와 특이값 분해에대해 요약을 하였는데, 이해가 되지 않는다면 공부하고 오는 것을 추천한다.

고유값과 고유벡터

정사각 행렬 \(A\)와 영벡터가 아닌 열벡터 \(v\), 스칼라 \(\lambda\)에 대해서

$$Av=\lambda v$$

를 만족할 때, \(\lambda\)를 고유값, \(v\)를 고유벡터라 한다.

*고유값과 고유벡터는 \(det(Av-\lambda I)=0\)임을 이용해 구할 수 있다.

 

고유값 분해(Eigen Value Decomposition, EVD)

고유값 행렬 \(\Lambda\), 고유 열벡터를 모아 만든 행렬 \(V\)에 대해서 \(AV=V\Lambda\) 이므로 \(V\)의 역행렬이 존재한다면

$$A=V\Lambda V^{-1}$$

만약 \(A\)가 \(A=A^T\)인 대칭행렬이라면 \(A^T = {(V^{-1})}^T \Lambda^T V^T \)이고, \(\Lambda\)가 대각행렬이니 \(\Lambda = \Lambda^T\)이다. 이를 계수비교하면 \(V^T = V^{-1}\)인 직교행렬임을 알 수 있고 아래와 같이 쓸 수 있다.

$$A=V\Lambda V^T = Q\Lambda Q^T$$

*대칭행렬인경우 \(V\)대신 \(Q\)로 쓰기도 한다.

 

특이값 분해(Singularity Value Decomposition, SVD)

고유값 분해를 임의의 \(m \times n\)행렬로 확장한 것으로, \(m \times n\) 행렬 \(A\),  \(m \times m\) Unitary Matrix \(U\), \(m \times n\) 특이값 행렬 \(S\), \(n \times n\) Unitary Matrix \(V\)에 대해 다음과 같이 행렬 \(A\)의 분해가 가능하다.

$$A=USV^T$$

\(AA^T\), \(A^T A\)는 대칭행렬이므로  항상 고유값 분해가 가능하다.(전치 시켜보면 대칭행렬의 정의를 만족한다.) 이를 고유값 분해하면

\begin{array}{cl} A A^T &= Q_1 S_1 Q_{1}^T \\ &= USV^T VS^T U^T\\ &=USS^T U^T \end{array}

따라서 \(Q_1 = U\) 임을 알 수 있다. 같은 방법으로

\begin{array}{cl} A^T A &= Q_2 S_2 Q_{2}^T\\ &=VS^T U^T USV^T\\ &=VS^T SV^T\end{array}

\(Q_2 = V\) 임을 알 수 있다.

또한 \(S^T S\), \(SS^T\)가 고유값행렬에 대응되므로, 고유값에 루트를 해주면 특이값이 되는 것을 알 수 있다. \(m>n\) 이면 \(S^T S\), \(m<n\) 이면 \(SS^T\)으로 부터 고유값을 얻으면 된다.

1.어떻게 차원을 축소해야하는가?: 공분산(Covariance)과 공분산 행렬

http://matrix.skku.ac.kr/math4ai-intro/W12/

 여러개의 변수가 존재하여도 모든 변수가 똑같은 정도로 데이터의 분포에 기여하지는 않을 것이다. (만약 그렇다면 차원을 축소할 의미가 없다.) 쉽게 차원을 줄일 수 있는 방법은 다차원 공간상의 점을 한 직선에 정사영 시키는 것이며, 이러한 과정은 단위백터에 내적하는 과정으로 이루어진다. 그렇다면 어떤 단위벡터로 내적을 해야하는가? 우리는 원래 데이터의 분포를 보기위해서 차원을 축소하는 것이므로 차원 축소후에도 분포가 잘 드러나야 한다. 이러한 분포를 수치로 나타내는 것이 분포도라고 하므로, 어떤 단위벡터랑 내적을 할 때, 분포도가 최대가 되는가라는 문제가 된다. 잠시 단위벡터는 미루어 두고, 데이터의 분포를 어떻게 볼지를 생각해보자.

 

 분포도라고 하면, 중학교때 배운 분산과 표준편차가 생각날 것이다. 그러나 이러한 분산과 표준편차는 변수가 1개에 대한 분포만 보여준다. 따라서 이러한 개념을 이변수로 확장한 공분산에 대해서 알아야할 필요가 있으며, 공분산을 구하는 식은 아래와 같다.

$$ Cov(X,Y) = \dfrac{\displaystyle\sum_i^N (x_i - \bar{x}) (y_i - \bar{y}) }{N}\ (x_i \in X, y_i \in Y) $$

여기서 \(X, Y\)는 확률변수이며, \(N\)은 원소의 숫자, 변수 위에 붙은 bar는 평균을 의미한다. 이 식에서 \(X=Y\)인 경우에는 우리가 알고 있는 분산식이 되므로, 공분산은 분산의 일반화된 개념임을 알 수 있다. 여기서 우리는 공분산을 편하게 계산하기 위해 데이터의 평균을 0으로 맞춰주는 전처리를 하였다고 하자. 그러면 공분산 식은 아래와 같이 변경된다.

$$ Cov(X,Y) = \dfrac{1}{N}\displaystyle\sum_i^N x_i y_i \ (x_i \in X, y_i \in Y) $$

이를 쉽게 계산하기 위해 행렬을 사용하자. 먼저 다음과 같이 Centering 된 데이터 행렬 \(D\)가 주어졌다고 생각하자

$$ \begin{equation} D = \begin{pmatrix} x_0 & y_0 \\ x_1 & y_1 \\ \vdots & \vdots \\ x_{N-1} & y_{N-1} \end{pmatrix} \end{equation} $$

여기서 \(x, y\)는 변수이며, 이러한 변수를 \(N\)개 수집한 것이다. 쉽게 비유하자면, 반 학생들의 키와 몸무게를 조사한 뒤, 데이터 행렬을 만들었다고 할 때, \(x\)는 키, \(y\)는 몸무게 이며, 학생은 0번부터 \(N-1\)번까지 \(N\)명이 있는 것이다. 그러면 공분산 행렬 \(\Sigma\)는 아래와 같다.

$$ \begin{equation} \Sigma = \dfrac{1}{N}D^T D = \begin{pmatrix} Var(X) & Cov(Y, X) \\ Cov(X, Y) & Var(Y) \end{pmatrix} \end{equation} $$

2.정사영했을 때 공분산이 최대가 되는 백터: 공분산행렬의 고유 벡터 

위에서 차원 축소후 분산이 최대가 될 때, 가장 차원 축소를 잘한 것이라고 하였다. 그러면 어느 직선 (단위 벡터)에 사영을 시켜야할까? 결과부터 예기하면 공분산 행렬의 고유벡터로 사영시킬 때, 가장 데이터의 분산이 커지는데 지금부터 이를 증명하고자 한다. 우선 \(N\)개의 \(n\)차원 데이터를 모아둔 데이터 행렬 \(D\)를 임의의 단위벡터 \(\vec{e}\)로 사영시킨다고 생각을 하자.

$$ D \in \mathbb{R}^{N \times n} $$

$$ \vec{e} \in \mathbb{R}^{N \times 1} $$

그럼 사영시킨 후, 데이터의 분산은 분산의 정의에 따라 다음과 같이 구할 수 있다.

$$ Var(D\vec{e}) = \dfrac{1}{N} \left\vert D\vec{e} - E( D\vec{e} ) \right\vert^2 $$

여기서 E는 기대값 (Expection)을 나타내는 기호이며, 이때 데이터는 Centering 되었으므로, 기대값또한 0이다. 그러므로 아래와 같이 식을 전개할 수 있다. 

$$ \begin{align} Var(D\vec{e}) &= \dfrac{1}{N} \left\vert D\vec{e} - E(D\vec{e} ) \right\vert^2\\ &= \dfrac{1}{N} \left\vert D\vec{e}) \right\vert^2\\ &= \dfrac{1}{N}( D\vec{e} )^T( D\vec{e} )\\ &= \dfrac{1}{N} \vec{e}^T D^T D \vec{e}\\ &= \vec{e}^T \Sigma \vec{e}\\ &= \dfrac{1}{N}\left( \sum_{i}^{N} x_i e_i \right)^2 \end{align} $$

 

이 분산값의 최대를 찾아야 하는데, \(\vec{e}\)가 단위벡터라는 제약이 있으므로 라그랑주 승수법(Lagrange Multiplier)를 사용하자

$$ \textrm{제약 조건}: |\vec{e}|^2 = 1 $$

$$ \dfrac{\partial}{\partial \vec{e}} Var(D\vec{e}) = \lambda \dfrac{\partial}{\partial \vec{e}} |\vec{e}|^2 $$

스칼라를 벡터로 미분하는 것은 좀더 일반화된 그래디언트(gradient) 연산이며, 아래와 같이 성분별로 미분하면 된다.

$$ \begin{align} \dfrac{\partial}{\partial \vec{e}} Var(D\vec{e}) &= \dfrac{2}{N}\left( \sum_{i}^{N} x_i e_i \right) ( x_1, x_2, \cdots, x_N )^T\\ &= \dfrac{2}{N} ( x_1, x_2, \cdots, x_N )^T \left( \sum_{i}^{N} x_i e_i \right)\\ &= \dfrac{2}{N} ( x_1, x_2, \cdots, x_N )^T (x_1, x_2, \cdots, x_N) (e_1, e_2, \cdots, e_N)^T\\ &= 2\left[ \dfrac{1}{N} ( x_1, x_2, \cdots, x_N ) (x_1, x_2, \cdots, x_N)^T \right]^T (e_1, e_2, \cdots, e_N)^T\\ &= 2\Sigma^T (e_1, e_2, \cdots, e_N)^T\\ &= 2\Sigma \vec{e}\ \left(\because \Sigma\ \textrm{is symmetric metrix},\ \Sigma^T = \left[\dfrac{1}{N}D^T D \right]^T = \dfrac{1}{N}D^T D = \Sigma \right) \end{align} $$

$$ \begin{align} \dfrac{\partial}{\partial \vec{e}} |\vec{e}|^2 &= \left[ \dfrac{\partial}{\partial e_1}\left( \sum_{i}^{N} e_i^2 \right),\ \dfrac{\partial}{\partial e_2}\left( \sum_{i}^{N} e_i^2 \right),\ \cdots,\ \dfrac{\partial}{\partial e_N}\left( \sum_{i}^{N} e_i^2 \right) \right]^T\\ &= \left[ 2e_1, 2e_2, \cdots, 2e_N \right]^T\\ &= 2\vec{e} \end{align} $$

대입하면 다음과 같다.

$$ 2\Sigma \vec{e} = 2\lambda \vec{e} $$

$$ \therefore \Sigma \vec{e} = \lambda \vec{e} $$

이는 고유값과 고유벡터의 정의이므로, 공분산 행렬의 고유벡터로 데이터를 사영시켰을 때 데이터의 분산이 최대가 되며, 이때의 분산은 고유값과 같다. 따라서 공분산 행렬을 만든 후 EVD를 통해 고유벡터와 고유값을 얻을 수 있으며, 공분산 행렬은 대칭 행렬 (Symetric Metrix)이므로 아래와 같이 분해된다.

$$ \Sigma = Q\Lambda Q^T $$

여기서 \(Q\)는 \(n\times n\)의 고유 열벡터 행렬이며, \(\Lambda\)는 고유값이 대각성분인 \(n\times n\)의 대각행렬이다. 따라서 데이터를 특정 직선으로 사영시킨 PC 점수 (Principal Component Score) 행렬 \(T\)는 아래와 같이 구할 수 있으며, \(T\)의 열이 PC 점수들이다.

$$\therefore T = DQ $$

3.SVD로 PC 점수 구하기: 더빠른 연산을 위해

데이터 행렬 \(D \in \mathbb{R}^{N\times n}\)를 SVD하면 다음과 같다.

$$ D = USV^T $$

여기서 U는 \(N \times N\)의 직교행렬이며, S는 주대각성분이 특이값인 \(N \times n\)인 유사 대각행렬, V는  \(n \times n\)의 직교행렬이다. (실수인 경우 Unitary Matrix는 직교행렬이다.) 이때

$$ \Sigma = \dfrac{1}{N}D^T D = \dfrac{1}{N} (USV^T)T(USV^T) = \dfrac{1}{N} VS^T U^T U S V^T = \dfrac{1}{N} V S^T S V^T $$

이므로 \( \Sigma = Q\Lambda Q^T \)에서 \(Q=V\)임을 알 수 있으며, 따라서 PC 점수 행렬 \(T\)는 아래와 같다.

$$\therefore T = DQ = USV^T V = US $$

또한 \( \Lambda = \dfrac{1}{N}S^T S \)이라는 것도 알 수 있으며, 따라서 공분산 행렬의 고유값 \(\lambda_i\)와 데이터 행렬의 특이값 \(\sigma_i \) 사이의 관계는 다음과 같다.

$$ \lambda_i = \dfrac{\sigma_i^2}{N} $$

따라서 공분산 행렬을 만든 뒤, EVD를 하지 않아도 데이터 행렬을 SVD하여 PCA를 수행할 수 있다. 

 

그러면 왜 SVD를 사용하는 것이 EVD하는 것보다 빠를까? 2에서 살펴보았듯이, \(n\)개의 차원이면 \(n\)개의 고유벡터가 나오며, \(n\)개의 PC가 도출된다. 그런데 실제 그래프를 그릴 때는 1~3개를 사용하므로 \(n\)개 고유벡터를 모두 구할 필요가 없다. 아래와 같이 SVD를 수행하면 모든 특이값을 구하지 않고, 특이값이 큰거 몇개만 계산할 수 있으므로 EVD를 사용하는 것보다 SVD를 사용하는 것이 효율적이다.

Full SVD

Thin SVD

Compact SVD

Truncated SVD

이미지 출처: https://darkpgmr.tistory.com/106

[선형대수학 #4] 특이값 분해(Singular Value Decomposition, SVD)의 활용

Thin SVD는 어짜피 0과 곱해질 부분을 제외하고 분해한 것이며, Compact SVD는 특이값중 0인 것까지 제외하고 분해한 것이다. 

 

4.PC 성분 고르기

2에서 살펴보았듯이, \(n\)개의 차원이면 \(n\)개의 고유벡터가 나오며, \(n\)개의 PC가 도출된다. 이제 여기에서 가장 잘 설명해주는 성분을 골라야 하는데, 당연히 분산이 큰 성분을 골라야 원 데이터의 분포를 잘 설명할 수 있다. 2에서 분산이 고유값과 같다고 하였으므로, 고유값이 큰 성분 1~3개를 골라 각 성분을 축으로 그래프를 그리면 될 것이다. 그런데 얼마나 차원을 축소하는 것이 허용될까? 이때 판단 기준이 되는 것이 해당 PC의 분산과 전체 분산에 대한 비율, 즉 \(\frac{\lambda}{\sum_{i=1}^{n} \lambda_i}\) 이며, 이를 그래프의 축에 몇 % explained라고 표시해준다.

 

 

5. 주의할점

만약 이와 같은 데이터가 주어졌을 때, Centering만 하는 것이 바람직할까? 다른 변수의 크기보다 칼로리의 크기가 훨씬 크므로 이대로 분석을 진행하면 칼로리의 영향이 지배적일 것이다. 따라서 이러한 경우는 데이터를 표준화 시켜주어야 하며, 이외의 경우에도 처음부터 데이터를 표준화시켜 분석을 하는 것이 편할 것이다. 

 

6.파이썬으로 PCA해보기: 성별에 따른 건강 데이터 분석

공공데이터 포털(https://www.data.go.kr/)에 있는 국민건강보험공단_건강검진정보를 사용하였다.

6.1 모듈 import 및 데이터 불러오기

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

#한 셀에 여러 입출력
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

#데이터 불러오기
origin = pd.read_csv("국민건강보험공단_건강검진정보_20211231.CSV", encoding="euc-kr")
origin.info()
origin

 

 

6.2 데이터 선별 (전처리)

다룰 데이터는 연속적인 데이터이므로, 음주 여부처럼 0, 1로 나오는 데이터나, 시도코드같은 데이터는 제외하여야 한다.

또 데이터가 너무 많으면 처리하는데 오래 걸리므로 적당히 신장, 체중, 식전혈당, 총콜레스테롤을 대상으로만 PCA를 진행할 것이다.

데이터 선택 및 결측값 제거

#데이터 선택 및 결측값 제거
data_selected = origin[["성별코드", "신장(5cm 단위)", "체중(5kg 단위)", "식전혈당(공복혈당)", "총 콜레스테롤"]]
data_selected = data_selected.dropna()
data_selected

데이터 표준화

#데이터 표준화 
data_std = pd.DataFrame(data_selected["성별코드"].reset_index(drop=True)) #drop=True로 하면 기존 인덱스를 버리고, 기본값(False)면 index Columns 생성함
for name in data_selected.columns[1:]:
    col = np.array(data_selected[name])
    col = (col - col.mean())/col.std()
    data_std[name] = pd.DataFrame(col)
data_std

남여 데이터 분리

다운 받은 페이지에서 같이 다운 받을 수 있는 국민건강정보데이터 건강검진정보 사용자 매뉴얼_20171027.hwp를 열어보면 아래와 같이 남자가 1, 여자가 2로 되어있음을 알 수 있다.

DataFrame의 행을 for로 순회하는 것은 매우 느리다. 반면 numpy ndarray로 성별 코드만 가져와 인덱스를 구한 뒤, drop 매서드를 사용하는 아래 방법은 빠르게 처리할 수 있다. 

#남여 분리
man_index = []
sex_code = np.array(data_std["성별코드"])
for i in range(len(data_std)):
    if sex_code[i] == 1:
        man_index.append(i)

data_woman = data_std.drop(man_index).drop("성별코드", axis=1)
data_man = data_std.drop(data_woman.index).drop("성별코드", axis=1)

data_man
data_woman

남자데이터

여자 데이터

PCA(Numpy)

우선 Numpy를 이용한 PCA 코드는 아래와 같다. 다만 Numpy의 경우 큰 행렬을 처리하려고 하면, 어마어마한 메모리를 사용하기 때문에 어쩔 수 없이 데이터를 슬라이싱하였다.

D_man = np.array(data_man[:10000].to_numpy())
U, s, VT = np.linalg.svd(D_man)
S = np.zeros([len(U), len(s)])
np.fill_diagonal(S, s)

T = U @ S

N = len(D_man.T)
Var_PC = np.array([x**2/N for x in s])
Var_PC_proportion = Var_PC / np.sum(Var_PC) * 100

var_table = pd.DataFrame(columns=[f"PC{i+1}" for i in range(len(s))])
var_table.loc["분산"] = Var_PC
var_table.loc["분산 비율"] = Var_PC_proportion
var_table

참고로 Numpy 공식 문서에 따르면 SVD를 진행한 후 특이값이 큰순서대로 정렬된다고 한다. 그러므로 항상 PC1이 가장 분산비율이 크다.

https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html

numpy.linalg.svd — NumPy v1.26 Manual

PC = T.T
plt.scatter(PC[0], PC[1], s=1, color="#78BAE5")

PCA(dask)

pip install dask

dask 라이브러리를 사용하면 적은 메모리를 사용하면서 큰 행렬을 다룰 수 있으며, CUDA를 사용하는 것도 가능하다. 남,여, 반복적으로 같은 작업을 해야하기 때문에 PCA 객체를 만들었다.

class PCA:
    def __init__(self, dataFrame:pd.DataFrame):
        if not type(dataFrame) is pd.DataFrame:
            raise ValueError("dataFrame must be pandas DataFrame")
        
        self.dataFrame = dataFrame
        
        #PC 구하기
        D = da.array(dataFrame.to_numpy())
        U, s, VT = da.linalg.svd(D)
        U = np.array(U)
        s = np.array(s)
        self.PC = (U @ np.diag(s)).T
        
        #분산 구하기
        N = len(self.PC)
        self.Var = np.array([x**2/N for x in s])
        self.Var_table = pd.DataFrame(columns=[f"PC{i+1}" for i in range(N)])
        self.Var_table.loc["분산"] = self.Var
        self.Var_table.loc["분산 (%)"] = self.Var / np.sum(self.Var) * 100
    
    def get_table(self):
        return self.Var_table

plt.figure(dpi=130)
plt.title("PCA")
plt.xlabel("PC1")
plt.ylabel("PC2")

plt.scatter(man.PC[0], man.PC[1], s=1, label="Man (71.4% explained)", color="#78BAE5", alpha=0.5)
plt.scatter(women.PC[0], women.PC[1], s=1, label="Women (68.8% explained)", color="#E7A46A", alpha=0.5)
plt.legend()

 

 

참고 자료

https://angeloyeo.github.io/2019/07/27/PCA.html

주성분 분석(PCA) - 공돌이의 수학정리노트 (Angelo's Math Notes)

http://matrix.skku.ac.kr/math4ai-intro/W12/

주성분 분석

https://ko.wikipedia.org/wiki/%EC%A3%BC%EC%84%B1%EB%B6%84_%EB%B6%84%EC%84%9D

주성분 분석 - 위키백과, 우리 모두의 백과사전

https://blog.dask.org/2020/05/13/large-svds

Large SVDs