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))

GPUを使って無線LANをクラックする話:Pyritの考古学/倫理学

拝啓。

自販機にたまにホット飲料を見かける季節となりましたが皆さまいかがおすごしでしょうか。

秋ですね。気まぐれと勢い 旅行で行った九州では彼岸花が咲いておりました。

今回は咲きかけギリギリのヒガンバナくんの写真です。

これぐらいのをお散歩中に見つけたら、だいたい翌日か翌々日ぐらいには咲きますので、翌日もぜひ見に行ってあげてくだせぇ。ヒガンバナくんは、写真で見るより実物のほうが絶対いいです。雨が降ってない事を祈ります。

あーでも。毒性があるので持ち帰って食べたりしちゃダメですよ。持ち帰るのは思い出だけ。山と一緒かな。

前回までのPyrit

じゃーまた、WiFiのパスワードでも、クラックしてきますか。

前回わかった事は、最初取ったベンチマークで「GPUはCPUより速くてすごーい!!」とか喜んでたけど、実はGPUは10%も使っていないし、CPUも全部は使いきってはいなかった…というなかなか衝撃的な事実でした。さらに言えば、そのベンチはGPUの分と称してCPUも使っている、公平かどうかは正直怪しいと言わざるを得ないものでありました。

言うなれば我々は…ただの幻を見て一喜一憂していたのです。

どうしてこうなった(AA略)

うーん。悲しい。

実際に何か手を動かす前に、どうしてこんなことになってしまったのか想像してみましょう。Pyritのクラスターを組んでいたみなさんは、実はほとんど寝てるだけのGPUを前に、ドヤっていたのか。まさか、そんな。

Pyritの考古学

Pyritが活発に開発されていたのは2009〜2010年ごろです。DeepLearningブームが始まったのが2012年とかなのでそれよりも前だし、その頃に生まれた子供は、今小学校で九九を覚えてる頃。それぐらい大昔です。そんな古代のソフトウェアを今動かそうというのですから、歴史を踏まえなければ、目の前の状況を理解することはできないでしょう。

まず、GPUをなんでこんなに遊ばせてしまっているのか。

ひょっとすると、前回使ったVisual Profilerが実はその頃なかったのかもしれません。だから使い切ってない事に誰も気づかなかったとか。が、しかし、NVIDIAの公式ページによると2008年から作ってるとの事なので(頑張ってるなNVIDIA)、この線は無い…と思いたい。つまり、「昔はGPUを全部使い切ってたけど、時間が経ったら状況が変わった」可能性が高い。

あの頃のGPU:GTX480が出た頃(買えたとは言ってない)

2010年といえば、GPUならGeForce GTX 295(2009年6月)とか、GTX 480(2010年3月)の頃です。…もう何も覚えてない…。まーとりあえず表でも書いてみますか。CUDAとは直接は関係ないですが、参考程度にPassMarkのベンチマークの数値も貼っときます。グラフィクスを描く時もCUDAコアが走るそうなので、無関係ではあるまい。

GPU発売時期CUDAコアクロックCUDAコア数PassMarkベンチ
GTX 2952009/061242MHz240×21052
GTX 4802010/031401MHz4804358
GTX 1080 Ti2017/021480MHz358414057

GTX480ではGPGPUへの最適化も進んだそうで、一気に性能が上がっていますが、曰く

最初のロットの歩留まりは、嘘か真かは不明だが2%ほどだったという。その後もステッピングを重ねてできる範囲での改良を施したものの、最終的に2010年3月のGeForce GTX 480発表時点で出荷可能な枚数は全世界を合わせて数千枚程でしかなかった。

ロードマップでわかる!当世プロセッサー事情 ― 第146回
GPU黒歴史 DX11への遅れが生んだ駄作 GeForce GTX 480

だそうで、2010年のゲーマーたちはだいたいGTX295を使っていたんじゃないかと思います。つまり、今の1080 Tiは当時に比べて10倍以上、逆に言えば、GTX295は今に比べたら1/10以下の文字通り「桁違い」の性能だったと思われます。…ところで、10%弱しか走ってなかったんですよね?

あの頃のCPU: Core i7が初めて出た頃

じゃあ一方CPUはどうなんだ。2010年にPCパーツショップに並んでいたのは、Core iシリーズの最初、Nehalem世代です。これももう覚えてないですねぇ。せいぜい登場時に「3と5と7って何?Core i7には7コア載ってるのか?」と思ったのを覚えている程度です(答え:載ってない)。

同じように表を書いてみましょう。

