Negative/Positive Thinking

2017-06-18

XBWを試す

はじめに

XBWをWaveletMatrixを使って、試しに実装してみた。

XBWとは

コード

XBWでの表のソートはそのまま文字列同士のソートをしている。
チェックは、コード中にあるように、適当にkey文字列とkey文字列じゃないのを生成してtrieと結果が一緒になるかだけ。
WaveletMatrixの方で理解のためにいくつか関数を書いているけど、rank()ぐらいしか使っていないので、それ以外はあまりVerifyできていない。
(g++はバージョン「5.4.0」、オプション「-std=gnu++1y -O2」で実行してる)

#include <vector>
#include <map>
#include <cstdint>
#include <algorithm>
#include <iostream>
#include <queue>
#include <random>

//完備辞書(Fully Indexable Dictionary)
// 【使いまわすときの注意】
// - 全部set()したら、最後にfinalize()を呼ぶこと
// - select(x)の実装で、xが要素数よりも多い場合-1を返す実装にしている
// - 32bit/64bit書き換えは、BIT_SIZE,BLOCK_TYPE,popcount,整数リテラルのサフィックスなどを書き換えること
class FID {
  static const int BIT_SIZE = 64;
  using BLOCK_TYPE = uint64_t;
  int size;
  int block_size;
  std::vector<BLOCK_TYPE> blocks;
  std::vector<int> s_rank;
public:
  //for BIT_SIZE == 32
  /*
  BLOCK_TYPE popcount(BLOCK_TYPE x){
    x = ((x & 0xaaaaaaaa) >> 1) + (x & 0x55555555);
    x = ((x & 0xcccccccc) >> 2) + (x & 0x33333333);
    x = ((x & 0xf0f0f0f0) >> 4) + (x & 0x0f0f0f0f);
    x = ((x & 0xff00ff00) >> 8) + (x & 0x00ff00ff);
    x = ((x & 0xffff0000) >> 16) + (x & 0x0000ffff);
    return x;
  }
   */
  //__builtin_popcount()
  
  //for BIT_SIZE == 64
  BLOCK_TYPE popcount(BLOCK_TYPE x){
    x = ((x & 0xaaaaaaaaaaaaaaaaULL) >> 1) + (x & 0x5555555555555555ULL);
    x = ((x & 0xccccccccccccccccULL) >> 2) + (x & 0x3333333333333333ULL);
    x = ((x & 0xf0f0f0f0f0f0f0f0ULL) >> 4) + (x & 0x0f0f0f0f0f0f0f0fULL);
    x = ((x & 0xff00ff00ff00ff00ULL) >> 8) + (x & 0x00ff00ff00ff00ffULL);
    x = ((x & 0xffff0000ffff0000ULL) >> 16) + (x & 0x0000ffff0000ffffULL);
    x = ((x & 0xffffffff00000000ULL) >> 32) + (x & 0x00000000ffffffffULL);
    return x;
  }
  //__builtin_popcountll()
public:
  FID(int size):
  size(size),
  block_size(((size + BIT_SIZE - 1) / BIT_SIZE) + 1),
  blocks(block_size, 0),
  s_rank(block_size, 0){}
  
  void set(int i){
    blocks[i/BIT_SIZE] |= 1ULL << (i%BIT_SIZE);
  }
  
  void finalize(){
    s_rank[0] = 0;
    for(int i=1; i<block_size; i++){
      s_rank[i] = s_rank[i-1] + popcount(blocks[i-1]);
    }
  }
  
  bool access(int i){
    return (blocks[i/BIT_SIZE] >> (i%BIT_SIZE)) & 1ULL;
  }

  //iより前のビットが立っている個数
  int rank(int i){
    BLOCK_TYPE mask = (1ULL << (i%BIT_SIZE)) - 1;
    return s_rank[i/BIT_SIZE] + popcount(mask & blocks[i/BIT_SIZE]);
  }

  //x番目にビットが立っている位置
  int select(int x){
    if(rank((block_size-1) * BIT_SIZE) <= x) return -1; //注意
    int lb = 0, ub = block_size-1;
    while(ub-lb>1){
      int m = (lb+ub)/2;
      if(s_rank[m]<=x) lb = m;
      else ub = m;
    }
    int lbb = lb*BIT_SIZE, ubb = (lb+1)*BIT_SIZE;
    while(ubb-lbb>1){
      int m = (lbb+ubb)/2;
      if(rank(m)<=x) lbb = m;
      else ubb = m;
    }
    return lbb;
  }
};

