自然言語処理研究会
id:nokunoさんが主宰する第2回自然言語処理勉強会@東京で"Latent Dirichlet Allocation入門"というタイトルで発表してきました。
内容としては機械学習ライブラリMalletに実装されているLDAのマルチスレッド実装クラスのParallelTopicModelで使われているトピックモデルの技術を紹介するという話でした。
本当は文章検索への応用とかの話もしたかったのですが準備に時間が足りず断念
[機械学習] PRML 14章の混合正規回帰モデルの実装
PRML 14.5.1の混合正規回帰モデルのEMアルゴリズムによるパラメータ推定を実装してみた。
一本の直線でフィッティングする代わりに
複数の直線でフィッティングするというモデルである。
これは混合正規分布の推定と同じくEMアルゴリズムで推定できる。詳しくはPRMLの14.5.1を参照。
package chapter14; import java.util.Random; public class Main { static double generate(double x){ if(Math.abs(x) >= 0.5){ return - 0.5 * x - 1; }else{ return 0.5 * x + 1; } } public static void main(String[] args) throws Exception{ Random rand = new Random(222); int N = 20; double Phi[][] = new double[N][2]; double target[] = new double[N]; for(int i = 0 ; i < N ; ++i){ double x = rand.nextDouble() * 2 - 1; double y = generate(x) + rand.nextGaussian() * 0.1; Phi[i][0] = 1; Phi[i][1] = x; target[i] = y; System.out.println(x+"\t"+y); } MixtureOfLL ml = new MixtureOfLL(2, Phi, target, rand.nextLong()); for(int iter = 0 ; iter < 30 ; ++iter){ System.out.println(iter); System.out.println(ml.W[0][1] + " * x + " + ml.W[0][0]); System.out.println(ml.W[1][1] + " * x + " + ml.W[1][0]); ml.iteration(); } } }
package chapter14; import java.util.Random; import org.apache.commons.math.linear.Array2DRowRealMatrix; import org.apache.commons.math.linear.LUDecompositionImpl; import org.apache.commons.math.linear.RealMatrix; import static java.lang.Math.*; public class MixtureOfLL { Random rand; int N , K; double gamma[][]; double pi[]; double betaInv; double W[][]; // design matrix double Phi[][]; Array2DRowRealMatrix PhiMatrix; double target[]; public MixtureOfLL(int classNum , double Phi[][] , double target[], long seed) { int dataNum = Phi.length; int featureNum = Phi[0].length; pi = new double[classNum]; W = new double[classNum][featureNum]; gamma = new double[dataNum][classNum]; this.Phi = Phi; this.target = target; this.N = dataNum; this.K = classNum; rand = new Random(seed); PhiMatrix = new Array2DRowRealMatrix(Phi); initParameter(); } void initParameter(){ betaInv = 1.0; double sum = 0.0; for(int k = 0 ; k < pi.length ; ++k){ pi[k] = rand.nextDouble(); sum += pi[k]; } double sinv = 1.0 / sum; for(int k = 0 ; k < pi.length ; ++k){ pi[k] *= sinv; } for(int k = 0 ; k < W.length ; ++k){ for(int f = 0 ; f < W[k].length ; ++f){ W[k][f] = rand.nextDouble() - 0.5; } } } double normal(double x , double mu , double sigma){ return 1.0 / sqrt(2.0 * PI * sigma * sigma) * exp( -(x-mu) * (x - mu)/(2*sigma * sigma)); } void optimizeBeta(){ double s = 0.0; for(int n = 0 ; n < gamma.length ; ++n){ double phi[] = Phi[n]; for(int k = 0 ; k < gamma[n].length ; ++k){ double mean = 0.0; for(int f = 0 ; f < phi.length ; ++f){ mean += W[k][f] * phi[f]; } s += gamma[n][k] * (target[n] - mean) * (target[n] - mean); } } betaInv = s / N; } void optimizePi(){ for(int n = 0 ; n < gamma.length ; ++n){ for(int k = 0 ; k < gamma[n].length ; ++k){ pi[k] += gamma[n][k]; } } for(int k = 0 ; k < pi.length ; ++k){ pi[k] /= N; } } void optimizeW(){ for(int k = 0 ; k < K ; ++k){ double R[][] = new double[N][N]; for(int n = 0 ; n < N ; ++n) R[n][n] = gamma[n][k]; Array2DRowRealMatrix Rmat = new Array2DRowRealMatrix(R); RealMatrix PhiR = PhiMatrix.transpose().multiply(Rmat); double[] prt = PhiR.operate(target); RealMatrix PRP = PhiR.multiply(PhiMatrix); // solve (PRP)W_k = prt for W_k LUDecompositionImpl luImpl = new LUDecompositionImpl(PRP); W[k] = luImpl.getSolver().solve(prt); } } void eStep(){ for(int n = 0 ; n < gamma.length ; ++n){ double sum = 0.0; double phi[] = Phi[n]; for(int k = 0 ; k < gamma[0].length ; ++k){ double mean = 0.0; for(int f = 0 ; f < phi.length ; ++f){ mean += W[k][f] * phi[f]; } gamma[n][k] = pi[k] * normal(target[n] , mean , betaInv); sum += gamma[n][k]; } double sinv = 1.0 / sum; for(int k = 0 ; k < gamma[0].length ; ++k){ gamma[n][k] *= sinv; } } } void mStep(){ optimizeBeta(); optimizePi(); optimizeW(); } void iteration(){ eStep(); mStep(); } }
[プログラミング]教師なし形態素解析
なんとなくid:nokunoさんが紹介してたNLTK Bookに乗ってる教師なし形態素解析を実装してみた。
http://d.hatena.ne.jp/nokuno/20100124/1264319152
切れ目の部分をflipした時にスコアをデータ全部見て再計算する代わりに、必要な部分だけ更新するようにしてみた。これでAlice's Adventures in Wonderlandぐらいのデータに対しては適応できるようになったけど、正直精度は微妙でした。
実行例
doyou see thekitty see thedoggy doyou like thekitty like thedoggy
プログラム
package segmentation; import java.util.*; public class Main { public static void main(String[] args) throws Exception{ String text = "doyouseethekittyseethedoggydoyoulikethekittylikethedoggy"; long ts = System.currentTimeMillis(); Segmentator seg = new Segmentator(text , 777); seg.run(1000000); long te = System.currentTimeMillis(); System.err.println("time = " + (te - ts) + " msec"); List<String> list = seg.getResult(); for(String l : list){ System.out.println(l); } } }
package segmentation; import java.util.*; public class Segmentator { Map<String, Integer> dict; int wordNum; int dictSize; String text; boolean[] z; Random rand; public Segmentator(String text, long seed) { this.text = text; int N = text.length(); z = new boolean[N - 1]; rand = new Random(seed); for (int i = 0; i < z.length; ++i) { z[i] = rand.nextBoolean(); } dict = new HashMap<String, Integer>(); init(); } int count(String word) { if (dict.containsKey(word)) { return dict.get(word); } else { return 0; } } void add(String word) { if (dict.containsKey(word)) { dict.put(word, dict.get(word) + 1); } else { dictSize += word.length(); dict.put(word, 1); } } void delete(String word) { int val = dict.get(word); if (val == 1) { dictSize -= word.length(); dict.remove(word); } else { dict.put(word, val - 1); } } void init() { wordNum = 1; for (int i = 0; i < z.length; ++i) { if (z[i]) { ++wordNum; } } int prev = 0; dictSize = 0; for (int i = 0; i < z.length; ++i) { if (z[i]) { String word = text.substring(prev, i + 1); add(word); prev = i + 1; } } String word = text.substring(prev); add(word); } void flip(int index) { int lstart = -1; for (int i = index - 1; i >= 0; --i) { if (z[i]) { lstart = i; break; } } String left = text.substring(lstart + 1, index + 1); int rend = z.length; for (int i = index + 1; i < z.length; ++i) { if (z[i]) { rend = i; break; } } String right = text.substring(index + 1, rend + 1); if (z[index]) { delete(left); delete(right); add(left + right); wordNum--; } else { delete(left + right); add(left); add(right); wordNum++; } z[index] = !z[index]; } boolean accept(double energy, double eneryNeighbor, double temperature) { if (eneryNeighbor < energy) { return true; } else { double r = Math.exp((energy - eneryNeighbor) * 0.5) * temperature; return rand.nextDouble() < r; } } void run(int iterMax) { double best = score(); boolean zBest[] = null; for (int iter = 0; iter < iterMax; ++iter) { int ind = rand.nextInt(z.length); double temp = (iterMax - iter) * 1.0 / iterMax; double current = score(); flip(ind); double next = score(); if (!accept(current, next, temp + 1e-8)) { flip(ind); } double score = score(); if (score < best) { best = score; zBest = z.clone(); } } if (zBest != null) { z = zBest.clone(); } } double score() { return dictSize + wordNum; } List<String> getResult() { List<String> list = new ArrayList<String>(); int prev = 0; for (int i = 0; i < z.length; ++i) { if (z[i]) { String word = text.substring(prev, i + 1); list.add(word); prev = i + 1; } } String word = text.substring(prev); list.add(word); return list; } }
Hadoopを使わずにWikipediaのテキスト処理を400倍高速化
タイトルは釣りです。id:mamorukさんの書いたHadoop で Wikipedia のテキスト処理を900倍高速化 - 武蔵野日記を読んで、そもそも1G程度のデータの単語頻度を数えるのに858分もかかるんだっけと思い、id:nokunoさんの資料を読んでみると単語頻度を求める際に
a b a a
みたいなデータを
a 3 b 1
に変形するのにsortしたファイルをuniq -cで処理するということをやっていた。これはあまり効率のよい方法ではなくて行数をNとしたときにO(N log N)の計算時間となる(文字列比較はO(1)でやれることにする)。
これに対して、単語の頻度をハッシュ表で保存すると理想的な条件の元ではO(N)の計算時間で頻度を求めることが出来、より高速に計算することが可能となることが期待される。
また、単語数をWとしたとき、C++のmapのような二分探索木を使ってもO(N log W)とsortしてuniqするのに比べ、より低い計算量を達成できる。
単語の頻度表を求めるコードをC++のunordered_map(gccのバージョンによっては使えない場合があります、その場合は適宜boost/unordered_mapなどに置き換えてください)を使って書くと以下のようになる
#include <unordered_map> #include <string> #include <iostream> using namespace std; int main(){ unordered_map<string , int> m; string line; while(getline(cin , line)) m[line]++; unordered_map<string , int>::iterator it; for(it = m.begin(); it != m.end(); ++it) cout << it->first << " " << it->second << endl; return 0; }
これを利用して日本語Wikideiaの生のunigramデータから単語頻度を求めると1分50秒で終了。コア数は1(CPUはCore i7 920, 2.67GHzでメモリが3G)でやってるので大体400倍ぐらい高速になっている計算になる。ちなみにunordered_mapをstd::mapに置き換えた場合4分49秒となった。
計算機を大量に並べてHadoopで分散処理というのも非常に価値のあることとは思いますが、会社によってはサーバを1台買うのにも非常に煩雑な手続きが必要で、しかも注文してから半年ぐらい来なかったり、開発環境のサーバが1台しかなかったりするので、こういった1台のサーバで処理を高速化するテクニックもまだまだ重要なのではないかと思いました。
あと、結局サーバを並べるといってもデータが2倍に増えたときにサーバが4倍必要になるとか言う状況だと、いつかは破綻するので計算量の見積りは大変重要です。
追記
id:n_shuyoさんに900倍を越えてほしかったといわれたので多少の高速化を試みてみた。データが1G程度しかないのでファイルの中身を全部メモリに読み込む+getlineを使わずに自力で改行をパースするという方針で書いた
#include <unordered_map> #include <map> #include <string> #include <iostream> #include <stdio.h> #include <sys/types.h> #include <sys/stat.h> #include <fcntl.h> #include <unistd.h> using namespace std; int main(){ struct stat st; char * fileName = "unigram_raw.txt"; if(stat(fileName , &st)){ perror("stat"); return 1; } int fd = open(fileName , O_RDONLY); char * buffer = new char[st.st_size]; if(read(fd , buffer , st.st_size) != st.st_size){ return 1; } int prev = 0; unordered_map<string , int> m; for(size_t i = 0 ; i < st.st_size ; ++i){ if(buffer[i] == '\n'){ string s(buffer + prev , buffer + i); m[s]++; prev = i + 1; } } unordered_map<string , int>::iterator it; for(it = m.begin(); it != m.end(); ++it) cout << it->first << " " << it->second << endl; return 0; }
これだと手元の環境では62秒で終わり、800倍以上(900は行かなかったorz)の高速化になっている。
追追記
bufferを全部持っている以上、string => intのテーブルを持つのではなく、(文字列の開始位置,文字列の終了位置)をキーにしたほうが効率が良さそうと思い実装。
struct myeq : std::binary_function<pair<int , int> , pair<int , int> , bool>{ char * buf; myeq(char * b) : buf(b) {} bool operator() (const pair<int , int> & x , const pair<int , int> & y) const{ char * s = buf + x.first; char * t = buf + y.first; char * se = buf + x.second; // *se == '\n' for( ; *s == *t && s != se ; ++s , ++t) ; // *s != *t or (*s == '\n' => *t == '\n') return *s == *t; } }; struct myhash : std::unary_function<pair<int , int> , size_t>{ char * buf; myhash(char * b) : buf(b) {} size_t operator()(const pair<int , int> & p) const{ size_t h = 0; for(int i = p.first ; i < p.second ; ++i){ h = h * 31 + buf[i]; } return h; } }; int main(){ struct stat st; char * fileName = "jawiki_unigram_raw.txt"; if(stat(fileName , &st)){ perror("stat"); return 1; } int fd = open(fileName , O_RDONLY); char * buffer = new char[st.st_size]; if(read(fd , buffer , st.st_size) != st.st_size){ return 1; } int prev = 0; myhash h(buffer); myeq e(buffer); unordered_map<pair<int , int> , int , myhash , myeq> m(3 , h , e); for(size_t i = 0 ; i < st.st_size ; ++i){ if(buffer[i] == '\n'){ pair<int , int> p(prev , i); unordered_map< pair<int , int> , int , myhash , myeq>::iterator it = m.find(p); if(it != m.end()){ it->second++; }else{ m.insert(pair<pair<int,int> , int>(p , 1)); } prev = i + 1; } } unordered_map<pair<int , int> , int>::iterator it; for(it = m.begin(); it != m.end(); ++it){ string s(buffer + it->first.first , buffer + it->first.second); cout << s << " " << it->second << endl; } return 0; }
これだと手元の環境では49秒で終わり、無事900倍達成できたのであった。めでたし。
Power Iteration Clustering
岡野原さんのtweetで紹介されていたPower Iteration Clusteringという文章分類の手法に関する論文[1,2]を読んでみた。
背景
n個のデータX={x_1,...,x_n}が与えられたときに各データ間の類似度s(x_i,x_j)を成分に持つ類似度行列Aを考える。
また次数行列としてAのi行目の値を合計したd_{ii} = \sum_j A_{ij}を対角成分にもつ対角行列をDとする。
このときW:=D^{-1} Aをnormalized affinity matrixと定義する。簡単のためWはフルランクであるとする。
この行列はすべての要素が1となる固有ベクトルをもち、この時固有値は1となる。実はこれが最大固有値である(行列Aの行和が1となること+Gershgorin circle theorem(en)より導かれる)。また、行列Wの固有値を1=λ_1>=...>=λ_nとすると行列L=I - Wの固有値は0=1-λ_1<=...<=1-λ_nとなり、Wの最大固有値をk個を求めることとLの最小固有値をk個求めることは等価となる。ここでLはグラフラプラシアンマトリックスと呼ばれるスペクトルクラスタリングと呼ばれる手法[4]で用いられている行列である。
Wが下の図(http://ranger.uta.edu/~chqding/Spectral/spectralA.pdfより引用)のようにクラス内ではXが同じクラス内では類似度が互いに高く、クラス間では類似度が低いというデータとなってる場合、Wの最大固有値のうち自明解を除き、上からk-1個とってきたときに対応する固有ベクトルが(1,1,1,1...,1,0,0,0...)と(0,0,0,0...,0,1,1,1...)のようなクラス所属を表すような値となっている(完全にクラス間の値が0となっている場合の話)。
上の話をもう少し定式化すると固有ベクトルをe_1,...,e_nとしたときに固有ベクトルの第a成分をk個並べたベクトルs_a =
と定義する。Wがk個のクラスタを持つときにspec(a,b)が小さければa,bが同じクラスに続し,逆も成り立つ[3]。
Power Iteration
行列の最大固有ベクトルを求める方法としてべき乗法(Power Iteration)という方法がある。これは適当なベクトルv^0=\sum_i c_i e_iに対してWをt回かけたベクトルが
となり、最大の固有ベクトルe_1が優越してくることを利用した方法である。
また、第2固有ベクトルを求める際にはv^0からe_1の成分を除外したものを初期値として同じ反復を繰り返すことにより原理的には求めることができ、以降の固有値も同様に求めることができるが、この方法は数値的に不安定であることが知られている(数値計算の途中で丸め誤差が入り取り除いたはずのe_1の成分が混じってくるため)。
Power Iteration Clustering
Power Iteration ClusteringはPower Iterationのv^tの反復途中の値について注目する方法である。今v^tの第a成分と第b成分の差を
と書く。このときpic(a,b)の値は
となる。Wがk個のクラスタに分けれるときk+1以降の固有値は速やかに0へ収束するため、
と近似できる。
ここでspec(a,b)が0に近いときはpic(a,b)も0に近いことはすぐに分かる。逆についてはe_i(a)-e_i(b)の値が大きくてもお互いにcancel outしてpic(a,b)が0に近くてもspec(a,b)が0に近いとは限らないが[1]ではそのようなことは"uncommon in practice"であると主張している。
上の議論により、同一クラスのデータについてはpic(a,b)が小さくなるようなクラスタリングを行えばよいことが分かるが、これはvの各成分の値についてK-meansを行えばよい(1次元データのクラスタリングとなる)。
文章分類への応用[2]
文章-頻度行列をFとしたとき類似度行列をA=F F^Tとすることによって、文章のクラスタリングを行なうことができる。ここでAの値を明示的に計算するとO(n^2)のメモリが必要となるが、Aとvの行列・ベクトル積が計算できればよいので必要なメモリの量はFの要素数に比例する量で抑えることができる。
まとめ
Power Iteration Clusteringは行列のスペクトル分解に基づいたPower Iterationを途中でうち切った値を用いることにより、高速+高精度のクラスタリングを行なうことができる手法である。
また、参考までにPICを実装してみたものをgithubで公開した。 http://github.com/tsubosaka/PIC
いくつかのデータについて実験してみたときkが多い(といっても5ぐらい)のときにあまり精度がよくないように思えた。
参考文献
[1] Frank Lin, William W. Cohen. Power Iteration Clustering, In Proc ICML 2010 to be appear
[2] Frank Lin, William W. Cohen. A Very Fast Method for Clustering Big Text Datasets,In Proc ECAI 2010 to be appear
[3] Marina Meil, Jianbo Shi. A Random Walks View of Spectral Segmentation, In Proc AISTATS 2001
[4] Ulrike von Luxburg, A tutorial on spectral clustering,Statistics and Computing, pp 396--416, 2007
[機械学習] PRML勉強会
PRML勉強会で11.4のスライスサンプリングについて発表してきました。
発表スライドは以下となります。
また、参考として以前のPRMLハッカソンで作成したスライスサンプリングを用いたLDAのコードをgithubにアップロードしました。http://github.com/tsubosaka/LDASlice/tree/master/src/lda/
[NLP][機械学習] 言語モデル覚え書き
この文章について
最近言語モデル方面にも少し興味があるので自分の知識を整理する意味で書いてみた。NLPは専門ではないので、おかしなことを書いてある可能性がありますがその場合はご指摘ください。
本文章ではn-gramモデル、単語の出現確率がn-1個前の単語のみに依存するモデルを考える。
問題
who is * という文が与えられたときに*にくる文字の確率を求めることを考える。この場合だと*には例えばheが当てはまるかもしれないが, isが入ることはまずなさそうに思える。このことは文法的にも説明ができると思うが、文法のルールを作るのは大変だし、文法的に正しい単語の中でどれが出やすいかということはできない。
一方で機械学習を使った言語モデルの文脈では文法的知識を余り持たず、与えられたコーパスから自動的に出やすい単語/表現を学習する方針をとる。
最尤推定
一番簡単なモデルとしては最尤推定を使うもので、たとえばコーパスの中にwho isの後にsheが6回、heが4回出ていたとすると
p(he|who is)=6/10, p(she|who is)=4/10
と推定する方法である。しかしながらこれだとコーパスに出てこなかった単語の確率は全て0になるという問題がある。
出現しなかった単語の確率が0になるのを防ぐだけであれば出現しなかった単語の出現回数に対して一律な数を加えることによって防げるが、これだとコーパス上でこの文脈で出現しなかった単語の確率は全ておなじということになり予測性能は悪い。
Dirichlet smoothing
そこでwho isの後に出現しなかった単語についてはより低次のn-gram、すなわちisの後に出現する単語の確率を考えるとよさそうに思える。
例えばxyzという列があったときに
と計算する、ここでc(・)は頻度を表す。
ノンパラを分かる人に説明するとこれはCRPとしてみなせて各テーブルに座っている人の数が単語頻度に相当して、集中度パラメータに応じた確率で新規テーブルを生成するのですがそこで選ばれる単語の確率分布がn-1gramの確率分布を使っているということです。
上の例を使うと各単語の出現確率は次のように書けます。
p(he|who is)=(6+\alpha p(he | is))/(10 + \alpha), p(she|who is)=(4 + \alpha p(she| is) /(10 + \alpha) p(w |who is)=(\alpha p(w | is))/ (10 + \alpha)
Kneser-Ney smoothing
上のDirichlet smoothingでは\alphaが小さいときに\alpha/(10+\alpha)というのがものすごく小さい値となるので頻度をある程度おさえるというのがKneser-Ney smootingです。例によってxyzという列が与えられた場合
となります。ここで|xy*|はxyに続く文字列の種類数です。追記(4/15):上の式のカッコが抜けてたので修正しました。
たとえば上の例ではwho isに続く文字種はhe/sheなので2となり、例えばd=0.5とすると
p(he|who is) =(6 - 0.5)/ 10 + p(he | is) / 10 p(she|who is)=(4 - 0.5)/ 10 + p(she | is) / 10 p(w |who is)= p(w | is) / 10
d |xy*| / c(xy*)の確率でn-1グラム文脈から単語が生起されることがわかり、これはDirichlet smoothingが\alpha / (c(xy*)+\alpha)の確率でn-1グラム文脈から単語を生起してたのと比較するとn-1グラム文脈の重みが強くなっていることが分かる。
Hierarchical Pitman-Yor Processes Language Model
これはちょっとうまく説明できないので雰囲気だけ述べると、上のKneser-Ney smoothingでは頻度1000の単語も頻度1の単語も同じディスカウント係数dを用いているがHPYLMでは文脈と単語に応じてディスカウント係数が違ってくる。これについては文脈(レストラン)において単語(料理)の頻度(客)が高い状況においてはその単語に関するテーブルの数が多くなり、その分ディスカウント係数が多くなるということが導かれる。
Further Topics
今who is *に来る単語に関しては2gram前、1gram前の単語しか考慮していないが、例えばwhere is *で出現する単語についても考えるとwhoとwhereは同じような機能を持つ単語のため精度が向上する可能性がある。これを表現するのに単語一つ一つに潜在語を考えるlwlmのようなモデルが役に立つ。
また、一段落前の文からwho is *の後に来る単語が一意に決まる可能性もあり、長距離の単語履歴を使えるモデルについてもそのうち調べようと思っている。