GPU発売時期クロック(ターボ)コア数(スレッド数)PassMarkベンチ(マルチ/シングル)
Core i7 9602009/103.2GHz(3.46GHz)4(8)5828/1389
Core i7 8700K2017/103.7GHz(4.7GHz)6(12)15969/2704
Xeon Silver
4116
2017/022.1GHz(3.0GHz)12(24)15131/1599

ターボクロックはカッコの中に一応書いてはいますが、ターボブーストは「一時的に加速できる」という機能ですので、ずっと負荷を掛け続けるベンチマークではベースクロックの方が重要です。

ベンチに使ってるサーバに乗ってるのが最後のXeon Silver 4116です。PassMarkのベンチマーク的にはなんとシングルコア性能自体は2009年のCPUである960と大して変わらないという結果になってます。もちろん周波数は960よりかなり低いので、「周波数あたりの性能は上がってる」と言えるのは間違いなく事実なんですけどね。

最近のサーバ用CPUトレンドはこのようにそこまでの性能ではない低い 周波数のコアをたくさん束ねて全体としての性能を確保する戦略らしく、結果としてCPU全体では10年前のCPUの三倍ぐらいの性能は出るようになっています。という所は観察できました。

時を経ても変わらないもの―Pythonはシングルスレッド

GPUもCPUも速くなってはいるけれど、一つだけ変わらない事があります。

それは、Pythonはシングルスレッドでしか動かないという事実です。PythonにはGIL(Global Interpreter Lock)という仕組みがあるので、いくらthreading.Thread()してスレッドを作っても、同時に動くPythonのスレッドはたった1つだけです。100コア積んでも5000兆コア積んでも1つだけ。

ちょっと待て。じゃあなんでPythonで書かれてるはずのPyritはマルチコア対応してるんだ?もっというとDeepLearningはみんなPythonでやってるんじゃ無いのか?シングルスレッドなのか?

もちろんそんな事あるわけなくて、PythonからC言語の拡張を呼んだ時は、C言語で処理しつつPythonの他のスレッドを起こす事ができます。その時はPy_BEGIN_ALLOW_THREADS;Py_END_ALLOW_THREADS;で囲めばマクロがよしなにスレッドの面倒を見てくれます。

実際にPyritでCPUでクラックする時もこんな感じになってます:

    Py_DECREF(passwd_seq);

    if (arraysize > 0)
    {
        Py_BEGIN_ALLOW_THREADS;
        i = 0;
        do
            i += finalize_pmk(&pmk_buffer[i]);
        while (i < arraysize);
        Py_END_ALLOW_THREADS;

        result = PyTuple_New(arraysize);
        for (i = 0; i < arraysize; i++)
            PyTuple_SetItem(result, i, PyString_FromStringAndSize((char*)pmk_buffer[i].e1, 32));
    } else {
        result = PyTuple_New(0);
    }

    PyMem_Free(pmk_buffer);

    return result;
}

これを使う事で1つ以上のCPUコアを使う事が可能になっていますが、しかし、他の部分は全部シングルスレッドでしか動かないことは記しておく価値があると思います。C言語のコードでも、PythonのオブジェクトからCの値に変換したり、その逆をしたりするところはシングルスレッドでやらないといけません。

CPUを全コア使いきれていない理由は、前回書いた通り、ここだと思います。パスワードを生成する部分がPythonで書かれているので、そこがボトルネックになっているんだと思います。

Pyritの倫理学:「公平なベンチマーク」ってなんじゃらほい

まぁ、だいたい状況は掴めました。次にルールを考えましょう。

「ベンチマーク」という「勝負」のルールを、です。

現在のPyritのbenchmarkコマンドのルールでは、「CPU」と「GPU」の勝負のはずが、「GPUのために準備するCPUのコード」も同時に別々のCPUのスレッドで走っています。

CPUのスレッドはランダムなタイミングで割り込まれる可能性がありますから、CPUがクラックしているコード上の「行間」でGPUの為に準備するスレッドのコードが実行されて、CPUの性能が見かけ上悪くなってる可能性はかなり高いです。わたしはこんなんじゃあ、公平なベンチマークとは到底言えない思います。

「公平」とはなんじゃらほい??

ほう、「公平」ねぇ。ところで、「公平」ってなんでしょう?試しに「公平な社会」とか書いて見ると、これがいかにも政治家さんや官僚さんの口から出てきそうな単語な事からも分かるように、「公平」とは人によって意見が別れうる、一筋縄ではいかない手強い単語なわけです。

まぁ、「公平な社会」が何なのかは置いておいて(わたしにはわからん)、今回は「公平なPyritのベンチマーク」とは何なのかのわたしなりの意見を書いておきます。