//ウェーブレット行列(Wavelet Matrix)
// 【使いまわすときの注意】
// - 全部set()したら、最後にfinalize()を呼ぶこと
// - 32bit/64bit書き換えは、BIT_SIZE,VAL_TYPEなどを書き換えること
class WaveletMatrix {
  static const int BIT_SIZE = 8;
  using VAL_TYPE = uint8_t;
  int size;
  std::vector<VAL_TYPE> v;
  std::vector<FID> matrix;
  std::vector<int> sep;

  struct mytuple {
    int b, s, e;
    mytuple(int b, int s, int e):b(b),s(s),e(e){}
    bool operator<(const mytuple& x) const {
      return e-s < x.e-x.s;
    }
  };
public:
  WaveletMatrix(int size):
  size(size),
  v(size, 0),
  matrix(BIT_SIZE, FID(size)),
  sep(BIT_SIZE, 0){}

  void set(int i, VAL_TYPE val){
    v[i] = val;
  }

  void finalize(){
    std::vector<VAL_TYPE> w(v.size(), 0);
    for(int b=BIT_SIZE-1; b>=0; b--){
      for(int i=0; i<size; i++){
        if((v[i] >> b) & 1ULL) matrix[b].set(i);
        else sep[b]++;
      }
      int b1=0, b2=sep[b];
      for(int i=0; i<size; i++){
        if((v[i] >> b) & 1ULL) w[b2++] = v[i];
        else w[b1++] = v[i];
      }
      for(int i=0; i<size; i++){
        v[i] = w[i];
      }
      matrix[b].finalize();
    }
  }

  //元の配列のi番目の要素
  VAL_TYPE access(int i){
    VAL_TYPE ret = 0;
    for(int b=BIT_SIZE-1; b>=0; b--){
      if(matrix[b].access(i)){
        i = sep[b] + matrix[b].rank(i);
        ret = (ret << 1) + 1ULL;
      }else{
        i = i - matrix[b].rank(i);
        ret = (ret << 1);
      }
    }
    return ret;
  }

  //[0,i)の範囲にxが何個存在するか
  int rank(int i, VAL_TYPE x){
    int lb = 0, ub = i;
    for(int b=BIT_SIZE-1; b>=0; b--){
      if((x >> b) & 1ULL){
        lb = matrix[b].rank(lb);
        ub = matrix[b].rank(ub);
        lb += sep[b];
        ub += sep[b];
      }else{
        lb = lb - matrix[b].rank(lb);
        ub = ub - matrix[b].rank(ub);                
      }
    }
    return ub - lb;
  }

  //i番目(0-index)のxが出現する位置
  int select(int i, VAL_TYPE x){
    int lb = 0, ub = size;
    while(ub-lb>1){
      int m = (lb+ub)/2;
      if(rank(m, x)<=i) lb = m;
      else ub = m;
    }
    return lb;
  }

  //[s,e)の範囲を切り出してソートしたときのn番目(0-index)の要素
  VAL_TYPE quantile(int s, int e, int n){
    for(int b=BIT_SIZE-1; b>=0; b--){
      int zn = (e - s) - (matrix[b].rank(e) - matrix[b].rank(s));
      if(zn <= n){
        s = matrix[b].rank(s);
        e = matrix[b].rank(e);
        s += sep[b];
        e += sep[b];
        n = n - zn;
      }else{
        s = s - matrix[b].rank(s);
        e = e - matrix[b].rank(e);                
      }      
    }
    return v[s];
  }

  //[s,e)の範囲で出現回数が多い数値順に、その数値と出現回数のTop-K
  std::vector<std::pair<VAL_TYPE,int>> top_k(int s, int e, int k){
    std::vector<std::pair<VAL_TYPE,int>> ret;
    std::priority_queue<mytuple> que;
    que.push(mytuple(BIT_SIZE-1,s,e));
    while(!que.empty()){
      mytuple q = que.top(); que.pop();
      int b = q.b, st = q.s, en = q.e;
      if(b < 0){
        ret.push_back(std::make_pair(v[st], en-st));
        if((int)ret.size() >= k) break;
      }else{
        int os = matrix[b].rank(st) + sep[b];
        int oe = matrix[b].rank(en) + sep[b];
        int zs = st - matrix[b].rank(st);
        int ze = en - matrix[b].rank(en);
        if(ze-zs > 0) que.push(mytuple(b-1,zs,ze));
        if(oe-os > 0) que.push(mytuple(b-1,os,oe));
      }
    }
    return ret;
  }

