SdAで遊ぶ

はじめに

Deepな話で、簡単に試せそうだったStacked Denoising AutoEncoderを試しに遊んでみる。
あんまり詳しく調べていないので、お遊びレベルという感じで・・・

注意:下記では「特徴抽出器」として利用する方法を試しています。通常は事前学習として行い、それを初期値に使い、普通にニューラルネットの学習を行うことを指すと思いますが、下記のような特徴抽出器的使い方もありみたいですね(ref.機械学習プロフェッショナルシリーズ「深層学習」pp.72-74)。

Stacked AutoEncoderとは

• BengioやVinsentによって提案や紹介
• AutoEncoderを何層も重ねたもの
• 各層の学習は、一つ前の隠れ層を入力にAutoEncoderを学習し、出力部分を捨てて次の層を学習する
  • Unsupervised layer-wise pre-training
• 層の最後の部分でロジスティック回帰などを利用する事で教師あり分類器などを行える

コード

liblinear形式のデータを、学習し、学習したSdA通した後の特徴に変換してliblinear形式で出力まで。
結果をliblinearで学習・評価してみる。

#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <sstream>
#include <cstdlib>

//行列計算用クラス
class Matrix {
  int m, n;
  std::vector<double> val;

  void splice(int i, int j){
    int N = val.size();
    if(i+(j-1)<N){
      std::vector<double> vl(N-j, 0.0);
      int vl_idx = 0, val_idx = 0;
      while(val_idx<N){
        if(val_idx<i||i+j-1<val_idx){
          vl[vl_idx] = val[val_idx];
          vl_idx++;
        }
        val_idx++;
      }
      val = vl;
    }
  }
  Matrix clone(){
    Matrix ret(m, n);
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        ret.setVal(i, j, val[n*i+j]);
      }
    }
    return ret;
  }

public:
  Matrix():m(0),n(0){}
  Matrix(int m, int n):m(m),n(n){
    if(m>0 && n>0){
      for(int i=0; i<m*n; i++){
        val.push_back(0.0);
      }
    }
  }
  Matrix(const Matrix& mat){
    m = mat.m;
    n = mat.n;
    val = mat.val;
  }

  int getM() const { return m; }
  int getN() const { return n; }
  double getVal(int i, int j) const { return val[n*i+j]; }
  void setVal(int i, int j, double x){ val[n*i+j] = x; }
  bool isSquare(){ return n==m; }

  Matrix add(const Matrix& mat){
    if(m == mat.m && n == mat.n){
      Matrix ret(m, n);
      for(int i=0; i<m; i++){
        for(int j=0; j<n; j++){
          ret.setVal(i, j, val[n*i+j] + mat.getVal(i,j));
        }
      }
      return ret;
    }
    return Matrix(-1, -1);
  }
  Matrix operator+(const Matrix& mat){ return add(mat); }

  Matrix sub(const Matrix& mat){
    if(m == mat.m && n == mat.n){
      Matrix ret(m, n);
      for(int i=0; i<m; i++){
        for(int j=0; j<n; j++){
          ret.setVal(i, j, val[n*i+j] - mat.getVal(i,j));
        }
      }
      return ret;
    }
    return Matrix(-1, -1);
  }
  Matrix operator-(const Matrix& mat){ return sub(mat); }

  Matrix prod(const Matrix& mat){
    if(n == 1 && mat.n == 1 && m == mat.m){
      Matrix ret(m, 1);
      for(int i=0; i<m; i++){
        ret.setVal(i, 0, getVal(i, 0) * mat.getVal(i, 0));
      }
      return ret;
    }
    else if(n == mat.m){
      Matrix ret(m, mat.n);
      for(int i=0; i<m; i++){
        for(int j=0; j<mat.n; j++){
          double d = 0.0;
          for(int k=0; k<n; k++){
            d += val[n*i+k] * mat.getVal(k,j);
          }
          ret.setVal(i, j, d);
        }
      }
      return ret;
    }
    return Matrix(-1, -1);
  }
  Matrix operator*(const Matrix& mat){ return prod(mat); }

  void time(double x){
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        val[n*i+j] *= x;
      }
    }
  }

  Matrix transpose(){
    Matrix ret(n, m);
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        ret.setVal(j, i, val[n*i+j]);
      }
    }
    return ret;
  }
  
  double cofactor(int i, int j){
    Matrix mat = clone();
    mat.splice(i*mat.n, mat.m);
    mat.m -= 1;
    for(int k=mat.m-1; k>=0; k--){
      mat.splice(k*mat.n+j, 1);
    }
    mat.n -= 1;
    return mat.det() * ( ((i+j+2)%2==0) ? 1 : -1);
  }
  double det(){
    if(isSquare()){
      if(m == 2){
        return val[0]*val[3]-val[1]*val[2];
      }else{
        double d = 0;
        for(int k=0; k<n; k++){
          d += val[k] * cofactor(0, k);
        }
        return d;
      }
    }else{
      return 0.0;
    }
  }
  Matrix cofactorMatrix(){
    Matrix mat(m, n);
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        mat.setVal(j, i, cofactor(i, j));
      }
    }
    return mat;
  }
  Matrix inverse(){
    if(isSquare()){
      double d = det();
      if(d != 0){
        Matrix mat;
        if(m>2){
          mat = cofactorMatrix();
        } else {
          mat = Matrix(2, 2);
          mat.setVal(0, 0, val[3]);
          mat.setVal(0, 1, -val[1]);
          mat.setVal(1, 0, -val[2]);
          mat.setVal(1, 1, val[0]);
        }
        mat.time(1 / d);
        return mat;
      }
    }else{
      return Matrix(-1, -1);
    }
    return Matrix(-1, -1);
  }
};
std::ostream& operator<<(std::ostream& os, const Matrix& mat){
  for(int i=0; i<mat.getM(); i++){
    for(int j=0; j<mat.getN(); j++){
      os << mat.getVal(i,j) << " ";
    }
    os << std::endl;
  }
  return os;
}



