머신 러닝 교과서/ 차원 축소를 사용한 데이터 압축

주성분 분석을 통한 비지도 차원 축소

  • 특성 선택과 마찬가지로 여러 특성 추출 기법을 사용하여 데이터셋의 특성 개수를 줄일 수 있다.
    • 특성 선택과 특성 추출의 차이는 원본 특성을 유지하느냐에 있다.
    • 순차 후진 선택 같은 특성 선택 알고리즘을 사용할 때는 원본 특성을 유지하지만, 특성 추출은 새로운 특성 공간으로 데이터를 변환하거나 투영한다.
    • 차원 축소 관점에서 보면 특성 추출은 대부분의 관련 있는 정보를 유지하면서 데이터를 압축하는 방법으로 이해할 수 있다.
  • 특성 추출이 저장 공간을 절약하거나 학습 알고리즘의 계산 효율성을 향상할 뿐만 아니라 차원의 저주 문제를 감소시켜 예측 성능을 향상하기도 한다. 특히 규제가 없는 모델로 작업할 때 그렇다.

주성분 분석의 주요 단계

  • PCA는 비지도 선형 변환 기법으로 주로 특성 추출과 차원 축소 용도로 많은 분야에서 널리 사용된다.
  • PCA는 특성 사이의 상관관계를 기반으로 하여 데이터에 있는 특성을 잡아낼 수 있다.
    • 요약하자면 PCA는 고차원 데이터에서 분산이 가장 큰 방향을 찾고 좀 더 작거나 같은 수의 차원을 갖는 새로운 부분 공간으로 이를 투영한다.
    • 새로운 부분 공간의 직교 좌표는 주어진 조건 하에서 분산이 최대인 방향으로 해석할 수 있다.
    • 새로운 특성 축은 아래 그림과 같이 서로 직각을 이룬다. 아래 그림에서 x_{1}, x_{2} 는 원본 특성 축이고 PC1, PC2 는 주성분이다.

  • PCA를 사용하여 차원을 축소하기 위해 d \times k 차원의 변환행렬 W 를 만든다.
    • 이 행렬로 샘플 벡터 x 를 새로운 k 차원의 특성 부분 공간으로 매핑한다.
    • 이 부분 공간은 원본 d 차원의 특성 공간보다 작은 차원을 가진다.

x = [x_{1}, x_{2}, ... , x_{d}], x \in \mathbb{R}^{d}

\downarrow xW, W \in \mathbb{R}^{d \times k}

z = [z_{1}, z_{2}, ... , z_{k}], z \in \mathbb{R}^{k}

  • 원본 d 차원 데이터를 새로운 k 차원의 부분 공간(일반적으로 k < d )으로 변환하여 만들어진 첫 번째 주성분이 가장 큰 분산을 가질 것이다.
    • 모든 주성분은 다른 주성분들과 상관관계가 있더라도 만들어진 주성분은 서로 직각을 이룰 것이다.
    • PCA 방향은 데이터 스케일에 매우 민감하다.
    • 특성의 스케일이 다르고 모든 특성의 중요도를 동일하게 취급하려면 PCA를 적용하기 전에 특성을 표준화 전처리해야 한다.
  • 차원 축소 PCA 알고리즘의 단계는 다음과 같다.
    1. d 차원 데이터셋을 표준화 전처리한다.
    2. 공분산 행렬(covariance matrix)을 만든다.
    3. 공분산 행렬을 고유 벡터(eigenvector)와 고윳값(eigenvalue)으로 분해한다.
    4. 고윳값을 내림차순으로 정렬하고 그에 해당하는 고유 벡터의 순위를 매긴다.
    5. 고윳값이 가장 큰 k 개의 고유 벡터를 선택한다. 여기서 k 는 새로운 특성 부분 공간의 차원이다. (k \leq d )
    6. 최상위 k 개의 고유벡터로 투영행렬(projection matrix) W 를 만든다.
    7. 투영 행렬 W 를 사용해서 d 차원 입력 데이터셋 X 를 새로운 k 차원의 특성 부분 공간으로 변환한다.

주성분 추출단계

  • 우선 PCA 처음 4단계를 처리한다.
    1. 데이터를 표준화 전처리 한다.
    2. 공분산 행렬을 구성한다.
    3. 공분산 행렬의 고윳값과 고유벡터를 구한다.
    4. 고윳값을 내림차순으로 정렬하여 고유 벡터의 순위를 매긴다.
  • Wine 데이터셋을 가져와서 70:30 비율로 훈련 세트와 테스트 세트를 나누고 표준화를 적용하여 단위 분산을 갖도록 한다.
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

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

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:,0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y, random_state=0)

sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)
  • 그 다음으로 공분산 행렬을 만드는 단계를 진행한다. 공분산 행렬은 d \times d 차원의 대칭 행렬로 특성 상호 간의 공분산을 저장한다. d 는 데이터셋에 있는 차원 개수이다.
    • 예컨대 전체 샘플에 대한 두 특성 x_{j} x_{k} 사이의 공분산은 다음 식으로 계산할 수 있다.

