• このエントリーをはてなブックマークに追加

今なら、継続入会で月額会員費が1ヶ月分無料!

記事 2件
  • 【機械学習】Chainer vs C++ AVX最終決戦

    2016-09-01 07:55  
    102pt
     さて 前回の結果に気を良くした僕は、調子に乗って5層オートエンコーダの高速化を試みるのである。 といっても、前回作ったプログラムがほぼ流用できるので、予め学習させたpklファイルをヘッダファイルに書き出すプログラムを流用して5層オートエンコーダをC++の配列として出力する。
    def dumpArray(ar):
        for i in range(len(ar)):
            e = ar[i]
            if isinstance(e,np.ndarray):
                print "{",
                dumpArray(e)
                print "}",
            else:
                print ("{}".format(e)),
            if i < len(ar)-1:
                print ","
    def dumpLayer(name,l):
        print "double ",
        print "{}_W".format(name),
        for d in l.W.data.shape:
            print "[{}]".format(d),
        print "={"
        dumpArray(l.W.data)
        print "};"
        print "double ",
        print "{}_b".format(name),
        for d in l.b.data.shape:
            print "[{}]".format(d),
        print "={"
        dumpArray(l.b.data)
        print "};"
    print "double in[]={"
    dumpArray(x_train[0])
    print "};"
    dumpLayer("l1",model.l1)
    dumpLayer("l2",model.l2)
    dumpLayer("l3",model.l3)
    dumpLayer("l4",model.l4)
    dumpLayer("l5",model.l5)


     こいつで得られたlayers.hをC++プログラムから読み込んで前回と同じようなことをするわけだ。 今回も余裕だろうと思って試すと
    $ python bench_ae_chainer.py 
    start
    Chainer:0.223642sec
    $ ./hoge 
    0.129455 sec
    31.683862 


     どっどっどっ どういうことだってばよ!? 前回は500倍くらいのスピードだったのに今回は2倍程度のスピードになってしまっている。 いろいろ考えると、まあ結論はひとつしかない。 要するに、第一層がデカすぎるのである。 前回の第一層は1,407バイト、つまり1キロバイト強しかないのに対し、今回の第一層は3,513,339バイト、つまり3.5メガバイトもある。 コンピュータの世界ではキャッシュに乗るか乗らないかというのが非常に重要だ。 最新世代のCore i7であっても、L1キャッシュは32KB、L2キャッシュでも256KB、L3キャッシュになってようやく、4〜8MBになる。 1.4キロバイトのデータは全てL1キャッシュに乗るが3.5メガバイトはL3キャッシュにやっと乗るかどうか、というところだ。 つまり前回の戦いは単にネイティブコードとPythonの戦いではなく、キャッシュの戦いでもあったわけだ。 まあ普通に考えて同じコンピュータで同じアルゴリズムを走らせてそんな数百倍の違いになることのほうがどうかしてるのだから、こういうからくりでまあ間違いないだろう。 さて、問題はそんなことではない。 この、キャッシュが全然効かない状況でどうやって処理速度を稼ぐのか。 C++のネイティブコードがPythonの二倍程度の速度でしかないなんていうことは常識的にあり得ないのである。 もしかするとnumpyの中でSIMD命令を使っているのではないか? そう思って調べてみると案の定、numpyの内部ではSIMD命令を使って高速化されていた。 SIMDとは、Single Instruction Multiple Data streamの略で、要するに単一の命令で複数のデータを一括処理しちゃおうという思想である。 むかしはSIMD命令を使うとすると全てインラインアセンブリでやっていたが、今はC++からでもイントリンジックという愉快な仕組みを導入すると好き放題使えるらしい。 僕のMacBookAirはAVX2までをサポートしているらしいので、最大256ビットのSIMDレジスタを使うことができる。 256ビットのSIMDレジスタは、32ビット単精度を8つか、64ビット倍精度を4つ同時に処理することができる。 やってやろうじゃん。 ということで元はこんな感じになっていたコードを
    Vector *apply(Vector *v){
    int offset=0;
    for(int j=0;j<h;j++){
    double acc=0.0;
    offset+=w;
    for(int i=0;i<w;i++){
    acc += m[offset+i] * (v->v[i]);
    }
    r->v[j]=acc;
    }
    return r;
    }

     こんな感じに書き換えてみた。
    Vector *apply(Vector *v){
    int offset=0;
    for(int j=0;j<h;j++){
    double acc=0.0;
    offset+=w;
    __m256d M,V,Tmp,Acc; 
    double *pm = m+offset;
    double *pv = v->v;
    Acc = _mm256_setzero_pd();
    int end = (w/4)*4;
    int i;
    for(i=0;i<end;i+=4){
    M = _mm256_load_pd(pm);
    V = _mm256_load_pd(pv);
    Tmp=_mm256_mul_pd(M,V);
    Acc=_mm256_add_pd(Tmp,Acc);
    pm+=4;
    pv+=4;
    }
    acc = Acc[0]+Acc[1]+Acc[2]+Acc[3];
    for(;i<w;i++){
    acc += *pm * *pv;
    pm++;
    pv++;
    }
    r->v[j]=acc;
    }
    return r;
    }

     地味に配列参照ではなくポインタにするなど地道な努力もかいま見える。
     SIMD命令を使う場合は複数のデータを一括で渡さなければならないのでポインタ演算の方が渡すのに都合がいい。 __m256dという型を宣言するとこれはAVXレジスタに割り当てられる。 今のAVXはレジスタが8本あるらしいのでまあさらに病的な改造を加えればもうちょっと速くなるのかもしれないが汎用処理で回してるであろうnumpyごときと戦うのにはこの程度で充分だろう(いやどうかな行列とベクトルの乗算だからもっと特化した計算になっていても不思議はないが)。

    $ ./avx
    0.075887 sec
    31.951800 

      ふむ。これでまあなんとかC++の面目が保てたか。 おっと忘れてた。コンパイラ先生にも少し頑張ってもらおうか。-O2発動!



    $ ./avx 
    0.022472 sec
    31.683862 


     来ましたね。 さらに、単精度にして8つ同時に計算したらどうなるか。


    $ ./avx 
    0.017531 sec
    0.000000 


     でましたね。 5回繰り返して平均を取った数字でグラフ化するとこんな感じ まあなんですか。10倍程度高速っていう感じになりました。 いやー久しぶりにCPUの中身いじった清々しさがあるなあ。まったく。 ニューラル・ネットワークの世界にやってきても、最適化野郎の出番はあるぜ、という美しい結論で終わっておきましょうか。 
  • 【機械学習】Chainerで作ったニューラルネットをC++に移植してベンチマークしてみる

    2016-08-29 20:40  
    102pt
     マシン語 それは男のロマンである。 異論は認めん。 CPUの能力を限界まで引き出すコード。 そう、マシンたちの歓喜の声が聞こえてくる。 もっとオレを回してくれ! マイクロコードの1クロックまで! レジスタの限界まで オレを、オレをもっと愛してくれ! コアのひとつひとつ、パイプラインの全てのステージをコードで満たしてくれ! オレはもっと回りたい。 もっと回したいんだ・・・ ・・・そうだ、いい子だ。 それがマシン語である。 マシン語に比べると最近はPythonだのRubyだの適当な擬似言語で満足する輩が多すぎる。 軟弱者め!!! とまあ意気込んでは見たものの・・・いきなりマシン語はつらたんなのでとりあえずネイティブコードで動かしてみることを今回の記事では目標にする。 題材はいつものアレ 9-8-8-6の中間層二層のニューラル・ネットワーク。 こいつをC++に移植してみる。 ウヒョー、ワクワクするぜっ! まず、C++の文法を思い出すところからはじめなければならない。 いやー、忘れたよ。 最初はJSONで出力してJSONで扱おうかと思ったのだが、ネットワークが充分小さいのでヘッダファイルにコンバートすればいいか、ということでL1とL2のW(重み)とb(バイアス)をlayers.hに出力してみる。
    double  l1_W [8] [9] ={
    { 0.367395937443 ,
    0.620269238949 ,
    -0.0110975131392 ,
    0.0458976365626 ,
    -2.03605866432 ,

     こんな感じ。 まず、C++の多次元配列の作り方からして忘れてたからね。 ずっとエラーが出てて「なんで?」と思ったら、配列の初期化を[]でやってたからだった。 どんだけ離れてたんだよ。 まあ慣れてくれば昔とった杵付。 なんとかなりまさあ とりあえずデータをC言語で扱おうと思ったのだが・・・Cはつらすぎる。 だめだ、おれのオブジェクトが思考を始めてしまった。 オーバーヘッドは爆増するが、ここはひとつクラスを作らせてつかぁさい!
    class Vector{
    public:
    double *v;
    int size;
    Vector(double *_v,int _size){
    v = _v;
    size = _size;
    v = new double[_size];
    for(int i=0;i<_size;i++){
    v[i]=_v[i];
    }
    };
    Vector(int _size){
    size = _size;
    v = new double[_size];
    };
    Vector *add(Vector *_v){
    for(int i=0;i<size;i++){
    v[i]+=_v->v[i];
    }
    return this;
    };
    Vector *sigmoid(){
    for(int i=0;i<size;i++){
    v[i]=_sigmoid(v[i]);
    }
    return this;
    };
    void dump(){
    printf("dim:%d\n",size);
    for(int i=0;i<size;i++){
    printf("%f,",v[i]);
    }
    puts("");
    }
    };

     テンソルのクラスを作るのはつらそうだったので軟弱なオレはもうベクトルのクラスと行列のクラスだけでなんとかすることにする。いんだよ細けぇことは。畳込みNNだってテンソルっつうても画像がたくさんあるだけだからな。学習しないでインファレンス側で使うだけならこれでいんだよ別に。 しかしクラスの定義とか参照とかポインタとか盛大に忘れていて焦った。 やはりプログラミング言語が進化するというのは大きな理由と意義がある。 「.」と「->」の使い分けとか、本当に無駄すぎる。まあしかし、マシン語を意識する場合はこれがわかんないと話しになんないわけだけどな。 思い出しながら作ったのでnewとかが盛大に無駄こいてる。 本来はパフォーマンスを再優先するならばせめてオブジェクトプールを使うべきだ。 まあSTLとか使うという手もあるだろう。 しかし、コンストラクタのオーバーロードとか作ったの何年ぶりだろう。 だがしかし、STLまで思い出すのは辛すぎるので今回はいいや。まあこれでパフォーマンスがうまく上がんなかったら別の手を考えるベエ。
    class Matrix{
    public:
    double *m;
    int w,h;
    Matrix(double *_m,int _h,int _w){
    int i,j;
    h=_h;
    w=_w;
    m = new double[w*h];
    for(j=0;j<h;j++){
    for(i=0;i<w;i++){
    m[j*w+i] = _m[j*w+i];
    }
    }
    };
    Vector *apply(Vector *v){
    Vector *r = new Vector(h); 
    for(int j=0;j<h;j++){
    double acc=0.0;
    for(int i=0;i<w;i++){
    acc += m[j*w+i] * (v->v[i]);
    }
    r->v[j]=acc;
    }
    return r;
    }
    };

     男の行列クラス! もう男らしい。行列に対してベクトルを乗算する機能しかない! 乗算する度にnewしてるとか昔の自分なら絶対許さないコードだが、いんだよ細けぇことは後で考えれば。 これでとりあえず下準備完了。 次にフォワード(順伝播)を作る。 シグモイド使うやつと使わない奴
    Vector *forward(Vector *x,Matrix *W,Vector *b){
    Vector *a = W->apply(x); 
    a->add(b);
    return a;
    }
    Vector *forwardSigmoid(Vector *x,Matrix *W,Vector *b){
    Vector *a = W->apply(x); 
    a->add(b);
    a->sigmoid();
    return a;
    }

     ちなみにC++の場合、シグモイドを自前で実装するとこうなる。
    double _sigmoid(double x) {
        return 1.0 / (1.0 + exp(-1.0 * x));
    }


     さあこれで準備完了。 いくぜ男のメイン関数
    int main(){
    /*テスト用データ*/
    double in[] = {0.1083984375,0.213452398777,0.774722874165,0.110441893339,0.212938994169,0.780655622482,
    0.112617187202,0.212380990386,0.785996675491};
    Vector *inV = new Vector(in,9);
    Matrix *l1W = new Matrix((double*)l1_W,8,9);
    Vector *l1b = new Vector(l1_b,8);
    Matrix *l2W = new Matrix((double*)l2_W,6,8);
    Vector *l2b = new Vector(l2_b,6);
    Vector *h;
    clock_t start=clock();
    double dummy=0.0;
    for(int i=0;i<1000;i++){
    h=forwardSigmoid(inV,l1W,l1b);
    h=forward(h,l2W,l2b);
    dummy+=h->v[0];
    }
    clock_t end=clock();
    double t=(end-start)/(double)CLOCKS_PER_SEC;
    printf("%f sec\n",t);
    printf("%f \n",dummy);
    h->dump();

     テストデータは前回Pythonの整数版の精度検証に使ったのと同じものを使用。 これが同じになればちゃんと動いてると推定できる。 1000回ループ回してるところでなぞの変数dummyを足している理由は、計算結果を利用しないとコンパイラが万が一炎の最適化をかけてしまったらベンチマークにならないからだ(表示してるのも同様の理由)。 さあ、実行してみよう。ちなみにPython版では最も高速だったFloat64で1000回ループに22ミリ秒掛かったんだよねー。 レッツゴー! ・・・? 1.6ミリ秒・・・だ・・・と・・・!!!!!?!??!?!?!?!?!?! 念のためもう一回実行してみる。  1ミリ秒 おいおい。 ということは、Python版の1/200の実行時間で終わってしまうということか?どういうことだってばよ。 今回の速度比較はこうなった。 これはなんの冗談なんだろう。本当に計算してんのかな?と思って、試しに浮動小数点の足し算だけのループをまわしてみると確かに高速化しているのでまあ何らかの計算はしているのだろう。試しにPythonのfloat64版と同じくらいの実行時間になるにはループがどのくらいあればいいのか調べてみた。これでもPython/float64よりはだいぶ高速だが、まあこの結果になるのにどのくらいのループをしたかというと
    for(int i=0;i<200000;i++){
    h=forwardSigmoid(inV,l1W,l1b);
    h=forward(h,l2W,l2b);
    dummy+=h->v[0];
    }


    20万回でした。つまり、どうやらC++で書き直すだけでCPUでの実行時間において、全く同じニューラル・ネットワークでPython版の200倍高速に、Chainerの500倍高速に計算できることがわかりました。まあ畳込み絡んできたりGPU絡んできたりすれば状況は違うんだろうけど。やはりChainerは便利で簡単に書ける分、代償も少なくないんだろうなあ。