Negative/Positive Thinking

2016-11-08

Inside-Outsideアルゴリズムを試す

はじめに

確率文脈自由文法での生成規則の適用確率の推定アルゴリズムで紹介されている「Inside-Outsideアルゴリズム」について、Webで検索してみても、最尤導出の構文木や内側確率の計算例などはあっても、外側確率や生成確率の推定などまで計算例を書いてくれているのはなさそうだったので、手計算確認用にプログラムを書いた。

Inside-Outsideアルゴリズムとは

  • 内側・外側アルゴリズム
  • 確率文脈自由文法の生成規則(チョムスキー標準形)の確率値をコーパスを求める際に、内側確率βと外側確率αを導入することで効率よく求められる
    • 隠れマルコフモデルにおける前向き・後ろ向きアルゴリズムに似た感じ
    • 内側確率 : 非終端記号Aから終端記号列w_i^jが生成される確率
    • 外側確率 : 導出中に出現するAについて、Aが支配しているw_i^j以外の両側の終端記号列w_1^{i-1}とw_{j+1}^Nが現れる確率
  • 内側確率と外側確率を使って各生成規則の適用回数を推定でき、これを用いて生成規則の確率値を求める、ということを繰り返す

詳しくは「北, 言語と計算4 確率的言語モデル, 東京大学出版会」の第5章を参照。

コード

内側確率と外側確率、適用回数の推定が確認できればよいと思って書いたため、だいたい愚直。
最尤導出と内側確率が上の書籍と同じ値であるのと、繰り返し最適化で対数尤度が収束しているようなので、たぶん大丈夫。。。
(雑な確認しかしていないので、何かありましたら教えていただければと思います)

#include <iostream>
#include <vector>
#include <map>
#include <string>
#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); 
}

double logsumexp(const std::vector<double>& lps){
  if(lps.size() == 0) return -1e300;
  double mx = lps[0];
  for(size_t i=0; i<lps.size(); i++){
    if(lps[i] > mx) mx = lps[i];
  }
  double sum = 0;
  for(size_t i=0; i<lps.size(); i++){
    sum += exp(lps[i] - mx);
  }
  return mx + log(sum);
}

//計算用のdp[i][j][A]形式の3次元dpテーブル
template<class T>
class DPTable {
  int N;
  std::vector< std::vector< std::map<std::string,T> > > dp;
public:
  DPTable(int N):dp(N, std::vector< std::map<std::string,T> >(N)){}

  bool exists(size_t i, size_t j, const std::string& S){ return dp[i][j].count(S) > 0; }

  T get(size_t i, size_t j, const std::string& S){ return dp[i][j][S]; }

  void set(size_t i, size_t j, const std::string& S, const T& val){ dp[i][j][S] = val; }

  std::vector<std::string> get_str_list(size_t i, size_t j){
    std::vector<std::string> ret;
    for(typename std::map<std::string,T>::iterator itr = dp[i][j].begin();
        itr != dp[i][j].end();
        ++itr){
      ret.push_back(itr->first);
    }
    return ret;
  }
};

//文法< A -> B C, p >
struct Grammar {
  std::string lhs; // A
  std::pair<std::string,std::string> rhs; // B C (lexicalならCは空)
  double lp; //確率値の対数
  Grammar(const std::string lhs,
          const std::pair<std::string,std::string> rhs,
          double lp):
    lhs(lhs),rhs(rhs),lp(lp){}
};
bool operator==(const Grammar& a, const Grammar& b){
  return a.lhs == b.lhs && a.rhs == b.rhs;
}
bool operator<(const Grammar& a, const Grammar& b){
  if(a.lhs == b.lhs) return a.rhs < b.rhs;
  return a.lhs < b.lhs;
}

//文法の管理
class Grammars {
  std::string start;
  std::map< std::pair<std::string,std::string>, std::vector<Grammar> > grammars;
public:
  Grammars(std::string start = "S"):start(start){}

  void set_start(const std::string& st){
    start = st;
  }

  std::string get_start(){
    return start;
  }

  void add(const Grammar& grm){
    grammars[grm.rhs].push_back(grm);
  }

  std::vector<Grammar> search(const std::pair<std::string,std::string>& rhs){
    return grammars[rhs];
  }

  //確率値を適当な値で埋める
  void fill_random(){
    std::map< std::string,std::vector<double> > sum;
    for(std::map< std::pair<std::string,std::string>, std::vector<Grammar> >::iterator itr = grammars.begin();
        itr != grammars.end();
        ++itr){
      for(size_t i=0; i<(itr->second).size(); i++){
        (itr->second)[i].lp = log(frand() * 0.09 + 0.01);//0.01〜0.1で適当に与える(次の正規化でΣp(A->*)=1となるように調整する)
        sum[(itr->second)[i].lhs].push_back((itr->second)[i].lp);
      }
    }
    //正規化
    std::map<std::string,double> norm;
    for(std::map< std::string,std::vector<double> >::iterator itr = sum.begin();
        itr != sum.end();
        ++itr){
      norm[itr->first] = logsumexp(itr->second);
    }
    for(std::map< std::pair<std::string,std::string>, std::vector<Grammar> >::iterator itr = grammars.begin();
        itr != grammars.end();
        ++itr){
      for(size_t i=0; i<(itr->second).size(); i++){
        (itr->second)[i].lp -= norm[(itr->second)[i].lhs];
      }
    }
  }

  std::vector<Grammar> get_all_grammar(){
    std::vector<Grammar> ret;
    for(std::map< std::pair<std::string,std::string>, std::vector<Grammar> >::iterator itr = grammars.begin();
        itr != grammars.end();
        ++itr){
      for(size_t i=0; i<(itr->second).size(); i++){
        ret.push_back((itr->second)[i]);
      }
    }
    return ret;
  }  

