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】

最適解は求められているけど、左の山が分散大きくても獲得関数の値が低いまま。
【α=50】

αを大きくする(大きくしすぎ？)と、左の山も探してくれてる。

x * cos(30.0 * x)

もっとガタガタしたもの。
【α=10】

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


http://www.slideshare.net/issei_sato/deep-learning-41310617
> BOツール公開予定

ω・)待機。

参考