[라이트 머신러닝] Session 13. 비지도 차원축소! PCA!

 이번 세션에 들어가기 전에, 저는 도통 이 책으로 이해가 안되서 따로 강의를 들어서 PCA에 대한 개념을 잡아왔으니 간단하게 먼저 설명을 드리고 들어가도록 하겠습니다. 부디 이 짧은 설명이 여러분이 이해하시는데 도움이 되길 바랍니다:)

 PCA는 입력 데이터의 구조는 최대한 보존하면서 차원을 감소시키는 기법입니다. 이때, 데이터의 분산이 데이터의 구조라고 할 수 있습니다. 이 데이터의 분산을 최대한으로 유지하는 저차원 벡터에 사영을 시키는 방식이죠.

 예를 들어, 2차원에서 1차원으로 차원을 축소시킨다고 생각해봅시다. 우리는 일단 평균을 0으로 맞춘 데이터가 필요하고, 이 데이터를 공분산 행렬로 만들어줍니다. 그리고 이 공분산 행렬의 고유벡터를 기준으로 PCA를 실행하죠.

 이렇게 되면 우리는 또 공분산 행렬은 뭔가, 하고 의문이 생깁니다. 공분산은 각 특성의 변동이 얼마나 닮았는가? 입니다. 그러니까 각 특성에서 전체 특성의 평균을 뺀 값의 닮은 정도를 이야기하는 것이죠. 앞에서 언급했듯이 우리는 먼저 평균을 0으로 만들어주었으니, 닮은 정도만 알면 되겠네요.

 닮은 정도는 행렬 간의 내적으로 알 수 있습니다. 입력 데이터 행렬 X와 X의 전치행렬을 점곱하면 내적 행렬이 생깁니다. 이 내적 행렬은 대칭행렬이고, 행렬 안에 든 데이터들은 거의 공분산과 유사합니다. 왜 유사하다고 표현했냐면, 공분산은 여기서 전체 샘플의 개수로 나눠주어야하기 때문입니다.

 만약에 샘플의 개수가 너무 많아지만 각각의 데이터 값이 커지기 때문에 이 문제를 해결하기 위해 샘플 수로 나눠주는 것입니다.

 자, 이 단계에서 PCA를 수행해준다고 할 때, 가장 데이터들이 넓게 퍼져 있어야 데이터의 원본 구조를 가장 잘 갖고 있다고 판단할 수 있겠죠? 이 넓게 퍼져있는 정도를 분산이라고 합니다. 아래 사진을 참조하면, 분홍색 선으로 표현된 부분이 가장 분산이 큰 주축 벡터입니다. 이 주축에 사영을 해줌으로써 PCA가 완료되는 것이죠.

 그렇다면 분산이 큰 정도는 어떻게 알 수 있을까요? 이것은 고유값과 고유 벡터를 통해서 구할 수 있습니다. 데이터의 고유값을 내림차순으로 정렬하면 처음 두 벡터가 가장 큰 분산을 가짐을 알 수 있습니다(2차원 기준). 가장 큰 고윳값을 가진 고유벡터가 주축이 되는 것이죠. 이 주축에 정사영 해주는 것이 바로 PCA입니다.

 아래 내용은 책 내용을 정리한 내용이므로, 이 앞선 부분을 간단하게 머리에 새겨두고 읽어주시길 바랍니다.

1. 주성분 분석의 주요 단계

비지도 선형 변환 기법인 PCA(Principal Component Analysis)는 고차원 데이터에서 분산이 큰 방향을 찾아서 더 작거나 같은 수의 차원을 갖는 부분공간으로 투영하는 것입니다. 특성 사이의 상관관계를 기반으로 데이터의 특성을 잡아낼 수 있죠. 이때 투영한 새로운 부분공간의 주성분, 직교좌표는 조건 하에서 분산이 최대인 방향입니다. 아래 그림을 볼까요?

