Pythonで次元圧縮する方法

こんにちは、Link-Uの町屋敷です。

今回は次元圧縮について書いていこうと思います。

データの次元数が多いとどうなるのか

次元の呪いという単語を機械学習では度々目にします。
入力するデータの次元数が多いとモデルに対して与えられる点が相対的に少なくなっていろいろ不都合が出るとか、単純に計算量が多くなってやばいといったもので、
計算が終わらないから次元圧縮するという流れになるんですが。
そもそも使用する頻度の高いSVMなどの学習機で実際どのくらい終わらなくなるのか計測したことがなかったので、
次元圧縮の話の前にまずやってみました。

今回使用するデータは生成データで2クラス分類問題を行います。

DIMENTION_MAX = 2000
SEED = 12345
np.random.seed(SEED)

n_features = 1000
label = []

#ラベルの生成
for i in range(int(n_features/2)):
    label.append(0)
    label.append(1)
    
label = np.array(label)

dim = 100
#次元数がMAXを超えるまでループ
while dim < DIMENTION_MAX:
    feature = []
    x_axis.append("{0}".format(dim))
    
    for i in range(dim):
        #データ生成用の変数を作る
        mu_a = np.random.normal(0,1)
        omega_a = np.random.rand()
        
        mu_b = np.random.normal(0,1)
        omega_b = np.random.rand()
        
        for j in range(int(n_features/2)):
            #上がラベル0、下がラベル1
            feature.append(np.random.normal(mu_a,omega_a) + np.random.normal(0,10*dim/DIMENTION_MAX))
            feature.append(np.random.normal(mu_b,omega_b) + np.random.normal(0,10*dim/DIMENTION_MAX))
            
    #特徴量を要素数*次元数の形に変換        
    feature = np.reshape(feature, (dim, n_features))
    feature = feature.T

具体的には上コードで正規分布を使って作ったデータです。

最初と最後の特徴量だけ取ってきてラベルごとに色を変えてプロットするとこんな感じ。

実データを使うとデータごとに違った結果にまずなりますので今回の結果はあくまで参考です。
次元数は100から始めて2000まで100ずつ増加させてそれぞれ1000個ずつデータを作ります。

次に各次元で生成した特徴量を学習データとテストデータに分割するんですが、

sklearn.model_selection.train_test_split

をそのまま使うと分割したデータのラベルの数が偏ってしまうので、今回はここから関数を拝借してきて使ってます。

それをSVM,ロジスティック回帰、ランダムフォレスト、ニューラルネットで学習して、次元数ごとに学習にかかった時間と学習器の良さを計るF1値をノートパソコンとサーバーで計測しました。

ただし、scikit-learnがGPUに対応していないので、ニューラルネットのみGPUを使っています。

結果は以下の通りで左側のy軸が計算にかかった時間、右側のy軸がF1の値でx軸は次元数です。

このように次元数が増えると時間がかかるだけでなく、精度も落ちてしまうアルゴリズムが出てきてしまいます。

特にSVMは学習に使うデータ量より多い次元を分類しようとすると精度が絶望的になります。

このような問題を回避するためには、データ量を増やす他に、次元を削減したり、データのすべての次元を使うのではなく有用だろうと思われるものを選んで学習に使う必要があります。

次元圧縮

今回は次元を圧縮する方法として、主成分分析(PCA)線形判別分析(LDA)を使います。

タスクが簡単なので線形な変換をするやつだけ。

文章トピック分類に用いるLDAではないので注意。

PCAには教師データが必要ありませんが、LDAには必要です。

LDAは教師データを使って最もよくクラスタを分類するように新しい軸を検索します。

一方でPCAには教師データが必要ありません。

PCAとLDAの違いはここがわかりやすい(英語)

PCA,LDAともにskikit-learnに実装されているので簡単に行うことが出来ます。

from sklearn.decomposition import PCA
    
    dcp = PCA(n_components=10) #n_componentsで削減する次元数を指定
    dcp_data = dcp.fit_transform(data)
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
    dcp = LinearDiscriminantAnalysis(n_components=10) #n_componentsで削減する次元数を指定
    dcp_data = dcp.fit_transform(data, label) #教師データとしてlabelを渡す

上に載せてあるデータにPCAをかけて2次元に圧縮すると結果はこんな感じになります。

LDAをかけて1次元に圧縮する(カラスが2つなので1次元になる)とこんな感じ。

PCAを使って圧縮時間も含めて処理時間とF1値を計算した結果がこちら。

比較用に圧縮してないときの結果も載せた。

F1は1が最高ですべてのテストデータを正解したことを示す。

全体的に処理速度が早くなり、F1の値が向上している。

ニューラルネットだけ時間が伸び気味なのは何故だろうか…

まとめ

今回は次元の圧縮で処理が早くなり、精度も向上することを示したが、

実データじゃない簡単タスクだからかもしれない。

isomapとかの非線形の次元圧縮方法や、t-sne、オートエンコーダも紹介していないので、

もっと難しいデータでいつか紹介したい。

使用したソースコード

#Python3系じゃないとバグります

import time
from matplotlib import pyplot as plt
import numpy as np
from sklearn.decomposition import PCA, FastICA
from sklearn.linear_model import LogisticRegression
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.mixture import GaussianMixture
from sklearn.svm import SVC
from sklearn.model_selection import StratifiedShuffleSplit
from sklearn.metrics import confusion_matrix, f1_score
from sklearn.manifold import TSNE

import keras
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten

import joblib