//乱数計算用関数
//xorshift
// 注意: longではなくint(32bit)にすべき
unsigned long xor128(){
  static unsigned long x=123456789, y=362436069, z=521288629, w=88675123;
  unsigned long t;
  t=(x^(x<<11));
  x=y; y=z; z=w;
  return w=(w^(w>>19))^(t^(t>>8));
}
//[0,1)の一様乱数
// 注意: int(32bit)にしたらint_maxで割ったほうがよい
double frand(){
  return xor128()%10000000000/static_cast<double>(10000000000); 
}


// binary-inputs denoising autoencoder (tied W, W' = W^T)
//  y = sigmoid(W x' + b)
//  z = sigmoid(W' y + b') ~ x
class DenoisingAutoEncoder {
  int num_visible, num_hidden;
  Matrix W, bvis, bhid;
  double p, eps;

  double sigmoid(double x){
    return 1.0 / (1.0 + exp(-x));
  }

public:
  DenoisingAutoEncoder(int num_visible = 1, int num_hidden = 1, double p = 0.0, double eps = 0.0):
    num_visible(num_visible),
    num_hidden(num_hidden),
    W(num_hidden, num_visible),
    bvis(num_visible, 1),
    bhid(num_hidden, 1),
    p(p),
    eps(eps)
  {
    //initial W
    for(int i=0; i<num_hidden; i++){
      for(int j=0; j<num_visible; j++){
        W.setVal(i, j, 1.0 - 2.0 * frand());
      }
    }
  }

