Search
Duplicate

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

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

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

주성분 분석의 주요 단계

PCA는 비지도 선형 변환 기법으로 주로 특성 추출과 차원 축소 용도로 많은 분야에서 널리 사용된다.
PCA는 특성 사이의 상관관계를 기반으로 하여 데이터에 있는 특성을 잡아낼 수 있다.
요약하자면 PCA는 고차원 데이터에서 분산이 가장 큰 방향을 찾고 좀 더 작거나 같은 수의 차원을 갖는 새로운 부분 공간으로 이를 투영한다.
새로운 부분 공간의 직교 좌표는 주어진 조건 하에서 분산이 최대인 방향으로 해석할 수 있다.
새로운 특성 축은 아래 그림과 같이 서로 직각을 이룬다. 아래 그림에서 x1,x2x1, x2는 원본 특성 축이고 PC1,PC2PC1, PC2는 주성분이다.
PCA를 사용하여 차원을 축소하기 위해 d×kd \times k차원의 변환행렬 WW를 만든다.
이 행렬로 샘플 벡터 xx를 새로운 kk차원의 특성 부분 공간으로 매핑한다.
이 부분 공간은 원본 dd차원의 특성 공간보다 작은 차원을 가진다.
x=[x1,x2,...,xd],xRdx = [x_{1}, x_{2}, ... , x_{d}], x \in \mathbb{R}^{d}
xW,WRd×k\downarrow xW, W \in \mathbb{R}^{d \times k}
z=[z1,z2,...,zk],zRkz = [z_{1}, z_{2}, ... , z_{k}], z \in \mathbb{R}^{k}
원본 dd차원 데이터를 새로운 kk차원의 부분 공간(일반적으로 k<dk < d)으로 변환하여 만들어진 첫 번째 주성분이 가장 큰 분산을 가질 것이다.
모든 주성분은 다른 주성분들과 상관관계가 있더라도 만들어진 주성분은 서로 직각을 이룰 것이다.
PCA 방향은 데이터 스케일에 매우 민감하다.
특성의 스케일이 다르고 모든 특성의 중요도를 동일하게 취급하려면 PCA를 적용하기 전에 특성을 표준화 전처리해야 한다.
차원 축소 PCA 알고리즘의 단계는 다음과 같다.
1.
dd차원 데이터셋을 표준화 전처리한다.
2.
공분산 행렬(covariance matrix)을 만든다.
3.
공분산 행렬을 고유 벡터(eigenvector)와 고윳값(eigenvalue)으로 분해한다.
4.
고윳값을 내림차순으로 정렬하고 그에 해당하는 고유 벡터의 순위를 매긴다.
5.
고윳값이 가장 큰 kk개의 고유 벡터를 선택한다. 여기서 kk는 새로운 특성 부분 공간의 차원이다. (kdk \leq d)
6.
최상위 kk개의 고유벡터로 투영행렬(projection matrix) WW를 만든다.
7.
투영 행렬 WW를 사용해서 dd차원 입력 데이터셋 XX를 새로운 kk차원의 특성 부분 공간으로 변환한다.

주성분 추출단계

우선 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)
Python
그 다음으로 공분산 행렬을 만드는 단계를 진행한다. 공분산 행렬은 d×dd \times d차원의 대칭 행렬로 특성 상호 간의 공분산을 저장한다. dd는 데이터셋에 있는 차원 개수이다.
예컨대 전체 샘플에 대한 두 특성 xjx_{j}와 xkx_{k}사이의 공분산은 다음 식으로 계산할 수 있다.
σjk=1ni=1n(xj(i)μj)(xk(i)μk)\sigma_{jk} = {1 \over n} \sum_{i = 1}^{n} (x_{j}^{(i)} - \mu_{j})(x_{k}^{(i)} - \mu_{k})
여기서 μj\mu_{j}와 μk\mu_{k}는 특성 jj와 kk의 샘플 평균이다.
데이터셋을 표준화 전처리했기 때문에 샘플 평균은 0이다.
두 특성 간 양의 공분산은 특성이 함께 증가하거나 감소하는 것을 나타낸다. 반면 음의 공분산은 특성이 반대 방향으로 달라진다는 것을 나타낸다.
예컨대 세 개의 특성으로 이루어진 공분산 행렬은 다음과 같이 쓸 수 있다. (여기서 Σ\Sigma는 대문자 표기일 뿐 합을 뜻하는 것이 아니므로 주의)
Σ=[σ12σ12σ13σ21σ22σ23σ31σ32σ32]\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개의 고유벡터와 고윳값을 얻을 수 있다.
다음 단계로 공분산 행렬의 고유벡터와 고윳값의 쌍을 구해 보자.
선형대수학에 따르면 고유벡터 vv는 다음 식을 만족한다.
Σv=λ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)
Python
Note) numpy.linalg.eig 함수는 대칭과 비대칭 정방 행렬을 모두 다룰 수 있지만 이따금 복소수 고윳값을 반환한다. 이와 비슷하게 에르미트 행렬을 분해하기 위해 구현된 numpy.linalg.eigh 함수는 공분산 행렬과 같은 대칭 행렬을 다룰 때 수치적으로 더 안정된 결과를 만든다. (numpy.linalg.eigh는 항상 실수 고윳값을 반환한다.)
Note) 사이킷런의 PCA 클래스는 직접 고윳값과 고유 벡터를 계산하는 대신 특이값 분해 (singular value decomposition) 방식을 이용하여 주성분을 구한다.