  void dump(){
    for(std::map< std::pair<std::string,std::string>, std::vector<Grammar> >::iterator itr = grammars.begin();
        itr != grammars.end();
        ++itr){
      for(size_t i=0; i<(itr->second).size(); i++){
        Grammar grm = (itr->second)[i];
        std::cout << grm.lhs << " -> " << grm.rhs.first << " " << grm.rhs.second << "\t" << exp(grm.lp) << std::endl;
      }
    }
  }
};



class PCFG {
public:
  //最尤導出計算用
  struct Info {
    std::string lhs;
    std::pair<std::string,std::string> rhs;
    std::pair< std::pair<int,int>, std::pair<int,int> > idx;
    double delta;
    std::pair<int,int> ans; //結果パース用
    Info():delta(-INF){}
  };
  //Inside-Outside algorithm計算用
  struct BetaAlpha {
    double beta;
    double alpha;
    BetaAlpha():beta(-INF),alpha(-INF){}
    BetaAlpha(double beta, double alpha):beta(beta),alpha(alpha){}
  };
private:
  static const double INF;

  //文法の管理
  Grammars grammars;

  /* 最尤導出計算用 */
  //dpテーブルから最尤な導出を配列形式にして返す
  std::pair< double,std::vector<Info> > get_best_tree(size_t N, DPTable<Info>& dp){
    std::vector<Info> ret;

    if(N==0 || !dp.exists(0,N-1,grammars.get_start())) return std::make_pair(-1,ret);

    //開始記号から構文木をたどる
    size_t i = 0;
    Info tmp = dp.get(0,N-1,grammars.get_start());
    double p = tmp.delta;
    ret.push_back(tmp);
    ret[i].ans = std::make_pair(ret.size(), ret.size()+1);
    ret.push_back(dp.get(tmp.idx.first.first,tmp.idx.first.second,tmp.rhs.first));
    ret.push_back(dp.get(tmp.idx.second.first,tmp.idx.second.second,tmp.rhs.second));
    i++;
    while(i < ret.size()){
      tmp = ret[i];
      if(tmp.rhs.second != ""){ //A -> B C
        ret[i].ans = std::make_pair(ret.size(), ret.size()+1);
        ret.push_back(dp.get(tmp.idx.first.first,tmp.idx.first.second,tmp.rhs.first)); //B
        ret.push_back(dp.get(tmp.idx.second.first,tmp.idx.second.second,tmp.rhs.second));//C
      } else { //A -> a
        ret[i].ans = std::make_pair(-1, -1); //a
      }
      i++;
    }
    return std::make_pair(p, ret);
  }

  /* Inside-Outsideアルゴリズム用 */
  //DPテーブルを作成
  void make_dp_table(DPTable<BetaAlpha>& dp, const std::vector<std::string>& text){
    std::vector<Grammar> grms = grammars.get_all_grammar();
    size_t N = text.size();
    //内側確率betaの計算
    {//対角要素を更新
      for(size_t i=0; i<N; i++){
        std::vector<Grammar> v = grammars.search(std::make_pair(text[i],""));
        for(size_t j=0; j<v.size(); j++){
          dp.set(i, i, v[j].lhs, BetaAlpha(v[j].lp, -INF));
        }
      }
    }
    {//残りの部分を更新
      for(size_t n=1; n<N; n++){
        for(size_t i=0; i<N-n; i++){
          std::map< std::string, std::vector<double> > memo;
          for(size_t k=0; k<grms.size(); k++){
            std::vector<double> v;
            for(size_t j=1; j<=n; j++){
              bool ok = true; //A->B CでBとCが両方0以上の値が存在しているものだけ計算
              double lp = 0;
              if(dp.exists(i,i+j-1,grms[k].rhs.first)){
                lp += dp.get(i,i+j-1,grms[k].rhs.first).beta;
              }else{
                ok = false;
              }
              if(dp.exists(i+j,i+n,grms[k].rhs.second)){
                lp += dp.get(i+j,i+n,grms[k].rhs.second).beta;
              }else{
                ok = false;
              }
              if(ok) v.push_back(lp);
            }
            if(v.size() > 0) memo[grms[k].lhs].push_back(grms[k].lp + logsumexp(v));
          }
          for(std::map< std::string, std::vector<double> >::iterator itr = memo.begin();
              itr != memo.end();
              ++itr){
            if((itr->second).size() > 0) dp.set(i,i+n,itr->first,BetaAlpha(logsumexp(itr->second), -INF));
          }
        }
      }
    }

    if(!dp.exists(0,N-1,grammars.get_start())) return; //構築に失敗したら終了
    
    //外側確率alphaの計算
    {//一番右上の外側確率を設定
      BetaAlpha ba = dp.get(0,N-1,grammars.get_start());
      dp.set(0,N-1,grammars.get_start(),BetaAlpha(ba.beta,0));
    }
    {//A -> B Cの形の生成規則に対し、Aの外側確率からBとCの外側確率を更新
      for(size_t n=N-1; n>=1; n--){
        for(size_t i=0; i<N-n; i++){
          for(size_t j=0; j<grms.size(); j++){
            if(!dp.exists(i,i+n,grms[j].lhs)) continue;
          
            std::map< std::pair< std::pair<size_t,size_t>,std::string >, std::vector<double> > memo;
            for(size_t k=1; k<=n; k++){
              //alpha[i+k][i+n][C] += P(A->BC) * alpha[i][i+n][A] * beta[i][i-1+k][B]
              if(dp.exists(i+k,i+n,grms[j].rhs.second) && dp.exists(i,i-1+k,grms[j].rhs.first)){
                double lp = grms[j].lp + dp.get(i,i+n,grms[j].lhs).alpha + dp.get(i,i-1+k,grms[j].rhs.first).beta;
                memo[std::make_pair(std::make_pair(i+k,i+n),grms[j].rhs.second)].push_back(lp);
              }
              //alpha[i][i+n-k][B] += P(A->BC) * alpha[i][i+n][A] * beta[i+n+1-k][i+n][C]
              if(dp.exists(i,i+n-k,grms[j].rhs.first) && dp.exists(i+n+1-k,i+n,grms[j].rhs.second)){
                double lp = grms[j].lp + dp.get(i,i+n,grms[j].lhs).alpha + dp.get(i+n+1-k,i+n,grms[j].rhs.second).beta;
                memo[std::make_pair(std::make_pair(i,i+n-k),grms[j].rhs.first)].push_back(lp);
              }
            }
            for(std::map< std::pair< std::pair<size_t,size_t>,std::string >, std::vector<double> >::iterator itr = memo.begin();
                itr != memo.end();
                ++itr){
              size_t l = (itr->first).first.first;
              size_t r = (itr->first).first.second;
              std::string A = (itr->first).second;
              if(!dp.exists(l,r,A)) continue;
              std::vector<double> v = itr->second;
              BetaAlpha ba = dp.get(l,r,A);
              v.push_back(ba.alpha);
              ba.alpha = logsumexp(v);
              dp.set(l,r,A,ba);
            }
          }
        }
      }
    }
  }
  