x1과 x2가 원본 특성 축이고, PC1과 PC2는 주성분입니다. 두 직교 좌표는 서로 직각을 이루죠. 이는 벡터의 성질에 의해서 두 직교좌표가 상관관계가 없다는 뜻입니다. PCA를 사용해서 차원을 축소하기 위해서는 원본차원 d * 새로운 특성 부분공간 k 차원의 변환행렬 W를 만듭니다. 이 부분 공간은 당연히 원본 d차원보다 작은 차원을 갖습니다.

원본 차원에서 새로운 차원의 부분 공간으로 변환하여 만들어진 첫 번째 주성분은 가장 큰 분산을 갖습니다. 모든 주성분은 다른 주성분과 상관관계가 없다는 제약 하에 가장 큰 분산을 가집니다. 이는 입력특성간의 상관관계와는 관련이 없습니다.

PCA 방향은 데이터의 규모에 민감합니다. 특성의 규모, 즉 스케일이 다르고 모든 특성의 중요도를 동일하게 취급하려면 특성 표준화 처리가 꼭 필요합니다. 이제 PCA 알고리즘을 사용하기 전에 해야할 목록을 정리해볼까요?

1. d 차원 데이터셋을 표준화합니다.

2. 공분산 행렬(covariance matrix)을 만듭니다.

3. 공분산 행렬을 고유 벡터(eigenvector)와 교윳값(eigenvalue)로 분해합니다.

4. 고윳값을 내림차순 정렬하고 그에 해당하는 벡터의 순위를 매깁니다.

5. 고윳값이 가장 큰 k개의 고유 벡터를 선택합니다. 여기서 k는 특성 부분 공간의 차원입니다.

6. 최상위 k개의 고유백터로 투영 행렬(projection matrix) W를 만듭니다.

7. 투영행렬 W를 이용해 d차원 입력 데이터셋 X를 새로운 k 차원의 특성 부분공간으로 변환합니다.

 

2. 주성분 추출 단계

이번에는 PCA의 처리과정 중 처음 4단계를 구현해보도록 하겠습니다. 물론 파이썬으로요! 먼저, 앞선 세션에서 사용한 wine 데이터를 로드해주고, 훈련세트와 테스트세트로 나누고 표준화를 적용해줍니다. 훈련세트가 70%, 테스트세트가 30% 입니다.

import pandas as pd
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/''machine-learning-databases/wine/wine.data', header = None)

from sklearn.model_selection import train_test_split
X, y = df.wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, Xtest, y_train, y_test = \
    train_test_split(X, y, test_size = 0.3, stratify = y, random_state = 0)

from sklearn.preprocessing import StandartScaler
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)

위 코드로 전처리 단계를 완료하고, 공분산 행렬을 만들어보도록 하겠습니다. 공분산 행렬은 주어진 벡터의 각 요소 쌍 사이의 공분산을 제공하는 정사각 행렬입니다. 공분산은 2개의 확률변수의 상관정도를 나타내는 값입니다. 우리가 만들 공분산 행렬은 d x d차원의 대칭 행렬로 특성간의 공분산을 저장합니다. 전체 샘플에 대한 두 특성 Xj와 Xk 사이의 공분산은 아래와 같이 계산 가능합니다.

여기서 각 특성에서 빼주는 값은 특성 j와 k의 샘플 평균입니다. 여기서는 데이터셋을 표준화 전처리했기 때문에 평균은 0입니다. 정규분포와 같죠. 공분산의 값이 양이면 특성이 함께 증가하거나 감소한다는 뜻이고, 음이면 반대 방향으로 달라진다는 뜻입니다. 세 개의 특성으로 이루어진 공분산 행렬은 아래와 같습니다. (아래 그림에서 시그마는 합의 기호가 아닙니다!)

공분산 행렬의 고유 벡터가 주성분을 표현합니다. 주성분은 최대 분산의 방향입니다. 여기에 대응되는 고윳값은 주성분의 크기입니다. Wine 데이터의 경우에는 13 * 13의 공분산 행렬에서 13개 벡터와 고윳값을 얻을 수 있습니다.