총분산과 설명된 분산

데이터셋 차원을 새로운 특성 부분 공간으로 압축하여 줄여야 하기에 가장 많은 정보(분산)를 가진 고유 벡터(주성분) 일부만 선택한다.
고윳값은 고유 벡터의 크기를 결정하므로 고윳값을 내림차순으로 정렬한다.
고윳값 순서에 따라 최상위 kk개의 고유벡터를 선택한다.
가장 정보가 많은 kk개의 고유 벡터를 선택하기 전에 고윳값의 설명된 분산(explained variance) 그래프로 그려보겠다.
고윳값 λj\lambda_{j}의 설명된 분산 비율은 전체 고윳값의 합에서 고윳값 λj\lambda_{j}의 비율이다.
λjdj=1λ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()
Python
랜덤 포레스트는 클래스 소속 정보를 사용하여 노드의 불순도를 계산하는 방면, 분산은 특성 축을 따라 값들이 퍼진 정도를 측정한다.

특성 변환

공분산 행렬을 고유 벡터와 고윳값 쌍으로 성공적으로 분해한 후 Wine 데이터 셋을 새로운 주성분 축으로 변환하는 3개의 단계는 다음과 같다.
고윳값이 가장 큰 kk개의 고유 벡터를 선택한다. 여기서 kk는 새로운 특성 부분 공간의 차원이다 (kdk \leq d)
최상위 kk개의 고유벡터로 투영 행렬 WW를 만든다.
투영행렬 WW를 사용해서 dd차원 입력 데이터셋 XX를 새로운 kk차원의 특성 부분 공간으로 변환한다.
좀 더 쉽게 설명하면 고윳값의 내림차순으로 고유 벡터를 정렬하고 선택된 고유 벡터로 투영 행렬을 구성한다. 이 투영 행렬을 사용하여 데이터를 저차원 부분 공간으로 변환한다.
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)
Python
앞의 코드를 이용하여 최상위 두 개의 고유 벡터로부터 13 x 2 차원의 투영행렬 WW를 만든다.
투영 행렬을 사용하면 샘플 xx (1 x 13 차원의 행 벡터)를 PCA 부분 공간(두 개의 주성분)을 투영하여 xx'를 얻을 수 있다. 두 개의 특성으로 구성된 2차원 샘플 벡터이다.
x=xWx' = xW
w = np.hstack((eigen_pairs[0][1][:, np.newaxis], eigen_pairs[1][1][:, np.newaxis]))
Python
비슷하게 124 x 13 차원의 훈련 데이터셋을 행렬 내적으로 두 개의 주성분으로 변환할 수 있다.
X=XWX' = XW
X_train_pca = X_train_std.dot(w)
Python
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()
Python
위 이미지에서 볼 수 있듯이 데이터가 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)
Python
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)
Python
테스트 세트에 적용한 결과는 아래와 같다.

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

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

산포 행렬 계산

