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

参考文献

2016-04-15

検索エンジンの日本語トークナイズメモ

はじめに

検索エンジンのトークナイズ処理の部分で行われている基本処理や工夫を少し調べてみたのでメモ。

トークナイズ処理

  • 「検索クエリ」に対してマッチする「ドキュメント」を高速に検索するためにインデクス(索引)を作成する
    • 本の最後の方にある「用語 - ページ」のような感じで、速く目的の用語が書いてあるページを調べられる
  • インデクスは、日本語の場合文字が連続しているため、「形態素」や「(文字)N-gram」などが使われる
文1「六本木ヒルズに行った」
文2「青山さんから電話があった」

【形態素でインデクスを作成する場合の例】
文1:「六本木ヒルズ」「に」「行く」「た」
文2:「青山」「さん」「から」「電話」「が」「あっ」「た」

【文字2-gram(bigram)でインデクスを作成する場合】
文1:「六本」「本木」「木ヒ」「ヒル」「ルズ」「ズに」「に行」「行っ」「った」
文2:「青山」「山さ」「さん」「んか」「から」「ら電」「電話」「話が」「があ」「あっ」「った」
  • 「転置インデクス」はこれを逆にしたもので、以下のように保持することで、そのインデクスを持つ文を高速に探せる
六本木ヒルズ:文1
青山:文2
た:文1、文2
...

トークナイズ処理における問題

  • 「六本木ヒルズ」というインデクスを作ってしまうと、検索クエリで「ヒルズ」が来た場合、転置インデクス内には完全マッチするものはないので、検索できず検索漏れになる
    • なので、できるだけ短めの単位でインデクスが作られている方がよかったりするが、短くすると意味が違ってしまうものまで拾ってしまう可能性がでてくる
  • また、検索クエリ側でもインデクスと同じ処理をし単位をそろえる必要があるが、形態素解析などを使うと「クエリの解析結果」と「文章中での解析結果」がちょっと変わってしまったりすることがある
  • カタカナ、アルファベット、数字、記号などは、日本語の形態素解析ではきちんと解析されない場合、それらを含む検索に失敗する可能性もある
    • 字種交じりだと切れ目がおかしくなりやすい

基本処理・工夫

形態素N-gramでのハイブリットインデクシング

http://www.slideshare.net/techblogyahoo/17lucenesolr-solrjp-apache-lucene-solrnbest

  • インデクスを「形態素」と「N-gram」を一緒に使うことで、「形態素」インデクスで意味の合う検索を抑えつつ、「N-gram」インデクスで検索漏れを防げる
  • 検索ノイズ(意味の合わないドキュメントが検索される)が大量に発生してしまう可能性がある
  • インデクスサイズが大きくなる

辞書登録による解析単位の調整

同上

  • 短い単位に解析されるよう辞書に登録
  • 登録するのが大変(保守も大変)

N-Bestによるトークン拡張

同上

  • 複数の解釈が可能な文などは、どの解析が正しいかを判断するのは難しい
  • 解析結果の上位N個(正確には、一番スコアが高い解析結果からのコスト差が閾値以下)を利用する

kuromojiのSEARCH MODE

http://atilika.org/

  • KuromojiにはSEARCH MODEというできるだけ短い単位で形態素解析するモードがついている
  • 処理内容としては、長い形態素の場合、形態素コストにペナルティを加算し出にくくする、ということをしている模様
    • code
    • 追記:solrだと短い形態素だけでなく長い形態素もインデクシングされるのなんでかなと思ったら、プログラムが違うようで、以下が実際の処理だった
    • JapaneseTokenizer.java
  • EXTENDED MODEというのは、さらに未知語を1文字ずつに解析させる処理が加わるそう

品詞や文字の正規化・除外・ストップワード

http://www.rondhuit.com/solr%E3%81%AE%E6%97%A5%E6%9C%AC%E8%AA%9E%E5%AF%BE%E5%BF%9C.html

  • 区切り文字、記号、助詞、助動詞などを検索したい場合は少ないので、ノイズとなりうるそれらの品詞でフィルタ
    • ただし、記号とかは意味のある場合が多いので、注意が必要そう
  • 全角、半角はあまり区別する必要がないことが多いので、どちらかに統一しておく
  • 動詞形容詞などの活用語は基本形に直す
  • 組文字(箸覆)、異体字とかも(下記の書籍参照)
  • 中黒(「・」)の処理、など
  • よく出てきて意味のないものは「ストップワード」として除外する

トークンの部分一致検索

山田,末永,「検索エンジン自作入門」,技術評論社

表記ゆれ同義語の考慮

同上

  • 表記が「サーバー」「サーバ」や「バイオリン」「ヴァイオリン」など揺れる可能性がある
  • これらを同一視できるような辞書を用意して統一しておく


参考