\sigma_{jk} = {1 \over n} \sum_{i = 1}^{n} (x_{j}^{(i)} - \mu_{j})(x_{k}^{(i)} - \mu_{k})

  • 여기서 \mu_{j} \mu_{k} 는 특성 j k 의 샘플 평균이다.
  • 데이터셋을 표준화 전처리했기 때문에 샘플 평균은 0이다.
  • 두 특성 간 양의 공분산은 특성이 함께 증가하거나 감소하는 것을 나타낸다. 반면 음의 공분산은 특성이 반대 방향으로 달라진다는 것을 나타낸다.
    • 예컨대 세 개의 특성으로 이루어진 공분산 행렬은 다음과 같이 쓸 수 있다. (여기서 \Sigma 는 대문자 표기일 뿐 합을 뜻하는 것이 아니므로 주의)

\Sigma = \left[ \begin{array}{rrr} \sigma_{1}^{2} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{2}^{2} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{3}^{2} \end{array} \right]

  • 공분산 행렬의 고유 벡터가 주성분(최대 분산의 방향)을 표현한다. 이에 대응되는 고윳값은 주성분의 크기이다.
    • Wine 데이터셋의 경우 13 x 13 차원의 공분산 행렬로부터 13개의 고유벡터와 고윳값을 얻을 수 있다.
  • 다음 단계로 공분산 행렬의 고유벡터와 고윳값의 쌍을 구해 보자.
    • 선형대수학에 따르면 고유벡터 v 는 다음 식을 만족한다.

\Sigma v = \lambda v

  • 여기서 \lambda 는 스케일을 담당하는 고윳값이다. 고윳벡터와 고윳값을 직접 계산하는 것은 복잡한 작업이기 때문에 Numpy의 linalg.eig 함수를 이용하여 계산한다.
import numpy as np

cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
  • Note) numpy.linalg.eig 함수는 대칭과 비대칭 정방 행렬을 모두 다룰 수 있지만 이따금 복소수 고윳값을 반환한다. 이와 비슷하게 에르미트 행렬을 분해하기 위해 구현된 numpy.linalg.eigh 함수는 공분산 행렬과 같은 대칭 행렬을 다룰 때 수치적으로 더 안정된 결과를 만든다. (numpy.linalg.eigh는 항상 실수 고윳값을 반환한다.)
  • Note) 사이킷런의 PCA 클래스는 직접 고윳값과 고유 벡터를 계산하는 대신 특이값 분해 (singular value decomposition) 방식을 이용하여 주성분을 구한다.

총분산과 설명된 분산

  • 데이터셋 차원을 새로운 특성 부분 공간으로 압축하여 줄여야 하기에 가장 많은 정보(분산)를 가진 고유 벡터(주성분) 일부만 선택한다.
    • 고윳값은 고유 벡터의 크기를 결정하므로 고윳값을 내림차순으로 정렬한다.
    • 고윳값 순서에 따라 최상위 k 개의 고유벡터를 선택한다.
  • 가장 정보가 많은 k 개의 고유 벡터를 선택하기 전에 고윳값의 설명된 분산(explained variance) 그래프로 그려보겠다.
    • 고윳값 \lambda_{j} 의 설명된 분산 비율은 전체 고윳값의 합에서 고윳값 \lambda_{j} 의 비율이다.

{\lambda_{j} \over \sum_{d}^{j=1} \lambda_{j}}

  • numpy의 cumsum 함수로 설명된 분산의 누적 합을 계산하고 matplotlib의 step 함수로 그래프를 그리면 다음과 같다.
import matplotlib.pyplot as plt

tot = sum(eigen_vals)

var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]

cum_var_exp = np.cumsum(var_exp)

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 variacne ratio')
plt.xlabel('Principal component index')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

  • 랜덤 포레스트는 클래스 소속 정보를 사용하여 노드의 불순도를 계산하는 방면, 분산은 특성 축을 따라 값들이 퍼진 정도를 측정한다.

특성 변환

  • 공분산 행렬을 고유 벡터와 고윳값 쌍으로 성공적으로 분해한 후 Wine 데이터 셋을 새로운 주성분 축으로 변환하는 3개의 단계는 다음과 같다.
    • 고윳값이 가장 큰 k 개의 고유 벡터를 선택한다. 여기서 k 는 새로운 특성 부분 공간의 차원이다 (k \leq d )
    • 최상위 k 개의 고유벡터로 투영 행렬 W 를 만든다.
    • 투영행렬 W 를 사용해서 d 차원 입력 데이터셋 X 를 새로운 k 차원의 특성 부분 공간으로 변환한다.
  • 좀 더 쉽게 설명하면 고윳값의 내림차순으로 고유 벡터를 정렬하고 선택된 고유 벡터로 투영 행렬을 구성한다. 이 투영 행렬을 사용하여 데이터를 저차원 부분 공간으로 변환한다.
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)
  • 앞의 코드를 이용하여 최상위 두 개의 고유 벡터로부터 13 x  2 차원의 투영행렬 W 를 만든다.
  • 투영 행렬을 사용하면 샘플 x (1 x 13 차원의 행 벡터)를 PCA 부분 공간(두 개의 주성분)을 투영하여 x' 를 얻을 수 있다. 두 개의 특성으로 구성된 2차원 샘플 벡터이다.