  //テキストの確率値を取得
  double get_sentence_prob(DPTable<BetaAlpha>& dp, const std::vector<std::string>& text){
    size_t N = text.size();
    return dp.get(0,N-1,grammars.get_start()).beta;
  }
  
  //コーパスの尤度を計算
  double calc_log_likelihood(const std::vector< std::vector<std::string> >& corpus){
    double ret = 0.0;
    for(size_t i=0; i<corpus.size(); i++){
      DPTable<BetaAlpha> dp(corpus[i].size());
      make_dp_table(dp, corpus[i]);
      if(!dp.exists(0,corpus[i].size()-1,grammars.get_start())){ //構文木が作成できていない場合はスキップ
        //std::cerr << "text[" << i << "] skipped" << std::endl;
        continue;
      }
      double lp = get_sentence_prob(dp, corpus[i]);
      //std::cerr << "text[" << i << "] = " << lp << std::endl;
      ret += lp;
    }
    return ret;
  }
    
  //ルールの適用確率の推定
  void update_grammar_prob(const std::vector< std::vector<std::string> >& corpus){
    std::vector<Grammar> grms = grammars.get_all_grammar();
    std::map< Grammar,std::vector<double> > cnt;
    //適用された回数をカウント
    for(size_t t=0; t<corpus.size(); t++){
      size_t N = corpus[t].size();
      DPTable<BetaAlpha> dp(N);
      make_dp_table(dp, corpus[t]);
      if(!dp.exists(0,N-1,grammars.get_start())) continue; //構文木が作成できていない場合はスキップ
      double lp_sentence = get_sentence_prob(dp, corpus[t]);
      for(size_t g=0; g<grms.size(); g++){
        std::vector<double> v;
        if(grms[g].rhs.second == ""){ //A -> a
          for(size_t i=0; i<N; i++){
            if(dp.exists(i,i,grms[g].lhs)) v.push_back(grms[g].lp - lp_sentence + dp.get(i,i,grms[g].lhs).alpha);
          }
        }else{ //A -> B C
          for(size_t n=1; n<=N-1; n++){
            for(size_t i=0; i<N-n; i++){
              for(size_t j=1; j<=n; j++){
                if(!dp.exists(i,i+n,grms[g].lhs)) continue;
                if(!dp.exists(i,i+j-1,grms[g].rhs.first)) continue;
                if(!dp.exists(i+j,i+n,grms[g].rhs.second)) continue;
                
                double lp = grms[g].lp - lp_sentence + 
                  dp.get(i,i+n,grms[g].lhs).alpha + 
                  dp.get(i,i+j-1,grms[g].rhs.first).beta +
                  dp.get(i+j,i+n,grms[g].rhs.second).beta;
                v.push_back(lp);
              }
            }
          }
        }
        cnt[grms[g]].push_back(logsumexp(v));
      }
    }
    //適用回数から確率を計算し、grammarsを更新
    std::map< std::string,std::vector<double> > sum;
    for(std::map< Grammar,std::vector<double> >::iterator itr = cnt.begin();
        itr != cnt.end();
        ++itr){
      sum[(itr->first).lhs].push_back(logsumexp(itr->second));
    }
    Grammars new_grms;
    new_grms.set_start(grammars.get_start());
    for(size_t i=0; i<grms.size(); i++){
      double lp = logsumexp(cnt[grms[i]]) - logsumexp(sum[grms[i].lhs]);
      new_grms.add(Grammar(grms[i].lhs, grms[i].rhs, lp));
    }
    grammars = new_grms;
  }

public:
  //文法の開始記号を登録
  void set_start(const std::string& st){
    grammars.set_start(st);
  }

  //文法を登録
  void add_grammar(const std::string& lhs, const std::pair<std::string,std::string>& rhs, double lp){
    grammars.add(Grammar(lhs,rhs,lp));
  }