#Quoted from https://stackoverflow.com/questions/35472712/how-to-split-data-on-balanced-training-set-and-test-set-on-sklearn
def get_safe_balanced_split(target, trainSize=0.8, getTestIndexes=True, shuffle=False, seed=None):
    classes, counts = np.unique(target, return_counts=True)
    nPerClass = float(len(target))*float(trainSize)/float(len(classes))
    if nPerClass > np.min(counts):
        print("Insufficient data to produce a balanced training data split.")
        print("Classes found %s"%classes)
        print("Classes count %s"%counts)
        ts = float(trainSize*np.min(counts)*len(classes)) / float(len(target))
        print("trainSize is reset from %s to %s"%(trainSize, ts))
        trainSize = ts
        nPerClass = float(len(target))*float(trainSize)/float(len(classes))
    # get number of classes
    nPerClass = int(nPerClass)
    print("Data splitting on %i classes and returning %i per class"%(len(classes),nPerClass ))
    # get indexes
    trainIndexes = []
    for c in classes:
        if seed is not None:
            np.random.seed(seed)
        cIdxs = np.where(target==c)[0]
        cIdxs = np.random.choice(cIdxs, nPerClass, replace=False)
        trainIndexes.extend(cIdxs)
    # get test indexes
    testIndexes = None
    if getTestIndexes:
        testIndexes = list(set(range(len(target))) - set(trainIndexes))
    # shuffle
    if shuffle:
        np.random.shuffle(trainIndexes)
        if testIndexes is not None:
            np.random.shuffle(testIndexes)
    # return indexes
    return trainIndexes, testIndexes
#Quote end

DIMENTION_MAX = 2000
SEED = 12345
np.random.seed(SEED)

n_features = 1000
label = []

#ラベルの生成
for i in range(int(n_features/2)):
    label.append(0)
    label.append(1)
durations = []
f1_values = []
label = np.array(label)
for clf_name in ['SVM', 'LogisticRegression', 'RandomForest', 'NeuralNet']:

    x_axis = []
    for decomp_name in ['LDA']:#['LDA', 'ICA']:
        dim = 1000
        #次元数がMAXを超えるまでループ
        while dim < DIMENTION_MAX:
            feature = []
            x_axis.append("{0}".format(dim))
            
            for i in range(dim):
                #データ生成用の変数を作る
                mu_a = np.random.normal(0,1)
                omega_a = np.random.rand()
                
                mu_b = np.random.normal(0,1)
                omega_b = np.random.rand()
                
                for j in range(int(n_features/2)):
                    #上がラベル0、下がラベル1
                    feature.append(np.random.normal(mu_a,omega_a) + np.random.normal(0,10*dim/DIMENTION_MAX))
                    feature.append(np.random.normal(mu_b,omega_b) + np.random.normal(0,10*dim/DIMENTION_MAX))
            
                    
            #特徴量を要素数*次元数の形に変換        
            feature = np.reshape(feature, (dim, n_features))
            feature = feature.T

            
            #バランス良く訓練データとテストデータに分割
            train_index, test_index = get_safe_balanced_split(label, trainSize=0.8, getTestIndexes=True, shuffle=True, seed=SEED)
           
            train_feature, test_feature = feature[train_index], feature[test_index]
            train_label, test_label = label[train_index], label[test_index]
            
            #学習器の選択
            if clf_name == 'SVM':
                clf = SVC(C = 1, gamma = 0.01)
            elif clf_name == 'LogisticRegression':
                clf = LogisticRegression(C = 1e5)
            elif clf_name == 'RandomForest':
                clf = ExtraTreesClassifier(n_estimators=1000,random_state=0)
            
            start = time.time()
            
            #次元圧縮法の選択
            if decomp_name == 'PCA':
                dec = PCA(n_components = 2)
                train_feature = dec.fit_transform(train_feature)    
                test_feature  = dec.transform(test_feature)
            elif decomp_name == 'LDA':
                dec = LinearDiscriminantAnalysis(n_components=1)
                train_feature = dec.fit_transform(train_feature, train_label)
                test_feature  = dec.transform(test_feature)
            elif decomp_name == 'ICA':
                dec = FastICA(n_components = 2)
                train_feature = dec.fit_transform(train_feature)
                test_feature  = dec.transform(test_feature)
            elif decomp_name == 'TSNE':
                dec = TSNE(n_components=2, perplexity=30)
                train_feature = dec.fit_transform(feature)
            
            if clf_name == 'NeuralNet':
                clf = Sequential()
                clf.add(Dense(64, activation='relu', input_shape=(len(train_feature[0]),)))
                clf.add(Dense(1, activation='sigmoid'))
                
                clf.compile(loss='binary_crossentropy',
                            optimizer=keras.optimizers.Adadelta(),
                            metrics=['accuracy'])
                
                clf.fit(train_feature, train_label,
                          batch_size=128,
                          epochs=100,
                          verbose=0)
                predict = np.round(clf.predict(test_feature))
            else:
                clf.fit(train_feature, train_label)
                predict = clf.predict(test_feature)
            
            
            
            duration = time.time() - start
            
            print("Dim {0}, Time: {1}".format(dim, duration))
            print(confusion_matrix(test_label, predict, [0,1]))
            
            f1_values.append(f1_score(test_label, predict, [0,1]))
            durations.append(duration)
            dim += 100
            
        #データの保存、保存したデータはjoblib.loadでロードできる 
        joblib.dump(durations, '{0}_{1}_Note_durations.pkl'.format(clf_name, decomp_name))
        joblib.dump(f1_values, '{0}_{1}_Note_f1values.pkl'.format(clf_name, decomp_name))