PCA에서 했기 때문에 데이터셋의 특성을 표준화하는 단계는 건너뛰고 바로 평균 벡터 계산을 진행한다.
평균 벡터를 사용하여 클래스 간의 산포 행렬과 클래스 내 산포 행렬을 구성한다.
평균 벡터 mim_{i}는 클래스 ii의 샘플에 대한 특성의 평균값 μm\mu_{m}을 저장한다.
mi=1nixDicxmm_{i} = {1 \over n_{i}} \sum_{x \in D_{i}}^{c} x_{m}
3개의 평균 벡터가 만들어진다.
mi=[μi,alcoholμi,malic acid...μi,proline]i{1,2,3}m_{i} = \left[ \begin{array}{rrrr} \mu_{i, \text{alcohol}} \\ \mu_{i, \text{malic acid}} \\ ... \\ \mu_{i, \text{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))
Python
평균 벡터를 사용하여 클래스 내 산포 행렬 SWS_{W}를 계산할 수 있다.
SW=i=1cSiS_{W} = \sum_{i=1}^{c} S_{i}
이 행렬은 개별 클래스 ii의 산포 행렬 SiS_{i}를 더하여 구한다.
Si=xDic(xmi)(xmi)TS_{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
Python
개별 산포 행렬 SiS_{i}를 산포 행렬 SWS_{W}로 모두 더하기 전에 스케일 조정을 해야 한다.
산포 행렬을 클래스 샘플 개수 ηi\eta_{i}로 나누면 사실 산포 행렬을 계산하는 것이 공분산 행렬 Σi\Sigma_{i}를 계산하는 것과 같아진다.
즉 공분산 행렬은 산포 행렬의 정규화 버전이다.
Σi=1niSi=1nixDic(xmi)(xmi)T\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
Python
클래스 내 산포 행렬(또는 공분산 행렬)을 계산한 후 다음 단계로 넘어가 클래스 간의 산포 행렬 SBS_{B}를 계산한다.
SB=i=1cni(mim)(mim)TS_{B} = \sum_{i = 1}^{c} n_{i} (m_{i} - m)(m_{i} - m)^{T}
여기서 mm은 모든 클래스의 샘플을 포함하여 계산된 전체 평균이다.
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)
Python

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

LDA의 남은 단계는 PCA와 유사하다. 공분산 행렬에 대한 고윳값 분해를 수행하는 대신 행렬 SW1SBS_{W}^{-1}S_{B}의 고윳값을 계산하면 된다.
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
Python
고유 벡터와 고윳값 쌍을 계산한 후 내림차순으로 고윳값을 정렬한다.
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)
Python
LDA에서 선형 판별 벡터는 최대 c1c-1개이다. cc는 클래스 레이블의 개수이다. 클래스 내 산포 행렬 SBS_{B}가 랭크(rank) 1 또는 그 이하인 cc개의 행렬을 합한 것이기 때문이다.
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()
Python
두 개의 판별 고유 벡터를 열로 쌓아서 변환 행렬 WW를 만든다.
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real, eigen_pairs[1][1][:, np.newaxis].real))
Python

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

이전 절에서 만든 변환 행렬 WW를 훈련 세트에 곱해서 데이터를 변환할 수 있다.
X=XWX' = 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()
Python

사이킷런의 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()
Python
테스트 결과는 다음과 같다.

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

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

커널 함수와 커널 트릭

