5장 - 차원 축소를 사용한 데이터 압축¶
머신 러닝 교과서 3판¶
아래 링크를 통해 이 노트북을 주피터 노트북 뷰어(nbviewer.jupyter.org)로 보거나 구글 코랩(colab.research.google.com)에서 실행할 수 있습니다.
 주피터 노트북 뷰어로 보기
|
 구글 코랩(Colab)에서 실행하기
|
목차¶
- 주성분 분석을 통한 비지도 차원 축소
- 주성분 분석의 주요 단계
- 주성분 추출 단계
- 총분산과 설명된 분산
- 특성 변환
- 사이킷런의 주성분 분석
- 선형 판별 분석을 통한 지도 방식의 데이터 압축
- 주성분 분석 vs 선형 판별 분석
- 선형 판별 분석의 내부 동작 방식
- 산포 행렬 계산
- 새로운 특성 부분 공간을 위해 선형 판별 벡터 선택
- 새로운 특성 공간으로 샘플 투영
- 사이킷런의 LDA
- 커널 PCA를 사용하여 비선형 매핑
- 커널 함수와 커널 트릭
- 파이썬으로 커널 PCA 구현
- 예제 1 - 반달 모양 구분하기
- 예제 2 - 동심원 분리하기
- 새로운 데이터 포인트 투영
- 사이킷런의 커널 PCA
- 요약
In [1]:
from IPython.display import Image
5.1 주성분 분석을 통한 비지도 차원 축소¶
5.1.1 주성분 분석의 주요 단계¶
In [2]:
Image(url='https://git.io/JtsvW', width=400)
Out[2]:
$\boldsymbol{x}\boldsymbol{W}=\boldsymbol{z}$
5.1.2 주성분 추출 단계¶
In [3]:
import pandas as pd
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/'
'machine-learning-databases/wine/wine.data',
header=None)
# UCI 머신 러닝 저장소에서 Wine 데이터셋을 다운로드할 수 없을 때
# 다음 주석을 해제하고 로컬 경로에서 데이터셋을 적재하세요:
# df_wine = pd.read_csv('wine.data', header=None)
df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash',
'Alcalinity of ash', 'Magnesium', 'Total phenols',
'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',
'Color intensity', 'Hue',
'OD280/OD315 of diluted wines', 'Proline']
df_wine.head()
Out[3]:
| Class label | Alcohol | Malic acid | Ash | Alcalinity of ash | Magnesium | Total phenols | Flavanoids | Nonflavanoid phenols | Proanthocyanins | Color intensity | Hue | OD280/OD315 of diluted wines | Proline | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 14.23 | 1.71 | 2.43 | 15.6 | 127 | 2.80 | 3.06 | 0.28 | 2.29 | 5.64 | 1.04 | 3.92 | 1065 |
| 1 | 1 | 13.20 | 1.78 | 2.14 | 11.2 | 100 | 2.65 | 2.76 | 0.26 | 1.28 | 4.38 | 1.05 | 3.40 | 1050 |
| 2 | 1 | 13.16 | 2.36 | 2.67 | 18.6 | 101 | 2.80 | 3.24 | 0.30 | 2.81 | 5.68 | 1.03 | 3.17 | 1185 |
| 3 | 1 | 14.37 | 1.95 | 2.50 | 16.8 | 113 | 3.85 | 3.49 | 0.24 | 2.18 | 7.80 | 0.86 | 3.45 | 1480 |
| 4 | 1 | 13.24 | 2.59 | 2.87 | 21.0 | 118 | 2.80 | 2.69 | 0.39 | 1.82 | 4.32 | 1.04 | 2.93 | 735 |
70%는 훈련 세트로 30%는 테스트 세트로 나눕니다.
In [4]:
from sklearn.model_selection import train_test_split
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)
데이터를 표준화합니다.
In [5]:
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)
노트
X_test_std = sc.fit_transform(X_test) 대신에 X_test_std = sc.transform(X_test)를 사용했습니다. 이 경우에 테스트 데이터셋의 평균과 표준편차가 훈련 데이터셋과 매우 비슷하기 때문에 큰 차이가 없습니다. 하지만 3장에서 보았듯이 데이터를 변환할 때 훈련 데이터셋에서 학습한 파라미터를 재사용하는 것이 올바른 방법입니다. 테스트 데이터셋은 "새로운 본 적 없는" 데이터를 의미하기 때문입니다.
초기에 fit_transform(X_test)를 사용했는데 이것은 모델 훈련에서 얻은 파라미터를 재사용하여 새로운 데이터를 표준화하지 않는 일반적인 실수입니다. 왜 이것이 문제가 되는지 간단한 예를 살펴 보겠습니다.
훈련 데이터셋에 1개의 특성("길이")을 가진 샘플 3개가 들어 있다고 가정해 보죠:
- train_1: 10 cm -> class_2
- train_2: 20 cm -> class_2
- train_3: 30 cm -> class_1
mean: 20, std.: 8.2
표준화를 한 후에 변환된 특성 값은 다음과 같습니다:
- train_std_1: -1.22 -> class_2
- train_std_2: 0 -> class_2
- train_std_3: 1.22 -> class_1
그다음 표준화된 길이가 0.6보다 작은 샘플을 class_2로 분류한다고 가정해 보죠(그 외에는 class_1). 지금까지는 좋습니다. 이제 레이블이 없는 3개의 포인트를 분류한다고 가정해 보죠:
- new_4: 5 cm -> class ?
- new_5: 6 cm -> class ?
- new_6: 7 cm -> class ?
훈련 데이터셋에 있는 표준화되기 전의 "길이" 값과 비교해 보면 직관적으로 이 샘플들은 class_2로 보입니다. 하지만 훈련 데이터셋에서 했던 것처럼 평균과 표준편차를 다시 계산하여 표준화하면 아마도 분류기가 샘플 4번과 5번만 class_2로 분류할 것입니다.
- new_std_4: -1.22 -> class 2
- new_std_5: 0 -> class 2
- new_std_6: 1.22 -> class 1
하지만 훈련 데이터셋의 표준화에 사용했던 파라미터를 사용하면 다음과 같은 값을 얻습니다:
- example5: -1.84 -> class 2
- example6: -1.71 -> class 2
- example7: -1.59 -> class 2
5 cm, 6 cm, 7 cm는 훈련 데이터셋에 있는 어떤 것보다도 작습니다. 따라서 훈련 데이터셋을 표준화한 값보다도 훨씬 작은 값으로 표준화되어야 합니다.
공분산 행렬의 고윳값 분해
In [6]:
import numpy as np
cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print('\n고윳값 \n%s' % eigen_vals)
고윳값 [4.84274532 2.41602459 1.54845825 0.96120438 0.84166161 0.6620634 0.51828472 0.34650377 0.3131368 0.10754642 0.21357215 0.15362835 0.1808613 ]
노트:
위에서 numpy.linalg.eig 함수를 사용해 대칭 공분산 행렬을 고윳값과 고유벡터로 분해했습니다.
>>> eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
이것이 잘못된 것은 아니지만 최적은 아닙니다. 에르미트(Hermetian) 행렬를 위해서 설계된 numpy.linalg.eigh를 사용하는 것이 더 좋습니다. 이 함수는 항상 실수 고윳값을 반환합니다. 수치적으로 약간 덜 안정적인 np.linalg.eig는 비대칭 정방행렬을 분해할 수 있지만 어떤 경우에 복소수 고윳값을 반환할 수 있습니다.
5.1.3 총분산과 설명된 분산¶
In [7]:
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)
In [8]:
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.savefig('images/05_02.png', dpi=300)
plt.show()
5.1.4 특성 변환¶
In [9]:
# (고윳값, 고유벡터) 튜플의 리스트를 만듭니다
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)
In [10]:
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
eigen_pairs[1][1][:, np.newaxis]))
print('투영 행렬 W:\n', w)
투영 행렬 W: [[-0.13724218 0.50303478] [ 0.24724326 0.16487119] [-0.02545159 0.24456476] [ 0.20694508 -0.11352904] [-0.15436582 0.28974518] [-0.39376952 0.05080104] [-0.41735106 -0.02287338] [ 0.30572896 0.09048885] [-0.30668347 0.00835233] [ 0.07554066 0.54977581] [-0.32613263 -0.20716433] [-0.36861022 -0.24902536] [-0.29669651 0.38022942]]
노트:
사용하는 Numpy와 LAPACK 버전에 따라 행렬 W의 부호가 바뀔 수 있습니다. 이는 문제가 아닙니다. $v$가 행렬 $\Sigma$의 고유벡터라면 다음을 얻을 수 있습니다.
$$\Sigma v = \lambda v,$$여기에서 $\lambda$는 고윳값입니다.
$$\Sigma \cdot (-v) = -\Sigma v = -\lambda v = \lambda \cdot (-v).$$이기 때문에 $-v$도 동일한 고윳값을 가진 고유벡터입니다.
In [11]:
X_train_std[0].dot(w)
Out[11]:
array([2.38299011, 0.45458499])
In [12]:
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.savefig('images/05_03.png', dpi=300)
plt.show()
5.1.5 사이킷런의 주성분 분석¶
노트
이어지는 네 개의 셀은 책에 없는 내용입니다. 사이킷런에서 앞의 PCA 구현 결과를 재현하기 위해 추가했습니다:
In [13]:
from sklearn.decomposition import PCA
pca = PCA()
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_
Out[13]:
array([0.36951469, 0.18434927, 0.11815159, 0.07334252, 0.06422108,
0.05051724, 0.03954654, 0.02643918, 0.02389319, 0.01629614,
0.01380021, 0.01172226, 0.00820609])
In [14]:
plt.bar(range(1, 14), pca.explained_variance_ratio_, alpha=0.5, align='center')
plt.step(range(1, 14), np.cumsum(pca.explained_variance_ratio_), where='mid')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.show()
In [15]:
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
In [16]:
plt.scatter(X_train_pca[:, 0], X_train_pca[:, 1])
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.show()
In [17]:
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.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,
color=cmap(idx),
edgecolor=None if idx==1 else 'black',
marker=markers[idx],
label=cl)
처음 두 개의 주성분을 사용하여 로지스틱 회귀 분류기를 훈련합니다.
In [18]:
from sklearn.linear_model import LogisticRegression
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
lr = LogisticRegression(random_state=1)
lr = lr.fit(X_train_pca, y_train)
In [19]:
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.savefig('images/05_04.png', dpi=300)
plt.show()
In [20]:
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.savefig('images/05_05.png', dpi=300)
plt.show()
In [21]:
pca = PCA(n_components=None)
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_
Out[21]:
array([0.36951469, 0.18434927, 0.11815159, 0.07334252, 0.06422108,
0.05051724, 0.03954654, 0.02643918, 0.02389319, 0.01629614,
0.01380021, 0.01172226, 0.00820609])
n_components에 (0, 1) 사이 실수를 입력하면 설명된 분산의 비율을 나타냅니다. 이 비율을 달성하기 위해 필요한 주성분 개수를 선택합니다.
In [22]:
pca = PCA(n_components=0.95)
pca.fit(X_train_std)
print('주성분 개수:', pca.n_components_)
print('설명된 분산 비율:', np.sum(pca.explained_variance_ratio_))
주성분 개수: 10 설명된 분산 비율: 0.966271440655874
n_components='mle'로 지정하면 토마스 민카(Thomas Minka)가 제안한 차원 선택 방식을 사용합니다(Minka, T. P. “Automatic choice of dimensionality for PCA”. In NIPS, pp. 598-604).
In [23]:
pca = PCA(n_components='mle')
pca.fit(X_train_std)
print('주성분 개수:', pca.n_components_)
print('설명된 분산 비율:', np.sum(pca.explained_variance_ratio_))
주성분 개수: 9 설명된 분산 비율: 0.949975302918623
PCA의 가장 큰 제약 사항 중 하나는 배치로만 실행되기 때문에 대용량 데이터셋을 처리하려면 많은 메모리가 필요합니다. IncrementalPCA를 사용하면 데이터셋의 일부를 사용하여 반복적으로 훈련할 수 있습니다.
partial_fit() 메서드는 네트워크나 로컬 파일 시스템으로부터 조금씩 데이터를 받아와 훈련할 수 있습니다. fit() 메서드는 numpy.memmap을 사용하여 로컬 파일로부터 데이터를 조금씩 읽어 올 수 있습니다. 한 번에 읽어 올 데이터 크기는 IncrementalPCA 클래스의 batch_size로 지정합니다. 기본값은 특성 개수의 5배입니다.
IncrementalPCA의 n_components 매개변수는 정수 값만 입력할 수 있습니다. 다음은 partial_fit() 메서드를 사용하여 앞의 PCA로 찾은 주성분의 결과와 비교하는 간단한 예입니다.
In [24]:
from sklearn.decomposition import IncrementalPCA
ipca = IncrementalPCA(n_components=9)
for batch in range(len(X_train_std)//25+1):
X_batch = X_train_std[batch*25:(batch+1)*25]
ipca.partial_fit(X_batch)
print('주성분 개수:', ipca.n_components_)
print('설명된 분산 비율:', np.sum(ipca.explained_variance_ratio_))
주성분 개수: 9 설명된 분산 비율: 0.9478392700446645
5.2 선형 판별 분석을 통한 지도방식의 데이터 압축¶
5.2.1 주성분 분석 vs 선형 판별 분석¶
In [25]:
Image(url='https://git.io/Jtsv8', width=400)
Out[25]:
5.2.2 선형 판별 분석의 내부 동작 방식¶
- 표준화 전처리
- 클래스별 평균 벡터
- 클래스 간 산포 행렬 $\boldsymbol S_B$, 클래스 내 산포 행렬 $\boldsymbol S_W$
- $\boldsymbol S_W^{-1}\boldsymbol S_B$ 행렬의 고윳값
- 고윳값을 내림차순 정렬
- 고윳값이 가장 큰 k개의 고유 벡터 선택
- 고유 벡터로 만든 변환 행렬로 데이터셋 투영
5.2.3 산포 행렬 계산¶
각 클래스에 대한 평균 벡터를 계산합니다:
In [26]:
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))
print('MV %s: %s\n' % (label, mean_vecs[label - 1]))
MV 1: [ 0.9066 -0.3497 0.3201 -0.7189 0.5056 0.8807 0.9589 -0.5516 0.5416 0.2338 0.5897 0.6563 1.2075] MV 2: [-0.8749 -0.2848 -0.3735 0.3157 -0.3848 -0.0433 0.0635 -0.0946 0.0703 -0.8286 0.3144 0.3608 -0.7253] MV 3: [ 0.1992 0.866 0.1682 0.4148 -0.0451 -1.0286 -1.2876 0.8287 -0.7795 0.9649 -1.209 -1.3622 -0.4013]
클래스 내 산포 행렬을 계산합니다:
$\boldsymbol S_W=\sum_{i=1}^c \boldsymbol S_i$
$\boldsymbol S_i=\sum_{x\in D_i}(\boldsymbol x-\boldsymbol m_i)^T(\boldsymbol x-\boldsymbol m_i)$
In [27]:
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 # 클래스 산포 행렬을 더합니다
print('클래스 내의 산포 행렬: %sx%s' % (S_W.shape[0], S_W.shape[1]))
클래스 내의 산포 행렬: 13x13
클래스가 균일하게 분포되어 있지 않기 때문에 공분산 행렬을 사용하는 것이 더 낫습니다:
In [28]:
print('클래스 레이블 분포: %s'
% np.bincount(y_train)[1:])
클래스 레이블 분포: [41 50 33]
$\sum_i=\dfrac{1}{n_i}\boldsymbol S_i=\dfrac{1}{n_i}\sum_{x\in D_i}(\boldsymbol x-\boldsymbol m_i)^T(\boldsymbol x-\boldsymbol m_i)$
In [29]:
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)
S_W += class_scatter
print('스케일 조정된 클래스 내의 산포 행렬: %sx%s' %
(S_W.shape[0], S_W.shape[1]))
스케일 조정된 클래스 내의 산포 행렬: 13x13
클래스 간 산포 행렬을 계산합니다:
$\boldsymbol S_B=\sum_{i=1}^c n_i(\boldsymbol m_i-\boldsymbol m)^T(\boldsymbol m_i-\boldsymbol m)$
In [30]:
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_std[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)
print('클래스 간의 산포 행렬: %sx%s' % (S_B.shape[0], S_B.shape[1]))
클래스 간의 산포 행렬: 13x13
5.2.4 새로운 특성 부분 공간을 위해 선형 판별 벡터 선택하기¶
행렬 $\boldsymbol S_W^{-1}\boldsymbol S_B$의 일반적인 고윳값 분해 문제를 풉니다:
In [31]:
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
고윳값의 역순으로 고유 벡터를 정렬합니다(판별 벡터의 개수 = 클래스 개수 - 1):
In [32]:
# (고윳값, 고유벡터) 튜플의 리스트를 만듭니다.
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)
# 고윳값의 역순으로 올바르게 정렬되었는지 확인합니다.
print('내림차순의 고윳값:\n')
for eigen_val in eigen_pairs:
print(eigen_val[0])
내림차순의 고윳값: 349.61780890599397 172.7615221897938 3.342838214841367e-14 2.842170943040401e-14 2.5545786180111422e-14 1.7533939180734234e-14 1.7533939180734234e-14 1.6579193995960903e-14 1.6579193995960903e-14 8.242524002707225e-15 8.242524002707225e-15 6.36835506006027e-15 2.974634375545734e-15
In [33]:
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.savefig('images/05_07.png', dpi=300)
plt.show()
변환 행렬 $\boldsymbol W$:
In [34]:
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,
eigen_pairs[1][1][:, np.newaxis].real))
print('행렬 W:\n', w)
행렬 W: [[-0.1481 -0.4092] [ 0.0908 -0.1577] [-0.0168 -0.3537] [ 0.1484 0.3223] [-0.0163 -0.0817] [ 0.1913 0.0842] [-0.7338 0.2823] [-0.075 -0.0102] [ 0.0018 0.0907] [ 0.294 -0.2152] [-0.0328 0.2747] [-0.3547 -0.0124] [-0.3915 -0.5958]]
5.2.5 새로운 특성 공간으로 샘플 투영하기¶
$\boldsymbol X'=\boldsymbol X \boldsymbol W$
In [35]:
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.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower right')
plt.tight_layout()
# plt.savefig('images/05_08.png', dpi=300)
plt.show()
5.2.6 사이킷런의 LDA¶
In [36]:
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
lda = LDA(n_components=2)
X_train_lda = lda.fit_transform(X_train_std, y_train)
In [37]:
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(random_state=1)
lr = lr.fit(X_train_lda, y_train)
plot_decision_regions(X_train_lda, y_train, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('images/05_09.png', dpi=300)
plt.show()
In [38]:
X_test_lda = lda.transform(X_test_std)
plot_decision_regions(X_test_lda, y_test, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('images/05_10.png', dpi=300)
plt.show()
사이킷런의 LDA 구현 방식
In [39]:
y_uniq, y_count = np.unique(y_train, return_counts=True)
priors = y_count / X_train_std.shape[0]
priors
Out[39]:
array([0.3306, 0.4032, 0.2661])
$m = \sum_{i=1}^c \frac{n_i}{n} m_i$
$S_W = \sum_{i=1}^c \frac{n_i}{n} S_i = \sum_{i=1}^c \frac{n_i}{n} \Sigma_i$
In [40]:
s_w = np.zeros((X_train_std.shape[1], X_train_std.shape[1]))
for i, label in enumerate(y_uniq):
# 1/(n-1)이 아니라 1/n로 나눈 공분산 행렬을 얻기 위해 bias=True로 지정합니다.
s_w += priors[i] * np.cov(X_train_std[y_train == label].T, bias=True)
$ S_B = \sum_{i=1}^{c}\frac{n_i}{n}(m_i-m)(m_i-m)^T $
In [41]:
s_b = np.zeros((X_train_std.shape[1], X_train_std.shape[1]))
for i, mean_vec in enumerate(mean_vecs):
n = X_train_std[y_train == i + 1].shape[0]
mean_vec = mean_vec.reshape(-1, 1)
s_b += priors[i] * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)
In [42]:
import scipy
ei_val, ei_vec = scipy.linalg.eigh(s_b, s_w)
ei_vec = ei_vec[:, np.argsort(ei_val)[::-1]]
In [43]:
lda_eigen = LDA(solver='eigen')
lda_eigen.fit(X_train_std, y_train)
Out[43]:
LinearDiscriminantAnalysis(solver='eigen')In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
LinearDiscriminantAnalysis(solver='eigen')
In [44]:
# 클래스 내의 산포 행렬은 covariance_ 속성에 저장되어 있습니다.
np.allclose(s_w, lda_eigen.covariance_)
Out[44]:
True
클래스 간 산포 행렬은 총 산포 행렬에서 클래스 내 산포 행렬을 빼서 구할 수 있다.
$\boldsymbol S_B = \boldsymbol S_T - \boldsymbol S_W$
In [45]:
Sb = np.cov(X_train_std.T, bias=True) - lda_eigen.covariance_
np.allclose(Sb, s_b)
Out[45]:
True
In [46]:
# 고유 벡터는 scalings_ 속성에 저장되어 있습니다.
np.allclose(lda_eigen.scalings_[:, :2], ei_vec[:, :2])
Out[46]:
True
In [47]:
np.allclose(lda_eigen.transform(X_test_std), np.dot(X_test_std, ei_vec[:, :2]))
Out[47]:
True
5.3 커널 PCA를 사용하여 비선형 매핑하기¶
In [48]:
Image(url='https://git.io/JtsvB', width=500)
Out[48]:
5.3.2 파이썬으로 커널 PCA 구현하기¶
In [49]:
from scipy.spatial.distance import pdist, squareform
from scipy.linalg import eigh
import numpy as np
from distutils.version import LooseVersion as Version
from scipy import __version__ as scipy_version
# scipy 2.0.0에서 삭제될 예정이므로 대신 numpy.exp를 사용합니다.
if scipy_version >= Version('1.4.1'):
from numpy import exp
else:
from scipy import exp
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]
투영된 데이터셋
"""
# 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
<ipython-input-49-ccbb9cb0be3c>:9: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
if scipy_version >= Version('1.4.1'):
예제 1: 반달 모양 구분하기¶
In [50]:
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.savefig('images/05_12.png', dpi=300)
plt.show()
In [51]:
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.savefig('images/05_13.png', dpi=300)
plt.show()
In [52]:
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.savefig('images/05_14.png', dpi=300)
plt.show()
예제 2: 동심원 분리하기¶
In [53]:
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.savefig('images/05_15.png', dpi=300)
plt.show()
In [54]:
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.savefig('images/05_16.png', dpi=300)
plt.show()
In [55]:
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.savefig('images/05_17.png', dpi=300)
plt.show()
5.3.3 새로운 데이터 포인트 투영하기¶
In [56]:
from scipy.spatial.distance import pdist, squareform
from numpy 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
반환할 주성분 개수
Returns
------------
alphas: {넘파이 ndarray}, shape = [n_samples, k_features]
투영된 데이터셋
lambdas: list
고윳값
"""
# 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
In [57]:
X, y = make_moons(n_samples=100, random_state=123)
alphas, lambdas = rbf_kernel_pca(X, gamma=15, n_components=1)
In [58]:
x_new = X[25]
x_new
Out[58]:
array([1.8713, 0.0093])
In [59]:
x_proj = alphas[25] # 원본 투영
x_proj
Out[59]:
array([-0.0788])
In [60]:
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
Out[60]:
array([-0.0788])
In [61]:
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.yticks([], [])
plt.legend(scatterpoints=1)
plt.tight_layout()
# plt.savefig('images/05_18.png', dpi=300)
plt.show()
5.3.4 사이킷런의 커널 PCA¶
In [62]:
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.savefig('images/05_19.png', dpi=300)
plt.show()
사이킷런의 매니폴드 알고리즘을 반달 모양 데이터셋과 동심원 데이터셋에 적용해 보겠습니다. 먼저 변환된 2차원 데이터셋을 그래프로 그리기 위한 간단한 함수를 정의합니다.
In [63]:
def plot_manifold(X, y, savefig_name):
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.xlabel('PC1')
plt.ylabel('PC2')
plt.tight_layout()
# plt.savefig(savefig_name, dpi=300)
plt.show()
지역 선형 임베딩(Locally Linear Embedding)은 이웃한 샘플 간의 거리를 유지하는 저차원 투영을 찾습니다. 지역 선형 임베딩을 구현한 사이킷런의 LocallyLinearEmbedding 클래스를 앞에서 적재한 반달 모양 데이터셋에 적용해 보겠습니다.
In [64]:
from sklearn.manifold import LocallyLinearEmbedding
lle = LocallyLinearEmbedding(n_components=2, random_state=1)
X_lle = lle.fit_transform(X)
plot_manifold(X_lle, y, 'images/05_lle_moon.png')
t-SNE(t-distributed Stochastic Neighbor Embedding)는 데이터 포인트 간의 유사도를 결합 확률(joint probability)로 변환하고, 저차원과 고차원의 확률 사이에서 쿨백-라이블러(Kullback-Leibler) 발산을 최소화합니다. t-SNE는 특히 고차원 데이터셋을 시각화하는데 뛰어난 성능을 냅니다. 사이킷런에는 TSNE 클래스에 구현되어 있습니다.
In [65]:
from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, random_state=1)
X_tsne = tsne.fit_transform(X)
plot_manifold(X_tsne, y, 'images/05_tsne_moon.png')
위와 비슷한 방식으로 KernelPCA, LocallyLinearEmbedding, TSNE를 동심원 데이터셋에 적용해 보겠습니다.
In [66]:
from sklearn.datasets import make_circles
X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)
scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_skernpca = scikit_kpca.fit_transform(X)
plot_manifold(X_skernpca, y, 'images/05_kpca_circles.png')
In [67]:
from sklearn.manifold import LocallyLinearEmbedding
lle = LocallyLinearEmbedding(n_components=2, random_state=1)
X_lle = lle.fit_transform(X)
plot_manifold(X_lle, y, 'images/05_lle_circles.png')
In [68]:
from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, random_state=1)
X_tsne = tsne.fit_transform(X)
plot_manifold(X_tsne, y, 'images/05_tsne_circles.png')