x' = xW

w = np.hstack((eigen_pairs[0][1][:, np.newaxis], eigen_pairs[1][1][:, np.newaxis]))
  • 비슷하게 124 x 13 차원의 훈련 데이터셋을 행렬 내적으로 두 개의 주성분으로 변환할 수 있다.

X' = XW

X_train_pca = X_train_std.dot(w)
  • 124 x 2 차원 행렬로 변환된 Wine 훈련 세트를 2차원 산점도로 시각화 하면 아래와 같다.
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.ylabel('PC 1')
plt.xlabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

  • 위 이미지에서 볼 수 있듯이 데이터가 y축 보다 x 축을 따라 더 넓게 퍼져 있다. 이는 이전에 만든 설명된 분산의 그래프와 동일한 결과이다.

사이킷런의 주성분 분석

  • 다음의 코드를 실행하면 두 개의 주성분 축으로 줄어든 훈련 데이터로 만든 결정 경계를 볼 수 있다.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

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.3, 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)

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.fit(X_train_pca, y_train)

plot_decision_regions(X_train_pca, y_train, classifier=lr)

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

  • 테스트 세트에 적용한 결과는 아래와 같다.

선형 판별 분석을 통한 지도 방식의 데이터 압축

  • 선형 판별 분석(Linear Discriminant Analysis, LDA)은 규제가 없는 모델에서 차원의 저주로 인한 과대적합 정도를 줄이고 계산 효율성을 높이기 위한 특성 추추르이 기법으로 사용할 수 있다.
  • LDA 이면에 있는 일반적인 개념은 PCA와 매우 비슷하다. PCA가 데이터셋에 있는 분산이 최대인 직교 성분 축을 찾으려고 하는 반면, LDA 목표는 클래스를 최적으로 구분할 수 있는 특성 부분 공간을 찾는 것이다.

주성분 분석 vs 선형 판별 분석

  • PCA와 LDA 모두 데이터셋의 차원 개수를 줄일 수 있는 선형 변환 기법이다. 전자는 비지도 학습 알고리즘이지만, 후자는 지도 학습 알고리즘이다.
    • 직관적으로 LDA가 PCA 보다 분류 작업에서 더 뛰어난 특성 추출 기법이라고 생각할 수 있다.
    • 마르티네스는 PCA를 통한 전처리가 특정 이미지 인식 작업에 더 뛰어난 분류 결과를 내는 경향이 있다고 보고 했다. 예컨대 각 클래스에 속한 샘플이 몇 개 되지 않았을 때이다.
  • 다음 그림은 이진 분류 문제를 위한 LDA 개념을 요약하여 나타낸다. 클래스1의 샘플은 동그라미고, 클래스 2의 샘플은 덧셈 기호이다.

  • x축(LD 1)으로 투영하는 선형 판별 벡터는 두 개의 정규 분포 클래스를 잘 구분한다. y 축(LD 2)으로 투영하는 선형 판별 벡터는 데이터셋에 있는 분산을 많이 잡아내지만, 클래스 판별 정보가 없기 떄문에 좋은 선형 판별 벡터가 되지 못한다.
  • LDA는 데이터가 정규분포라고 가정한다. 또 클래스가 동일한 공분산 행렬을 가지고 샘플은 서로 통계적으로 독립적이라고 가정한다.
    • 하나 이상의 가정이 조금 위반되더라도 여전히 LDA는 차원 축소를 상당히 잘 수행한다.

선형 판별 분석의 내부 동작 방식

  • LDA 수행에 필요한 단계는 다음과 같다.
    1. d 차원의 데이터셋을 표준화 전처리한다 (d 는 특성 개수)
    2. 각 클래스에 대해 d 차원의 평균 벡터를 계산한다.
    3. 클래스 간의 산포 행렬(scatter matrix) S_{B} 와 클래스 내 산포 행렬 S_{W} 를 구성한다.
    4. S_{W}^{-1} S_{B} 행렬의 고유 벡터와 고윳값을 계산한다.
    5. 고윳값을 내림차순으로 정렬하여 고유 벡터의 순서를 매긴다.
    6. 고윳값이 가장 큰 k 개의 고유 벡터를 선택하여 d \times k 차원의 변환 행렬 W 를 구성한다. 이 행렬의 열이 고유벡터이다.
    7. 변환 행렬 W 를 사용하여 새로운 특성 부분 공간으로 투영한다.
  • 여기서 볼 수 있듯 LDA는 행렬을 고윳값과 고유 벡터로 분해하여 새로운 저차원 특성 공간을 구성한다는 점에서 PCA와 매우 닮았다.

