第5章：次元削減データを圧縮する

今回は第5章の内容を解説します。
全ての内容は書いてないので、本気で勉強したい方は下記を購入してもらえればなと。

[第3版]Python機械学習プログラミング 達人データサイエンティストによる理論と実践 (impress top gear)

前章では特徴量選択の様々な手法を用いてデータセットの次元を削減する方法を取り上げました。特徴量選択に変わる次元削減の1つの方法は「特徴量抽出」です。本章ではデータセットの情報を要約するのに役立つ基本的な手法の3つを解説します。

・教師なしデータ圧縮での主成分分析（PCA）
・クラスの分離を最大化する次元削減方としての線型判別分析（LDA）
・カーネル主成分分析（KPCA）による非線形次元削減

5.1 主成分分析による教師なし次元削減

特徴量選択と同様に、特徴量抽出を利用すればデータセット内の特徴量の個数を減らす事ができます。逐次後退選択（SBS）といった特徴量選択アルゴリズムを使った時には、元の特徴量は変換される事なくそのままの形で維持されてました。これに対し特徴量抽出ではデータが新しい特徴量空間に変換または射影される。次元削減ではデータに含まれる情報の大部分を維持することを目標としたデータ圧縮の手段として、特徴量抽出を捉える事ができる。実際には、特徴量抽出は学習アルゴリズムの記憶域や計算黄チルの改善だけでなく、特に生息化されていないモデルを扱う時に、次元の呪いを減らすことで予測性能を向上させる目的でも利用できる。

5.1.1 主成分分析の主要なステップ

ここでは主成分分析（PCA）について説明します。PCAは様々な分野にわたって使われている教師なし線形変換法であり、最もよく用いられるタスクは特徴量抽出と次元削減です。PCAは特徴量同士の相関関係に基づいてデータからパターンを抽出するのに役立ちます。

5.1.2 主成分を抽出する

PCAの最初の4つの手順に取り組みます。

1.データを標準化する
2.共分散行列を作成する
3.共分散行列の固有値と固有ベクトルを取得する
4.固有値を光順でソートすることで固有ベクトルをランク付けする

まずWineデータセットを読み込みます。

import pandas as pd
df_wine = pd.read_csv(
    'https://archive.ics.uci.edu/ml/machine-learning-databases/wine/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()

次にWineデータセットを処理して訓練データセット（データの70%を使用）とテストデータセット（30%）に分割して、分散が1となるように標準化する。

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
# 2列目以降のデータをXに、1列目のデータをyに格納
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)

Wineデータセットの共分散行列の固有対を取得します。

import numpy as np
cov_mat = np.cov(X_train_std.T)  # 共分散行列を作成
eigen_vals, eigen_vacs = np.linalg.eig(cov_mat) # 固有値と固有ベクトルを計算
print('\nEigenvalues \n%s' % eigen_vals)

標準化されたデータセットの共分散行列の計算にはNumpyのcov関数を使いました。また、linalg.eig関数を使って固有分解を実行し、13個の固有値からなるベクトル（eigen_vals）と13*13次元の行列として格納された対応する固有ベクトルを生成してます。

5.1.3 全分散と説明分散

ここではデータセットを新しい特徴量部分空間に圧縮するという方法でデータセットの次元を削減します。そこで、データセットに含まれる大半の情報（分散）を含んでいる固有ベクトル（主成分）だけを選択します。そこで固有ベクトルの大きい順に見てみましょう。その前に、固有値の分散説明率をプロットします。

分散説明率の累積和はNumpyのcumsum関数と使って計算できます。計算した値のプロットにはmatplotlibの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()

上記のグラフから1つの主成分だけでなく40%近くを近くを占めてるのが分かります。最初の2つの主成分を合わせると分散の60%近くになります。

PCAは教師なしのデータ圧縮法であることを再認識します。クラスラベル情報は使ってません。ランダムフォレストはデータの所属情報を使ってノードの不純度を計算しますが、PCAで計算される分散は特徴量軸に沿った値の散らばりを計測します。

5.1.4 特徴量変換

共分散行列を固有対（固有ベクトルと固有値）に上手く分解できたところで、最後の3つの手順に進みWIneデータセットを新しい主成分軸に変換します。手順は次の通りです。

1.最も大きいk個の固有値に対応するk個の固有ベクトルを計算する。この場 合の新しいkは新しい特徴量部分空間の次元数を表す(k=<d)
2.上位k個の固有ベクトルから射影行列Wを作成する。
3.射影行列Wを使ってd次元の入力データセっとXを変換し、新しいk次元の特徴量部空間を取得する。

