小規模データで単語の数を数えてみた (2)

研究

大規模データで単語の数を数える - ny23の日記 でみたように，単純に高頻度の item が欲しい場合には，小規模データで単語の数を数えてみた (1) - ny23の日記 で使った sketch-based なアルゴリズムよりは，counter-based なアルゴリズムの方が（キーを陽に保存するので）使い勝手が良い．というわけで，space saving アルゴリズムを実装してみた．カウンタのデータ構造には，原著論文

に書いてある stream summary をほぼ定義通り素直に実装して利用．

// GNU GPL version 2 copyright@ny23
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <tr1/unordered_set>

// easygoing implementation of stream-summary data structure
template <typename T>
class connector { // base class for an element in a linked list
public:
  T * prev;
  T * next;
  void pop () {   // pop self from the linked list
    prev->next = next;
    next->prev = prev;
  }
  void push (T * prev_) { // push self into the linked list (next to prev_)
    prev = prev_;
    next = prev_->next;
    prev->next = next->prev = static_cast <T*> (this);
  }
  bool is_singleton () const { return prev == this && next == this; }
private:
  connector (const connector&);
  connector& operator= (const connector&);
protected:
  connector (T * prev_ = 0, T * next_ = 0) : prev (prev_), next (next_) {}
  ~connector () {}
};

class bucket_t;
class item_t : public connector <item_t> {
public:
  char * key; // keep public for coming item
  item_t () :
    connector <item_t> (), key (0), _capacity (0), _delta (0), _parent (0) {}
  ~item_t () { delete [] key; }
  bucket_t * get_parent ()            const { return _parent; }
  void       set_parent (bucket_t * bucket) { _parent = bucket; }
  size_t     get_delta  ()            const { return _delta; }
  void reset_item (char* new_key, size_t new_capacity, size_t delta) {
    if (_capacity < new_capacity) { // avoid memory reallocation
      delete [] key;
      key = new char[new_capacity];
      _capacity = new_capacity;
    }
    std::strcpy (key, new_key);
    _delta = delta;
  }
private:
  size_t     _capacity; // trade space for fast update with larger epsilon
  size_t     _delta;
  bucket_t * _parent;
};

class bucket_t : public connector <bucket_t> {
public:
  bucket_t () :
    connector <bucket_t> (this, this), _count (0), _child (0) {};
  void pop_item  (item_t * item) {
    item->pop ();
    if (_child == item) _child = item->next;
  }
  void push_item (item_t * item) {
    item->push (_child);
    item->set_parent (this);
  }
  size_t   get_count ()        const { return _count; }
  void     set_count (size_t count)  { _count = count; }
  item_t * get_child ()        const { return _child; }
  void     set_child (item_t * item) {
    _child = item->next = item->prev = item; // singleton
    item->set_parent (this);
  }
private:
  size_t   _count;
  item_t * _child;
};

// eco-friendly recycling plant for bucket objects
class bucket_pool {
public:
  void recycle_bucket (bucket_t * bucket) { bucket->next = _M; _M = bucket; }
  bucket_t * create_bucket () {
    if (! _M) return new bucket_t; // You can't get blood out of a stone.
    bucket_t * bucket = _M;
    _M = bucket->next;
    return bucket;
  }
  bucket_pool () : _M (0) {}
  ~bucket_pool () {
    while (_M) { bucket_t * tmp = _M->next; std::swap (_M, tmp); delete tmp; }
  }
private:
  bucket_t * _M;
  bucket_pool (const bucket_pool&);
  bucket_pool& operator= (const bucket_pool&);
};

struct eq_ptr {
  bool operator () (const item_t * const a, const item_t * const b) const {
    for (size_t i = 0; a->key[i] == b->key[i]; ++i)
      if (! a->key[i]) return true;
    return false;
  }
};

struct hashfun_ptr { // FNV-1
  unsigned int operator () (const item_t * const a) const {
    unsigned int h = 0x811c9dc5;
    for (char const *p = a->key; *p; ++p)
      { h *= 0x01000193; h ^= *p; }
    return h;
  }
};

