Negative/Positive Thinking

RBMで遊ぶ

機械学習プログラミング

はじめに

深イイ学習とかで使われているらしいRestricted Boltzmann Machineでちょっと遊んでみる。

Restricted Boltzmann Machineとは

制約Boltzmann Machine
各層内のノード間の結合がないようなBoltzmann Machine
- この制約によって、学習時の隠れ層の各ノードが独立に計算できる

RBMの学習は、「Contrastive Divergence Learning」で行うことができる
- http://www.cs.toronto.edu/~hinton/absps/cdmiguel.pdf
- http://dl.dropbox.com/u/2048288/RestrictedBoltzmannMachine.pdf

コード

できたモデルからどうやって生成するのかよくわからなかったので、適当にサンプリングしてみる。

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

using namespace std;

//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_maxぐらいで割った方がよい
double frand(){
  return xor128()%1000000/static_cast<double>(1000000); 
}



//binary Restricted Boltzmann Machine
class RBM {
  static const double alpha = 0.01; //学習率

  int N_visible, N_hidden; //ノード数
  std::vector<int> state_hidden; //隠れ素子の状態
  std::vector<double> bias_visible, bias_hidden; //バイアス
  std::vector< std::vector<double> > weights; //ウェイトパラメータ
  
  double sigmoid(double x){
    return 1.0 / (1.0 + exp(-x));
  }

  //内部状態から可視素子v_iの状態{0,1}を決定する
  int rand01_for_visible_state(int i){
    double sum = bias_visible[i];
    for(size_t j=0; j<N_hidden; j++){
      sum += weights[i][j] * state_hidden[j];
    }
    if(frand() < sigmoid(sum)) return 1;
    return 0;
  }

  //可視素子の状態と内部状態から隠れ素子h_jの状態{0,1}を決定する
  int rand01_for_hidden_state(const std::vector<int> &state_visible, int j){
    double sum = bias_hidden[j];
    for(size_t i=0; i<N_visible; i++){
      sum += state_visible[i] * weights[i][j]; 
    }
    if(frand() < sigmoid(sum)) return 1;
    return 0;
  }

  //Reconstructし、勾配を求めるための計算をする
  double calc(int n,
              const std::vector<int>& data, 
              std::vector< std::vector<double> >& CD_pos,
              std::vector< std::vector<double> >& CD_neg,
              std::vector<double>& bias_visible_pos,
              std::vector<double>& bias_visible_neg,
              std::vector<double>& bias_hidden_pos,
              std::vector<double>& bias_hidden_neg
              ){
    //dataからn回reconstractしたものを利用する
    std::vector<int> model = data;
    for(int t=0; t<n; t++){
      for(size_t j=0; j<N_hidden; j++){
        state_hidden[j] = rand01_for_hidden_state(model, j);
      }
      for(size_t i=0; i<N_visible; i++){
        model[i] = rand01_for_visible_state(i);
      }
    }

    //bias勾配計算用
    for(size_t i=0; i<N_visible; i++){
      bias_visible_pos[i] += data[i];
      bias_visible_neg[i] += model[i];
    }
    for(size_t j=0; j<N_hidden; j++){
      double sum;
      sum = bias_hidden[j];
      for(size_t k=0; k<N_visible; k++){
        sum += data[k] * weights[k][j];
      }
      bias_hidden_pos[j] += sigmoid(sum);

      sum = bias_hidden[j];
      for(size_t k=0; k<N_visible; k++){
        sum += model[k] * weights[k][j];
      }
      bias_hidden_neg[j] += sigmoid(sum);
    }
    
    //weights勾配計算用
    for(size_t i=0; i<N_visible; i++){
      for(size_t j=0; j<N_hidden; j++){
        double sum;
        //CD_posの計算
        sum = bias_hidden[j];
        for(size_t k=0; k<N_visible; k++){
          sum += data[k] * weights[k][j];
        }
        CD_pos[i][j] += data[i] * sigmoid(sum);
    
        //CD_negの計算
        sum = bias_hidden[j];
        for(size_t k=0; k<N_visible; k++){
          sum += model[k] * weights[k][j];
        }
        CD_neg[i][j] += model[i] * sigmoid(sum);
      }
    }
    
    //誤差計算
    double error = 0;
    for(size_t i=0; i<N_visible; i++){
      error += (data[i]-model[i]) * (data[i]-model[i]);
    }
    return error;
  }

public:
  RBM(int n_visible, int n_hidden):
    N_visible(n_visible),
    N_hidden(n_hidden),
    state_hidden(n_hidden),
    bias_visible(n_visible, 1.0),
    bias_hidden(n_hidden, 1.0),
    weights(n_visible, std::vector<double>(n_hidden, 0))
  {
    //weightsは適当な値で初期化しておく
    for(size_t i=0; i<n_visible; i++){
      for(size_t j=0; j<n_hidden; j++){
        weights[i][j] = frand() * 0.2;
      }
    }
  }

