[Python] PCA (Principal Component Analysis)

들어가기에 앞서..

PCA를 이해하기 위해서는 선형대수학 개념이 필수적입니다. 본 실습에서 흐름 자체를 따라갈 수는 있어도 완전히 이해하기 위해서는 고윳값, 고유벡터, 내적, 직교 등의 개념을 숙지해야합니다. 아래 코드를 진행하다가 이해가지 않는 점이 있다면 따로 공부하는 것을 추천드립니다. 바로 아래 링크에 설명이 아주 잘되어 있습니다. 추천드려요!

주성분 분석(PCA) - 공돌이의 수학정리노트

 

데이터 준비

실습에 사용할 데이터입니다. 머신러닝 분야에서 매우 자주 사용되는 붓꽃(iris) 데이터셋입니다.

UCI Machine Learning Repository: Iris Data Set

 

데이터 불러오기

이번 실습에서는 pandas를 활용하여 데이터를 불러오도록 하겠습니다. 이전 KNN 실습과 동일한 데이터입니다.

import pandas as pd

data = pd.read_csv("iris.data", names=['sepal length', 'sepal width', 'petal length', 'petal width', 'label'])
data

 

Tensor로 변환, 중심 옮기기

불러온 데이터를 torch tensor 형태로 변환해줍니다. data_x는 붓꽃의 꽃받침, 꽃잎의 정보를 tensor 형태로 저장하고, data_y는 붓꽃 품종에 대한 라벨을 저장합니다.

import torch

data_x = torch.from_numpy(data.iloc[:, :4].values)
data_y = data.iloc[:, 4]

 

축을 옮기면 분산을 구하기 쉬워집니다. data_x에서 각 행마다 평균을 빼주어 축을 옮겨줍니다.

# 중심 옮기기
X = data_x - data_x.mean(dim=0)

 

중심을 옮긴 데이터를 시각화합니다.

import matplotlib.pyplot as plt

plt.scatter(X[:,0], X[:,1])
plt.show()

 

주성분 찾기: findPC

본격적으로 주성분을 찾는 함수를 구현해보도록 하겠습니다. 목적은 분산을 최대로 하는 w벡터를 찾는 것입니다. 분산이 클수록 각 데이터의 특성 차이가 크게 반영된다는 것이고, 이는 데이터의 성분을 잘 분리한다는 것이기 때문입니다. 즉, 분산을 최대로 하는 w벡터가 찾고자하는 주성분 벡터가 되는 것입니다.

 

w벡터를 이용해 분산을 구해주고, 매 epoch 마다 분산을 backward 시켜줍니다. 이때 -를 붙힌 이유는 backward를 이용한 학습 방식 자체가 최소가 되는 값을 찾아가는 것이기 때문입니다. 음수가 되면 분산의 절댓값이 커질수록 더 작아지기 때문입니다. 반환해줄 때는 다시 양수로 변환해주면 됩니다.

def findPC(X):
  w = torch.randn(X.shape[1])
  w_size = (w * w).sum() ** 0.5
  w /= w_size

  lr = 0.1

  for epoch in range(100):
    w.requires_grad_(True)
    var = -((X * w).sum(dim=1) ** 2).mean()

    var.backward()

    with torch.no_grad():
      w = w - lr * w.grad
      w /= ((w * w).sum() ** 0.5)
      # print(f'{w}, {var}')

  w.requires_grad_(False)
  return w, -var.item()


w, v = findPC(X)
print(w, v)

tensor([ 0.3616, -0.0823, 0.8566, 0.3588]) 4.1966750403008435

 

주성분 찾기: PCA

k개의 주성분을 구하는 함수입니다. 앞서 구현했던 findPC 함수를 활용합니다. k번 반복하면서 이전에 찾았던 주성분을 제거하고 새로운 주성분을 탐색합니다. 찾은 주성분과 분산은 따로 리스트를 만들어 저장했다가 반환해줍니다.