  void train(const std::vector<int>& in){
    //make input x
    Matrix x(num_visible, 1);
    for(int i=0; i<in.size(); i++){
      x.setVal(i, 0, in[i]);
    }
    //noise addition
    Matrix noisex(num_visible, 1);
    for(int i=0; i<in.size(); i++){
      double val = in[i];
      if(frand() < p){
        if(val < 0.5) val = 1.0;
        else val = 0.0;
      }
      noisex.setVal(i, 0, val);
    }
    //vector one
    Matrix one(num_hidden, 1);
    for(int i=0; i<num_hidden; i++){
      one.setVal(i, 0, 1.0);
    }
    //make y
    Matrix y = W * noisex + bhid;
    for(int i=0; i<num_hidden; i++){
      y.setVal(i, 0, sigmoid(y.getVal(i, 0)));
    }
    //make z
    Matrix z = W.transpose() * y + bvis;
    for(int i=0; i<num_visible; i++){
      z.setVal(i, 0, sigmoid(z.getVal(i, 0)));
    }

    //make grad bhid
    Matrix gradbhid = W * (x-z) * y * (one-y);
    //make grad bvis
    Matrix gradbvis = x - z;
    //make grad W
    Matrix gradW = (gradbhid * noisex.transpose()) + ((gradbvis * y.transpose()).transpose());

    //update W
    for(int i=0; i<num_hidden; i++){
      for(int j=0; j<num_visible; j++){
        W.setVal(i, j, W.getVal(i, j) + eps * gradW.getVal(i, j));
      }
    }

    //update b
    for(int i=0; i<num_hidden; i++){
      bhid.setVal(i, 0, bhid.getVal(i, 0) + eps * gradbhid.getVal(i, 0));
    }

    //update b'
    for(int i=0; i<num_visible; i++){
      bvis.setVal(i, 0, bvis.getVal(i, 0) + eps * gradbvis.getVal(i, 0));
    }
  }
  
  //隠れ層の値を取得
  std::vector<double> getHiddenValues(const std::vector<int>& in){
    std::vector<double> ret;
    
    //make input x
    Matrix x(num_visible, 1);
    for(int i=0; i<in.size(); i++){
      x.setVal(i, 0, in[i]);
    }
    //make y
    Matrix y = W * x + bhid;
    for(int i=0; i<num_hidden; i++){
      ret.push_back( sigmoid(y.getVal(i, 0)) );
    }
    
    return ret;
  }

  //出力層の値を取得
  std::vector<double> getOutputValues(const std::vector<int>& in){
    std::vector<double> ret;

    //make input x
    Matrix x(num_visible, 1);
    for(int i=0; i<in.size(); i++){
      x.setVal(i, 0, in[i]);
    }
    //make y
    Matrix y = W * x + bhid;
    for(int i=0; i<num_hidden; i++){
      y.setVal(i, 0, sigmoid(y.getVal(i, 0)));
    }
    //make z
    Matrix z = W.transpose() * y + bvis;
    for(int i=0; i<num_visible; i++){
      ret.push_back( sigmoid(z.getVal(i, 0)) );
    }

    return ret;
  }

  //重み行列Wを出力
  void dumpW(){
    std::cout << "Weight W" << std::endl;
    std::cout << W << std::endl;
  }
  //バイアスbとb'を出力
  void dumpb(){
    std::cout << "bhid" << std::endl;
    std::cout << bhid << std::endl;

    std::cout << "bvis" << std::endl;
    std::cout << bvis << std::endl;
  }
};


//SdA
class StackedAutoEncoders {
  int epoch;
  DenoisingAutoEncoder* das[10];
  std::vector<int> hN;

  StackedAutoEncoders(const StackedAutoEncoders&);
public:
  StackedAutoEncoders(int dim_input, const std::vector<int>& hN_input, double noiseP, double eps, int epoch_input):
    hN(hN_input),
    epoch(epoch_input)
  {
    //10層まで
    if(hN.size() > 10) throw "too deep";
    for(int i=0; i<10; i++) das[i] = NULL;

    //各層のDenoisingAutoEncoderを用意する
    for(int i=0; i<hN.size(); i++){
      if(i==0){
        das[i] = new DenoisingAutoEncoder(dim_input, hN[i], noiseP, eps);
      }else{
        das[i] = new DenoisingAutoEncoder(hN[i-1], hN[i], noiseP, eps);
      }
    }   
  }
  ~StackedAutoEncoders(){
    for(int i=0; i<hN.size(); i++){
      if(das[i] != NULL){
        delete das[i];
        das[i] = NULL;
      }
    }
  }