簡単に説明しますと、固有値の大きいものから順に固有対を並び替え、選択された固有ベクトルから射影行列を生成する。そして、この射影行列を使ってデータをより低い次元の部分空間に変換します。まず、固有値の大きいものから順に固有対を並び替えます。

# (固有値、固有ベクトル)のタプルのリストを作成
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vacs[:,i]) for i in range(len(eigen_vals))]
# (固有値、固有ベクトル)のタプルを大きいものから順に並び替え
eigen_pairs.sort(key=lambda k : k[0], reverse=True)

次に、最も大きい2つの固有値に対応する2つの固有ベクトルを集めます。それにより、このデータセットにおける分散の約60%を捉える事ができます。ここでは説明のために固有値を2つだけ選択してますが、このデータは後ほど第1主成分と第2主成分による2次元の散布図を使ってプロットします。実務では計算効率と分類器の性能バランスを見ながら主成分の個数を決定する必要があります。

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

このコードを実行すると、上位2つの固有ベクトルから13*2次元の射影行列Wが作成されます。

この射影行列を使うことでデータ点x（13次元の行ベクトルとして表される）をPCAの部分空間（主成分1と主成分2）に変換すると、2つの新しい特徴量からなる2次元のベクトルを生成できます。以下のコードではdot関数を使って行列の内積を求めてます。

同様に行列のない席を計算することにより124*13次元の訓練データセット全体を2つの主成分に変換できます。

X_train_pca = X_train_std.dot(w)
X_train_pca

最後に124*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.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

この散布図では説明のためにクラスラベルの情報を符号化しましたが、PCAがクラスラベルの情報を使わない教師なしのデータ圧縮法であることに注意しましょう。

5.1.5 sklearnの主成分分析

前項のアプローチは主成分分析（PCA）の内部の仕組みについて理解するのに役立ちました。ここではsklearnで実装されているPCAクラスを使う方法についてみていきます。PCAは変換器クラス（transform）の1つであり、訓練データを使ってモデルを適合させた上で、訓練データとテストデータを同じモデルパラメータに基づいて変換します。まずはロジスティック回帰を使って変換後のデータ点を分類します。決定領域をプロットするための関数を作成します。

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))
    # 各特徴量を1次元配列に変換して予測を実行
    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), edgecolors='black', marker=markers[idx], label=cl)

from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
# 主成分数を指定して、PCAインスタンスを生成
pca = PCA(n_components=2)
# ロジスティック回帰のインスタンスを生成
lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
# 次元削減
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
# 削減したデータセットでロジスティック回帰を適合
lr.fit(X_train_pca, y_train)
# 決定境界をプロット
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

5.2 線形判別分析による教師ありデータ圧縮

線形判別分析（LDA）は特徴量抽出手法の1つです。LDAは計算効率を高め、正則化されてないモデルで「次元の呪い」による過学習を抑制するために利用できます。LDAの基本的な考え方は主成分分析と似てますが、PCAがデータセットにおいて分散が最も大きい直交成分軸を見つけ出そうとするのに対し、LDAはクラスの分離を最適化する特徴量部分空間を見つけ出そうとします。LDAとPCAの類似性を詳しく取り上げ、LDAのアプローチ方法も見ていきます。

5.2.1 主成分分析と線形判別分析

LDAとPCAはどちらも線形変換法であり、データセットの次元数を減らすために利用できます。PCAが教師なしのアルゴリズムに対し、LDAは教師ありのアルゴリズムです。このため、分類タスクの特徴量抽出ほ手段としてはLDAの方が優れてます。

5.2.2 線形判別分析の仕組み

下記が大まかな流れです。

1.d次元のデータセットを標準化する（dは特徴量の個数）。
2.クラスごとにd次元の平均ベクトルを計算する。
3.平均ベクトルを使ってクラス間変動行列とクラス内変動行列を生成する。
4.行列の固有ベクトルと対応する固有値を計算する。
5.固有値を降順でソートすることで、対応する固有ベクトルをランク付けす    る。
6.d * k次元の変換行列を生成するために最も大きいk個の固有値に対応するk
        個の固有ベクトルを選択する。固有ベクトルは行列の列である。
7.変換行列を使ってデータ点を新しい特徴量部分空間に射影する。

