読者です 読者をやめる 読者になる 読者になる

kivantium活動日記

プログラムを使っていろいろやります

C++による多層パーセプトロンの実装

注意:このページの実装は古いものです。更新はkivantium/libNN · GitHubで行っています

前回の記事でPylearn2というDeep learningライブラリを導入しましたが、導入したもののどうやって使えばいいのかが全く分からず使用を断念しました。他にもCaffeなどいくつか実装があるようですが、どれもよく分かりませんでした。
「ライブラリがないなら作ればいいじゃない」と脳内マリー・アントワネットが語りかけてきたので、自分で実装するという茨の道を選ぶことにしました。

というわけでDeep learning実装の第一歩としてニューラルネットワークの基本になる多層パーセプトロンの実装を行いました。多層パーセプトロンご注文は機械学習ですか?で使った分類器です。OpenCVに用意されているので車輪の再発明になりますが、勉強としてはちょうどいいでしょう。

多層パーセプトロンとは

多層パーセプトロンは神経細胞のメカニズムを参考に考案されたニューラルネットワークの一種です。図のようなユニットと呼ばれる単位が結びついて層のような構造をしています。
f:id:kivantium:20141221234412p:plain:w300
一番最初の層を入力層、最後の層を出力層、中間の層を隠れ層と呼びます。

入力層には判別したいデータをベクトル表示したものの各成分が与えられます。データをベクトル表示という考え方が馴染みにくいかもしれませんが、8x8のグレースケール画像から各ピクセルの値を使って64次元ベクトルを作るというように考えてみてください。

隠れ層・出力層の入力は一つ下の層の値に重み付けして足し合わせた和によって決定されます。
入力層の1番目のユニットが1つ目の隠れ層の1番目に情報を伝えるときの値を、入力層の2番目のユニットが隠れ層の1番目のユニットに情報を伝えるときの値を……、というように重み付けの値に名前をつけています。

隠れ層の出力は、入力値に対する非線形関数によって決められることが多いです。tanhを使うこともありますが、今回はシグモイド関数

を使います。
シグモイド関数には

という性質があります。

バイアスとして入力層と隠れ層には常に1を出力するユニットを作ることが多いようです。

出力層は入力をそのまま出力します。

出力層では正解とするクラスを1、正解ではないクラスを0となるようにすることで分類を行います。(クラス5に分類するなら出力層の5番目のユニットが1、その他のユニットが0であればよい)実際には正解の出力をぴったり1にするのは不可能なので、最大の数字を出力したクラスに分類するようにしています。

多層パーセプトロンの学習は、訓練データに対して正しいクラス分類を返すような重みづけの設定を行うことを目標としています。それを行うためには誤差逆伝播法という、うまい重み付けを見つけ出すアルゴリズムを使います。

誤差逆伝播

証明は省略しますが、K層からなる多層パーセプトロン誤差逆伝播法はこのようにして行います。(入門 パターン認識と機械学習を参考にしました)

  • 初期化
    • それぞれの重みの初期値を乱数で与える
  • 順伝播

学習パターンx(ベクトル)に対し、i番目の層の入力ξiを次のように定める。


出力は次のように決める。(出力層ではf(x)=xとしている)

  • 誤差の計算
    • 次の式に従って層ごとの誤差を計算する。

k=K(出力層)のとき (y_jは正しい分類結果を表す)

k≠Kのとき

  • 重みの更新

次の式に従って重みを更新する (ηは学習定数と呼ばれる定数)

これを素直に実装したのが次のクラスです。

C++による実装例

#include <string>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <cmath>

class mlp {
private:
    // number of each layer
    int in, hid, out;
    // input layer
    float *xi1, *xi2, *xi3;
    // output layer
    float *o1, *o2, *o3;
    // error value
    float *d2, *d3;
    // wait
    float *w1, *w2;

    //sigmoid function
    float sigmoid(float x){
        return 1/(1+exp(-x));
    }

    //derivative of sigmoid function
    float d_sigmoid(float x){
        return (1-sigmoid(x))*sigmoid(x);
    }

