Negative/Positive Thinking

2016-10-08

FastBDTでの高速化

はじめに

勾配ブースティング木の高速化はどうすればいいだろうと調べていたら、arxivで流れているのを見かけたのでメモ。

Stochastic Gradient Boosted Decision Tree(SGBDT)

高速化

メモリアクセス、CPUキャッシュ、前処理、並列化を考慮。

メモリアクセスパターン
  • 学習データのメモリ配置について最適化
  • 一般的に以下の2つのレイアウト方式が考えられる
    • array of structs(構造体の配列) : x_1, y_1, z_1, x_2, y_2, z_2, x_3, y_3, z_3, x_4, y_4, z_4
    • struct of arrays(配列の構造体) : x_1, x_2, x_3, x_4, y_1, y_2, y_3, y_4, z_1, z_2, z_3, z_4
  • FastBDTは「array of structs」を利用

【正例と負例のメモリマップ】

  • array of structsを利用する場合、そのまま並べてしまうと、正例と負例の判定で条件分岐が必要になってしまう
  • そこで、正例と負例のデータを分けて、保存や処理を行う

レイヤー毎の複数ノードの処理】

  • ノード単位で別々にノードを処理していく方法だと、使用する学習データへのアクセスがとびとびになってしまい、メモリレイアウト方式によらずランダムにジャンプしてアクセスしてしまう
    • (DFS的な感じでのノードの処理?)
  • 木のレイヤー単位で、そのレイヤーすべてのノードを並列処理することができると、連続的なメモリアクセスにできる

【確率的サブサンプリング

  • SGBDTだと、学習データからサンプリングしたデータを使用するため、使うデータと使わないデータを選ぶ必要がある
  • array of structsだと、各データポイントの範囲内であれば、条件分岐なしでデータへ連続的にアクセスできる

前処理
  • 決定木では、「値それ自身」というより「値同士の大小関係」しか使わない
  • そこで、入力素性についてequal-frequency binningをして整数値(インデクス)にマッピング
    • equal-frequency binning : 入力の値をbinningするとき、各ビンが同じデータ数になるようにbinningする(ソートしてm個ずつビンにいれる)
    • 速度面だけでなく、ピークやヘビーテイルを含む入力の分布を一様分布にマッピングすることで分離性能も改善しうる
  • また、木の学習が終わったら、整数値から小数値へ逆に変換すれば、実行時に変換のオーバーヘッドなく利用できる

並列化
  • タスクレベル(マルチプロセッサマルチコア、マルチスレッドなど)の並列化を適用できる
  • 命令レベル(命令パイプライン化、multiple execution unit、vectorizationなど)の並列化の適用は難しい
    • アルゴリズム固有の条件分岐が大量にあるため
    • ただ、メモリレイアウトと前処理から、条件分岐をインデクス操作に置き換えることで高速化はできる(branch locality)
      • 「if(X) a++; else b++;」ではなく、「arr[X]++;」のような書き方

参考

2016-10-05

勾配ブースティング木を試す

はじめに

ここしばらくサボってしまっているので、(のんびりと)いろいろ勉強していきたい。
コンテスト系などでもよい成績を残している勾配ブースティング木について試してみた。

勾配ブースティングとは

  • Gradient Boosting
  • 加法的モデルH(x)=Σρ_t * h_t(x)を、前向き段階的に各ρ,hを学習・追加していく
    • h_tは弱学習器、ρ_tは弱学習器tの貢献度
  • 雑に言うと、tでは(t-1)まででの予測の弱点を補うような弱学習器h_tを学習・追加する
  • 弱点は「予測結果と正解との差」が考えられるが、差ではなく(負の)勾配として考える
    • なんらかの損失関数を定めることで、その勾配を弱学習器h_tに学習させていく(勾配降下している感じ)
    • ρ_tは学習率に相当
  • 勾配ブースティング木は、弱学習器として決定木・回帰木を用いたもの
    • 多重加法的回帰木(multiple additive regression trees)とも

損失関数と負の勾配
  • Square Loss
    • loss = 1/2 * { y_i - H(x_i) }^2
    • -grad = y_i - H(x_i)
      • 「正解と予測値の差」を学習していく
  • Absolute Loss
    • loss = |y_i - H(x_i)|
    • -grad = sign(y_i - H(x_i))
  • Huber Loss
    • loss = 1/2 * { y_i - H(x_i) }^2 (|y+i - H(x_i)|<=δの場合), δ * (|y_i - H(x_i)| - δ/2) (otherwise)
    • -grad = y_i - H(x_i) (|y_i - H(x_i)|<=δの場合), δ * sign(y_i - H(x_i)) (otherwise)
      • ノイズなどで大きく外れやすいデータに対してロバスト

ライブラリ

実用で使う場合は、以下のような有用ライブラリがあるので、そちらを使うのがよい。

コード

毎度ながら、アルゴリズムを追いたいので、工夫や効率化などはしてない。
回帰木は、木の深さを指定し、そこまでノードを増やす(作れない場合はそこまで)。
変数の探索は、全データ全変数を探すので重い。

#include <iostream>
#include <vector>
#include <fstream>

//回帰木
class RT {
  struct NODE {
    int ln, rn; //左のノードID、右のノードID
    int j; //分割変数
    double th; //分割閾値
    std::vector<int> ids; //割り当てられた学習データのID
    double ybar; //平均値
    NODE():ln(-1),rn(-1),j(-1),th(0),ybar(0){}
  };

  int N; //木のノード数
  int depth; //木の深さ
  std::vector<NODE> nodes; //木のノード情報

  const double INF;
  
