[라이트 머신러닝] Session 15. 커널 PCA를 이용한 비선형 매핑

C. 커널 PCA를 사용한 비선형 매핑

여태까지 많은 머신 러닝 알고리즘은 입력 데이터가 선형적으로 구분이 가능하다는 가정을 합니다. 다른 알고리즘들-아달린, 로지스틱 회귀, SVM-은 선형적으로 완벽하게 분리되지 않는 이유를 잡음때문이라고 이야기합니다.

실전에서는 더 자주 비선형 문제들을 맞닥뜨립니다. 이 경우에 항상 PCA나 LDA와 같은 차원 축소 기법이 최선이라고는 말할 수 없겠죠. 이제부터 알아볼 것은 PCA의 커널화 버전인 KPCA입니다.

 

 

1. 커널 함수와 커널 트릭

앞선 세션에서 커널 SVM에 대해 이야기한 것을 떠올려보면, 비선형 문제를 풀기 위해 고차원 공간으로 데이터를 투영해 풀었습니다. k 고차원 부분 공간에 있는 샘플을 변환하기 위해 비선형 매핑 함수를 정의합니다.

 

 

이 함수를 d차원 보다 더 큰 k차원으로 매핑하기 위해 원본 특성의 비선형 조합을 만듭니다. 2차원의 특성벡터를 매핑 가능한 3차원 공간은 아래와 같습니다.

 

 

중간 정리를 하면, 커널 PCA로 비선형 매핑을 수행해서 고차원 공간으로 변환하고, 표준 PCA로 샘플이 선형 분류기로 구분될 수 있는 저차원으로 데이터를 투영시키는 것입니다.

이 방식의 단점은 계산 비용이 매우 비싸다는 점입니다. 이 부분에서 커널 트릭(kernel trick)이 등장합니다. 커널 트릭을 사용하면 원본 특성 공간에서 두 고차원 특성 벡터의 유사도를 계산할 수 있습니다. 커널 트릭의 역할을 알아보기 전에 다시 PCA 방식의 수식을 생각해볼까요?

특성 k와 j의 공분산은 아래와 같습니다.

 

 

여기서 우리는 특성 평균 뮤를 0에 맞추었으므로 아래와 같이 간단히 쓸 수 있습니다.

 

 

이 식은 두 특성 간의 공분산입니다. 이제 공분산 행렬의 일반식으로 만들어줍니다.

 

 

베른하르트 슐코프는 이 방식을 일반화해 파이를 통한 비선형 특성 조합으로 원본 특성공간의 샘플간 점곱을 대체했습니다. 그 식이 아래와 같 습니다.

 

 

그리고 나서 고유벡터를 얻기 위해 아래 식을 풀어줍니다.

 

 

여기서 람다와 v는 공분산 행렬의 고윳값과 고유 벡터입니다. a는 커널(유사도) 행렬 K의 고유벡터를 추출함으로써 구할 수 있습니다.

 

 

커널 SVM을 생각해보면 커널 트릭을 사용해 샘플 x끼리의 파이 함수 점곱을 커널 함수 K로 바꿔 고유 벡터를 계산할 필요가 없었습니다.

커널 PCA는 그냥 PCA처럼 투영행렬을 구성한 게 아니라 각각 성분에 이미 투영된 샘플입니다. 기본적으로 커널은 두 벡터 사이의 점곱을 계산할 수 있는 함수입니다. 즉 유사도를 측정할 수 있는 함수인 것이죠.

가장 널리 사용되는 커널은 아래와 같습니다.

1. 다항 터널 - 세타는 임계값이고 P는 사용자가 지정한 거듭제곱

 

 

2. 하이퍼볼릭 탄젠트 커널

 

 

3. 방사 기저 함수 혹은 가우시안 커널

 

 

지금까지 배운 것을 요약하면 커널 PCA를 구현하기 위해 다음 세 단계를 정의할 수 있습니다.

1. 커널 행렬 K를 계산하고 샘플의 모든 쌍에 대해 구합니다.

 

 

 

2. 아래 식을 사용해 커널 행렬 K를 중앙에 맞춰줍니다. 여기사 In은 모든 값이 1/n인 n * n 차원 행렬입니다.

 

 

3. 고윳값 크기대로 내림차순 정렬해서 중앙에 맞춘 커널 행렬 중 최상위 k개 벡터를 고릅니다. 여기서 고유벡터는 주성분 축은 아닙니다.이미 이 축에 투영된 샘플입니다.

여기서 왜 커널 행렬을 중앙에 맞추었는지 궁금하시죠? 앞서 우리는 이미 표준화 전처리된 데이터를 다룬다고 가정했고, 공분산 행렬 구성 후 비선형 특성 조합의 점곱을 파이를 이용한 비선형 특성 조합으로 점곱 대체할 때 모든 특성 평균이 0임을 기억하시죠?

하지만 새로은 특성 공간을 계산하지 않기 때문에 이 특성 공간이 중앙인지 확신할 수가 없기 때문입니다.

2. 파이썬으로 커널 PCA 구현하기