  //文に対して最尤な導出を計算
  std::pair< double, std::vector<Info> > calc_best_tree(const std::vector<std::string>& text){
    size_t N = text.size();
    DPTable<Info> dp(N);
    //対角要素を計算
    for(size_t i=0; i<N; i++){
      std::vector<Grammar> v = grammars.search(std::make_pair(text[i],""));
      for(size_t j=0; j<v.size(); j++){
        Info info;
        info.lhs = v[j].lhs;
        info.rhs = v[j].rhs;
        info.idx.first = std::make_pair(i,i);
        info.idx.second = std::make_pair(-1,-1);
        info.delta = v[j].lp;
        
        if(dp.get(i, i, v[j].lhs).delta < v[j].lp){
          dp.set(i, i, v[j].lhs, info);
        }
      }
    }
    //三角行列の2番目以降の対角線上の要素を計算
    for(size_t n=1; n<N; n++){
      for(size_t i=0; i<N-n; i++){
        for(size_t j=1; j<=n; j++){
          std::vector<std::string> v1 = dp.get_str_list(i, i+j-1);
          std::vector<std::string> v2 = dp.get_str_list(i+j, i+n);
          for(size_t v1i=0; v1i<v1.size(); v1i++){
            for(size_t v2i=0; v2i<v2.size(); v2i++){
              std::vector<Grammar> v = grammars.search(std::make_pair(v1[v1i],v2[v2i]));
              for(size_t k=0; k<v.size(); k++){
                if(dp.get(i, i+n, v[k].lhs).delta <
                   v[k].lp + dp.get(i,i+j-1,v1[v1i]).delta + dp.get(i+j,i+n,v2[v2i]).delta){
                  Info info;
                  info.lhs = v[k].lhs;
                  info.rhs = v[k].rhs;
                  info.idx.first = std::make_pair(i,i+j-1);
                  info.idx.second = std::make_pair(i+j,i+n);
                  info.delta = v[k].lp + dp.get(i,i+j-1,v1[v1i]).delta + dp.get(i+j,i+n,v2[v2i]).delta;
                  
                  dp.set(i,i+n,v[k].lhs,info);
                }
              }
            }
          }
        }
      }
    }
    return get_best_tree(N, dp);
  }

  //Inside-Outside algorithmで文法の確率値を推定
  void estimate(const std::vector< std::vector<std::string> >& corpus, size_t ITER = 100000, double EPS = 1e-6){
    std::cerr << "Start Training..." << std::endl;
    //文法の確率を適当な値で埋める
    grammars.fill_random();
    
    double LL_prev = -INF, LL_now;
    for(size_t iter=1; iter<=ITER; iter++){
      //確率値の更新
      update_grammar_prob(corpus);
      //尤度の計算
      LL_now = calc_log_likelihood(corpus);

      std::cerr << "ITER = " << iter << "\tLL : " << LL_now << std::endl;
      if(fabs(LL_now - LL_prev) < EPS) break;
      LL_prev = LL_now;
    }

    //推定後の文法情報をダンプ
    std::cout << "[New Grammar]" << std::endl;
    grammars.dump();
    std::cout << std::endl;
  }

  //文を解析し、結果を出力、文の確率を返す
  double dump(const std::vector<std::string>& text){
    size_t N = text.size();
    std::pair< double, std::vector<Info> > best_tree = calc_best_tree(text);
    std::vector<Info> ret = best_tree.second;

    DPTable<BetaAlpha> dp(N);
    make_dp_table(dp, text);
    double lp_sentence = 0.0;
    if(dp.exists(0,N-1,grammars.get_start())){
      lp_sentence = get_sentence_prob(dp, text);
    }
    
    std::cout << "Text : ";
    for(size_t i=0; i<text.size(); i++){
      if(i!=0) std::cout << " ";
      std::cout << text[i];
    }
    std::cout << std::endl;
    std::cout << "P(W,T_best) : " << exp(best_tree.first) << std::endl;
    std::cout << "P(W) : " << exp(lp_sentence) << std::endl;
    
    for(size_t i=0; i<best_tree.second.size(); i++){
      std::cout << i << "\t";
      std::cout << ret[i].ans.first << "\t" << ret[i].ans.second << "\t";
      std::cout << ret[i].lhs << " -> " << ret[i].rhs.first << " " << ret[i].rhs.second << std::endl;
    }
    std::cout << std::endl;
    return best_tree.first;
  }
  
};
const double PCFG::INF = 1e100;

int main(){
  PCFG pcfg;

  //文法の登録(書籍のp.128図5.1)
  pcfg.set_start("S"); //開始記号
  pcfg.add_grammar("S", std::make_pair("N","V"), log(0.4));
  pcfg.add_grammar("S", std::make_pair("S","PP"), log(0.5));
  pcfg.add_grammar("S", std::make_pair("V","N"), log(0.1));
  pcfg.add_grammar("V", std::make_pair("V","N"), log(0.4));
  pcfg.add_grammar("PP", std::make_pair("P","N"), log(1.0));
  pcfg.add_grammar("N", std::make_pair("N","PP"), log(0.1));
  pcfg.add_grammar("N", std::make_pair("I",""), log(0.4));
  pcfg.add_grammar("N", std::make_pair("Maria",""), log(0.3));
  pcfg.add_grammar("N", std::make_pair("pizza",""), log(0.2));
  pcfg.add_grammar("V", std::make_pair("eat",""), log(0.6));
  pcfg.add_grammar("P", std::make_pair("with",""), log(1.0));

  //コーパス
  std::vector< std::vector<std::string> > corpus;
  { //I eat pizza with Maria
    std::vector<std::string> text;
    text.push_back("I");
    text.push_back("eat");
    text.push_back("pizza");
    text.push_back("with");
    text.push_back("Maria");
    corpus.push_back(text);
  }
  { //Maria eats pizza
    std::vector<std::string> text;
    text.push_back("Maria");
    text.push_back("eat");
    text.push_back("pizza");
    corpus.push_back(text);
  }

  //コーパス全体の対数尤度の合計
  double ll_sum;
  
  //指定した確率値での解析結果
  ll_sum = 0;
  for(size_t i=0; i<corpus.size(); i++){
    ll_sum += pcfg.dump(corpus[i]);
  }
  std::cout << "LLsum = " << ll_sum << std::endl;
  std::cout << std::endl;
  
  //適用確率を推定
  pcfg.estimate(corpus);

  //推定後の確率値での解析結果
  ll_sum = 0;
  for(size_t i=0; i<corpus.size(); i++){
    ll_sum += pcfg.dump(corpus[i]);
  }
  std::cout << "LLsum = " << ll_sum << std::endl;
  std::cout << std::endl;
  
  return 0;
}