  void train(const std::vector< std::vector<int> >& inputs){
    std::vector< std::vector<int> > now(inputs.size(), std::vector<int>(inputs[0].size()));
    for(int i=0; i<inputs.size(); i++){
      for(int j=0; j<inputs[i].size(); j++){
        now[i][j] = inputs[i][j];
      }
    }

    //greedy & layer-wise pre-training
    for(int d=0; d<hN.size(); d++){
      //層dについて、入力nowで学習する
      for(int t=0; t<epoch; t++){
        std::cout << "epoch" << t << std::endl;
        for(int i=0; i<now.size(); i++){
          std::cout << i << " " << now[i].size() << std::flush;
          das[d]->train(now[i]);
          std::cout << "             \r";
        }
        std::cout << std::endl;

        if(t%10==0){
          //エラー率の確認
          int errSum = 0;
          for(int i=0; i<now.size(); i++){
            std::vector<double> res = das[d]->getOutputValues(now[i]);
            for(int j=0; j<now[i].size(); j++){
              if( now[i][j] != (res[j]>0.5?1:0) ) errSum++;
            }
          }
          std::cout << "layer[" << d << "]" << " : Err = " << (100.0 * errSum / (now.size() * now[0].size())) << "%" << std::endl;
        }

      }
      
      //学習結果から隠れ層を次の入力にする
      for(int i=0; i<now.size(); i++){
        std::vector<double> resH = das[d]->getHiddenValues(now[i]);
        std::vector<int> newIn;
        for(int j=0; j<resH.size(); j++){
          newIn.push_back( (resH[j]>0.5?1:0) ); //01入力のDAEなので、0.5より大きければ1,そうでなければ0と決定的に決めたものを使用する
        }
        now[i] = newIn;
      }
    }
  }

  //SdAを通した後の特徴の出力
  std::vector<int> output(const std::vector<int>& input){
    std::vector<int> res = input;

    for(int d=0; d<hN.size(); d++){
      std::vector<double> resH = das[d]->getHiddenValues(res);
      res.clear();
      for(int j=0; j<resH.size(); j++){
        res.push_back( (resH[j]>0.5?1:0) ); //同上
      }
    }
    return res;
  }

};


int main(){

  std::string trainFilename, testFilename;
  int InputDim, dataNum, Epoch;
  double noiseP, epsilon;
  int hNsize;

  //設定の読み込み
  std::cin >> trainFilename >> testFilename;
  std::cin >> InputDim >> dataNum;
  std::cin >> noiseP >> epsilon >> Epoch;
  std::cin >> hNsize;
  std::vector<int> hN;
  for(int i=0; i<hNsize; i++){
    int tmp;
    std::cin >> tmp;
    hN.push_back(tmp);
  }


  std::ifstream trainfile(trainFilename.c_str());
  std::string line, feat;
  //inputs
  //liblinear形式だけど、素性値は無視して、素性がでているものは1、でていないものは0として利用
  //「+1 1:0 2:1」などでも1:1 2:1として利用するので注意
  std::vector<int> in_y;
  std::vector< std::vector<int> > in(dataNum, std::vector<int>(InputDim, 0));
  for(int i=0; i<dataNum; i++){
    std::getline(trainfile, line);
    int t;
    std::stringstream ss(line);
    ss >> t;
    if(t>0) t = 1;
    else t = 0;
    in_y.push_back(t);

    while(ss >> feat){
      std::string::size_type idx = feat.find(":");
      int id = atoi(feat.substr(0, idx).c_str());
      int value = atoi(feat.substr(idx+1).c_str());

      in[i][id-1] = value;
    }
  }

  //train
  StackedAutoEncoders sae(InputDim, hN, noiseP, epsilon, Epoch);
  sae.train(in);

  //学習データの出力
  std::ofstream trainout((trainFilename + ".sda").c_str());
  for(int i=0; i<dataNum; i++){
    std::vector<int> res = sae.output(in[i]);

    trainout << in_y[i];
    for(int j=0; j<res.size(); j++){
      if(res[j] == 1){
        trainout << " " << j+1 << ":1";
      }
    }
    trainout << std::endl;
  }

  //評価データの出力
  std::ifstream testfile(testFilename.c_str());
  std::ofstream testout((testFilename + ".sda").c_str());
  while(std::getline(testfile, line)){
    int t;
    std::stringstream ss(line);
    ss >> t;
    if(t>0) t = 1;
    else t = 0;
    testout << t;
    
    std::vector<int> tmp(InputDim, 0);
    while(ss >> feat){
      std::string::size_type idx = feat.find(":");
      int id = atoi(feat.substr(0, idx).c_str());
      int value = atoi(feat.substr(idx+1).c_str());

      tmp[id-1] = value;
    }
    std::vector<int> res = sae.output(tmp);
    for(int i=0; i<res.size(); i++){
      if(res[i] == 1){
        testout << " " << i+1 << ":1";
      }
    }
    testout << std::endl;
  }

  return 0;
}

