[機械学習] PRML Hackathon #2

PRML Hackathonに行ってきました。
今回のHackathonでは昨日書いたスライスサンプリングという手法をLDAの推論に使ってみて通常のGibbs samplerと比較してみました。

結果としてはサンプリング速度が2-3倍程度高速になり、手法の有効性を確かめることができました。また、perplexityの値は
Gibbs samplerよりも少し悪い結果となりました。(誤解を招く表現でした、詳しくは下に追記しました)

追記:
perplexityの値が悪くなったというと変分ベイズのように近似が入って性能が悪くなる印象を与えますがそうではないです。
slice samplerはgibbs samplerと同様に十分な回数反復すれば正しい確率分布からのサンプルを取ることができる方法です。
実験で言えることとしては単に同一iterationの条件ではgibbs samplerの方がperplexityが低くなっているということだけです。また、すべてのトークンに関して一回ずつサンプリングを行なうのを 1 iterationとしていますが、この 1 iterationの速度もslice samplerの方が数倍高速に行えるため、実時間で見た時のperplexityの収束速度についてもslice samplerの方が高速です。

[機械学習] スライスサンプリング

持橋さんのlwlmに関する記事を読んで、スライスサンプリング[1,2]というのが有用そうということで調べてみたのでメモ。

スライスサンプリング概要

今確率分布f(x)が
f(x) \propto l(x) \pi(x)
の形で与えられており、このf(x)からxをサンプリングすることを考える。ここでl(x)は非負の関数、\pi(x)は確率密度関数とする。

サンプリングを以下の手順で行なう

  1. 初期値x_0を適当に決定する
  2. u_t = Uniform(0 , l(x_t))で一様分布からサンプリング
  3. X_t = { x | u_t < l(x)}とする
  4. x_{t+1}をX_tに含まれるxから\pi(x)に従ってサンプリングする
  5. 2へ戻る

ここでu_tの値は捨てて、{x_t}だけ取り出すとf(x)に従うxをサンプリングできる。

何が嬉しいのか

スライスサンプリングの話は以前から聞いたことがあったのですが、連続の場合だと4の部分が簡単にできそうではなくあんまりありがたみが分からなかったのですが、離散の場合ですと結構使える方法となっている。

例えばlwlmの実装では隠れ変数h_iを
p(h_i |o_i, H) \propto p(h_i|H) p(h_{i+1}|h_i,H) ... p(h_{i+n}|h_i,H) p(o_i | h_i)
に従ってサンプリングしている。ここで便宜上h_i以外の隠れ変数の集合をHと書いた。lwlmにおいてはh_iの取りうる値が単語数と同じだけあり、上の値を計算してそれに従ってサンプリングするのは非常に効率が良くない。

そこで
l(h_i) =  p(h_i|H) , \pi(h_i)=p(h_{i+1}|h_i,H) ... p(h_{i+n}|h_i,H) p(o_i | h_i)
とすると上のスライスサンプリングが使えて、u = Uniform(0 , l(h_i))以上の潜在変数のみに関して\pi(h_i)を計算して、その確率に従いサンプリングすると候補数を少なくすることができる。u以上の値を列挙するには文脈上での潜在変数の頻度順*1にデータをソートしておくことにより効率的に行なうことができる。

離散分布においてこの手法は他にも文献[3]などで用いられている。

参考文献

  • [1] Damien, P. and Wakefield, J. and Walker, S., Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables, Journal of the Royal Statistical Society. Series B, Statistical Methodology, pp. 331--344, 1999
  • [2] Neal, R.M., Slice sampling, Annals of Statistics, vol 31, pp. 705--741, 2003
  • [3] Van Gael, J. and Saatci, Y. and Teh, Y.W. and Ghahramani, Z., Beam sampling for the infinite hidden Markov model, In Proceedings of the 25th international conference on Machine learning, pp. 1088--1095, 2008

*1:実際にはもう少し工夫をしているようです

JavaでPDFから文章を抽出