산포 행렬 계산

  • PCA에서 했기 때문에 데이터셋의 특성을 표준화하는 단계는 건너뛰고 바로 평균 벡터 계산을 진행한다.
    • 평균 벡터를 사용하여 클래스 간의 산포 행렬과 클래스 내 산포 행렬을 구성한다.
    • 평균 벡터 m_{i} 는 클래스 i 의 샘플에 대한 특성의 평균값 \mu_{m} 을 저장한다.

m_{i} = {1 \ over n_{i}} \sum_{x \in D_{i}}^{c} x_{m}

  • 3개의 평균 벡터가 만들어진다.

m_{i} = \left[ \begin{array}{rrrr} \mu_{i, alcohol} \\ \mu_{i, malic acid} \\ ... \\ \mu_{i, proline} \end{array} \right] i \in \{1 , 2, 3 \}

np.set_printoptions(precision=4)
mean_vecs = []

for label in range(1, 4):
   mean_vecs.append(np.mean(X_train_std[y_train==label], axis=0))
  • 평균 벡터를 사용하여 클래스 내 산포 팽렬 S_{W} 를 계산할 수 있다.

S_{W} = \sum_{i=1}^{c} S_{i}

  • 이 행렬은 개별 클래스 i 의 산포 행렬 S_{i} 를 더하여 구한다.

S_{i} = \sum_{x \in D_{i}}^{c} (x - m_{i})(x - m_{i})^{T}

d = 13
S_W = np.zeros((d,d))

for label, mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.zeros((d, d))

   for row in X_train_std[y_train == label]:
       row, mv = row.reshape(d, 1), mv.reshape(d, 1)
       class_scatter += (row - mv).dot((row - mv).T)   

    S_W += class_scatter
  • 개별 산포 행렬 S_{i} 를 산포 행렬 S_{W} 로 모두 더하기 전에 스케일 조정을 해야 한다.
    • 산포 행렬을 클래스 샘플 개수 \eta_{i} 로 나누면 사실 산포 행렬을 계산하는 것이 공분산 행렬 \Sigma_{i} 를 계산하는 것과 같아진다.
    • 즉 공분산 행렬은 산포 행렬의 정규화 버전이다.

\Sigma_{i} = {1 \over n_{i}} S_{i} = {1 \over n_{i}} \sum_{x \in D_{i}}^{c} (x - m_{i})(x - m_{i})^{T}

d = 13
S_W = np.zeros((d,d))

for label, mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.cov(X_train_std[y_train==label].T, bias=True)
    S_W += class_scatter
  • 클래스 내 산포 행렬(또는 공분산 행렬)을 계산한 후 다음 단계로 넘어가 클래스 간의 산포 행렬 S_{B} 를 계산한다.

S_{B} = \sum_{i = 1}^{c} n_{i} (m_{i} - m)(m_{i} - m)^{T}

  • 여기서 m 은 모든 클래스의 샘플을 포함하여 계산된 전체 평균이다.
mean_overall = np.mean(X_train_std, axis=0)
mean_overall = mean_overall.reshape(d, 1)
d = 13
S_B = np.zeros((d, d))

for i, mean_vec in enumerate(mean_vecs):
   n = X_train[y_train == i +1, :].shape[0]
    mean_vec = mean_vec.reshape(d, 1)
    S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)

새로운 특성 부분 공간을 위해 선형 판별 벡터 선택

  • LDA의 남은 단계는 PCA와 유사하다. 공분산 행렬에 대한 고윳값 분해를 수행하는 대신 행렬 S_{W}^{-1} S_{B} 의 고윳값을 계산하면 된다.
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
  • 고유 벡터와 고윳값 쌍을 계산한 후 내림차순으로 고윳값을 정렬한다.
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
  • LDA에서 선형 판별 벡터는 최대 c-1 개이다. c 는 클래스 레이블의 개수이다. 클래스 내 산포 행렬 S_{B} 가 랭크(rank) 1 또는 그 이하인 c 개의 행렬을 합한 것이기 때문이다.
    • 0이 아닌 고윳값이 두 개만 있는 것을 볼 수 있다.
  • 선형 판별 벡터(고유 벡터)로 잡은 클래스 판별 정보가 얼마나 많은지 측정하기 위해 선형 판별 벡터를 그려보면 아래 그림과 같다.
tot = sum(eigen_vals.real)
discr = [(i/tot) for i in sorted(eigen_vals.real, reverse=True)]
cum_discr = np.cumsum(discr)

plt.bar(range(1, 14), discr, alpha=0.5, align='center', label='individual discriminability')
plt.step(range(1, 14), cum_discr, where='mid', label='cumulative "discriminability"')
plt.ylabel('"discriminability" ratio')
plt.xlabel('Linear Discriminants')
plt.ylim([-0.1, 1.1])
plt.legend(loc='best')
plt.tight_layout()
plt.show()

  • 두 개의 판별 고유 벡터를 열로 쌓아서 변환 행렬 W 를 만든다.
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real, eigen_pairs[1][1][:, np.newaxis].real))