結果

実行結果。
導出木の配列の部分は、以下の形式。

木のノード番号(0は根ノード)   左の子ノードの番号   右の子ノードの番号  適用規則

Text : I eat pizza with Maria
P(W,T_best) : 0.001152     //書籍のp.128 P(W,T_2)と同じ
P(W) : 0.0013824           //書籍のp.129 P(W)と同じ
0       1       2       S -> S PP
1       3       4       S -> N V
2       5       6       PP -> P N
3       -1      -1      N -> I
4       7       8       V -> V N
5       -1      -1      P -> with
6       -1      -1      N -> Maria
7       -1      -1      V -> eat
8       -1      -1      N -> pizza

Text : Maria eat pizza
P(W,T_best) : 0.00576
P(W) : 0.00576
0       1       2       S -> N V
1       -1      -1      N -> Maria
2       3       4       V -> V N
3       -1      -1      V -> eat
4       -1      -1      N -> pizza

LLsum = -11.9231

Start Training...
ITER = 1        LL : -11.9597
ITER = 2        LL : -11.7674
ITER = 3        LL : -11.7563
ITER = 4        LL : -11.7553
ITER = 5        LL : -11.7552
ITER = 6        LL : -11.7552
ITER = 7        LL : -11.7552
ITER = 8        LL : -11.7552
[New Grammar]
N -> I  0.532391
N -> Maria      0.355308
N -> N PP       1.2904e-08
S -> N V        0.666667
PP -> P N       1
S -> S PP       0.333333
S -> V N        0
V -> V N        0.5
V -> eat        0.5
N -> pizza      0.112301
P -> with       1

Text : I eat pizza with Maria
P(W,T_best) : 0.00118018
P(W) : 0.00118018
0       1       2       S -> S PP
1       3       4       S -> N V
2       5       6       PP -> P N
3       -1      -1      N -> I
4       7       8       V -> V N
5       -1      -1      P -> with
6       -1      -1      N -> Maria
7       -1      -1      V -> eat
8       -1      -1      N -> pizza

Text : Maria eat pizza
P(W,T_best) : 0.00665024
P(W) : 0.00665024
0       1       2       S -> N V
1       -1      -1      N -> Maria
2       3       4       V -> V N
3       -1      -1      V -> eat
4       -1      -1      N -> pizza

LLsum = -11.7552

ランダムな確率から(toy)コーパスを使って推定しなおしても構文木は同じ結果になっている。
1つ目の文の構文木は以下。
f:id:jetbead:20161108011708p:image

上記の推定では、コーパスに出現してない生成規則「N -> N PP」の確率値が小さくなって、推定規則を使ったコーパスの対数尤度が元の対数尤度より大きくすることができている模様(-11.9231→-11.7552)。
また、初期値に与える確率値を変えると結果も変わることも確認できる。(EMアルゴリズム的に)

参考

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万回回してもまだ収束しきっていないようにも見える。。。

参考文献

2016-01-19

ラティスのNbestを求める

はじめに

形態素解析とかの解析時に使うラティス(形態素候補をグラフにしたもの)のうち、1番ベストな解だけが欲しい場合が多いが、2番目以降の解も欲しい場合がある。
N番目までの解を効率よく求める方法があり、使いたいケースが出てきたのに書いてみる。

Nbestとは

  • 解析したときに、スコアが良い順に、第一候補(1best)の解だけでなく、第二候補、第三候補・・・の解が考えられるとき、第N候補までのことをNbestという
  • Nbest解

前向きDP後ろ向きA*アルゴリズム
  • http://ci.nii.ac.jp/naid/110002725063
  • ラティスにおいて、効率よくNbestを求める方法が提案されている
  • 最初に、1bestを求める方法と同じようにスタートノードからあるノードまでの最小コストhをすべてのノードについて求めておく
    • 1bestの時は、ゴールノードからスタートノードに向かって、最小コストが選ばれるノードを逆に辿る事で求められた
    • Nbestの時は、「A*アルゴリズム」の考え方を利用して求める
  • ゴールノードからスタートノードに向かって、今考えているノードから(最適とは限らない、これまで取ったパスによって変わる)ゴールまでのコストgを求めていく
  • 前側のコスト合計hと後ろ側のコスト合計gを使うと、スタートから今考えているノードを通って(今まで通ってきた最適とは限らない)ゴールまでのパスの合計値fが求められる
  • fを使って小さい順に考えるノードを増やしながら、スタートノードにたどり着くパスを見つけると、Nbestの解を順次求められる

書いてみたけどわかりにくい。
徳永さんの「日本語入力を支える技術」の中で丁寧に紹介されている。

コード

スタートノード(BOS)とゴールノード(EOS)が1つずつあり、必ずスタートノードからゴールノードにパスを持つようなDAGであれば大丈夫だと思うので、それで書いてみる。
(雑なので、パスのノードを別に持ったり、トポロジカル順序を求めてる無駄がある・・・)

#include <iostream>
#include <vector>
#include <queue>
#include <algorithm>
#include <map>
#include <set>

#define INF 99999999