  //分割関数
  double split_value(int j,
                     double th,
                     int id,
                     const std::vector< std::vector<double> >& x,
                     const std::vector<double>& y){
    //分割したときのそれぞれのyの平均値
    double lybar = 0.0, rybar = 0.0;
    int lnum = 0, rnum = 0;
    for(size_t i=0; i<nodes[id].ids.size(); i++){
      int ii = nodes[id].ids[i];
      if(x[ii][j] < th){
        lybar += y[ii];
        lnum++;
      }else{
        rybar += y[ii];
        rnum++;
      }
    }
    if(lnum == 0 || rnum == 0) return INF;
    lybar /= lnum;
    rybar /= rnum;

    //それぞれの二乗和を求める
    double lcost = 0.0, rcost = 0.0;
    for(size_t i=0; i<nodes[id].ids.size(); i++){
      int ii = nodes[id].ids[i];
      if(x[ii][j] < th){
        lcost += (y[ii] - lybar) * (y[ii] - lybar);
      }else{
        rcost += (y[ii] - rybar) * (y[ii] - rybar);
      }
    }
    //二乗和の合計を最小化の基準にする
    return lcost + rcost;
  }

  //ノードidにおける最適な分割基準を見つける
  std::pair<int, double> search_best_split(int id,
                 const std::vector< std::vector<double> >& x,
                 const std::vector<double>& y){
    std::pair<int, double> best(-1, INF);
    double best_value = INF;

    int dim = x[0].size();
    for(int j=0; j<dim; j++){
      //閾値候補
      std::vector<double> v;
      for(size_t i=0; i<nodes[id].ids.size(); i++){
        v.push_back(x[nodes[id].ids[i]][j]);
      }
      //変数jと閾値v[i]で分割したときの分割関数の値を計算
      for(size_t i=0; i<v.size(); i++){
        double res = split_value(j, v[i], id, x, y);
        if(best_value > res){
          best_value = res;
          best.first = j;
          best.second = v[i];
        }
      }
    }
    return best;
  }

  //木を成長させる
  void grow_tree(int id,
                 const std::vector< std::vector<double> >& x,
                 const std::vector<double>& y){
    //ノードに割り当てられた学習データのyの平均値
    nodes[id].ybar = 0.0;
    for(size_t i=0; i<nodes[id].ids.size(); i++){
      nodes[id].ybar += y[nodes[id].ids[i]];
    }
    if(nodes[id].ids.size() == 0) return;
    nodes[id].ybar /= nodes[id].ids.size();
    
    //子ノードを準備
    if(2*id + 1 >= N) return;
    nodes[id].ln = 2*id + 1;
    nodes[id].rn = 2*id + 2;

    //最適な分割変数、分割閾値を探索
    std::pair<int, double> best = search_best_split(id, x, y);
    nodes[id].j = best.first;
    nodes[id].th = best.second;
    if(best.first == -1) return;

    //子ノードに学習データを分類する
    for(size_t i=0; i<nodes[id].ids.size(); i++){
      if(x[nodes[id].ids[i]][nodes[id].j] < nodes[id].th){ //閾値未満なら左ノードへ
        nodes[nodes[id].ln].ids.push_back(nodes[id].ids[i]);
      }else{ //閾値以上なら右ノードへ
        nodes[nodes[id].rn].ids.push_back(nodes[id].ids[i]);
      }
    }
    std::vector<int>(nodes[nodes[id].ln].ids).swap(nodes[nodes[id].ln].ids);
    std::vector<int>(nodes[nodes[id].rn].ids).swap(nodes[nodes[id].rn].ids);
  }

  
public:
  RT(int depth_):depth(depth_),INF(1e23){
    N = (1<<depth) - 1;
  }

  void fit(const std::vector< std::vector<double> >& x, const std::vector<double>& y){
    //木の初期化
    nodes.clear();
    nodes.resize(N);
    for(int i=0; i<static_cast<int>(y.size()); i++){
      nodes[0].ids.push_back(i);
    }
    std::vector<int>(nodes[0].ids).swap(nodes[0].ids);

    //木を成長させる
    for(int i=0; i<N; i++){
      grow_tree(i, x, y);
    }

    //木の各ノードの不要な情報を削除
    for(int i=0; i<N; i++) nodes[i].ids.clear();
  }

  double predict(const std::vector<double>& x){
    if(N == 0) return 0;
    int id = 0;
    while(true){
      if(nodes[id].j == -1) return nodes[id].ybar;
      
      if(x[nodes[id].j] < nodes[id].th) id = nodes[id].ln;
      else id = nodes[id].rn;
    }
  }
};

//勾配ブースティング木
class GradientBoostingTree {
  static const int treedepth = 3;
  int iter;
  std::vector<RT> trees;
  double rho; //学習率

  //勾配を計算
  double calc_ngrad(double y, double f){
    //Square loss
    return y - f;

    //Absolute loss
    //return (y-f>0?1:-1);

    //Huber loss
    //double delta = 1.0;
    //if(fabs(y-f) < delta) return y - f;
    //else return (delta * (y-f>0?1:-1));
  }
  
public:
  GradientBoostingTree(int iter, double rho = 1.0):
  iter(iter),trees(iter, RT(treedepth)),rho(rho){
  }
  
  void fit(const std::vector< std::vector<double> >& x, const std::vector<double>& y){
    //最初に普通にtree[0]をfit
    trees[0].fit(x, y);

    //現在までの計算結果をfに保存
    std::vector<double> f(y.size());
    for(size_t i=0; i<y.size(); i++){
      f[i] = trees[0].predict(x[i]);
    }

    //(iter-1)回勾配降下する    
    for(int i=1; i<iter; i++){
      //負の勾配を計算
      std::vector<double> ngrad(y.size(), 0.0);
      for(size_t j=1; j<ngrad.size(); j++){
        ngrad[j] = calc_ngrad(y[j], f[j]);
      }
      //tree[i]でfit
      trees[i].fit(x, ngrad);
      //fit後の予測値で値を更新
      for(size_t j=0; j<y.size(); j++){
        f[j] += rho * trees[i].predict(x[j]);
      }
    }
  }