プログラム上からPDFの文章を取り出したいと思うことがあったので、方法を調べてみた。
PDFBoxというツールを使うと結構いい感じに抽出できた。
以下に簡単なサンプルプログラムを示す。

import java.io.*;

import org.apache.pdfbox.pdfparser.PDFParser;
import org.apache.pdfbox.pdmodel.PDDocument;
import org.apache.pdfbox.util.PDFTextStripper;

public class ExtractPDF {
  private static String extractText(String filePath)
      throws FileNotFoundException, IOException {
    FileInputStream pdfStream = new FileInputStream(filePath);
    PDFParser parser = new PDFParser(pdfStream);
    parser.parse();
    PDDocument pdf = parser.getPDDocument();
    PDFTextStripper stripper = new PDFTextStripper();
    ByteArrayOutputStream out = new ByteArrayOutputStream();
    stripper.writeText(pdf, new BufferedWriter(new OutputStreamWriter(out)));
    return out.toString();
  }
  public static void main(String[] args) throws Exception {
    String pdfFile = "test.pdf";
    String data = extractText(pdfFile);
    System.out.println(data);
  }
}

手元の適当な論文に対する出力結果

Online EM for Unsupervised Models
Percy Liang Dan Klein
Computer Science Division, EECS Department
University of California at Berkeley
Berkeley, CA 94720
fpliang,kleing@cs.berkeley.edu
Abstract
The (batch) EM algorithm plays an important
role in unsupervised induction, but it some-
(中略)
EM sEM‘ EM sEM‘
POS 56:2  1:36 58:8  0:73; 1:41  6:01  6:09
DOC 41:2  1:97 51:4  0:97; 2:82  7:93  7:88
SEG(en) 80:5  0:0 81:0  0:0; 0:42  4:1  4:1
SEG(ch) 78:2  0:0 77:2  0:0; 0:04  7:26  7:27
ALIGN 79:0  0:14 78:8  0:14; 0:25  5:04  5:11
Table 2: Mean and standard deviation over different ran-
dom seeds. For EM and sEM, the first number after  
(以下略)

参考

Luceneと併用して全文検索を行なうことも可能なよう

残念ながら日本語のPDFには対応してないようですが、パッチを当てれば読める場合もあるようです

PRML勉強会

PRML勉強会の発表資料公開しました。
自分の担当は10.1の変分推論のところでした。次回は10.2以降からでPRMLで計算が一番多いところなので予習を頑張らないとなといった感じです。

[プログラミング] Google Sparsehashを使うときの注意点

持橋さんの書かれたgoogle-sparsehashと自作のsplay-treeとの速度比較をした結果の記事を読んで、さすがに速度に200倍近くの差がでるのはおかしいだろうということで原因を探ってみた。

結論としてはGoogle Sparsehashを使うときに__gnu_cxx::hashを使わない方がよいということが分かった。

時間の測定に用いられているコードは概ね以下のコードと同じである。

#include <iostream>
#include <google/sparse_hash_map>
#include <cstdio>
#include <cstring>
#include <ext/hash_map>

using namespace std;
using google::sparse_hash_map;

typedef __gnu_cxx::hash<const char *> HashFunc;

struct eqstr
{
  bool operator() (const char *a, const char *b) const
  {
    return (a == b) || (a && b && strcmp(a,b) == 0);
  }
};

int read_word (FILE *fp, const char *buf)
{
  return fscanf(fp, "%s", buf);
}

int main (int argc, char *argv[])
{
  FILE *fp;
  char s[BUFSIZ];
  sparse_hash_map <const char *,int, HashFunc ,eqstr> lexicon;
  sparse_hash_map <const char *,int, HashFunc ,eqstr>::iterator it;

  if ((fp = fopen(argv[1], "r")) == NULL){
    perror("fopen");
    exit(1);
  }

  while (read_word(fp, s) != EOF){
    if ((it = lexicon.find(s)) == lexicon.end()){
      lexicon[strdup(s)] = 1;
    } else {
      it->second++;   // count++;
    }
  }
  fclose(fp);
  return 0;
}