  //[s,e)の範囲でx<=c<yを満たすような数値cの合計出現数
  int rangefreq(int s, int e, VAL_TYPE x, VAL_TYPE y){
    int ret = 0;
    std::queue<std::pair<mytuple,VAL_TYPE>> que;
    que.push(std::make_pair(mytuple(BIT_SIZE-1,s,e),0));
    while(!que.empty()){
      std::pair<mytuple,VAL_TYPE> q = que.front(); que.pop();
      int b = q.first.b, st = q.first.s, en = q.first.e;
      VAL_TYPE mn = q.second;
      VAL_TYPE mx = q.second | ((b>=0)?0:((-1ULL) >> (BIT_SIZE - 1 - b)));
      if(x <= mn && mx < y){
        ret += en-st;
      }
      else if(mx < x || y <= mn){
        continue;
      }
      else {
        if(b < 0) continue;
        int os = matrix[b].rank(st) + sep[b];
        int oe = matrix[b].rank(en) + sep[b];
        int zs = st - matrix[b].rank(st);
        int ze = en - matrix[b].rank(en);
        if(ze-zs > 0) que.push(std::make_pair(mytuple(b-1,zs,ze), q.second));
        if(oe-os > 0) que.push(std::make_pair(mytuple(b-1,os,oe), q.second | (1ULL << b)));
      }
    }
    return ret;
  }
};

//XBW
class XBW {
  using VAL_TYPE = uint8_t;
  const char LAST_CHAR = (char)(0xff);
  
  struct Trie {
    bool flg;
    std::string rpp;
    std::map<char,Trie> next;
    Trie(){ flg = false; }
    void insert(const std::string &str){
      Trie *r = this;
      for(size_t i=0; i<str.length(); i++){
        r = &(r->next[str[i]]);
      }
      r->flg = true;
    }
    bool find(const std::string &str){
      Trie *r = this;
      for(size_t i=0; i<str.length(); i++){
        if(r->next.count(str[i]) == 0) return false;
        r = &(r->next[str[i]]);
      }
      return r->flg;
    }
  };
  struct ST {
    std::string children;
    std::string rpp;
    ST(std::string children, std::string rpp):children(children),rpp(rpp){}
    bool operator<(const ST& x) const { return rpp < x.rpp; }
  };

  Trie root;
  int xbw_size;
  std::string xbw_str;

  WaveletMatrix wm;
  FID fid;
  std::map<char,int> C;

  void build(std::vector<ST>& v){
    //XBWのサイズと文字の出現数のカウント
    std::map<char,int> cnt;
    xbw_size = 0;
    for(size_t i=0; i<v.size(); i++){
      xbw_size += v[i].children.length();
      for(size_t j=0; j<v[i].children.length(); j++){
        cnt[v[i].children[j]]++;
      }
    }

    //構築
    wm = WaveletMatrix(xbw_size);
    fid = FID(xbw_size);
    int idx = 0;
    for(size_t i=0; i<v.size(); i++){
      fid.set(idx);
      for(size_t j=0; j<v[i].children.length(); j++){
        wm.set(idx, (VAL_TYPE)(v[i].children[j]));
        idx++;
      }
    }
    wm.finalize();
    fid.finalize();

    C[(char)(0)] = 1;
    for(int i=1; i<256; i++){
      C[(char)(i)] = C[(char)(i-1)] + cnt[(char)(i-1)];
    }

    //trieの削除
    //root.next.clear();
  }

  int rank(int i, VAL_TYPE x){
    int pos = fid.select(i);
    return wm.rank(((pos<0)?xbw_size:pos),x);
  }  
public:
  XBW():wm(1),fid(1){}

  void add(const std::string& key){
    root.insert(key);
  }

  void finalize(){
    //chilren, reverse prefix pathの表の作成
    std::vector<ST> v;
    std::queue<Trie*> que;
    que.push(&root);
    while(!que.empty()){
      Trie *r = que.front(); que.pop();
      std::string children;
      std::string rpp = r->rpp;
      for(std::map<char,Trie>::iterator it=(r->next).begin(); it!=(r->next).end(); ++it){
        children += it->first;
        (it->second).rpp = rpp + it->first;
        que.push(&(it->second));
      }
      if(r->flg){
        children += LAST_CHAR;
      }
      std::reverse(rpp.begin(), rpp.end());
      v.push_back(ST(children, rpp));
    }

    std::sort(v.begin(), v.end());

    //XBW文字列の作成
    for(size_t i=0; i<v.size(); i++){
      if(i>0) xbw_str += "|";
      for(size_t j=0; j<v[i].children.length(); j++){
        if(v[i].children[j] == LAST_CHAR) xbw_str += "__LAST__"; //表示の都合上
        else xbw_str += v[i].children[j];
      }
    }

    build(v);
  }