  double predict(const std::vector<double>& x){
    double ret = trees[0].predict(x);
    for(int i=1; i<iter; i++){
      ret += rho * trees[i].predict(x);
    }
    return ret;
  }
};



int main(int argc, char** argv){
  std::vector< std::vector<double> > train_x, test_x;
  std::vector<double> train_y, test_y;

  {//学習データの読み込み
    std::ifstream trainfs(argv[1]);
    int dim;
    trainfs >> dim;
    double x, y;
    while(trainfs >> y){
      train_y.push_back(y);
      std::vector<double> tmp;
      for(int i=0; i<dim; i++){
        trainfs >> x;
        tmp.push_back(x);
      }
      train_x.push_back(tmp);
    }
  }
  std::cerr << "load file : " << argv[2] << std::endl;
  {//テストデータの読み込み
    std::ifstream testfs(argv[2]);
    int dim;
    testfs >> dim;
    double x, y;
    while(testfs >> y){
      test_y.push_back(y);
      std::vector<double> tmp;
      for(int i=0; i<dim; i++){
        testfs >> x;
        tmp.push_back(x);
      }
      test_x.push_back(tmp);
    }
  }
  std::cerr << "done..." << std::endl;

  //学習
  GradientBoostingTree rt(100, 0.5);

  std::cerr << "fitting..." << std::endl;
  rt.fit(train_x, train_y);
  std::cerr << "done..." << std::endl;

  //予測
  std::cout << "#test_y\tpredict" << std::endl;
  for(size_t i=0; i<test_y.size(); i++){
    std::cout << test_y[i] << "\t" << rt.predict(test_x[i]) << std::endl;
  }
  
  return 0;
}

実験

データ

https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/regression.html#eunite2001
libsvmのregressionデータの「eunite2001」の学習データとテストデータを使ってみてみる。

データの形式

入力の形式は、

次元
y_t x_1 x_2 x_3 ...
...

のように変換して入力。

$ cat ./data/eunite2001.in
16
818.0 0 0 1 0 0 0 0 0 0 0.859223 0.645631 0.589806 0.711165 0.808252 0.759709 0.808252
803.0 0 0 0 1 0 0 0 0 0 0.859223 0.859223 0.645631 0.589806 0.711165 0.808252 0.759709
804.0 0 0 0 0 1 0 0 0 0 0.822816 0.859223 0.859223 0.645631 0.589806 0.711165 0.808252
746.0 0 0 0 0 0 1 0 0 0 0.825243 0.822816 0.859223 0.859223 0.645631 0.589806 0.711165
730.0 0 0 0 0 0 0 0 0 0 0.684466 0.825243 0.822816 0.859223 0.859223 0.645631 0.589806
...

実験

「eunite2001.t」の結果でいくつか比較。
RMSE(Root Mean Squared Error = sqrt( 1/N * Σ(正解値_i - 予測値_i)^2 ))を見てみる。


iterρdepthRMSE
501.0630.8698
1001.0630.8651
1000.5627.9962
1000.5322.2205

それっぽい値がでているのでたぶん大丈夫。
iter(反復回数)よりもρ(学習率)やdepth(回帰木の深さ)の調整の方が効果が出ているよう。

参考文献

2016-06-15

t-SNEで遊ぶ

はじめに

最近よく見かける「t-SNE」という非線形次元圧縮手法を試してみた。

t-SNEとは

  • t-Distributed Stochastic Neighbor Embedding
  • SNEと呼ばれる次元圧縮手法の問題点を改善した手法
    • SNEは、「各点間の"ユークリッド距離"を、類似度に相当する"条件付き確率"に変換して低次元にマッピングする手法」のこと
    • 各点について、高次元での確率分布と低次元での確率分布が一致するように低次元での点を探している
    • 確率分布の違いは「カルバックライブラー情報量」で、これを損失関数にして勾配法で低次元の点を求める
  • 低次元での分布に「自由度1のt-分布」を利用する
  • さらに、高速化を行った「Barnes-Hut t-SNE」という手法ではO(n log n)で計算できる

詳しい説明は、元論文や紹介記事がたくさんある。

また、実用で使う場合は、以下に各種実装へのリンクがまとめられているので、そちらを使うのがよい。
https://lvdmaaten.github.io/tsne/

コード

論文にある「Algorithm 1: Simple version of t-Distributed Stochastic Neighbor Embedding」をそのまま実装したつもり。
論文で紹介されている「early compression」「early exaggeration」はいれていない。

#include <iostream>
#include <vector>
#include <cmath>

// 注意: 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));
}
// 注意: int_maxぐらいで割るべき
double frand(){
  return xor128()%1000000/static_cast<double>(1000000); 
}
static const double PI = 3.14159265358979323846264338;
double normal_rand(double mu, double sigma2){
  double sigma = sqrt(sigma2);
  double u1 = frand(), u2 = frand();
  double z1 = sqrt(-2*log(u1)) * cos(2*PI*u2);
  //double z2 = sqrt(-2*log(u1)) * sin(2*PI*u2);
  return mu + sigma*z1;
}

typedef std::vector<double> P;
static const double EPS = 1e-5;

class tSNE {
  int fromDim, toDim;
  double perp;
  int T;
  double eta;

  double dist2(const P& xi, const P& xj){
    double ret = 0.0;
    for(int i=0; i<xi.size(); i++){
      ret += (xi[i]-xj[i]) * (xi[i]-xj[i]);
    }
    return ret;
  }
  
public:
  tSNE(int fromDim, int toDim = 2, double perp = 5.0, int T = 100000, double eta = 100):
  fromDim(fromDim), toDim(toDim), perp(perp), T(T), eta(eta) { }
  