struct Node {
  int id;
  int w;
  //id,w
  std::map<int,int> forward;
  std::map<int,int> backward;

  int h;
};

struct Path {
  int id;
  int g;
  int f;
  std::vector<int> rev; //現在idからゴールノードまでの逆順パス(注意:無駄にメモリを使っている)
};
bool operator>(const Path& a, const Path& b){
  return a.f > b.f;
}

//スタートとゴールとなるノードがそれぞれ1つずつあり、ゴールノード以外で終わるようなノードが無いDAG
class stDAG {
  int n; //ノード数
  int s, t; //スタートノード番号、ゴールノード番号
  std::vector<Node> nodes; //ノード情報

  //トポロジカル順序を求めるため
  bool topological_sort(std::vector<int>& order){
    std::vector<int> color(n);
    for(int i=0; i<n; i++){
      if(!color[i] && !topological_visit(i, order, color)) return false;
    }
    std::reverse(order.begin(), order.end());
    return true;
  }
  bool topological_visit(int v, std::vector<int>& order, std::vector<int>& color){
    color[v] = 1;
    for(std::map<int,int>::iterator itr = nodes[v].forward.begin();
        itr != nodes[v].forward.end();
        ++itr){
      if(color[itr->first] == 2) continue;
      if(color[itr->first] == 1) return false;
      if(!topological_visit(itr->first, order, color)) return false;
    }
    order.push_back(v);
    color[v] = 2;
    return true;
  }

public:
  //ノード数、スタートノードID、ゴールノードID
  stDAG(int n, int s, int t):
    n(n),
    nodes(n),
    s(s),
    t(t)
  {
    for(int i=0; i<n; i++){
      nodes[i].id = i;
      nodes[i].h = INF;
    }
  }

  //ノード情報の追加
  void add_node(int id, int w){
    nodes[id].w = w;
  }
  //エッジ情報の追加
  void add_edge(int id, int next_id, int w){
    nodes[id].forward[next_id] = w;
    nodes[next_id].backward[id] = w;
  }
  
  //前向きDP
  void forward_dp(){
    std::vector<int> order;
    if(!topological_sort(order)){
      std::cerr << "topological sort error" << std::endl;
      return;
    }
    if(nodes[order[0]].id != s || nodes[order[n-1]].id != t){
      std::cerr << "s,t id error" << std::endl;
      return;
    }

    nodes[order[0]].h = 0;
    for(int i=0; i<n; i++){
      int id = order[i];
      //std::cout << id << "\t" << nodes[id].h << std::endl;
      for(std::map<int,int>::iterator itr = nodes[id].forward.begin();
          itr != nodes[id].forward.end();
          ++itr){
        int next_id = itr->first;
        int next_w = nodes[next_id].w;
        int edge_w = itr->second;
        nodes[next_id].h = std::min(nodes[next_id].h, nodes[id].h + edge_w + next_w);
      }      
    }
  }

  //後向きA*
  void backward_astar(int N){
    int no = 1;
    std::priority_queue<Path, std::vector<Path>, std::greater<Path> > que;

    Path path;
    path.id = t;
    path.g = 0;
    path.f = nodes[t].h;
    path.rev.push_back(t);
    que.push(path);

    while(!que.empty()){
      path = que.top(); que.pop();
      int id = path.id;
      int bestf = path.f;

      if(id == s){
        std::cout << no << "best: ";
        for(int i=path.rev.size()-1; i>=0; i--){
          std::cout << path.rev[i] << " ";
        }
        std::cout << " => " << path.f + nodes[s].w << std::endl;
        if(no == N) return;

        no++;
      }

      for(std::map<int,int>::iterator itr = nodes[id].backward.begin();
          itr != nodes[id].backward.end();
          ++itr){
        int gg = nodes[id].w + itr->second + path.g;
        int ff = gg + nodes[itr->first].h;

        Path new_path;
        new_path.id = itr->first;
        new_path.g = gg;
        new_path.f = ff;
        new_path.rev = path.rev;
        new_path.rev.push_back(itr->first);
        que.push(new_path);
      }
    }
  }
};


int main(){
  //graph1
  stDAG dag(6, 0, 5);
  dag.add_node(0, 0);
  dag.add_node(1, 1);
  dag.add_node(2, 2);
  dag.add_node(3, 3);
  dag.add_node(4, 4);
  dag.add_node(5, 5);
  dag.add_edge(0, 1, 1);
  dag.add_edge(0, 2, 2);
  dag.add_edge(1, 3, 2);
  dag.add_edge(1, 4, 1);
  dag.add_edge(2, 3, 1);
  dag.add_edge(2, 4, 2);
  dag.add_edge(3, 5, 2);
  dag.add_edge(4, 5, 1);
  dag.forward_dp();
  dag.backward_astar(4);

  //graph2
  /*
  stDAG dag(7, 0, 6);
  dag.add_node(0, 3);
  dag.add_node(1, 1);
  dag.add_node(2, 2);
  dag.add_node(3, 1);
  dag.add_node(4, 2);
  dag.add_node(5, 3);
  dag.add_node(6, 2);
  dag.add_edge(0, 1, 2);
  dag.add_edge(0, 2, 1);
  dag.add_edge(1, 3, 3);
  dag.add_edge(2, 3, 1);
  dag.add_edge(3, 4, 1);
  dag.add_edge(3, 5, 2);
  dag.add_edge(4, 6, 4);
  dag.add_edge(5, 6, 1);
  dag.forward_dp();
  dag.backward_astar(4);
  */

  return 0;
}

結果


グラフ1

f:id:jetbead:20160119014826p:image