  std::string get_xbw_string(){
    return xbw_str;
  }
  
  bool trie_find(const std::string& key){
    return root.find(key);
  }
  
  bool find(const std::string& key){
    int r = 0;
    for(size_t i=0; i<key.length(); i++){
      if(rank(r+1,(VAL_TYPE)(key[i])) - rank(r,(VAL_TYPE)(key[i])) == 0) return false;
      r = C[key[i]] + rank(r,(VAL_TYPE)(key[i]));
    }
    if(rank(r+1,(VAL_TYPE)(LAST_CHAR)) - rank(r,(VAL_TYPE)(LAST_CHAR)) == 0) return false;
    return true;
  }
};

int main(){
  std::mt19937 rnd{ std::random_device()() };
  std::map<std::string,int> keys, no_keys;
  int max_size = 500000; //keyとno_keysの要素数
  int max_len = 40; //文字列の長さの最大値
  
  int turn = 0;
  while(keys.size() < max_size || no_keys.size() < max_size){
    //generate key string

    int len = rnd() % max_len + 1;
    std::string key = "";
    for(int j=0; j<len; j++){
      key += (char)(' ' + rnd()%95);
    }

    if(keys.count(key) != 0 || no_keys.count(key) != 0) continue;
    
    if(turn == 0 && keys.size() < max_size){
      keys[key] = 1;
      turn = 1 - turn;
    }
    else if(turn == 1 && no_keys.size() < max_size){
      no_keys[key] = 1;
      turn = 1 - turn;
    }
  }
  std::cout << "key generated..." << std::endl;
  
  XBW xbw;
  for(const auto& x : keys){
    xbw.add(x.first);
  }
  xbw.finalize();
  std::cout << "XBW built..." << std::endl;
  
  //std::cout << "XBW = " << xbw.get_xbw_string() << std::endl;
  
  bool error = false;
  for(const auto& x : keys){
    if(xbw.trie_find(x.first) != xbw.find(x.first)){
      std::cout << "error : " << x.first << std::endl;
      error = true;
    }
  }
  for(const auto& x : no_keys){
    if(xbw.trie_find(x.first) != xbw.find(x.first)){
      std::cout << "error : " << x.first << std::endl;
      error = true;
    }
  }
  if(!error) std::cout << "no error" << std::endl;
  
  return 0;
}


確認のために、解説ページのTrieからXBWを出力してみる。
main()の内容を以下のように変更すると確認できる。

int main(){
  std::vector<std::string> v{"to","tea","ten","i","in","inn","we"};
  XBW xbw;
  for(const auto& x : v){
    xbw.add(x);
  }
  xbw.finalize();
  std::cout << xbw.get_xbw_string() << std::endl;
  return 0;
}

結果。

itw|__LAST__|an|__LAST__|n__LAST__|__LAST__|n__LAST__|__LAST__|__LAST__|eo|e

参考

2017-02-23

Kneser-Ney smoothingで遊ぶ

はじめに

100-nlp-papersで紹介されてた一番最初の論文に、クナイザーネイスムージングのスッキリな実装が載っていたので書いてみる。

Joshua Goodman: A bit of progress in language modeling, MSR Technical Report, 2001.

Kneser-Ney smoothingとは

  • 言語モデルのスムージング(平滑化)手法一種で、高い言語モデル性能を実現している
    • ニューラル言語モデルでも比較によく使われる
  • アイデアとしては「(n-1)-gramが出現した文脈での異なり数」を使うこと
    • 頻度を使うと、高頻度なn-gramではその(n-1)-gramも多くなってしまうため、特定文脈でしかでないような(n-1)-gramに対しても高い確率値ことになっていて、歪んだ結果になってしまう
      • 「San Francisco」の頻度が多いと「Francisco」の頻度も高くなるが、P(Francisco|on)とかはあまり出現しないので低くなってほしいところ、「Franciscoの頻度」を使って確率値を推定すると高くなってしまう
    • 頻度ではなく、異なり数で(n-1)-gramの確率を推定することで、補正する
  • 上のレポートでは、Interpolatedな補間方法での実装例を紹介している
    • back-offな方法も考えらえる
    • discount(割引値)パラメータn-gramごとに分けた方法は「modified Kneser-Ney smoothing」と呼ばれている