    //caluculate forward propagation of input x
    void forward(float *x){
        //calculation of input layer
        for(int j=0; j<in; j++){
            xi1[j] = x[j];
            xi1[in] = 1;
            o1[j] = xi1[j];
        }
        o1[in] = 1;

        //caluculation of hidden layer
        for(int j=0; j<hid; j++){
            xi2[j] = 0;
            for(int i=0; i<in+1; i++){
                xi2[j] += w1[i*hid+j]*o1[i];
            }
            o2[j] = sigmoid(xi2[j]);
        }
        o2[hid] = 1;

        //caluculation of output layer
        for(int j=0; j<out; j++){
            xi3[j] = 0;
            for(int i=0; i<hid+1; i++){
                xi3[j] += w2[i*out+j]*o2[i];
            }
            o3[j] = xi3[j];
        }
    }

public:
    mlp(const std::string filename){
        std::ifstream ifs(filename.c_str());
        if(!ifs){
            std::cout << "File does not exist" << std::endl;
            exit(1);
        }
        std::string str;
        ifs >> str;
        if(str!="LIBNN_MLP_0"){
            std::cout << "File type error!" << std::endl;
            exit(1);
        }
        ifs >> in >> hid >> out;
        // allocate inputs
        xi1 = new float[in+1];
        xi2 = new float[hid+1];
        xi3 = new float[out];

        // allocate outputs
        o1 = new float[in+1];
        o2 = new float[hid+1];
        o3 = new float[out];

        // allocate errors
        d2 = new float[hid+1];
        d3 = new float[out];

        // allocate memories
        w1 = new float[(in+1)*hid];
        w2 = new float[(hid+1)*out];

        for(int i=0; i<(in+1)*hid; i++)  ifs >> w1[i];
        for(int i=0; i<(hid+1)*out; i++) ifs >> w2[i]; 
    }
    /* constructor
     * n_in:  number of input layer
     * n_hid: number of hidden layer
     * n_out: number of output layer */
    mlp(int n_in, int n_hid, int n_out){
        // unit number
        in = n_in;
        hid = n_hid;
        out = n_out;

        // allocate inputs
        xi1 = new float[in+1];
        xi2 = new float[hid+1];
        xi3 = new float[out];

        // allocate outputs
        o1 = new float[in+1];
        o2 = new float[hid+1];
        o3 = new float[out];

        // allocate errors
        d2 = new float[hid+1];
        d3 = new float[out];

        // allocate memories
        w1 = new float[(in+1)*hid];
        w2 = new float[(hid+1)*out];

        // initialize wait
        // range depends on the research of Y. Bengio et al. (2010)
        srand ((unsigned) (time(0)));
        float range = std::sqrt(6)/std::sqrt(in+hid+2);
        std::srand ((unsigned) (std::time(0)));
        for(int i=0; i<(in+1)*hid; i++) w1[i] = (float) 2*range*std::rand()/RAND_MAX-range;
        for(int i=0; i<(hid+1)*out; i++) w2[i] = (float) 2*range*std::rand()/RAND_MAX-range;
    }

    // deconstructor
    ~mlp(){
        // delete arrays
        delete[] xi1, xi2, xi3;
        delete[] o1, o2, o3;
        delete[] d2, d3;
        delete[] w1, w2; 
    }

    /* train: multi layer perceptron
     * x: train data(number of elements is in*N)
     * t: correct label(number of elements is N)
     * N: data size
     * repeat: repeat times
     * eta: learning rate */
    void train(float x[], int t[], float N, int repeat=1000, float eta=0.1){
        for(int times=0; times<repeat; times++){
            for(int sample=0; sample<N; sample++){
                // forward propagation
                forward(x+sample*in);

                // calculate the error of output layer
                for(int j=0; j<out; j++){
                    if(t[sample] == j) d3[j] = o3[j]-1;
                    else d3[j] = o3[j];
                }
                // update the wait of output layer
                for(int i=0; i<hid+1; i++){
                    for(int j=0; j<out; j++){
                        w2[i*out+j] -= eta*d3[j]*o2[i];
                    }
                }
                // calculate the error of hidden layer
                for(int j=0; j<hid+1; j++){
                    float tmp = 0;
                    for(int l=0; l<out; l++){
                        tmp += w2[j*out+l]*d3[l];
                    }
                    d2[j] = tmp * d_sigmoid(xi2[j]);
                }
                // update the wait of hidden layer
                for(int i=0; i<in+1; i++){
                    for(int j=0; j<hid; j++){
                        w1[i*hid+j] -= eta*d2[j]*o1[i];
                    }
                }
            }
        }
    }

    // return most probable label to the input x
    int predict(float x[]){
        // forward propagation
        forward(x);
        // biggest output means most probable label
        float max = o3[0];
        int ans = 0;
        for(int i=1; i<out; i++){
            if(o3[i] > max){
                max = o3[i];
                ans = i;
            }
        }
        return ans;
    }

     // save mlp to csv file
    void save(const std::string filename){
        std::ofstream ofs(filename.c_str());
        ofs << "LIBNN_MLP_0" << std::endl;
        ofs << in << " " << hid << " " << out << std::endl;
        for(int i=0; i<(in+1)*hid-1; i++) ofs << w1[i] << " ";
        ofs << w1[(in+1)*hid-1] << std::endl;
        for(int i=0; i<(hid+1)*out-1; i++) ofs << w2[i] << " "; 
        ofs << w2[(hid+1)*out-1] << std::endl;
        ofs.close();
    }
};

このクラスをmlp.hという名前で保存します。

このクラスを使うとニューラルネットワーク練習問題にのせた3クラス分類問題のq2.csvは次のようなコードで分類することができます。

#include <iostream>
#include <cstdio>
#include "mlp.h"
using namespace std;

int main(void){
    // サンプルの数
    const int sample = 150;
    // 入力ベクトルの次元
    const int size = 2;
    // ラベルの種類
    const int label = 3;
    
    // 入力ユニット2つ, 隠れユニット3つ, 出力ユニット3つの多層パーセプトロンを作る
    mlp net(size,3,label);
  
    // 訓練データ
    float x[size*sample];
    // ラベルデータ
    int t[sample];
    
    // CSVの読み込み
    FILE *fp = fopen("sample.csv", "r");
    if(fp==NULL) return -1;
    for(int i=0; i<sample; i++){
        // ラベルの読み込み
        fscanf(fp, "%d,", t+i);
        // データの読み込み
        for(int j=0; j<size; j++) fscanf(fp, "%f,", x+size*i+j);
    }
   
    // 多層パーセプトロンの訓練
    // 繰り返しの回数
    const int repeat = 500;
    net.train(x, t, sample, repeat);
    // 結果の表示
    for(int i=0; i<sample; i++){
        cout << net.predict(x+size*i) << endl;
    }
    
    return 0;
}

というわけで車輪の再発明でしたが、多層パーセプトロンを自前で実装することでブラックボックスを一つ減らすことができるようになりました。次はDeep learningの実装を行いたいところですが、いつになるやら……。