$ ./a.out
1best: 0 1 4 5  => 13
2best: 0 1 3 5  => 14
3best: 0 2 3 5  => 15
4best: 0 2 4 5  => 16

グラフ2

f:id:jetbead:20160119014827p:image

$ ./a.out
1best: 0 2 3 5 6  => 16
2best: 0 2 3 4 6  => 17
3best: 0 1 3 5 6  => 18
4best: 0 1 3 4 6  => 19

とりあえず合ってそうなので、大丈夫そう。

参考

2015-10-14

Minimal Acyclic Subsequential Transducerで遊ぶ

はじめに

https://pycon.jp/2015/ja/proposals/vote/11/
Pycon2015で発表された「Pythonで作って学ぶ形態素解析」で紹介されていた辞書データ構造の「Minimal Acyclic Subsequential Transducer」について、勉強のために書いてみた。

Minimal Acyclic Subsequential Transducerとは

  • Finite State Transducerの一種
  • Transducerにおいて、initial stateが一つで、同じ入力ラベルを共有する同じ状態からのの遷移が2つ以上なく、各最終状態での最終出力文字列が高々p個のとき、p-subsequentialで、pが整数ならfinitely subsequentialというらしい
  • minimal(状態数が最少)、Acyclic(サイクルが無い)
    • Cyril Allauzen, Mehryar Mohri, Finitely Subsequential Transducersとp-Subsequentiable Transducersあたりを読む
  • 「Lucene/Solrで使われているFST」はこれの事
    • Kagome(Golang), Janome(python)でも採用
    • 辞書サイズがコンパクトになるので内包してもそこまで大きくならない


わかりやすい他の方の参考資料

構築方法

OpenFstとか使うイメージでいたけど、上記で紹介されている通り、一時的な状態を作らずに、ソート済みの入力から直接FSTを構築する方法が提案されている。
詳しい手順も上記の資料(qiita)で紹介されいているので、省略。

コード

論文にできるだけ従って、入力・出力は文字列、出力は一つのみ、で実装。
問題がある部分は、状態の探索(member関数)で線形に等価な状態を探索している点や、状態の等価性判定(equal関数)が下記の方法でもよいのか怪しい点(Def.2-3?)、など。
他の実装は全然見てないので、間違ってたら後で修正する。

#include <iostream>
#include <vector>
#include <algorithm>
#include <list>
#include <cstdio>
#include <map>
#include <queue>

//入力文字列の最大長
#define MAX_WORD_SIZE 100

struct State {
  State* next[0x100];
  std::string output[0x100];
  std::string state_out;
  bool final;
  int _num;

  State(){
    clear();
  }
  State(State* s){
    _num = -1;
    state_out = s->state_out;
    final = s->final;
    for(size_t i=0; i<0x100; i++){
      next[i] = s->next[i];
      output[i] = s->output[i];
    }
  }
  ~State(){}
  void clear(){
    _num = -1;
    state_out = "";
    final = false;
    for(size_t i=0; i<0x100; i++){
      next[i] = (State*)0;
      output[i] = "";
    }
  }

  //状態の等価チェック
  // 与えられた状態sに対して、遷移状態と出力記号などが完全に一致するかを確認
  bool equal(State* s){
    if(final && s->final){
      if(state_out == s->state_out) return true;
      return false;
    }

    std::queue< std::pair< std::pair<State*,std::string>, std::pair<State*,std::string> > > que;
    for(size_t i=0; i<0x100; i++){
      que.push(std::make_pair(std::make_pair(next[i],""), std::make_pair(s->next[i],"")));
    }

    while(!que.empty()){
      std::pair<std::pair<State*,std::string>, std::pair<State*,std::string> > p = que.front(); que.pop();
      if(p.first.second != p.second.second) return false;
      if(p.first.first == NULL && p.second.first == NULL) continue;
      if(p.first.first == NULL || p.second.first == NULL) return false;
      if(p.first.first->state_out != p.second.first->state_out) return false;

      for(size_t i=0; i<0x100; i++){
        if(p.first.first->next[i] == NULL && p.second.first->next[i] == NULL) continue;
        que.push(std::make_pair(std::make_pair(p.first.first->next[i],p.first.second + (char)i), 
                                std::make_pair(p.second.first->next[i],p.second.second + (char)i)));
      }
    }
    return true;
  }
};

struct Dictionary {
  Dictionary(){}
  ~Dictionary(){
    for(std::list<State*>::iterator itr = states.begin();
        itr != states.end();
        ++itr){
      delete *itr;
    }
    states.clear();
    for(size_t i=0; i<MAX_WORD_SIZE; i++){
      if(TempState[i]){
        delete TempState[i];
        TempState[i] = NULL;
      }
    }
  }

  //辞書の状態に存在する状態かどうかチェック
  // [注意] 線形に全部等価チェックしているので非常に重い
  State* member(State* s){
    for(std::list<State*>::iterator itr = states.begin();
        itr != states.end();
        ++itr){
      State* p = *itr;
      if(s->equal(p)) return p;
    }
    return NULL;
  }

  void insert(State* s){
    states.push_back(s);
  }

  State* copy_state(State* s){
    return new State(s);
  }

  State* find_minimized(State* s){
    State* r = member(s);
    if(r == NULL){
      r = copy_state(s);
      insert(r);
    }
    return r;
  }
  void set_transition(State* s, char c, State* t){
    s->next[c] = t;
  }