UNKの扱い

コード

UNKのためにちょっと修正した。あんまりちゃんとチェックできていないけど、それっぽい数値を返しているのでおそらく大丈夫。

#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <unordered_map>

class InterpolatedKneserNeyTrigram {
  const std::string delim = "\t";
  double discount; //(0,1)の範囲で指定する
  std::unordered_map<std::string,int> TD, TN, TZ, BD, BN, BZ, UN;
  int UD;
public:
  InterpolatedKneserNeyTrigram():discount(0.1),UD(0){}
  InterpolatedKneserNeyTrigram(double d):discount(d),UD(0){}

  //ファイルに書き出し
  void save(const std::string& filename){
    std::ofstream fout(filename);
    if(!fout){ std::cerr << "cannot open file" << std::endl; return; }
    fout << discount << std::endl;
    fout << TD.size() << std::endl;
    std::unordered_map<std::string,int>::iterator it;
    for(it=TD.begin(); it!=TD.end(); ++it) fout << it->first << std::endl << it->second << std::endl;
    fout << TN.size() << std::endl;
    for(it=TN.begin(); it!=TN.end(); ++it) fout << it->first << std::endl << it->second << std::endl;
    fout << TZ.size() << std::endl;
    for(it=TZ.begin(); it!=TZ.end(); ++it) fout << it->first << std::endl << it->second << std::endl;
    fout << BD.size() << std::endl;
    for(it=BD.begin(); it!=BD.end(); ++it) fout << it->first << std::endl << it->second << std::endl;
    fout << BN.size() << std::endl;
    for(it=BN.begin(); it!=BN.end(); ++it) fout << it->first << std::endl << it->second << std::endl;
    fout << BZ.size() << std::endl;
    for(it=BZ.begin(); it!=BZ.end(); ++it) fout << it->first << std::endl << it->second << std::endl;
    fout << UN.size() << std::endl;
    for(it=UN.begin(); it!=UN.end(); ++it) fout << it->first << std::endl << it->second << std::endl;
    fout << UD << std::endl;    
    fout.close();
  }
  //ファイルから読み込み
  void load(const std::string& filename){
    std::ifstream fin(filename);
    if(!fin){ std::cerr << "cannot open file" << std::endl; return; }
    fin >> discount;
    int td, tn, tz, bd, bn, bz, un;
    std::string s;
    int c;
    fin >> td;
    for(int i=0; i<td; i++){ getline(fin,s); getline(fin, s); fin >> c; TD[s] = c; }
    fin >> tn;
    for(int i=0; i<tn; i++){ getline(fin,s); getline(fin, s); fin >> c; TN[s] = c; }
    fin >> tz;
    for(int i=0; i<tz; i++){ getline(fin,s); getline(fin, s); fin >> c; TZ[s] = c; }
    fin >> bd;
    for(int i=0; i<bd; i++){ getline(fin,s); getline(fin, s); fin >> c; BD[s] = c; }
    fin >> bn;
    for(int i=0; i<bn; i++){ getline(fin,s); getline(fin, s); fin >> c; BN[s] = c; }
    fin >> bz;
    for(int i=0; i<bz; i++){ getline(fin,s); getline(fin, s); fin >> c; BZ[s] = c; }
    fin >> un;
    for(int i=0; i<un; i++){ getline(fin,s); getline(fin, s); fin >> c; UN[s] = c; }
    fin >> UD;
    fin.close();
  }

  void set_discount(double d){ discount = d; }
  double get_discount() const { return discount; }
  
  void add_sentence(const std::vector<std::string>& sentence){
    std::string w2 = "", w1 = "";

    for(size_t i=0; i<sentence.size(); i++){
      std::string w0 = sentence[i];
      TD[ w2 + delim + w1 ]++;
      if(TN[ w2 + delim + w1 + delim + w0 ]++ == 0){
        TZ[ w2 + delim + w1 ]++;

        BD[ w1 ]++;
        if(BN[ w1 + delim + w0 ]++ == 0){
          BZ[ w1 ]++;

          UD++;
          UN[ w0 ]++;
        }
      }
      w2 = w1;
      w1 = w0;
    }
  }