새로운 특성 공간으로 샘플 투영

  • 이전 절에서 만든 변환 행렬 W 를 훈련 세트에 곱해서 데이터를 변환할 수 있다.

X' = XW

X_train_lda = 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_lda[y_train==l, 0], X_train_lda[y_train==l, 1] * -1, c=c, label=l, marker=m)

plt.ylabel('LD 1')
plt.xlabel('LD 2')
plt.legend(loc='lower right')
plt.tight_layout()
plt.show()

사이킷런의 LDA

  • 사이킷런에 구현된 LDA 클래스를 살펴보겠다.
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

lda = LDA(n_components=2)
X_train_lda = lda.fit_transform(X_train_std, y_train)

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

plot_decision_regions(X_train_lda, y_train, classifier=lr)

plt.ylabel('LD 1')
plt.xlabel('LD 2')
plt.legend(loc='lower right')
plt.tight_layout()
plt.show()

  • 테스트 결과는 다음과 같다.

커널 PCA를 사용하여 비선형 매핑

  • 많은 머신 러닝 알고리즘은 입력 데이터가 선형적으로 구분 가능하다는 가정을 한다. 그래서 지금까지 다루었던 알고리즘들은 선형적으로 완벽하게 분리되지 못한 이유를 잡음 때문이라고 가정한다.
  • 그러나 실제 애플리케이션에서는 비선형 문제를 더 자주 맞닥뜨리게 될 것이다.
    • 이런 비선형 문제를 다룰 때 PCA와 LDA 같은 차원 축소를 위한 선형 변환 기법은 최선의 선택이 아니다.
    • 이절에서는 PCA의 커널화 버전 또는 KPCA를 다루겠다.
    • 커널 PCA를 사용하여 선형적으로 구분되지 않는 데이터를 선형 분류기에 적합한 새로운 저차원 부분 공간으로 변환하는 방법을 살펴보겠다.

커널 함수와 커널 트릭

  • 커널 SVM에 관해 배운 것을 떠올려 보면 비선형 문제를 해결하기 위해 클래스가 선형으로 구분되는 새로운 고차원 특성 공간으로 투영할 수 있었다.
    • k 고차원 부분 공간에 있는 샘플 x \in \mathbb{R}^{d} 를 변환하기 위해 비선형 매핑 함수 \phi 를 정의한다.

\phi : \mathbb{R}^{d} \to \mathbb{R}^{k} (k >> d)

  • \phi 함수를 d 차원의 원본 데이터셋에서 더 큰 k 차원의 특성 공간으로 매핑하기 위해 원본 특성의 비선형 조합을 만드는 함수로 생각할 수 있다.
    • 예컨대 2차원 (d=2) 특성 벡터 x \in \mathbb{R}^{d} 가 있으면 (x d 개의 특성으로 구성된 열 벡터) 매핑 가능한 3D 공간은 다음과 같다.

x = [x_{1}, x_{2}]^{T}

\downarrow \phi

z = [x_{1}^{2}, \sqrt{2x_{1}x_{2}}, x_{2}^{2}]^{T}

  • 다른 말로 하면 커널 PCA를 통한 비선형 매핑을 수행하여 데이터를 고차원 공간으로 변환한다.
    • 그다음 고차원 공간에 표준 PCA를 사용하여 샘플이 선형 분류기로 구분될 수 있는 저차원 공간으로 데이터를 투영한다 (샘플이 이 입력 공간에서 잘 구분될 수 있다고 가정).
    • 이 방식의 단점은 계산 비용이 매우 비싸다는 것이다. 여기에 커널 트릭(kernel trick)이 등장한다.
    • 커널 트릭을 사용하면 원본 특성 공간에서 두 고차원 특성 벡터의 유사도를 계산할 수 있다.
  • 커널 트릭에 대해 알아보기 전에 표준 PCA 방식을 다시 살펴보자
    • 두 개의 특성 k, j 사이의 공분산은 다음과 같이 계산한다.

\sigma_{jk} = {1 \over n} \sum_{i = 1}^{n} (x_{j}^{(i)} - \mu_{j})(x_{k}^{(i)} - \mu_{k})

  • \mu_{j} = 0, \mu_{k} = 0 처럼 특성 평균을 0에 맞추었으므로 이 식은 다음과 같이 간단히 쓸 수 있다.

\sigma_{jk} = {1 \over n} \sum_{i = 1}^{n} x_{j}^{(i)} x_{k}^{(i)}

  • 이 식은 두 특성 간의 공분산을 의미한다. 공분산 행렬 \Sigma 를 계산하는 일반식으로 써보자.

\Sigma = {1 \over n} \sum_{i = 1}^{n} x^{(i)} x^{(i)^{T}}

  • 베른하르트 슐코프(Bernhard Scholkopf)는 이 방식을 일반화 하여 \phi 를 통한 비선형 특성 조합으로 원본 특성 공간의 샘플 사이의 내적을 대체했다.

\Sigma = {1 \over n} \sum_{i = 1}^{n} \phi (x^{(i)}) \phi(x^{(i)})^{T}

  • 이 공분산 행렬에서 고유 벡터(주성분)를 얻기 위해서는 다음 식을 풀어야 한다.