int main (int argc, char** argv) {

  typedef std::tr1::unordered_set <item_t *, hashfun_ptr, eq_ptr> hash_t;
  
  if (argc < 3)
    { std::fprintf (stderr, "Usage: %s file epsilon\n", argv[0]); std::exit (1); }
  
  FILE * fp = std::fopen (argv[1], "r");
  if (! fp)
    { std::fprintf (stderr, "no such file: %s\n", argv[1]); std::exit (1); }
  char buf[1 << 18]; std::setvbuf (fp, &buf[0], _IOFBF, 1 << 18);

  size_t k = static_cast <size_t> (std::ceil (1.0 / std::strtod (argv[2], NULL)));

  // monitoring streams
  hash_t hash;
  bucket_pool bucket_pool; // trade space for faster update
  bucket_t * root        = new bucket_t;  // dummy
  item_t *   coming_item = new item_t;
  size_t n = 0;
#ifdef __APPLE__
  size_t  read = 0;
  while ((coming_item->key = fgetln (fp, &read)) != NULL) {
#else
  ssize_t read = 0;
  size_t  size = 0;
  while ((read = getline (&coming_item->key, &size, fp)) != -1) {
#endif
    * (coming_item->key + read - 1) = '\0';
    item_t *   item   = 0;     // item to be updated
    bucket_t * bucket = root;  //      its bucket
    hash_t::iterator it = hash.find (coming_item);
    if (it == hash.end ()) { // not monitored
      if (hash.size () == k) { // replace with
        bucket = bucket->next;
        item = bucket->get_child ();  // the least significant item
        hash.erase (hash.find (item));
      } else { // new item
        item = new item_t;
      }
      item->reset_item (coming_item->key, read, bucket->get_count ());
      hash.insert (item);
    } else {  // monitored
      item  = *it;
      bucket  = item->get_parent ();
    }
    const size_t new_count  = bucket->get_count () + 1;
    if (bucket->next->get_count () == new_count) { // push to the sibling
      if (bucket != root) {
        if (item->is_singleton ()) { // pop singleton
          bucket->pop ();
          bucket_t * tmp = bucket->prev;
          std::swap (bucket, tmp);
          bucket_pool.recycle_bucket (tmp);
        } else {
          bucket->pop_item (item);
        }
      }
      bucket->next->push_item (item); // move bucket to bucket->next
    } else if (item->is_singleton ()) { // reuse singleton
      bucket->set_count (new_count);
    } else { // prepare a new bucket
      if (bucket != root) bucket->pop_item (item);
      bucket_t * new_bucket = bucket_pool.create_bucket ();
      new_bucket->set_count (new_count);
      new_bucket->set_child (item);
      new_bucket->push (bucket);
    }
    ++n;
  }
  std::fclose (fp);
  std::fprintf (stderr, "%ld tokens, %ld space.\n", n, k);
  
  // print
  for (bucket_t * bucket = root->prev; bucket != root; ) {
    item_t * item = bucket->get_child ();
    do {
      std::fprintf (stdout, "%s %ld %ld\n",
                    item->key, bucket->get_count (), item->get_delta ());
      item_t * tmp = item; item = item->next; delete tmp;
    } while (item != bucket->get_child ());
    bucket_t * tmp = bucket; bucket = bucket->prev; delete tmp;
  }
  delete root;
  delete coming_item;
  return 0;
}

stream summary は，頻度情報を保持する bucket に，同じ頻度を共有する item（誤差情報つき）の集合を対応させたデータ構造．アルゴリズム自体は単純なのだけど，実装上は更新時の場合分けがやや煩雑*1．正直なところ，c++ のプログラムとしてはクラス設計がイマイチな気もするけど（取り敢えず動いてるからいいか）．

> run ./space_saving unigram_raw.txt 0.01

実行してみる（カウンタのサイズ k を越えるところは NA と表示）．

// compiled with gcc 4.6.0 20101113 (-O2 -march=core2 -m64)
---------------------------------------------------------------------
ε                 |   0.01  |   0.001 |  0.0001 | 0.00001 |  exact
k                 |     100 |    1000 |   10000 |  100000 |  1594517
---------------------------------------------------------------------
Nε                | 2790250 |  279025 |   27902 |    2790 |
---------------------------------------------------------------------
memory [MiB]      |   0.8   |   1.0   |   2.6   |   13.7  |  106.5  
memory [MiB](key) |  <0.01  |   0.06  |   0.52  |    3.46 |   28.22
memory [MiB](key)*|  <0.01  |   0.07  |   0.59  |    4.20 |   39.19
time   [s]        |  37.5   |  34.4   |  29.9   |   41.9  |   28.9
---------------------------------------------------------------------
Rank 1-5
の　　　　　　　　　　|       ->|       ->|       ->|       ->| 13406476
、　　　　　　　　　　|       ->|       ->|       ->|       ->| 12466310
。　　　　　　　　　　|       ->|       ->|       ->|       ->|  9253307
に　　　　　　　　　　|       ->|       ->|       ->|       ->|  8296578
は　　　　　　　　　　|       ->|       ->|       ->|       ->|  7891964
Rank 10^1+1-5
と　　　　　　　　　　|       ->|       ->|       ->|       ->|  4144898
し　　　　　　　　　　|  3688148|       ->|       ->|       ->|  3688134
（　　　　　　　　　　|  2942772|       ->|       ->|       ->|  2942768
）　　　　　　　　　　|  2852888|       ->|       ->|       ->|  2852884
年　　　　　　　　　　|  2833216|       ->|       ->|       ->|  2814393
Rank 10^2+1-5
において　　　　　　　|       NA|       ->|       ->|       ->|   178891
名　　　　　　　　　　|       NA|   178428|       ->|       ->|   178165
上　　　　　　　　　　|       NA|       ->|       ->|       ->|   178124
町　　　　　　　　　　|       NA|   180929|       ->|       ->|   177584
性　　　　　　　　　　|       NA|       ->|       ->|       ->|   177481
Rank 10^3+1-5
2010　　　　　　　　 |       NA|       NA|       ->|       ->|    24094
含め　　　　　　　　　|       NA|       NA|       ->|       ->|    24074
特定　　　　　　　　　|       NA|       NA|       ->|       ->|    24072
値　　　　　　　　　　|       NA|       NA|       ->|       ->|    24035
飛行　　　　　　　　　|       NA|       NA|    24029|       ->|    24022
Rank 10^4+1-5
武蔵野　　　　　　　　|       NA|       NA|       NA|       ->|     1818
文字通り　　　　　　　|       NA|       NA|       NA|       ->|     1818
基部　　　　　　　　　|       NA|       NA|       NA|       ->|     1818
前部　　　　　　　　　|       NA|       NA|       NA|       ->|     1818
優等　　　　　　　　　|       NA|       NA|       NA|       ->|     1818
--------------------------------------------------------------------

memory [MB] (key) は単語の字面の保存に使用するメモリのサイズを std::strlen で測ったもの．key に対して使用メモリ量が大きく見えるのは，メモリ確保が 16 bytes 単位だからかもしれない（* が 16 bytes 単位のメモリ確保を考慮した補正値）．この実装では，メモリの動的な確保をなるべく避けるためにやや富豪的にメモリを使っているので，更新時間を追求するのでなければ，使用メモリはもう少し減らすことができる．
例によって誤差は Nε より大幅に抑えられている．このデータでは，頻度が top-k の item が欲しい時は，その10倍ぐらいのサイズのカウンタを用意すれば (ε = 1/10k)，十分という感じ（なお，表には載せていないが，item ごとに見積もり誤差も得られる）．
lossy counting の結果も載せようかと思ったけど，疲れたのでやめた．結果だけ言うと，このデータでは頻度分布が power law に従っている関係で，ε が小さくなるにつれカウンタの最大数 kmax が 1/ε よりも大きく抑えられた（例えば，ε=0.01 -> kmax=128 なのに対し，ε=0.00001 -> kmax=21677; その分，低頻度の単語については頻度の誤差が大きくなるが）*2．space saving で lossy counting のカウンタ最大数と揃えるように ε を設定すると，精度は大差ない感じ．この手の skew なデータでは，lossy counting でも space saving でもどちらでも使いたい方を使えば良いのではないかな．*3．

*1:観測された item がカウンタに登録されているか否か，更新対象の item が bucket に所属する唯一の item か否か，更新語の頻度を持つ bucket が既にあるかどうか，などの場合分けがあり，これらの組み合わせによってデータ構造に対して必要な操作が変わる．

*2:skew でないデータが入力される場合，lossy counting が必要とするカウンタのサイズは O(1/ε log Nε)．

*3:大規模データで単語の数を数える - ny23の日記 のサーベイ論文で用いられている lossy counting の実装は，1/ε 単位でデータを区切って qsort というデータの skewness を無視した超富豪的実装になっているので，結果はあまり参考にしないほうが良いと思う（skew なデータを相手にする場合，ハッシュを使った方が，実装も平易だし，更新速度の点でも使用メモリの点でも有利で良いような気がする）．