  //学習
  //  data:学習データ  epoch:学習の繰り返し回数  n:GibbsSamplingの回数
  void train(const std::vector< std::vector<int> >& data, int epoch, int n){
    //Contrastive Divergence Training
    for(size_t t=0; t<epoch; t++){
      //勾配計算用
      std::vector< std::vector<double> > CD_pos(N_visible, std::vector<double>(N_hidden,0));
      std::vector< std::vector<double> > CD_neg(N_visible, std::vector<double>(N_hidden,0));
      std::vector<double> bias_visible_pos(N_visible,0), bias_visible_neg(N_visible,0);
      std::vector<double> bias_hidden_pos(N_hidden,0), bias_hidden_neg(N_hidden,0);
      //誤差確認用
      double error = 0.0;

      //データセットに対して、CDやbaisの勾配を求めるための値を計算
      for(size_t i=0; i<data.size(); i++){
        error += calc(n, data[i],
                      CD_pos, CD_neg,
                      bias_visible_pos, bias_visible_neg,
                      bias_hidden_pos, bias_hidden_neg
                      );
      }

      //CD_pos合計、CD_neg合計をデータセット数で割って平均値を求める
      for(size_t i=0; i<N_visible; i++){
        for(size_t j=0; j<N_hidden; j++){
          CD_pos[i][j] /= data.size();
          CD_neg[i][j] /= data.size();
        }
      }
      
      //bias_*_pos合計, bias_*_neg合計をデータセット数で割って平均値を求める
      for(size_t i=0; i<N_visible; i++){
        bias_visible_pos[i] /= data.size();
        bias_visible_neg[i] /= data.size();
      }
      for(size_t j=0; j<N_hidden; j++){
        bias_hidden_pos[j] /= data.size();
        bias_hidden_neg[j] /= data.size();        
      }
      
      //weightsを更新(最急降下法: W' = W + alpha * CD = W + alpha * (CD_pos - CD_neg))
      for(size_t i=0; i<N_visible; i++){
        for(size_t j=0; j<N_hidden; j++){
          weights[i][j] += alpha * (CD_pos[i][j] - CD_neg[i][j]);
        }
      }
      
      //baisを更新(weightsと同じ)
      for(size_t i=0; i<N_visible; i++){
        bias_visible[i] += alpha * (bias_visible_pos[i] - bias_visible_neg[i]);
      }
      for(size_t j=0; j<N_hidden; j++){
        bias_hidden[j] += alpha * (bias_hidden_pos[j] - bias_hidden_neg[j]);
      }
      
      //std::cout << t << " : " << error/data.size() << std::endl;
    }
  }

  //与えられた可視素子から隠れ素子の状態を更新する
  void set_hidden_node(const std::vector<int>& data){
    //隠れ素子を更新しておく
    for(size_t j=0; j<N_hidden; j++){
      state_hidden[j] = rand01_for_hidden_state(data, j);
    }
  }

  //現在の隠れ素子の状態から、p(v|h)を用いて可視素子の状態を求め、さらに隠れ素子の状態を更新しておく
  std::vector<int> get_visible_node(){
    std::vector<int> ret(N_visible);
    //可視素子を更新
    for(size_t i=0; i<N_visible; i++){
      ret[i] = rand01_for_visible_state(i);
    }
    //隠れ素子を更新
    for(size_t j=0; j<N_hidden; j++){
      state_hidden[j] = rand01_for_hidden_state(ret, j);
    }
    return ret;
  }
};



int main(){

  int N_visible, N_hidden; //可視層のノード数、隠れ層のノード数
  cin >> N_visible >> N_hidden;

  //RBMのインスタンス生成
  RBM rbm(N_visible, N_hidden);

  //学習データ
  int N_testcase;  
  cin >> N_testcase;
  std::vector< std::vector<int> > train_data;
  for(size_t t=0; t<N_testcase; t++){
    std::vector<int> dat(N_visible);
    for(size_t i=0; i<N_visible; i++){
      cin >> dat[i];
    }
    train_data.push_back(dat);
  }

  //学習
  rbm.train(train_data, 3000, 1);


  rbm.set_hidden_node(train_data[0]); //最初の可視ノードは学習データの最初のやつを指定しておく
  for(int t=0; t<500; t++){
    //サンプリング
    std::vector<int> ret = rbm.get_visible_node();
    //表示
    for(size_t i=0; i<ret.size(); i++){
      std::cout << ret[i] << " ";
    }
    std::cout << std::endl;
  }
        
  return 0;
}

結果

適当な例

入力

出力

まぁまぁできてるっぽい。

さらに適当な例

入力

出力

え、こういう入力でも学習できるの

参考文献