試し

とりあえず、いつも使ってるa9aを利用。
http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/

a9a

• 設定ファイル

# liblinear形式の学習ファイル liblinear形式のテストファイル
# 素性数 学習ファイルの行数
# ノイズの割合 学習率 各層のIteration回数
# 層の数
# 各層の隠れ層のノード数
# ...
$ cat a9a.test
a9a a9a.t
123 32561
0.1 0.001 50
9
110
100
90
80
70
60
50
40
30

• 実行

$ g++ -O2 sda.cc
$ ./a.out < a9a.test
...
$ ../liblinear/liblinear-1.93/train -s 0 -c 0.1 a9a.sda
$ ../liblinear/liblinear-1.93/predict a9a.t.sda a9a.sda.model out
Accuracy = 77.8515% (12675/16281)


# 参考
$ ../liblinear/liblinear-1.93/train -s 0 -c 0.1 a9a
$ ../liblinear/liblinear-1.93/predict a9a.t a9a.model out 
Accuracy = 85.0132% (13841/16281)

うあああ、がっつり下がった。。。

けど、123次元→30次元に落とした特徴でも77.8%ぐらい出てるので、まあまあ方向はあっていそうかも。
パラメータ未調整というのはおおいにあるけど、Iteration回数が十分に収束しないのに次の層の学習に進んでいるのでErrorが蓄積してしまっているっぽくも見える。(3層ぐらいだとそんなにAcc落ちずに次元だけ落ちる)
あと、時間結構かかる・・・

a9a(2014/1/16追記)

ちゃんと収束がしてくれる程度にIteration回数を確保＆学習率をちょっと大きめに取ってやってみる。
9層はおおすぎるかなと思うので、5層に減らす。

• 設定ファイル

a9a a9a.t
123 32561
0.1 0.005 500
5
110
100
90
80
70

• 実験

$ ../liblinear/liblinear-1.93/predict a9a.t.sda a9a.sda.model out
Accuracy = 84.2884% (13723/16281)

結構Accを落とさず次元減らせた。
いい感じの表現を事前学習できたかも。

画像とかで試した方がよさそうなので、やってみたい。

参考

• 人工知能学会誌、DeepLearning連載解説
• http://deeplearning.net/tutorial/SdA.html
• http://www.slideshare.net/alembert2000/learning-deep-architectures-for-ai-3-deep-learning
• http://cl.naist.jp/~kevinduh/a/deep2014/140116-ResearchSeminar.pdf
• http://blog.yusugomori.com/post/43636583955/c-c-deep-learning-deep-belief-nets-stacked
• http://www.sakurai.comp.ae.keio.ac.jp/classes/infosem-class/2012/15DeepLearning.pdf