5.2.3 変動行列を計算する

Wineデータセットの特徴量はPCAに関する最初の節ですでに標準化してます。このため、最初の手順をとばして平均ベクトルの計算に進みます。

np.set_printoptions(precision=4)
mean_vesc = []
for label in range(1, 4):
    mean_vesc.append(np.mean(X_train_std[y_train==label], axis=0))
    print('MV %s: %s\n' %(label, mean_vesc[label-1]))

d = 13
s_w = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vesc):
    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('Within-class scatter matrix: %sx%s' %(s_w.shape[0], s_w.shape[1]))

変換行列を計算する時には、訓練データセットにおいてクラスラベルが一様に分布していることが前提になります。ですが、クラスラベルの個数を出力してみるとこの前提を見たいしてないことが分かります。

print('Class label distribution : %s' % np.bincount(y_train[1:]))

スケーリングされた内変動行列を計算するコードは次のようになります。

d = 13
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vesc):
    class_scatterb = np.cov(X_train_std[y_train == label].T)
    S_W += class_scatter

print('Scaled within-class scatter matrix : %sx%s' % (S_W.shape[0], S_W.shape[1]))

スケーリングされた内変動行列（共分散行列）を計算した後は、次の手順であるクラス間変動行列の計算に進みます。

mean_overall = np.mean(X_train_std, axis=0)
d = 13 # 特徴量の個数
S_B = np.zeros((d, d))
for i, mean_vec in enumerate(mean_vesc):
    n = X_train_std[y_train == i + 1, :].shape[0]
    mean_vec = mean_vec.reshape(d, 1)
    mean_overall = mean_overall.reshape(d, 1)
    S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)
    
print('Between-class scatter matrix: %sx%s' % (S_B.shape[0], S_B.shape[1]))

5.2.4 新しい特徴量部分空間の線型判別を選択する

LDAの残りの手順はPCAの手順と似ています。ただし、共分散行列で固有分解を実行するのではなく、行列を一般化された固有値問題を解きます。

# inv関数で逆行列、dot関数で行列積、eig関数で固有値を計算
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)
print('Eigenvalues in descending order\n')
for eigen_val in eigen_pairs:
    print(eigen_val[0])

線型判別（固有ベクトル）によってクラスの判別情報がどの程度捕捉されるのか計測されるため、PCAの節で作成した分散説明のグラフと同様に、固有値を減らしながら線型判別をプロットします。

# 固有値の有数部の総和を求める
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 Discriminaants')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

グラフからわかるように、最初の2つの線型判別だけでWine訓練データセットないの有益な情報を100%補足してます。次に最も判別力のある2つの固有ベクトルを列方向に並べて、変換行列Wを作成します。

# 2つの固有ベクトルから変換行列を作成
import numpy as np
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,
                            eigen_pairs[1][1][:, np.newaxis].real))
print('Matrix W:\n', w)

5.2.6 sklearnによる線型判別分析

sklearnによるLDAクラスをを調べて、ロジスティック回帰の分類器がLDAによって変換された低次元の訓練データセットをどのように処理するか見てみよう。

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
# 次元数を指定してLDAのインスタンスを生成
ida = LinearDiscriminantAnalysis(n_components=2)
X_train_ida = ida.fit_transform(X_train_std, y_train)

lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
lr = lr.fit(X_train_ida, y_train)
plot_decision_regions(X_train_ida, y_train, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

テストデータセットでの結果を見てみます。

X_test_ida = ida.transform(X_test_std)
plot_decision_regions(X_test_ida, y_test, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

ロジスティック回帰の分類器はWineデータセットの元の13この特徴量ではなく、2次元の特徴量部分だけを使うことにより、テストデータセットのデータ点を完全に分類できることが分かります。

5.3.4 sklearnのカーネル主成分分析

線型判別できないデータに対しては「カーネル主成分分析」を実装します。

from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_moons
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)

最初の2つの主成分に変換された半月形データをプロットしてみます。

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

まとめ

今回は特徴量抽出のための基本的な次元削減方法を実装しました。標準のPCA、LDA、カーネルPCAです。「PCA」は教師なしのデータに対して適用し、「LDA」は教師ありのデータに対して適用し、クラスの分類を最大化しようとします。「カーネルPCA」は非線形の特徴量抽出が実装できます。codeで書いてみましたが、まだまだ理解は浅いです。日々勉強、というより日々作業ですね。実際のデータに対して適用し、数こなそうかなと。