\Sigma v = \lambda v

\Rightarrow {1 \over n} \sum_{i = 1}^{n} \phi (x^{(i)}) \phi(x^{(i)})^{T} v = \lambda v

\Rightarrow v = {1 \over n \lambda} \sum_{i = 1}^{n} \phi (x^{(i)}) \phi(x^{(i)})^{T} = \sum_{i=1}^{n} a^{(i)} \phi(x^{(i)})

  • 여기서 \lambda v 는 공분산 행렬 \Sigma 의 고윳값과 고유 벡터이다. a 는 커널(유사도) 행렬 K 의 고유 벡터를 추출함으로써 구할 수 있다.
  • 커널 SVM을 사용하여 비선형 문제 풀기를 떠올리면 커널 트릭을 사용하여 샘플 x 끼리의 \phi 함수 내적을 커널 함수 K 로 바꿈으로써 고유 벡터를 명싲거으로 계산할 필요가 없었다.

\kappa (x^{(i)}, x^{(j)}) = \phi (x^{(i)})^{T} \phi (x^{(j)})

  • 다른 말로 하면 커널 PCA로 얻은 것은 표준 PCA 방식에서처럼 투영 행렬을 구성한 것이 아니고 각각의 성분에 이미 투영된 샘플이다.
    • 기본적으로 커널 함수는 두 벡터 사이의 내적을 계산할 수 있는 함수이다. 즉, 유사도를 측정할 수 있는 함수이다.
  • 가장 널리 사용되는 커널은 다음과 같다.
  • 다항 커널
    • \kappa (x^{(i)}, x^{(j)}) = (x^{(i)T} x^{(j)} + \theta)^{P}
    • 여기서 \theta 는 임계 값이고 P 는 사용자가 지정한 거듭제곱이다.
  • 하이퍼볼릭 탄젠트(hyperbolic tangent) (시그모이드(sigmoid)) 커널
    • \kappa (x^{(i)}, x^{(j)}) = tanh (\eta x^{(i)T} x^{(j)} + \theta)
  • 방사 기저 함수(Radial Basis Function, RBF) 또는 가우시안 커널
    • \kappa (x^{(i)}, x^{(j)}) = \exp(- {\|x^{(i)} - x^{(j)}\|^{2} \over 2 \sigma^{2}})
    • 변수 \gamma = {1 \over 2 \sigma^{2}} 을 도입하여 종종 다음과 같이 쓴다.
    • \kappa (x^{(i)}, x^{(j)}) = \exp(- \gamma \|x^{(i)} - x^{(j)}\|^{2})
  • 지금까지 배운 것을 요약하면 RBF 커널 PCA를 구현하기 위해 다음 3 단계를 정의할 수 있다.
    1. 커널 (유사도) 행렬 K 를 다음 식으로 계산한다.
      • \kappa (x^{(i)}, x^{(j)}) = \exp(- \gamma \|x^{(i)} - x^{(j)}\|^{2})
      • 샘플의 모든 쌍에 대해 구한다.
        • K = \left[ \begin{array}{rrrr} \kappa (x^{(1)}, x^{(1)}) & \kappa (x^{(1)}, x^{(2)}) & ... & \kappa (x^{(1)}, x^{(n)}) \\  \kappa (x^{(2)}, x^{(1)}) & \kappa (x^{(2)}, x^{(2)}) & ... & \kappa (x^{(2)}, x^{(n)}) \\ ... & ... & ... & ... \\ \kappa (x^{(n)}, x^{(1)}) & \kappa (x^{(n)}, x^{(2)}) & ... & \kappa (x^{(n)}, x^{(n)}) \end{array} \right] 
        • 100개의 훈련 샘플이 담긴 데이터셋이라면 각 싸으이 유사도를 담은 대칭 커널 행렬은 100 x 100 차원이 된다.
    2. 다음 식을 사용하여 커널 행렬 K 를 중앙에 맞춘다.
      • K' = K - 1_{n} K - K 1_{n} + 1_{n} K 1_{n}
      • 여기서 1_{n} 은 모든 값이 {1 \over n} n \times n 차원 행렬이다. (커널 행렬과 같은 차원)
    3. 고윳값 크기대로 내림차순으로 정렬하여 중앙에 맞춘 커널 행렬에서 최상위 k 개의 고유 벡터를 고른다. 표준 PCA와 다르게 고유 벡터는 주성분 축이 아니며, 이미 이 축에 투영된 샘플이다.
  • 위 단계의 2번째에서 왜 커널 행렬을 중앙에 맞추었는지 궁금할 수 있다.
    • 우리는 앞서 표준화된 전처리된 데이터를 다룬다고 가정했다. 공분산 행렬을 구성하고 비선형 특성 조합으로 내적을 \phi 를 사용한 비선형 특성 조합으로 내적을 대체할 때 사용한 모든 특성의 평균이 0이다.
    • 반면 새로운 특성 공간을 명시적으로 계산하지 않기 때문에 이 특성 공간이 중앙에 맞추어져 있는지 보장할 수 없다.
    • 이것이 새로운 두 번째 단계에서 커널 행렬의 중앙을 맞추는 것이 필요한 이유이다.

