Negative/Positive Thinking

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(回帰木の深さ)の調整の方が効果が出ているよう。

参考文献

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/jetbead/20161005/1475597438