커널 SVM에 관해 배운 것을 떠올려 보면 비선형 문제를 해결하기 위해 클래스가 선형으로 구분되는 새로운 고차원 특성 공간으로 투영할 수 있었다.
kk 고차원 부분 공간에 있는 샘플 xRdx \in \mathbb{R}^{d}를 변환하기 위해 비선형 매핑 함수 ϕ\phi 를 정의한다.
ϕ:RdRk(k>d)\phi : \mathbb{R}^{d} \to \mathbb{R}^{k} (k > d)
ϕ\phi 함수를 dd 차원의 원본 데이터셋에서 더 큰 kk 차원의 특성 공간으로 매핑하기 위해 원본 특성의 비선형 조합을 만드는 함수로 생각할 수 있다.
예컨대 2차원 (d=2) 특성 벡터 xRdx \in \mathbb{R}^{d} 가 있으면 (xx는 dd개의 특성으로 구성된 열 벡터) 매핑 가능한 3D 공간은 다음과 같다.
x=[x1,x2]Tx = [x_{1}, x_{2}]^{T}
ϕ\downarrow \phi
z=[x12,2x1x2,x22]Tz = [x_{1}^{2}, \sqrt{2x_{1}x_{2}}, x_{2}^{2}]^{T}
다른 말로 하면 커널 PCA를 통한 비선형 매핑을 수행하여 데이터를 고차원 공간으로 변환한다.
그다음 고차원 공간에 표준 PCA를 사용하여 샘플이 선형 분류기로 구분될 수 있는 저차원 공간으로 데이터를 투영한다 (샘플이 이 입력 공간에서 잘 구분될 수 있다고 가정).
이 방식의 단점은 계산 비용이 매우 비싸다는 것이다. 여기에 커널 트릭(kernel trick)이 등장한다.
커널 트릭을 사용하면 원본 특성 공간에서 두 고차원 특성 벡터의 유사도를 계산할 수 있다.
커널 트릭에 대해 알아보기 전에 표준 PCA 방식을 다시 살펴보자
두 개의 특성 k,jk, j 사이의 공분산은 다음과 같이 계산한다.
σjk=1ni=1n(xj(i)μj)(xk(i)μk)\sigma_{jk} = {1 \over n} \sum_{i = 1}^{n} (x_{j}^{(i)} - \mu_{j})(x_{k}^{(i)} - \mu_{k})
μj=0,μk=0\mu_{j} = 0, \mu_{k} = 0처럼 특성 평균을 0에 맞추었으므로 이 식은 다음과 같이 간단히 쓸 수 있다.
σjk=1ni=1nxj(i)xk(i)\sigma_{jk} = {1 \over n} \sum_{i = 1}^{n} x_{j}^{(i)} x_{k}^{(i)}
이 식은 두 특성 간의 공분산을 의미한다. 공분산 행렬 Σ\Sigma 를 계산하는 일반식으로 써보자.
Σ=1ni=1nx(i)x(i)T\Sigma = {1 \over n} \sum_{i = 1}^{n} x^{(i)} x^{(i)^{T}}
베른하르트 슐코프(Bernhard Scholkopf)는 이 방식을 일반화 하여 ϕ\phi를 통한 비선형 특성 조합으로 원본 특성 공간의 샘플 사이의 내적을 대체했다.
Σ=1ni=1nϕ(x(i))ϕ(x(i))T\Sigma = {1 \over n} \sum_{i = 1}^{n} \phi (x^{(i)}) \phi(x^{(i)})^{T}
이 공분산 행렬에서 고유 벡터(주성분)를 얻기 위해서는 다음 식을 풀어야 한다.
Σv=λv\Sigma v = \lambda v
1ni=1nϕ(x(i))ϕ(x(i))Tv=λv\Rightarrow {1 \over n} \sum_{i = 1}^{n} \phi (x^{(i)}) \phi(x^{(i)})^{T} v = \lambda v
v=1nλi=1nϕ(x(i))ϕ(x(i))T=i=1na(i)ϕ(x(i))\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와 vv는 공분산 행렬 Σ\Sigma의 고윳값과 고유 벡터이다. aa는 커널(유사도) 행렬 KK의 고유 벡터를 추출함으로써 구할 수 있다.
커널 SVM을 사용하여 비선형 문제 풀기를 떠올리면 커널 트릭을 사용하여 샘플 xx 끼리의 ϕ\phi 함수 내적을 커널 함수 KK로 바꿈으로써 고유 벡터를 명시적으로 계산할 필요가 없었다.
κ(x(i),x(j))=ϕ(x(i))Tϕ(x(j))\kappa (x^{(i)}, x^{(j)}) = \phi (x^{(i)})^{T} \phi (x^{(j)})
다른 말로 하면 커널 PCA로 얻은 것은 표준 PCA 방식에서처럼 투영 행렬을 구성한 것이 아니고 각각의 성분에 이미 투영된 샘플이다.
기본적으로 커널 함수는 두 벡터 사이의 내적을 계산할 수 있는 함수이다. 즉, 유사도를 측정할 수 있는 함수이다.
가장 널리 사용되는 커널은 다음과 같다.
다항 커널
κ(x(i),x(j))=(x(i)Tx(j)+θ)P\kappa (x^{(i)}, x^{(j)}) = (x^{(i)T} x^{(j)} + \theta)^{P}
여기서 θ\theta는 임계 값이고 PP는 사용자가 지정한 거듭제곱이다.
하이퍼볼릭 탄젠트(hyperbolic tangent) (시그모이드(sigmoid)) 커널
κ(x(i),x(j))=tanh(ηx(i)Tx(j)+θ)\kappa (x^{(i)}, x^{(j)}) = tanh (\eta x^{(i)T} x^{(j)} + \theta)
방사 기저 함수(Radial Basis Function, RBF) 또는 가우시안 커널
κ(x(i),x(j))=exp(x(i)x(j)22σ2)\kappa (x^{(i)}, x^{(j)}) = \exp(- {\|x^{(i)} - x^{(j)}\|^{2} \over 2 \sigma^{2}})
변수 γ=12σ2\gamma = {1 \over 2 \sigma^{2}} 을 도입하여 종종 다음과 같이 쓴다.
κ(x(i),x(j))=exp(γx(i)x(j)2)\kappa (x^{(i)}, x^{(j)}) = \exp(- \gamma \|x^{(i)} - x^{(j)}\|^{2})
지금까지 배운 것을 요약하면 RBF 커널 PCA를 구현하기 위해 다음 3 단계를 정의할 수 있다.
1.
커널 (유사도) 행렬 KK를 다음 식으로 계산한다.
κ(x(i),x(j))=exp(γx(i)x(j)2)\kappa (x^{(i)}, x^{(j)}) = \exp(- \gamma \|x^{(i)} - x^{(j)}\|^{2})
샘플의 모든 쌍에 대해 구한다.
K=[κ(x(1),x(1))κ(x(1),x(2))...κ(x(1),x(n)) κ(x(2),x(1))κ(x(2),x(2))...κ(x(2),x(n))............κ(x(n),x(1))κ(x(n),x(2))...κ(x(n),x(n))]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.
다음 식을 사용하여 커널 행렬 KK를 중앙에 맞춘다.
K=K1nKK1n+1nK1nK' = K - 1_{n} K - K 1_{n} + 1_{n} K 1_{n}
여기서 1n1_{n}은 모든 값이 1n{1 \over n}인  n×nn \times n차원 행렬이다. (커널 행렬과 같은 차원)
3.
고윳값 크기대로 내림차순으로 정렬하여 중앙에 맞춘 커널 행렬에서 최상위 kk개의 고유 벡터를 고른다. 표준 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
Python