  std::vector<P> reduce(const std::vector<P>& X){
    int N = X.size();
    std::vector<P> Y(N, P(toDim,0));
    
    //compute pairwise affinities p(j|i) with perplexity perp
    std::vector< std::vector<double> > condP(N, std::vector<double>(N,0));
    for(int i=0; i<N; i++){
      //binary search
      double lb = 0.0, ub = 1e10;
      for(int bi=0; bi<200; bi++){
        double sigma2 = (lb+ub)/2.0;
        
        for(int j=0; j<N; j++){
          if(i!=j){
            condP[j][i] = exp( -dist2(X[i],X[j]) / (2.0 * sigma2) );
          }else{
            condP[j][i] = 0.0;
          }
        }
        double sum = 0.0;      
        for(int k=0; k<N; k++){
          if(i != k) sum += condP[k][i];
        }
        for(int j=0; j<N; j++){
          condP[j][i] /= sum;
        }
        //compute perplexity
        double HPi = 0.0;
        for(int j=0; j<N; j++){
          if(i!=j) HPi -= condP[j][i] * log(condP[j][i]);
        }
        double perpPi = pow(2.0, HPi);
        if(perpPi < perp) lb = sigma2;
        else ub = sigma2;
      }
      std::cerr << "sigma^2 of " << i << " = " << lb << std::endl;
    }

    //set p(i,j)
    std::vector< std::vector<double> > jointP(N, std::vector<double>(N,0));
    for(int i=0; i<N; i++){
      for(int j=0; j<N; j++){
        if(i!=j) jointP[i][j] = (condP[j][i] + condP[i][j]) / (2.0 * N);
      }
    }
    
    //sample initial solution Y^(0)
    for(int i=0; i<N; i++){
      for(int j=0; j<toDim; j++){
        Y[i][j] = normal_rand(0.0, 0.001);
      }
    }
    
    //iteration
    std::vector<P> diff(N, P(toDim, 0));
    double alpha = 0.5;
    
    for(int t=1; t<=T; t++){
      std::vector<P> newY = Y;
      //compute low-dimensional affinities q(i,j)
      std::vector< std::vector<double> > jointQ(N, std::vector<double>(N,0));
      double sum = 0.0;
      for(int k=0; k<N; k++){
        for(int l=0; l<N; l++){
          if(k!=l) sum += 1.0 / (1.0 + dist2(Y[k], Y[l]));
        }
      }
      for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
          if(i!=j){
            jointQ[i][j] = 1.0 / (1.0 + dist2(Y[i], Y[j]));
            jointQ[i][j] /= sum;
          }
        }
      }

      //(compute cost C = KL(P||Q))
      double cost = 0.0;
      for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
          if(i!=j){
            cost += jointP[i][j] * log( jointP[i][j] / jointQ[i][j] );
          }
        }
      }
      
      //compute gradient
      std::vector<P> delta(N, P(toDim, 0));
      for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
          if(i!=j){
            double a = (jointP[i][j] - jointQ[i][j]) * (1.0 / (1.0 + dist2(Y[i], Y[j])));
            for(int k=0; k<toDim; k++){
              delta[i][k] += 4.0 * a * (Y[i][k] - Y[j][k]);
            }
          }
        }
      }

      //set Y^(t)
      for(int i=0; i<N; i++){
        for(int j=0; j<toDim; j++){
          newY[i][j] = Y[i][j] - eta * delta[i][j] + alpha * diff[i][j];
          diff[i][j] = newY[i][j] - Y[i][j];
        }
      }

      Y = newY;
      
      if(t%100==0){
        std::cerr << "t = " << t << " : cost = " << cost << std::endl;
      }
      
      if(t >= 250) alpha = 0.8;
    }
    return Y;
  }
};


int main(){
  int N, fromDim, toDim;
  std::cin >> N;
  std::cin >> fromDim >> toDim;

  //input
  std::vector<P> X;
  for(int i=0; i<N; i++){
    P p(fromDim, 0);
    for(int j=0; j<fromDim; j++){
      std::cin >> p[j];
    }
    X.push_back(p);
  }

  //compute
  tSNE tsne(fromDim, toDim, 5.0, 100000, 100.0);
  std::vector<P> Y = tsne.reduce(X);

  //output
  for(int i=0; i<N; i++){
    for(int j=0; j<toDim; j++){
      if(j!=0) std::cout << "\t";
      std::cout << Y[i][j];
    }
    std::cout << std::endl;
  }
  return 0;
}

実験

MNISTの「test」のデータから500件分をプロット。各要素は0〜255だったようなので、256で割った値を入力として利用。
入力の形式は、

データ点数
高次元での次元数 低次元での次元数
高次元でのデータ点1(タブかスペース区切り、高次元での次元数分)
高次元でのデータ点2(タブかスペース区切り、高次元での次元数分)
...

にしている。

$ cat in.txt
500
784 2
0       0       0       0       0       0       0       0 ...
0       0       0       0       0       0       0       0 ...
...

最終状態は以下。一応収束してくれている。
f:id:jetbead:20160615022815p:image:w640


ステップ5000ずつの低次元マッピングの状態をプロットしてみる。mnistの数字(0〜9)ごとに色や記号の形をわけている。
f:id:jetbead:20160615021706g:image
最初のバラバラな状態からすぐに塊ができて、それらが微調整されていっている。
効率化テクニックを入れていないせいか、どんどん大きくなっているし、10万回回してもまだ収束しきっていないようにも見える。。。

参考文献

2015-07-20

GP-MIで遊ぶ