파이썬으로 커널 PCA 구현

  • RBF 커널 PCA 코드
from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh
import numpy as np

def rbf_kernel_pca(X, gamma, n_components):
    """
   RBF 커널 PCA 구현

    매개변수
   ----------
   X: {넘파이 ndarray}, shape = [n_samples, n_features]

    gamma: float
     RBF 커널 튜닝 매개변수
    n_components: int
     변환한 주성분 개수   

    반환값
   -----------
   X_pc: {넘파이 ndarray}, shape = [n_samples, k_features]
      투영된 데이터셋
    """

   # M x N 차원의 데이터셋에서 샘플 간의 유클리디안 거리의 제곱을 계산
    sq_dists = pdist(X, 'sqeuclidean')

   # 샘플 간의 거리를 정방 대칭 행렬로 반환
    mat_sq_dists = squareform(sq_dists)

    # 커널 행렬을 계산
    K = exp(-gamma * mat_sq_dists)

   # 커널 행렬을 중앙에 맞춘다
    N = K.shape[0]
   one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

   # 중앙에 맞춰진 커널 행렬의 고윳값과 고유 벡터를 구한다.
    # scipy.linalg.eigh 함수는 오름차순으로 반환한다.
    eigvals, eigvecs = eigh(K)
    eigvals, eigvecs = eigvals[::-1], eigvecs[:, ::-1]

   # 최상위 k개의 고유 벡터를 선택한다(투영 결과)
    X_pc = np.column_stack([eigvecs[:, i] for i in range(n_components)])

    return X_pc

예제 1

  • 반달 모양을 띤 100개의 샘플로 구성된 2차원 데이터셋을 구성
from sklearn.datasets import make_moons
import matplotlib.pyplot as plt

X, y = make_moons(n_samples=100, random_state=123)

plt.scatter(X[y==0, 0], X[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y==1, 0], X[y==1, 1], color='blue', marker='o', alpha=0.5)
plt.show()

  • PCA의 주성분에 데이터셋을 투영한다.
scikit_pca = PCA(n_components=2)

X_spca = scikit_pca.fit_transform(X)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))

ax[0].scatter(X_spca[y==0, 0], X_spca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y==1, 0], X_spca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_spca[y==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y==1, 0], np.zeros((50,1))-0.02, color='blue', marker='^', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')

plt.show()

  • 이 데이터에 앞서 작성한 rbf_kernel_pca를 적용해 보자
X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))

ax[0].scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_kpca[y==1, 0], X_kpca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_kpca[y==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_kpca[y==1, 0], np.zeros((50,1))-0.02, color='blue', marker='^', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')

plt.tight_layout()
plt.show()

  • 이제 두 클래스는 선형적으로 구분이 잘 되므로 선형 분류기를 위한 훈련 데이터로 적합하다.
  • 아쉽지만 보편적인 \gamma 파라미터 값은 없다. 주어진 문제에 적합한 \gamma 를 찾으려면 실험이 필요하다.

예제 2

  • 동심원 데이터 셋을 구성
from sklearn.datasets import make_circles

X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)

plt.scatter(X[y==0, 0], X[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y==1, 0], X[y==1, 1], color='blue', marker='o', alpha=0.5)
plt.tight_layout()
plt.show()

  • 기본 PCA 적용
scikit_pca = PCA(n_components=2)

X_spca = scikit_pca.fit_transform(X)

fit, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))

ax[0].scatter(X_spca[y==0, 0], X_spca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y==1, 0], X_spca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_spca[y==0, 0], np.zeros((500,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y==1, 0], np.zeros((500,1))-0.02, color='blue', marker='^', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')

plt.tight_layout()
plt.show()

  • 적절한 gamma를 부여하여 RBF 커널 PCA를 구현
X_kpca = RbfKernelPCA.rbf_kernel_pca(X, gamma=15, n_components=2)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))