예제 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()
Python
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()
Python
이 데이터에 앞서 작성한 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()
Python
이제 두 클래스는 선형적으로 구분이 잘 되므로 선형 분류기를 위한 훈련 데이터로 적합하다.
아쉽지만 보편적인 γ\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()
Python
기본 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()
Python
적절한 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()
Python

새로운 데이터 포인트 투영

앞선 예제는 하나의 데이터셋을 새로운 특성에 투영했는데, 실전에서는 하나 이상의 변환해야 할 데이터셋이 있다.
장의 서두에서 보았던 기본 PCA 방법을 보면 변환 행렬과 입력 샘플 사이의 내적을 계산해서 데이터를 투영했다.
변환 행렬의 열은 공분산 행렬에서 얻은 최상위 kk개의 고유 벡터(v)(v)이다.
커널 PCA 이면의 아이디어로 돌아가 보면 중심을 맞춘 커널 행렬의 고유 벡터(a)(a)를 구했다.
즉 샘플은 이미 주성분 축 vv에 투영되어 있다.
새로운 샘플 xx'를 주성분 축에 투영하려면 다음을 계산해야 한다.
ϕ(x)Tv\phi (x')^{T} v
다행히 커널 트릭을 사용하여 명시적으로 ϕ(x)Tv\phi (x')^{T} v를 계산할 필요가 없다.
기본 PCA와 다르게 커널 PCA는 메모리 기반 방법이다.
즉, 새로운 샘플을 투영하기 위해 매번 원본 훈련 세트를 재사용해야 한다.
훈련 세트에 있는 ii번째 새로운 샘플과 새로운 샘플 xx' 사이 RBF 커널(유사도)을 계산해야 한다.
ϕ(x)Tv=ia(i)ϕ(x)Tϕ(x(i))=ia(i)κ(x,x(i))\phi (x')^{T} v = \sum_{i} a^{(i)} \phi (x')^{T} \phi (x^{(i)}) = \sum_{i} a^{(i)} \kappa (x', x^{(i)})
커널 행렬 KK의 고유 벡터 aa와 고윳값 λ\lambda는 다음 식을 만족한다.
Ka=λaKa = \lambda a
새로운 샘플과 훈련 세트의 샘플 간 유사도를 계산한 후 고윳값으로 고유 벡터 aa를 정규화 해야 한다.
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
Python
수정된 커널 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개의 고유 벡터를 선택한다 alphas = np.column_stack([eigvecs[:, i] for i in range(n_components)]) # 고유 벡터에 상응하는 고윳값을 선택한다. lambdas = [eigvals[i] for i in range(n_components)] return alphas, lambdas
Python

사이킷런의 커널 PCA

편리하게도 사이킷런은 sklearn.decomposition 모듈 아래 커널 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()
Python