다음으로 고유 벡터와 고윳값 쌍을 구해보도록 하겠습니다. 고유벡터 v는 아래 식을 만족합니다.

위 식의 람다는 스케일을 담당하는 고윳값입니다. 고유 벡터와 고윳값을 계산하기는 귀찮기 때문에 넘파이 linalg.eig 함수를 사용합니다. 코드는 아래와 같습니다.

import numpy as np
cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print('\nEigenvalues \n%s' % eigen_vals)

numpy.cov 함수를 사용해서 훈련 데이터의 공분산 행렬을 계산해주고, linalg.eig 함수로 고윳값을 계산합니다. 이를 통해서 13개 고윳값이 들어있는 벡터인 eigen_vals와 대응하는 고유벡터가 저장된 행렬을 얻을 수 있습니다.

3. 총분산과 설명된 분산

우리의 목표는 데이터셋 차원을 새로운 특성 부분 공간으로 압축해야하기 때문에 가장 많은 정보를 가진 고유 벡터의 일부만 선택합니다. 이때 선택하는 고유벡터는 고윳값 순서에 따라 최상위 k개 입니다. 일단, 소제목에서 알 수 있는 설명된 분산이 무엇인지부터 알아볼까요?

설명된 분산(explained variance)은 통계에서 주어진 데이터의 분산을 설명하는 비율을 측정합니다. 그 중에서, 설명된 분산의 비율은 전체 고윳값 중에서 원하는 고윳값의 비율입니다.

람다 j가 여기에서 원하는 고윳값이 됩니다. 식으로 나타내면 위와 같겠죠? 이제 넘파이 cumsum 함수로 설명된 분산 누적 합을 계산하고, matplot lib의 step함수로 그래프를 그려보도록 하겠습니다.

tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)

...

import matplotlib.pyplot as plt


plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',
        label='individual explained variance')
plt.step(range(1, 14), cum_var_exp, where='mid',
         label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal component index')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

이 코드를 사용하면 결과 그래프에서 설명된 분산의 비율을 확인할 수 있습니다.

그려진 그래프에서 첫 번째 고윳값이 전체 분산의 40%를 담당하고, 두 번째 고윳값까지가 60%를 차지한다는 것을 알 수 있습니다.

아마도 여러분은 뭔가 랜덤 포레스트의 특성 중요도랑 비슷하지 않은가 하는 생각이 드실 겁니다. 하지만 우리가 지금 이야기하고 있는 이 PCA라는 친구는 비지도 학습입니다! 또한 랜덤 포레스트에서는 노드의 불순도를 계산하지만 분산은 값들이 퍼진 정도를 측정합니다.

4. 특성 변환

이제 우리가 해야할 일은 세 가지입니다.

1. 고윳값의 내림차순으로 고유벡터를 정렬하기.
2. 선택된 고유벡터로 투영 행렬 구성하기.
3. 투영행렬을 사용해 데이터를 저차원 부분 공간으로 변환하기.

먼저, 고유벡터와 고윳값 쌍을 내림차순으로 정렬합니다.

# (고윳값, 고유벡터) 튜플의 리스트를 만듭니다
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
               for i in range(len(eigen_vals))]

# 높은 값에서 낮은 값으로 (고윳값, 고유벡터) 튜플을 정렬합니다
eigen_pairs.sort(key=lambda k: k[0], reverse=True)

그리고 가장 큰 두 개의 고윳값과 고유벡터 세트를 선택합니다. 앞선 그래프에서 확인하신 것과 같이 이 두 세트로 총 분산의 60% 를 잡아낼 수 있습니다. 여기 예제에서는 2차원 산점도를 위해 2 개를 선택했지만, 실전에서는 계산 효율과 모델 성능 사이의 중간점을 찾아서 개수를 선택해주면 됩니다.