CPUとGPUは明らかに対等な立場ではありません。というのも、CPUだけでパスワードのクラックはできるけれど、GPUだけではできないからです。CPUがGPUの処理しやすい問題の形式に整えてから、GPUが一気に処理をする。そういう、協力関係で計算を進めるのがGPGPUです。 こういうのを一般的かつかっちょいい言葉では「ヘテロジニアス・コンピューティング」といいます。

そういう意味では、同じ勝負でも「サッカー」とはだいぶ雰囲気が違います。サッカーみたいに、お互いが対称なルール、例えば「GPUめ、CPUに頼らずクラックしてみろ。CPUもお前には頼らないから平等だよな?」みたいなルールでは、単にGPUが何もできなくて降参して終わりです。

そもそも勝負のルールってどうして必要なんだ?

サッカーという勝負のルールは、選手や観客が楽しく、かつ安全に遊ぶためのルールですが(たぶん)、じゃあ、Pyritのベンチマークという勝負のルールは何の為のルールでしょう?

まぁ結局は同じ、「面白いかどうか」が基準でいいと思うんですが、それをもうちょっと具体的に言えば、「GPUとCPU同士の本気の勝負が見たい」に尽きると思います。もうちょっと大人じみた言葉で言えば、「GPUの得意な処理をCPUからGPUにオフロードした時に、そこがCPUと比べてどれぐらい速くなるか?」の限界が見たい。

もちろん、一部だけ高速化できても全体で速くなるとは限りません。全体の10%がGPUで10倍速くなっても、1.1倍ぐらいしか速くなりません(アムダールの法則)。例えば、GPUがCPUに比べてそこまで得意ではない処理もGPUに持って言って「並列化」する事で、その部分はそこまで高速化しなくとも、システム全体としては高速化する可能性はあります。…が、個人的にそこはあんまり面白くないなぁ…GPUの最高性能、つまり「ガチ」が見たいなぁ…と思ったので独断により今後はこの上に書いた方針で行きたいと思います。Pyrit全体よりも局所的な情報の方が、他のシステムについて見積もる時に便利そうだし。まぁこのいかにもな理由の方は後付けなんですが。

流石にこれだけではGPUに甘々な感じがするので、今邪魔されているCPUにも救いの手を差し伸べましょう。CPUが邪魔されている原因は、結局CPUとGPUを同時に走らせてクラックしていることです。これは単にベンチマークをとる時に、別々に走らせれば邪魔されずにすみます。

ベンチマークのルールは与えられるものじゃない、考えるものだ

というわけで考えた「面白い勝負」のためのルールをまとめると:

  • GPUの一番得意な処理をCPUからGPUにオフロードする
    • そこが全く同じ処理をするCPUと比べ、どれぐらい速くなるかを測定
  • CPUとGPUは別々に走らせて、お互いに邪魔しあう事を防ぐ

そろそろ長くなってきたので、もっと具体的な話は次回で。

本当に今のルールは「公平」ではないのか?

一応擁護しておくと、今のPyritのベンチマークのルールも「公平」だと言えなくもありません。

というのも、今のPyritのベンチマークは会社でよくある「成果主義」そのものだ、と捉えることも出来るからです。会社というのも協力関係がお互いにあるわけですが、「成果主義」の会社ではお給料を決めたりする時はその上で「そいつが、会社の利益に、どれくらい貢献したか?」みたいなのを(定量的にじゃないにしても)測ります。Pyritの今のベンチマークは、それと似たような事を測っている、と捉えることは可能だと思います。

で、この「成果主義」をどう思うか、具体的にどうするか、色々意見があることと、上でわたしがベンチマークのルールをどうするか悩んでいるのと、対応がつくわけです。もちろん、人とコンピュータの世界では前提が色々違うので全く同じというわけではありませんが、「一筋縄ではいかないねぇ〜」「意見が別れるだろうね〜」という点では全く同じです。

まとめ

Pyritが想像以上に大昔のソフトウェアだったので、歴史を調べました。

Pyritが盛んに開発されていた頃のGPUはCUDAへの最適化が進む一歩手前で今とは桁違いの性能差がありそうなこと。

シングルコア性能は今のサーバのXeonは昔のハイエンドCore i7と同じぐらいの性能だけれどマルチコアにすることで性能を稼いでいること、Pythonは相変わらずシングルスレッドなこと。そんなことがわかりました。

そして、わたしの興味に応じた「ぼくのかんがえたベンチマークのルール」も考えました。意外と人の仔の世界とも通じるものが有るんじゃのぅ。

さ〜て次回のPyritは?

このルールに応じてPyritを書き換えて、測定しなおしましょう!

山の上の石に日が暮れるまで座って色々考えてると今回みたいな記事になりまして。
まぁ、たまにはこういうのをじっくり考えるのも大事なんじゃないかなぁ、って最近思うんですよ。