これを例えば以下のperlスクリプトで生成した1から10万までの数字がランダムに1回ずつ現れるテキストで実行時間を測ると2.4秒かかり、50万件の場合は36秒かかる

use List::Util 'shuffle';
print "$_\n" for(shuffle((1 .. 100000)));

データを5倍にしたときに実行時間がリニアに増えないため、大量のハッシュの衝突が起きてそうだということが予想できる。
実際にハッシュ関数をtypedefしている部分を

typedef tr1::hash<const char *> HashFunc;

にして時間を測ったところ10万件の場合0.149秒、50万件の場合0.844秒かかった。

なぜこれだけかかるのかについてですが、 __gnu_cxx::hashハッシュ関数はext/hash_fun.hの中をみると

  inline size_t
  __stl_hash_string(const char* __s)
  {
    unsigned long __h = 0;
    for ( ; *__s; ++__s)
      __h = 5 * __h + *__s;
    return size_t(__h);
  }

となっており、別の方の実験でも非常に衝突が起こりやすいことが報告されており、あまりよくないハッシュ関数となっている。

一方でtr1::hashはfnv hashといわれる、高速で比較的ばらつきやすいハッシュ関数を利用しておりこれが実行時間の短縮につながっていると思われる。

gcc 4系列を使っている場合は./configure時にtr1::hashがデフォルトのハッシュとして設定されるが、そうではない場合 __gnu_cxx::hashが選ばれ、ハッシュ関数が適切でないためsparse_hash_mapの性能が低下することがありえます。これを回避する方法としては自分でハッシュ関数を定義するというもので例えば各種マップ実装の性能比較のテストで用いられている

struct stringhash {                      // hash function for string
  size_t operator()(const string& str) const {
    const char *ptr = str.data();
    int size = str.size();
    size_t idx = 19780211;
    while(size--){
      idx = idx * 37 + *(uint8_t *)ptr++;
    }
    return idx;
  }
};

を用いるという方法があります。

いくつかの検索結果を見たところ、ハッシュ関数が適切でないことからsparse_hash_mapが思ったような性能が出ないという記事がいくつか見受けられたので、もし性能がでないと思われたらハッシュ関数を変えることを検討してみたらどうでしょうか。

[機械学習] クラスタリングにおけるコサイン類似度に関する性質の証明

bayonやCLUTOが爆速な理由 - download_takeshi’s diaryを読んで、すぐには成り立つかどうか分からなかったので証明してみた。

上の記事で述べられていることはクラスタ中のベクトルとその中心ベクトルのコサイン類似度の和と、クラスタ中のベクトルを全て足したベクトルのノルムが一致するというである。
ただしここでクラスタ中の要素ベクトルはすべて大きさ1の規格化されたベクトルであるとする。

証明

クラスタ内に含まれるベクトルをx_1,\dots,x_N とする。
このとき全ベクトルを足しこんだ複合ベクトルを
 X_s = \sum_i x_i
とする。またこのクラスタのセントロイドは
 C   = \frac{1}{N} \sum_i x_i = \frac{X_s}{N}
となる。このときセントロイドと各ベクトルx_iとのコサイン類似度は
[tex: s_i = \frac{}{||C|| ||x_i||} = \frac{}{||{C}||}]
となる。ここで ||x_i|| = 1と正規化されていることを用いた。この類似度の合計は
[tex: S = \sum_i s_i = \frac{1}{||C||} \sum_i = \frac{1}{||C||} ]
となり、ここで X_s,Cの定義より
[tex: S = \frac{1}{||C||} = \frac{N}{||X_s||} <\frac{1}{N} X_s,X_s> = \frac{1}{||X_s||} = ||X_s||]
となり、複合ベクトルのノルムがコサイン類似度の合計に等しいことが示せた (Q.E.D.)。

[機械学習] PRML勉強会発表会資料

3/6の第11回PRML読書会と2/9の第10回PRML読書会の資料を今さらですが、SlideShareにアップしたのでリンク貼っておきます。