ax[0].scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_kpca[y==1, 0], X_kpca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_kpca[y==0, 0], np.zeros((500,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_kpca[y==1, 0], np.zeros((500,1))-0.02, color='blue', marker='^', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')

plt.tight_layout()
plt.show()

새로운 데이터 포인트 투영

  • 앞선 예제는 하나의 데이터셋을 새로운 특성에 투영했는데, 실전에서는 하나 이상의 변환해야 할 데이터셋이 있다.
  • 장의 서두에서 보았던 기본 PCA 방법을 보면 변환 행렬과 입력 샘플 사이의 내적을 계산해서 데이터를 투영했다.
    • 변환 행렬의 열은 공분산 행렬에서 얻은 최상위 k 개의 고유 벡터(v) 이다.
  • 커널 PCA 이면의 아이디어로 돌아가 보면 중심을 맞춘 커널 행렬의 고유 벡터(a) 를 구했다.
    • 즉 샘플은 이미 주성분 축 v 에 투영되어 있다.
    • 새로운 샘플 x' 를 주성분 축에 투영하려면 다음을 계산해야 한다.

\phi (x')^{T} v

  • 다행히 커널 트릭을 사용하여 명시적으로 \phi (x')^{T} v 를 계산할 필요가 없다.
    • 기본 PCA와 다르게 커널 PCA는 메모리 기반 방법이다.
    • 즉, 새로운 샘플을 투영하기 위해 매번 원본 훈련 세트를 재사용해야 한다.
    • 훈련 세트에 있는 i 번째 새로운 샘플과 새로운 샘플 x' 사이 RBF 커널(유사도)을 계산해야 한다.

\phi (x')^{T} v = \sum_{i} a^{(i)} \phi (x')^{T} \phi (x^{(i)}) = \sum_{i} a^{(i)} \kappa (x', x^{(i)})

  • 커널 행렬 K 의 고유 벡터 a 와 고윳값 \lambda 는 다음 식을 만족한다.

Ka = \lambda a

  • 새로운 샘플과 훈련 세트의 샘플 간 유사도를 계산한 후 고윳값으로 고유 벡터 a 를 정규화 해야 한다.
from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh
import numpy as np

def rbf_kernel_pca(X, gamma, n_components):
    """
   RBF 커널 PCA 구현

    매개변수
   ----------
   X: {넘파이 ndarray}, shape = [n_samples, n_features]

    gamma: float
     RBF 커널 튜닝 매개변수
    n_components: int
     변환한 주성분 개수   

    반환값
   -----------
   X_pc: {넘파이 ndarray}, shape = [n_samples, k_features]
      투영된 데이터셋
    """

   # M x N 차원의 데이터셋에서 샘플 간의 유클리디안 거리의 제곱을 계산
    sq_dists = pdist(X, 'sqeuclidean')

   # 샘플 간의 거리를 정방 대칭 행렬로 반환
    mat_sq_dists = squareform(sq_dists)

    # 커널 행렬을 계산
    K = exp(-gamma * mat_sq_dists)

   # 커널 행렬을 중앙에 맞춘다
    N = K.shape[0]
   one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

   # 중앙에 맞춰진 커널 행렬의 고윳값과 고유 벡터를 구한다.
    # scipy.linalg.eigh 함수는 오름차순으로 반환한다.
    eigvals, eigvecs = eigh(K)
   eigvals, eigvecs = eigvals[::-1], eigvecs[:, ::-1]

   # 최상위 k개의 고유 벡터를 선택한다
alphas = np.column_stack([eigvecs[:, i] for i in range(n_components)])

# 고유 벡터에 상응하는 고윳값을 선택한다.
lambdas = [eigvals[i] for i in range(n_components)]

return alphas, lambdas
  • 수정된 커널 PCA를 이용하여 반달 데이터 셋을 테스트 해 보자
X, y = make_moons(n_samples=100, random_state=123)
alphas, lambdas = RbfKernelPCA2.rbf_kernel_pca(X, gamma=15, n_components=1)

x_new = X[25]
x_proj = alphas[25]

def project_x(x_new, X, gamma, alphas, lambdas):
    pair_dist = np.array([np.sum((x_new - row)**2) for row in X])
   k = np.exp(-gamma * pair_dist)
    return k.dot(alphas/lambdas)

x_reproj = project_x(x_new, X, gamma=15, alphas=alphas, lambdas=lambdas)

plt.scatter(alphas[y==0, 0], np.zeros((50)), color='red', marker='^', alpha=0.5)
plt.scatter(alphas[y==1, 0], np.zeros((50)), color='blue', marker='o', alpha=0.5)
plt.scatter(x_proj, 0, color='black', label='original projection of point X[25]', marker='x', s=100)
plt.scatter(x_proj, 0, color='green', label='remapped point X[25]', marker='x', s=500)
plt.legend(scatterpoints=1)
plt.tight_layout()
plt.show()

사이킷런의 커널 PCA

  • 편리하게도 사이킷런은 sklearn.decomposition 모듈 아래 커널 PCA 클래스를 구현해 두었으므로 아래와 같이 사용하면 된다.
from sklearn.datasets import make_moons
from sklearn.decomposition import KernelPCA
import matplotlib.pyplot as plt

X, y = make_moons(n_samples=100, random_state=123)

scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)

X_skernpca = scikit_kpca.fit_transform(X)

plt.scatter(X_skernpca[y==0, 0], X_skernpca[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X_skernpca[y==1, 0], X_skernpca[y==1, 1], color='blue', marker='o', alpha=0.5)
plt.xlabel('PC1')
plt.xlabel('PC2')
plt.tight_layout()
plt.show()

[ssba]

The author

Player가 아니라 Creator가 되려는 자/ suyeongpark@abyne.com

댓글 남기기

이 사이트는 스팸을 줄이는 아키스밋을 사용합니다. 댓글이 어떻게 처리되는지 알아보십시오.