def PCA(k, X):

  vars = []
  pcs = []

  Z = X - X.mean(dim=0)

  for i in range(k):
    w, v = findPC(Z)
    pcs.append(w)
    vars.append(v)

    # Z = Z - w * Z dot b
    Z = Z - w * (Z * w).sum(dim=1).unsqueeze(1)
  
  return pcs, vars


pcs, vars = PCA(4, X)
print(pcs)
print(vars)

[tensor([ 0.3616, -0.0823, 0.8566, 0.3588]), tensor([ 0.6492, 0.7339, -0.1867, -0.0706]), tensor([-0.6580, 0.6302, 0.3982, 0.1064]), tensor([ 0.1106, -0.4004, 0.1975, -0.8879])]

[4.196674962590516, 0.24058053835753687, 0.062016943866016895, 0.03248544594386099]

 

시각화

그래프를 통해 시각화해봅니다. 각 주성분 별 분산을 확인할 수 있습니다.

plt.plot(vars, "-o")
plt.show()

 

변환된 그래프 그리기

앞서 구했던 주성분들이 축이 되도록 변환해줍니다.

# @ = matmal
# 실습할 때 함수로 k 변경해가며 해보기

X = torch.tensor(X, dtype=torch.float)
XX = X @ torch.stack(pcs[:2]).T

 

구했던 주성분들 중 분산이 높은 순서대로 PC1과 PC2를 축으로 하여 그래프를 그립니다.

species = {a: i for i, a in enumerate(set(data['label']))}

plt.scatter(XX[:,0], XX[:,1], c=[species[x] for x in data['label']])
plt.show()

 

sklearn

PCA 함수를 직접 구현하는 것은 꽤 어렵습니다. 사이킷런 모델을 활용하면 좀 더 쉽게 구현이 가능합니다. x를 표준화하는 것도 함수 하나만 적용하면 됩니다.

from sklearn.preprocessing import StandardScaler

x = data.drop(['label'], axis=1).values # 독립변인들의 value값만 추출 (위의 data_x에 해당)
y = data['label'].values # 종속변인 추출 (위의 data_y에 해당)

x = StandardScaler().fit_transform(x) # x객체에 x를 표준화한 데이터를 저장

features = ['sepal length', 'sepal width', 'petal length', 'petal width']
pd.DataFrame(x, columns=features).head()

 

사이킷런의 PCA 함수를 이용해 주성분 2개를 구합니다.

from sklearn.decomposition import PCA

pca = PCA(n_components=2) # 주성분 개수
principalComponents = pca.fit_transform(x)
principalDf = pd.DataFrame(data=principalComponents, columns = ['pc1', 'pc2'])

principalDf.head()

 

다음과 같은 코드를 통해 주성분이 차지하고 있는 분산 비율을 확인할 수 있습니다. pca.explained_variance_로 분산값도 얻을 수 있습니다.

pca.explained_variance_ratio_

array([0.72770452, 0.23030523])

 

두 주성분의 분산 비율을 합하면 전체의 약 96% 가량임을 확인할 수 있습니다.

# 두 개의 분산이 전체의 약 96% 설명 가능
sum(pca.explained_variance_ratio_)

0.9580097536148197

 

시각화하기 전에 데이터의 라벨을 데이터프레임에 합해줍니다.

finalDf = pd.concat([principalDf, data[['label']]], axis = 1)

 

PC1과 PC2를 축으로 하여 데이터를 시각화합니다.

plt.xlabel('pc1')
plt.ylabel('pc2')
plt.title('PCA')

labels = ['Iris-setosa', 'Iris-versicolor', 'Iris-virginica']
colors = ['r', 'g', 'b']
for label, color in zip(labels, colors):
    indicesToKeep = finalDf['label'] == label
    plt.scatter(finalDf.loc[indicesToKeep, 'pc1']
               , finalDf.loc[indicesToKeep, 'pc2']
               , c = color
               , s = 20)
plt.legend(labels)
plt.grid()