はじめに

http://live.nicovideo.jp/watch/lv228162988
先週のNL研のニコ生で、ベイズ的最適化についての招待講演を見ていて「SEは滑らかすぎる」という発言がよくわからなかったので、GP-MIを試してみる。


Contal et al., Gaussian Process Optimization with Mutual Information
http://arxiv.org/abs/1311.4825
http://econtal.perso.math.cnrs.fr/presentations/slides_icml14.pdf

コード

カーネルは、通常(?)のSquared Exponential kernel。
逆行列求めるところが重いのでEigen使う版も一緒に。
いつものごとく雑。

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

static const double PI = 3.14159265358979323846264338;

//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)の一様乱数
double frand(){
  return xor128()%1000000/static_cast<double>(1000000); 
}
//正規乱数
double normal_rand(double mu, double sigma2){
  double sigma = sqrt(sigma2);
  double u1 = frand(), u2 = frand();
  double z1 = sqrt(-2*log(u1)) * cos(2*PI*u2);
  //double z2 = sqrt(-2*log(u1)) * sin(2*PI*u2);
  return mu + sigma*z1;
}

#ifdef USE_EIGEN
#include <Eigen/Core>
#include <Eigen/LU>
using namespace Eigen;
#else
//行列計算用
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;
}
#endif



class GP_MI {  
  double d; //xの次元
  double alpha; //活用探索のトレードオフ調整用パラメータ
  double noise_sigma2; //ノイズの分散
  int ITER; //サンプリング回数
  double gamma;

  bool flg_obs; //推定,観測を繰り返しているかどうかの確認用フラグ

  //観測データ
  std::vector< std::vector<double> > obs_x;
  std::vector<double> obs_y;

  //各x_iの探索範囲
  std::vector<double> x_min, x_max;

  //kernel(x,x')の値を求める
  double calc_kernel(const std::vector<double>& x1, const std::vector<double>& x2){
    //Squared exponential kernel
    double sigma2 = 1.0;
    double l2 = 0.01; //比較的ガタガタを試すので小さめ

    double z = 0.0;
    for(size_t i=0; i<x1.size(); i++){
      z += pow(((x1[i] - x2[i])), 2.0) / (2.0 * l2);
    }
    return sigma2 * exp(-z);
  }

public:
  GP_MI(int d, double alpha, double noise_sigma2, int ITER):
    d(d),
    alpha(alpha),
    noise_sigma2(noise_sigma2),
    gamma(0.0),
    x_min(d, 0.0),
    x_max(d, 1.0),
    ITER(ITER),
    flg_obs(false)
  { }

  //xの探索範囲を指定
  void set_x_range(int di, double mn, double mx){
    x_min[di] = mn;
    x_max[di] = mx;
  }

  //現在までの観測結果から推定される関数をもとに、一番探した方がよいxを返す
  std::vector<double> get_best_x(){
    //前回の観測が行われたかチェック
    if(!flg_obs){
      std::cerr << "please observe and set" << std::endl;
      return std::vector<double>();
    }
    flg_obs = false;
    
    //前計算C_T, Y_T
#ifdef USE_EIGEN
    MatrixXd C_T = MatrixXd::Zero(obs_y.size(), obs_y.size());
    MatrixXd Y_T = MatrixXd::Zero(obs_y.size(), 1);
    for(size_t i=0; i<obs_y.size(); i++){
      for(size_t j=0; j<obs_y.size(); j++){
        C_T(i,j) = calc_kernel(obs_x[i], obs_x[j]) + (i==j?noise_sigma2:0.0);
      }
    }
    for(size_t i=0; i<obs_y.size(); i++){
      Y_T(i,0) = obs_y[i];
    }
    MatrixXd C_T_inv = C_T.inverse();
#else
    Matrix C_T(obs_y.size(), obs_y.size());
    Matrix Y_T(obs_y.size(), 1);
    for(size_t i=0; i<obs_y.size(); i++){
      for(size_t j=0; j<obs_y.size(); j++){
        C_T.setVal(i,j, calc_kernel(obs_x[i], obs_x[j]) + (i==j?noise_sigma2:0.0) );
      }
    }
    for(size_t i=0; i<obs_y.size(); i++){
      Y_T.setVal(i,0, obs_y[i]);
    }
    Matrix C_T_inv = C_T.inverse();
#endif

    //ITER回サンプリングしてargmaxを求める
    double best_a = -1.0e10;
    std::vector<double> best_x;
    double best_sigma2 = 0;

    for(size_t iter=0; iter<ITER; iter++){
      //適当なxを生成
      std::vector<double> x(d);
      for(size_t i=0; i<d; i++) x[i] = (x_max[i]-x_min[i])*frand() + x_min[i];
      
#ifdef USE_EIGEN
      MatrixXd kT = MatrixXd::Zero(obs_y.size(), 1);
      for(size_t i=0; i<obs_y.size(); i++){
        kT(i,0) = calc_kernel(obs_x[i], x);
      }

      double kxx = calc_kernel(x, x);

      MatrixXd kCY = kT.transpose() * ( C_T_inv * Y_T );
      MatrixXd kCk = kT.transpose() * ( C_T_inv * kT );

      double mu = kCY(0,0);
      double sigma2 = kxx - kCk(0,0);
#else
      Matrix kT(obs_y.size(), 1);
      for(size_t i=0; i<obs_y.size(); i++){
        kT.setVal(i,0, calc_kernel(obs_x[i], x));
      }

      double kxx = calc_kernel(x, x);

      Matrix kCY = kT.transpose() * ( C_T_inv * Y_T );
      Matrix kCk = kT.transpose() * ( C_T_inv * kT );

      double mu = kCY.getVal(0,0);
      double sigma2 = kxx - kCk.getVal(0,0);
#endif

      //獲得関数の値を求める
      double phi = sqrt(alpha) * (sqrt(sigma2 + gamma) - sqrt(gamma));      
      if(best_a < mu + phi){
        best_a = mu + phi;
        best_x = x;
        best_sigma2 = sigma2;
      }
    }

    gamma += best_sigma2;
    return best_x;
  }