  double prob(const std::vector<std::string>& sentence){
    std::string w2 = "", w1 = "";
    double ret = 0;

    for(size_t i=0; i<sentence.size(); i++){
      std::string w0 = sentence[i];
      double prob = 0;

      //そのままだとUN[w0]==0のときprob==0になるため、1/|V|を使うように変更
      double uniform = 1.0 / UN.size();
      double unigram = 0.0;
      if(UN.count( w0 ) > 0){
        unigram = (UN[ w0 ] - discount) / (double)UD;
      }
      unigram += discount * uniform;
      if(BD.count( w1 ) > 0){
        double bigram = 0;
        if(BN.count( w1 + delim + w0 ) > 0){
          bigram = (BN[ w1 + delim + w0 ] - discount) / BD[ w1 ];
        }
        bigram += BZ[ w1 ] * discount / BD[ w1 ] * unigram;

        if(TD.count( w2 + delim + w1 ) > 0){
          double trigram = 0;
          if(TN.count( w2 + delim + w1 + delim + w0 ) > 0){
            trigram = (TN[ w2 + delim + w1 + delim + w0 ] - discount) / TD[ w2 + delim + w1 ];
          }
          trigram += TZ[ w2 + delim + w1 ] * discount / TD[ w2 + delim + w1 ] * bigram;
          prob = trigram;
        }else{
          prob = bigram;
        }
      }else{
        prob = unigram;
      }
      ret += log(prob);      
      w2 = w1;
      w1 = w0;
    }
    return ret;
  }
};


int main(){
  InterpolatedKneserNeyTrigram lm;  
  std::vector< std::vector<std::string> > train_v, valid_v;
  
  {//ファイルの読み込み
    std::ifstream trainfs("train.txt");
    std::ifstream validfs("valid.txt");
    std::string w;
    std::vector<std::string> tmp;
    while(trainfs >> w){
      tmp.push_back(w);
      if(w == "EOS"){
        train_v.push_back(tmp);
        tmp.clear();
      }
    }
    
    tmp.clear();
    while(validfs >> w){
      tmp.push_back(w);
      if(w == "EOS"){
        valid_v.push_back(tmp);
        tmp.clear();
      }
    }
  }
  
  {//学習用の文を全部入れる
    for(size_t i=0; i<train_v.size(); i++){
      lm.add_sentence(train_v[i]);
    }
  }
  
  {//よさそうなdを探す
    double best = log(0), best_d = 0;
    double prec = 0.001;
    for(double d=prec; d<1; d+=prec){
      lm.set_discount(d);
      double logq = 0.0;
      for(size_t i=0; i<valid_v.size(); i++){
        logq += lm.prob(valid_v[i]);
      }
      std::cerr << d << "\t" << logq << std::endl;
      if(best < logq){
        best = logq;
        best_d = d;
      }
    }
    lm.set_discount(best_d);
    std::cerr << "best: " << best << " (d = " << best_d << ")" << std::endl;
  }

  lm.save("lm.data");

  return 0;
}

実験

データの準備

「坊ちゃん」の言語モデルを作ってみる。
青空文庫から「坊ちゃん」のテキストを取得し、「≪≫」などで囲まれた部分を削除したものを用意。
全部で470行で、10行ごとをdiscount係数確認用にする。

さらにそれを、mecab+ipadicで1行1単語にした以下のようなテキストを準備する。
(学習用train.txt 424行分、確認用valid.txt 46行分)

親譲り
の
無鉄砲
で
小
供
の
時
から
損
...
と
答え
た
。
EOS
親類
の
もの
...

(1行分の終わりには「EOS」を含む)

最適なパラメータの探索

確認用のデータで一番尤度が高くなるパラメータを採用する。

f:id:jetbead:20170222021821p:image
最適なのは、discount=0.897のとき、対数尤度が-29116ぐらい。


事例

いくつかの例で確率を見てみる。

s=「親譲り の 無鉄砲 だ EOS」 → log(P(s)) = -23.7238
s=「親譲り の ブレイブハート だ EOS」 → log(P(s)) = -33.8758
s=「吾輩 は 猫 だ EOS」 → log(P(s)) = -36.8097
s=「無鉄砲 な フレンズ だ EOS」 → log(P(s)) = -38.3098


参考

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-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

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

トークンの部分一致検索

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

  • (長めの)形態素の部分文字列を高速に検索できるようにしておくことで、検索ヒット数が少ない場合の検索再現率を改善
  • 形態素(トークン)に対して、半無限文字列(i文字目から最後の文字までの部分文字列)を前方一致で検索
    • suffix array的な感じ。パトレシアトライなどを使っている模様

表記ゆれ、同義語の考慮

同上

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


参考

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

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

参考