w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis]))
print('투영 행렬 W:\n', w)

위 코드를 실행하면 13 * 2의 투영 행렬 W를 만들 수 있습니다. 이 투영행렬을 사용하면 샘플 x(1 * 13의 행 백터)를 PCA 부분공간으로 투영해 x'를 얻을 수 있습니다. 2개 특성으로 이루어진 2차원 샘플 벡터가 되죠.

식은 요렇게 됩니다. 비슷하게 전체 124 * 13차원의 훈련 데이터를 행렬 곱으로 두 개의 주성분으로 변환할 수 있습니다. 위 식을 대문자로 바꿔서 표현할 수 있습니다.

X_train_std[0].dot(w)

X_train_pca = X_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train == l, 0], 
                X_train_pca[y_train == l, 1], 
                c=c, label=l, marker=m)

plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

위 코드를 이용하면 124 * 2차원의 행렬로 변환한 Wine 훈련 세트를 2차원 산점도로 시각화합니다.

이 그래프에서 볼 수 있듯이 데이터가 y축(두 번째 고윳값)보다 x축(첫 번째 고윳값)을 따라서 더 넓게 펼쳐져 있습니다.

5. 사이킷런의 주성분 분석

이전 단계들에서 PCA의 원리를 이해해 볼 수 있었다면 이젠 사이킷런을 이용해볼 차례입니다. 일단 먼저 결정 경계를 그리는 함수는 초반 세션에서 만들었던 plot_decision_regions를 이용하겠습니다.

def plot_decision_regions(X, y, classifier, resolution=0.02):

    # 마커와 컬러맵을 준비합니다
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # 결정 경계를 그립니다
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # 클래스 샘플을 표시합니다
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], 
                    y=X[y == cl, 1],
                    alpha=0.6, 
                    c=cmap.colors[idx],
                    edgecolor='black',
                    marker=markers[idx], 
                    label=cl)

그리고 훈련 데이터를 사용해 모델을 훈련하고 같은 모델 파라미터로 훈련과 테스트 데이터를 변환합니다. 사이킷런 PCA 클래스를 Wine데이터셋의 훈련 세트에 적용하고 로지스틱 회귀로 분류해보겠습니다.

from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)

lr = LogisticRegression(solver='liblinear', multi_class='auto')
lr = lr.fit(X_train_pca, y_train)

#그리기!
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

코드를 실행하면 두 개의 주성분 축으로 줄어든 결정 경계를 볼 수 있습니다.

사이킷런의 PCA 투영과 우리가 앞에서 쭉 구현해본 PCA는 서로 거울처럼 뒤집혀있는 것을 확인할 수 있습니다. 틀린 것은 아니고, 그저 계산 방법에 따라 고유벡터가 양수나 음수 값을 가질 수 있기 때문입니다. 필요하다면 데이터에 -1을 곱해줄 수 있습니다.

고유벡터는 일반적으로 단위 길이가 1이 되도록 정규화되어있습니다. 이제 테스트 세트를 변환하고 결정경계를 그려볼까요?

plot_decision_regions(X_test_pca, y_test, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

만약 전체 주성분의 설명된 분산의 비율을 알고 싶다면 PCA 객체를 생성할 때 n_components를 None으로 해두고 explained_variance_ratio_ 속성에서 모든 주성분의 설명된 분산 비율을 알 수 있습니다.

pca = PCA(n_components=None)
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_

---

자, 여기까지가 PCA에 대한 설명이었습니다. 다음 세션에서는 LDA를 들고 돌아올텐데요, 위에 수많은 수식에 대해서 너무 신경쓰실 필요는 없습니다. PCA가 어떤 거고, 어떻게 수행되는지에 대해서만 인지하면 됩니다. 실제로 아주 잦게 쓰이는 기법도 아니고, 관심있는 분야가 아니면 아는 사람도 적은 부분입니다. 그럼 다음 세션에서 뵙겠습니다!