  //get_best_x()で得られたxに対する実際の関数の値を、観測結果として登録
  void set_observation(std::vector<double>& x, double y){
    obs_x.push_back(x);
    obs_y.push_back(y);
    flg_obs = true;
  }

};



//ブラックボックスな関数
double func(const std::vector<double>& x){
  return x[0] * sin(10.0 * x[0]);
}

int main(){
  //xの次元数  alpha  noise_sigma2  サンプリング回数
  GP_MI gpmi(1, 10.0, 0.0001, 500);

  //各次元の範囲
  gpmi.set_x_range(0, 0.0, 1.0);

  //最初に何点か求めておく
  std::vector<double> x(1);
  x[0] = 0;
  gpmi.set_observation(x, func(x));
  x[0] = 0.5;
  gpmi.set_observation(x, func(x));
  x[0] = 1.0;
  gpmi.set_observation(x, func(x));

  //自動で探索点を選ぶ
  for(int i=0; i<30; i++){
    std::vector<double> best_x = gpmi.get_best_x();
    std::cout << i << " : best_x = " << best_x[0] << std::endl;
    gpmi.set_observation(best_x, func(best_x));
  }

  return 0;
}

結果

最初に0,0.5,1.0の3点だけ求めた状態からスタート。
αがlog(2/δ), 0<δ<1とかだけど、小さいと隣の山に探しに行ってくれないので、大きいけど、α=10と50で試してみる。
グラフは、赤がμとσ^2、緑が獲得関数μ+φ。

x * sin(10.0 * x)

【α=10】
f:id:jetbead:20150720123215g:image
最適解は求められているけど、左の山が分散大きくても獲得関数の値が低いまま。
【α=50】
f:id:jetbead:20150720123719g:image
αを大きくする(大きくしすぎ?)と、左の山も探してくれてる。

x * cos(30.0 * x)

もっとガタガタしたもの。
【α=10】
f:id:jetbead:20150720123915g:image
うまくいっていない(本当の最適点はx=0.84ぐらい)。
動画で「滑らかすぎる」という発言があったけど、滑らかすぎて、関数の形の推定がうまくいかない=適切な探索点をみつけられない、ということか(?)。
ただ、ARDなSEカーネルやMaternカーネルの実装、最適なパラメータ探索(MCMC)がちょっとよくわかってない・・・


http://www.slideshare.net/issei_sato/deep-learning-41310617
> BOツール公開予定
|ω・)待機。

参考

2015-07-12

Elman netを試す

はじめに

プロフェッショナルな「深層学習」本で紹介されているRNNの一種のElman netを書いてみる。

Recurrent Neural Network(RNN)とは

  • 再帰型ニューラルネット
  • ネットワーク内部に内部有向閉路を持つようなニューラルネットの総称
    • Feedforwardの時は、入力層から出力層へ一方方向
  • この構造のおかげで、時系列や言語モデル、系列ラベリングなど前の状態を持つような問題に対して考慮できる
  • いろんな種類がある(以下はwikipediaから)
    • Fully Recurrent network
      • Hopfield network
      • Boltzmann machine
    • Simple recurrent network
      • Elman net
      • Jordan net
    • Echo state network
    • Long short term memory(LSTM) network
    • Bi-directional RNN
    • Hierarchical RNN
    • Stochastic neural networks

Simple recurrent network
  • 言葉のとおり、シンプルな構造のRNN
  • 代表的には「Elman net」と「Jordan net」がある
  • 3層(入力層、隠れ層、出力層)のFeedforward Neural netについて、どこからどこへ閉路を作るかで異なる
    • Elman net : 隠れ層から隠れ層への辺がある
    • Jordan net : 出力層から隠れ層への辺がある

コード

「深層学習」本では、Elman netのBPTT法(Backpropagation through time、時間でネットワークを分けて普通に逆誤差伝播)での学習が紹介されていると思われるので、これを書いてみる。

#include <iostream>
#include <vector>
#include <cstdio>
#include <algorithm>
#include <cmath>

static const double PI = 3.14159265358979323846264338;

//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)の一様乱数
double frand(){
  return xor128()%1000000/static_cast<double>(1000000); 
}
//正規乱数
double normal_rand(double mu, double sigma2){
  double sigma = sqrt(sigma2);
  double u1 = frand(), u2 = frand();
  double z1 = sqrt(-2*log(u1)) * cos(2*PI*u2);
  //double z2 = sqrt(-2*log(u1)) * sin(2*PI*u2);
  return mu + sigma*z1;
}


//隠れ層が1層の単純なRNN
// BPTT(Backpropagation through time)法で学習
class SimpleRNN {
  double eps; //学習率
  int Nin, Nhid, Nout; //各層のユニット数(バイアスを除く)
  std::vector< std::vector<double> > W; //時刻t-1での隠れ層から時刻tでの隠れ層へ
  std::vector< std::vector<double> > Win; //時刻tでの入力層から隠れ層へ
  std::vector< std::vector<double> > Wout; //時刻tでの隠れ層から出力層へ