사이파이와 넘파이 헬퍼 함수를 사용해 커널 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):
    # MxN 차원의 데이터셋에서 샘플 간의 유클리디안 거리의 제곱을 계산합니다.
    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

차원 축소에 RBF 커널 PCA를 사용하는 단점은 사전에 람다 매개변수를 지정해야 하는데, 적절한 값을 찾기 위해 실험이 필요하다는 점입니다. 이 다음 세션들에서 소개할 그리드 서치와 같은 튜닝 알고리즘을 사용하는 것이 최선입니다.

이제 예제 두 가지를 사용해 실험하도록 하겠습니다. 첫 번째, 반달모양 구분하기 입니다. 먼저 데이터 셋을 만들어줍니다.

import matplotlib.pyplot as plt
from sklearn.datasets import make_moons

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.tight_layout()
plt.show()

 

 

확실히 비선형 데이터이죠? 이제 커널 PCA로 펼쳐서 선형 분류기에 적합하게 만들어 줍니다. 먼저, 기본 PCA에 투영하면 어떻게 보이는지부터 확인해보겠습니다.

from sklearn.decomposition import 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='o', 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 함수 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='o', 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()

 

 

이렇게 예쁜 데이터를 얻을 수 있습니다.

두 번째 예제는 원을 분리하는 것입니다. 일단 데이터를 확인해보겠습니다.

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)

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((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='o', 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는 선형 분류기에 적합한 결과를 낼 수 없습니다. 이제 적절한 람다 값을 주고 RBF 커널 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((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='o', 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()

 

 

3. 새로운 데이터 포인트 투영

이제까지는 하나의 데이터셋을 새로운 특성에 투영했지만, 실전에서는 변환해야할 데이터셋이 하나 이상입니다. 예를 들면 훈련과 테스트 데이터처럼 말이죠.

PCA에서는 변환행렬 W와 입력 샘플의 점곱으로 데이터를 투영했습니다. 그렇다면 커널 PCA는 어떨까요? 커널 PCA는 중심을 맞춘 커널 행렬의 고유벡터를 구했습니다. 따라서 샘플은 이미 주성분 축인 v에 들어있습니다. 새로운 샘플 x'를 투영하기 위해서는 아래를 계산해야합니다.

 

 

다행히 우리는 커널 트릭을 사용하므로 이것을 계산할 필요는 없습니다. 일반 PCA와 커널 PCA의 다른 점은 메모리 기반이라는 것입니다. 즉 새로운 샘플을 투영하기 위해서는 매번 원본 훈련 세트를 재사용해야합니다.

훈련 세트의 i번째 샘플과 새로운 샘플 x'사이의 RBF 유사도를 계산해야합니다.

 

 

커널 행렬 K의 고유벡터 a와 고윳값 람다는 아래 관계를 만족합니다.

 

 

from scipy import exp
from scipy.linalg import eigh
import numpy as np

def rbf_kernel_pca(X, gamma, n_components):
    # MxN 차원의 데이터셋에서 샘플 간의 유클리디안 거리의 제곱을 계산합니다.
    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 구현으로 1차원 부분공간에 투영하겠습니다.

X, y = make_moons(n_samples=100, random_state=123)
alphas, lambdas = rbf_kernel_pca(X, gamma=15, n_components=1)

새로운 샘플을 투영하기 전에 26번째 포인트가 새로운 데이터 포인트 x'라고 가정합니다.

x_new = X[25]
x_new
#array([1.8713, 0.0093])

x_proj = alphas[25] # 원본 투영
x_proj
#array([0.0788])


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)
x_reproj
#array([0.0788])

마지막으로 투영한 것을 그래프로 그려보겠습니다.

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='^', s=100)
plt.scatter(x_reproj, 0, color='green',
            label='remapped point X[25]', marker='x', s=500)
plt.legend(scatterpoints=1)

plt.tight_layout()
plt.show()

 

 

이렇게 새로운 데이터 포인트 x'가 올바르게 투영되었습니다.

4. 사이킷런의 커널 PCA

사이킷런은 sklearn.decomposition 모듈 아래에 커널 PCA를 구현해두었습니다. kernel 매개변수로 커널 종류를 지정합니다. 그리고 우리가 구현한 커널 PCA와 비교를 위해 변환된 반달 데이터를 두 개 주성분에 그려보도록 하겠습니다.

from sklearn.decomposition import KernelPCA

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.ylabel('PC2')
plt.tight_layout()
plt.show()

 

 

위 그래프를 통해 사이킷런의 kernelPCA 결과는 직접 구현한 것과 같음을 알 수 있습니다.

 

---

 여기까지 커널 PCA를 알아보았습니다. 사실 저는 아직도 계산 비용을 줄이는 커널 트릭에 대해서 제대로 이해가 되지 않았는데요ㅜㅜ 이 부분에 대해서는 이후에 공부해서 어떻게 작동하는지 추가로 업로드(수정) 하도록 하겠습니다.(2020.02.20) 다음 세션에서는 모델을 평가하는 방법에 대해서 다루도록 하겠습니다. 다음 세션에서 봬요!