  //辞書登録。inは辞書順ソートされている必要がある
  void create(const std::vector<std::string>& in, const std::vector<std::string>& out){
    for(size_t i=0; i<MAX_WORD_SIZE; i++) TempState[i] = new State;
    PreviousWord = "";
    TempState[0]->clear();

    std::string CurrentWord = "";
    std::string CurrentOutput = "";
    for(size_t t=0; t<in.size(); t++){
      CurrentWord = in[t];
      CurrentOutput = out[t];
      int i, j;
      i = 1;
      while(i<=CurrentWord.length() && i<=PreviousWord.length() && (CurrentWord[i-1] == PreviousWord[i-1])) i++;
      int PrefixLengthPlus1 = i;


      for(i=PreviousWord.size(); i>=PrefixLengthPlus1; i--){
        set_transition(TempState[i-1], PreviousWord[i-1], find_minimized(TempState[i]));
      }

      for(i=PrefixLengthPlus1; i<=CurrentWord.length(); i++){
        TempState[i]->clear();
        set_transition(TempState[i-1], CurrentWord[i-1], TempState[i]);
      }

      if(in[t] != PreviousWord){
        TempState[CurrentWord.length()]->final = true;
      }


      for(j=1; j<=PrefixLengthPlus1-1; j++){
        if(TempState[j-1] == NULL) continue;
        std::string outputStr = TempState[j-1]->output[CurrentWord[j-1]];
        std::string CommonPrefix = "";
        for(int k=0; k<std::min(outputStr.length(), CurrentOutput.length()); k++){
          if(outputStr[k] != CurrentOutput[k]) break;
          CommonPrefix += outputStr[k];
        }
        std::string WordSuffix = outputStr.substr(CommonPrefix.length());

        TempState[j-1]->output[CurrentWord[j-1]] = CommonPrefix;
        for(size_t c=0; c<0x100; c++){
          if(TempState[j]->next[c] != NULL){
            TempState[j]->output[c] = WordSuffix + TempState[j]->output[c];
          }
        }
        CurrentOutput = CurrentOutput.substr(CommonPrefix.length());
      }
      if(in[t] == PreviousWord){
        TempState[CurrentWord.length()]->state_out += CurrentOutput;
      }else{
        TempState[PrefixLengthPlus1-1]->output[CurrentWord[PrefixLengthPlus1-1]] = CurrentOutput;
      }
      PreviousWord = CurrentWord;
    }
    for(int i=CurrentWord.length(); i>=1; i--){
      TempState[i-1]->next[PreviousWord[i-1]] = find_minimized(TempState[i]);
      InitialState = find_minimized(TempState[0]);
    }
  }

  //検索
  std::string search(const std::string& query){
    std::string ret = "";
    State* p = InitialState;
    ret += p->state_out;
    for(int i=0; i<query.length(); i++){
      if(!p->next[query[i]]) return "";
      ret += p->output[query[i]];
      p = p->next[query[i]];
      ret += p->state_out;
    }
    if(p->final) return ret;
    return "";
  }


  //内部状態を出力
  void dump(){
    std::queue<State*> que;
    que.push(InitialState);
    int state_id = 0;
    while(!que.empty()){
      State* p = que.front(); que.pop();
      if(p->_num >= 0) continue;
      p->_num = state_id;
      for(int i=0; i<0x100; i++){
        if(p->next[i] == NULL) continue;
        que.push(p->next[i]);
      }
      state_id++;
    }

    std::cout << "size : " << states.size() << std::endl;
    std::map<int,int> memo;
    que.push(InitialState);
    while(!que.empty()){
      State* p = que.front(); que.pop();
      if(memo.count(p->_num) > 0) continue;
      memo[p->_num] = 1;

      std::cout << "node " << p->_num << " : [" << p->state_out << "] " << (p->final?"final":"") << std::endl;
      for(int i=0; i<0x100; i++){
        if(p->next[i] == NULL) continue;
        std::cout << "  " << (char)i << "->" << p->next[i]->_num << " [" << p->output[i] << "]" << std::endl;
        que.push(p->next[i]);
      }
    }
  }

private:
  std::list<State*> states;
  State* TempState[MAX_WORD_SIZE];
  std::string PreviousWord;
  State* InitialState;
};


int main(){
  Dictionary dic;

  std::vector<std::string> in, out;
  in.push_back("apr"); out.push_back("30");
  in.push_back("aug"); out.push_back("31");
  in.push_back("dec"); out.push_back("31");
  in.push_back("feb"); out.push_back("28");
  in.push_back("jan"); out.push_back("31");
  in.push_back("jul"); out.push_back("31");

  dic.create(in, out);

  dic.dump();

  std::cout << "-------" << std::endl;

  std::cout << "apr : " << dic.search("apr") << std::endl;
  std::cout << "aug : " << dic.search("aug") << std::endl;
  std::cout << "dec : " << dic.search("dec") << std::endl;
  std::cout << "feb : " << dic.search("feb") << std::endl;
  std::cout << "jan : " << dic.search("jan") << std::endl;
  std::cout << "jul : " << dic.search("jul") << std::endl;
  std::cout << "abc : " << dic.search("abc") << std::endl;

  return 0;
}

結果

size : 12
node 0 : [] 
  a->1 [3]
  d->2 [31]
  f->3 [28]
  j->4 [31]
node 1 : [] 
  p->5 [0]
  u->6 [1]
node 2 : [] 
  e->7 []
node 3 : [] 
  e->8 []
node 4 : [] 
  a->9 []
  u->10 []
node 5 : [] 
  r->11 []
node 6 : [] 
  g->11 []
node 7 : [] 
  c->11 []
node 8 : [] 
  b->11 []
node 9 : [] 
  n->11 []
node 10 : [] 
  l->11 []
node 11 : [] final
-------
apr : 30
aug : 31
dec : 31
feb : 28
jan : 31
jul : 31
abc : 

ランダムなアルファベット文字列でやってみても問題ないので一応大丈夫そう。
メモリ効率や高速化など実用レベルにするなら結構賢く書かないと大変そう。

参考