  std::vector< std::vector<double> > u, v; //各時刻での隠れ層の入力, 出力層の入力
  std::vector< std::vector<double> > z, y; //各時刻での隠れ層の出力, 出力層の出力
  std::vector< std::vector<double> > delta, delta_out; //各時刻での各ユニットのδ値

public:
  SimpleRNN(int Nin, int Nhid, int Nout, double eps):
    eps(eps),
    Nin(Nin),
    Nhid(Nhid),
    Nout(Nout),
    W(Nhid+1, std::vector<double>(Nhid, 0)),
    Win(Nin+1, std::vector<double>(Nhid, 0)),
    Wout(Nhid+1, std::vector<double>(Nout, 0))
  {
    for(int i=0; i<Nhid+1; i++){
      for(int j=0; j<Nhid; j++){
        W[i][j] = normal_rand(0.0, 0.1);
      }
    }
    for(int i=0; i<Nin+1; i++){
      for(int j=0; j<Nhid; j++){
        Win[i][j] = normal_rand(0.0, 0.1);
      }
    }
    for(int i=0; i<Nhid+1; i++){
      for(int j=0; j<Nout; j++){
        Wout[i][j] = normal_rand(0.0, 0.1);
      }
    }
  }

  std::vector<int> forward_propagation(const std::vector< std::vector<double> >& in){
    int T = in.size();
    
    u = std::vector< std::vector<double> >(T, std::vector<double>(Nhid, 0));
    v = std::vector< std::vector<double> >(T, std::vector<double>(Nout, 0));
    z = std::vector< std::vector<double> >(T, std::vector<double>(Nhid+1, 0));
    y = std::vector< std::vector<double> >(T, std::vector<double>(Nout, 0));

    //各時刻tでの値を求める
    std::vector<int> ret;

    for(int t=0; t<T; t++){
      for(int i=0; i<Nhid; i++){
        u[t][i] = 0.0;
      }
      //入力層の出力->隠れ層の入力
      for(int i=0; i<Nin; i++){
        for(int j=0; j<Nhid; j++){
          u[t][j] += in[t][i] * Win[i][j];
        }
      }
      for(int i=0; i<Nhid; i++){ //バイアス
        u[t][i] += 1.0 * Win[Nin][i];
      }
      //時刻t-1での隠れ層の出力->時刻tでの隠れ層の入力
      for(int i=0; i<Nhid+1; i++){
        if(t!=0){
          for(int j=0; j<Nhid; j++){
            u[t][j] += z[t-1][i] * W[i][j];
          }
        }
        if(t==0 && i==Nhid){ //バイアス
          for(int j=0; j<Nhid; j++){
            u[t][j] += 1.0 * W[Nin][j];
          }
        }
      }
      //時刻tでの隠れ層の出力
      for(int i=0; i<Nhid; i++){
        z[t][i] = 1.0 / (1.0 + exp(-u[t][i]));
      }
      z[t][Nhid] = 1.0;
      //時刻tでの出力層の入力
      for(int i=0; i<Nhid+1; i++){
        for(int j=0; j<Nout; j++){
          v[t][j] += z[t][i] * Wout[i][j];
        }
      }
      //時刻tでの出力層の出力
      double Z = 0;
      double mx = -1.0;
      int mx_i = -1;
      for(int i=0; i<Nout; i++){
        Z += exp(v[t][i]);
      }
      for(int i=0; i<Nout; i++){
        y[t][i] = exp(v[t][i]) / Z;

        if(mx < y[t][i]){
          mx = y[t][i];
          mx_i = i;
        }
      }
      ret.push_back(mx_i);
    }
    return ret;
  }

  double back_propagation(const std::vector< std::vector<double> >& in, const std::vector<int>& out){
    double err = 0.0;
    int T = in.size();
    std::vector<int> res = forward_propagation(in);

    delta = std::vector< std::vector<double> >(T, std::vector<double>(Nhid, 0));
    delta_out = std::vector< std::vector<double> >(T, std::vector<double>(Nout, 0));

    //時刻の逆方向からdelta値と重みの更新を連続的にしていく
    for(int t=T-1; t>=0; t--){
      //時刻tでの出力層のdelta値
      for(int i=0; i<Nout; i++){
        if(out[t] == i){
          delta_out[t][i] = y[t][i] - 1.0;
          err += - 1.0 * log(y[t][i]);
        }else{
          delta_out[t][i] = y[t][i] - 0.0;
        }
      }
      //時刻tでの隠れ層の出力層からのdelta値
      for(int i=0; i<Nhid; i++){
        for(int j=0; j<Nout; j++){
          double fu = 1.0 / (1.0 + exp(-u[t][i]));
          double fdash = fu * (1.0 - fu);
          delta[t][i] += Wout[i][j] * delta_out[t][j] * fdash;
        }
      }
      //時刻tでの隠れ層の時刻t+1での隠れ層からのdelta値
      for(int i=0; i<Nhid; i++){
        for(int j=0; j<Nhid; j++){
          if(t+1>=T) continue; //時刻T以降の場合はdelta=0として扱う

          double fu = 1.0 / (1.0 + exp(-u[t][i]));
          double fdash = fu * (1.0 - fu);
          delta[t][i] += W[i][j] * delta[t+1][j] * fdash;
        }
      }
      //delta値を使って各重みの微分を計算し、SGDで重みを更新
      for(int j=0; j<Nhid+1; j++){ //隠れ層->出力層
        for(int k=0; k<Nout; k++){
          Wout[j][k] -= eps * (delta_out[t][k] * z[t][j]);
        }
      }
      for(int pj=0; pj<Nhid+1; pj++){ //隠れ層->隠れ層
        for(int j=0; j<Nhid; j++){
          if(t-1>=0){
            W[pj][j] -= eps * (delta[t][j] * z[t-1][pj]);
          }else{ //時刻t<0の場合は、バイアスだけ1でそれ以外の出力は0として扱う
            if(pj==Nhid){ //バイアス
              W[pj][j] -= eps * (delta[t][j] * 1.0);
            }
          }
        }
      }
      for(int i=0; i<Nin+1; i++){ //入力層->隠れ層
        for(int j=0; j<Nhid; j++){
          if(i<Nin){
            Win[i][j] -= eps * (delta[t][j] * in[t][i]);
          }else{ //バイアス
            Win[i][j] -= eps * (delta[t][j] * 1.0);
          }
        }
      }

    }
    return err;
  }

};

int main(){
  const int TERM = 500; //TERM個単位でerrを確認
  const double finish_err = 0.01; //errがfinish_err未満になったら学習終了

  int Nin, Nhid, Nout;
  double eps;
  int trainN, testN;

  std::cin >> Nin >> Nhid >> Nout;
  std::cin >> eps;
  
  /// 学習 ////////////////////
  std::cin >> trainN;
  std::vector< std::vector< std::vector<double> > > in;
  std::vector< std::vector<int> > out;
  for(int i=0; i<trainN; i++){
    int t;
    std::cin >> t;

    std::vector< std::vector<double> > in_one;
    std::vector<int> out_one;
    for(int j=0; j<t; j++){
      double val;
      int outval;
      std::cin >> outval;
      std::vector<double> v;
      for(int k=0; k<Nin; k++){
        std::cin >> val;
        v.push_back(val);
      }
      in_one.push_back(v);
      out_one.push_back(outval);
    }

    in.push_back(in_one);
    out.push_back(out_one);
  }
  //データシャッフル用
  std::vector<int> rnd;
  for(int i=0; i<trainN; i++){
    rnd.push_back(i);
  }
  std::random_shuffle(rnd.begin(), rnd.end());

  //学習ループ
  SimpleRNN rnn(Nin, Nhid, Nout, eps);
  int iter = 0;
  double err = 0.0;
  while(true){
    //TERM個単位で確認
    if(iter!=0 && iter%TERM == 0){
      err /= TERM;
      std::cerr << "err = " << err << std::endl;
      if(err < finish_err) break;
      err = 0;
    }

    err += rnn.back_propagation(in[rnd[iter%trainN]], out[rnd[iter%trainN]]);
    iter++;
  }

  /// テスト ////////////////////
  std::cin >> testN;
  for(int i=0; i<testN; i++){
    int t;
    std::cin >> t;

    std::vector< std::vector<double> > in_one;
    std::vector<int> out_one;
    for(int j=0; j<t; j++){
      double val;
      int outval;
      std::cin >> outval;
      std::vector<double> v;
      for(int k=0; k<Nin; k++){
        std::cin >> val;
        v.push_back(val);
      }
      in_one.push_back(v);
      out_one.push_back(outval);
    }

    //結果の出力
    std::vector<int> res = rnn.forward_propagation(in_one);
    std::cout << "ref:";
    for(int j=0; j<out_one.size(); j++){
      std::cout << " " << out_one[j];
    }
    std::cout << std::endl;

    std::cout << "out:";
    for(int j=0; j<res.size(); j++){
      std::cout << " " << res[j];
    }
    std::cout << std::endl;

    std::cout << std::endl;
  }

  return 0;
}

入力ファイルの形式

[入力層のユニット数]  [隠れ層のユニット数]  [出力層のユニット数]
[学習率]
[学習インスタンス数]
[ケース1の数]
[出力1]  [次元1の値]  [次元2の値] ...
[出力2]  [次元1の値]  [次元2の値] ...
[出力3]  [次元1の値]  [次元2の値] ...
...
[テストインスタンス数]
[ケース1の数]
[出力1]  [次元1の値]  [次元2の値] ...
[出力2]  [次元1の値]  [次元2の値] ...
[出力3]  [次元1の値]  [次元2の値] ...
...

結果

おもちゃなケースで確認だけ。
http://d.hatena.ne.jp/jetbead/20111125/1322236626 でのサンプルを今回の形式に変換したもの。

  • 入力
10 10 3
0.01
6
3
2  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4
2  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
1  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3
2  0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
2  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4
2  0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
0  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
1  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3
2  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
2  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
4
2  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
6
3
2  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4
2  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
1  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3
2  0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
2  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4
2  0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
0  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
1  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3
2  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
2  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
4
2  0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
2  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
  • 出力
$ ./a.out < in
err = 2.20794
err = 1.84924
err = 1.54379
err = 1.22335
err = 0.931081
err = 0.681177
err = 0.493641
err = 0.360733
err = 0.265542
err = 0.198631
err = 0.153081
err = 0.121373
err = 0.0987013
err = 0.0824887
err = 0.0702189
err = 0.0606353
err = 0.0533114
err = 0.0473444
err = 0.0423609
err = 0.0383958
err = 0.034996
err = 0.0320187
err = 0.0295968
err = 0.0274426
err = 0.0254882
err = 0.0238809
err = 0.0224107
err = 0.0210389
err = 0.0199064
err = 0.0188471
err = 0.0178354
err = 0.0170008
err = 0.0162053
err = 0.0154302
err = 0.0147935
err = 0.0141766
err = 0.0135647
err = 0.0130652
err = 0.0125743
err = 0.0120795
err = 0.0116788
err = 0.0112797
err = 0.0108714
err = 0.010544
err = 0.0102138
err = 0.00987132
ref: 2 2 2
out: 2 2 2

ref: 2 2 0 1
out: 2 2 0 1

ref: 2 2 2
out: 2 2 2

ref: 2 2 0 1
out: 2 2 0 1

ref: 2 2 2
out: 2 2 2

ref: 2 2 2 2
out: 2 2 2 2

ちゃんと収束して学習データはちゃんとラベルつけられている。
しかし、展開する時間tが長くなるほどネットワークが深く(deep)なってしまい、勾配消失(または発散)問題が起きてしまってうまく学習できない。

参考