« March 2009 | Main | May 2009 »

2009.04.25

部分文字列の話

ここしばらく、部分文字列の統計量を利用した機械学習やデータマイニングをやっている。そこの話からちょっと抜粋。

長さnの文字列T[1,...,n]が与えられた時、T中に出現する部分文字列T[i...j] (1≦i≦j≦n)の数はn個の中からiとjの2箇所を選ぶのでO(n^2)個ある。例えば、n=10^6(1MB)だったら、部分文字列の数は約10^12個(1T)と非常に大きい。

しかし、これらの部分文字列の出現位置は同じである場合が多い。例えばT="abracadabra"であれば、"abra"と"abr"の出現場所は1番目と8番目であり、全く同じである。

では出現位置(部分文字列の左端を出現位置とする)が全く同じであるような部分文字列をまとめてグループにした場合、グループの数はいくつになるのだろうか。

これは接尾辞木(wikipedia 授業の資料)を知っているなら簡単に説明できる。

Tに対しT[i,...,n]、つまり頭の文字を削った文字列、をTの接尾辞と呼ぶ。この接尾辞から構築したtrieに対し、子が一つであるような節点を文字を残したまま取り除いて枝にしてしまったものが接尾辞木である(正確には接尾辞木のスケルトン。実際にはsuffix linkとかつける)。各葉には、接尾辞の出現位置を記録しておく。

この接尾辞木において、節点に対応する部分文字列の出現位置は、その節点の子孫の葉の集合に対応する。枝の途中に対応する部分文字列の出現位置も、その先の節点の出現位置と同じであり、同じ枝に対応する部分文字列の出現位置は同じである。

この接尾辞木の節点の数は、接尾辞を一個ずつ加えて接尾辞木を構築する過程を考えると、各追加で、高々一つしか節点は増えないので、高々n-1個である。

一回しか出現していない部分文字列のグループは、接尾辞木中の葉に隣接する枝に対応し、二回以上出現している部分文字列のグループは節点と節点の間の枝に対応している。葉の数は接尾辞の数と同じなのでn個。よって、グループの数は(葉の数)+(節点の数) = 2n-1、二回以上出現しているグループの数はn-1個と分かる。

最初の問いに答えると、部分文字列の数はn^2なのに、それらをまとめると高々n個しかない。ものすごく少ないことがわかる。

さらに、これらの極大部分文字列は左側に伸ばした場合の冗長も考えるとさらにまとめることができる。(たとえば先ほどの例はbraとabraの出現位置もおなじ)。

こうすると例えば、さっきの例のTにおける極大部分文字列で二回以上出現しているものはaとabraのみである。

では、実際にこれらのグループに対応する部分文字列(極大部分文字列)を取り出すにはどうすればいいのか。接尾辞木を構築するのが一番直感的だが、接尾辞木のサイズは非常に大きい。

現在では接尾辞配列と共通接頭辞配列を組み合わせて、線形時間で求めるのが一番手っ取り早い。
接尾辞配列の線形時間はSAISとか使えばできて、線形時間での共通接頭辞配列構築と、それを利用した極大部分文字列の列挙は[pdf]中の論文の疑似コードをそのまま書けばできる。左側の冗長性に関してはBWT配列に関するrankを使えば全部トータルでも線形時間で極大な部分文字列を求めることができる(ちょっと面倒なのでやらなくてもいいが)。

これでテキスト長が10億文字ぐらいになっても10分程度で全ての極大な部分文字列を抜き出すことができる。

出現位置がまとめられたらいろいろ使いどころはある。例えば、私が今やっている全ての部分文字列を利用した文書分類では、全ての文書をつなげた文字列から極大部分文字列を抽出し、そこから各部分文字列に関する統計量(尤度の微分とか)を効率良く計算することができている。

ちょっと毛色の違う例として、例えば、アクセスログのパターンを解析をする場合を考えてみる。ログデータとして、アクセスしたサイトの記録がセッション毎に残されているとする。各セッション中に含まれるサイト数の平均が約百で、サイトの総種類数が1万、セッションが1000万回分記録されているとする。ここから10回以上出現しているアクセスパターン(サイト列)を、そのパターンの長さを問わずに抽出しろというお仕事があったとする。

見積もろうとすると、まず全通りのパターンを力づくで記録するのは、長さ2のアクセスパターンで1万×1万=1億通りのパターンがでる可能があるので長さ10とかはまず無理である。
ちょっと考えて、実際に出現したパターンをハッシュとかで記録して求めようとすると、各セッション毎に一万個(百×百)のアクセスパターンを記録する必要があり、合計で1000億個のアクセスパターンを記録する必要がある。とてもメモリに載らないので10台ぐらい必要、もしくはDBとかでガリガリで1日ぐらいかかるのではないかということになる。まず長さを切ってしまうだろう。

この場合のスマートな回答は、全てのセッションをつなげた長い配列を作り、そこから極大部分文字列を抽出し、10回以上出現する極大部分文字列のみを報告すればよい(極大部分文字列に含まれる部分文字列は全て頻出)。このときセッションを跨いだパターンを抽出しないようにセッション間には特殊文字をいれておけばよい。これらの特殊文字は全て違うものにしといてもよい。こうすると文字種類数も数千万となるが今の接尾辞配列構築アルゴリズムの計算量はアルファベット数には依存しない。1台のマシンでしかも10分程度で漏れなくパターンが抽出できますよといえば、驚かれるかもしれない
(実際できる見積もりを、そのまま言っていいかどうかは、また考える必要があるが)。

| | Comments (137) | TrackBack (0)

2009.04.16

オンラインEMアルゴリズム

EMアルゴリズム(Expectation Maximizationアルゴリズム、期待値最大化法、以下EMと呼ぶ)は、データに観測できない隠れ変数(潜在変数)がある場合のパラメータ推定を行う時に有用な手法である。

EMは何それという人のために簡単な説明を下の方に書いたので読んでみてください。

EMのきちんとした説明なら持橋さんによる解説「自然言語処理のための変分ベイズ法」や「計算統計 I―確率計算の新しい手法 統計科学のフロンティア 11」が丁寧でわかりやすい。

EMは教師無学習では中心的な手法であり、何か観測できない変数を含めた確率モデルを作ってその確率モデルの尤度を最大化するという枠組みで、観測できなかった変数はなんだったのかを推定する場合に用いられる。
例えば自然言語処理に限っていえば文書や単語クラスタリングから、文法推定、形態素解析、機械翻訳における単語アライメントなどで使われる。

EMは大量のデータを使えば使うほど精度があがり、しかもデータは教師無でよいので事実上使い放題である。結果として計算量は大きくなる。

EMはMap Reduceの枠組みでの大規模並列化しやすいのが知られているが[pdf:chu nips2006] [pdf:Wolfe icml2008]、それでも重い。
で、高速化の話がちらほらでていたが、教師有のオンライン学習のように、データ毎に更新してもうまくいくという研究報告がでてきた。

"Online EM for Unsupervised Models", [pdf] P. Liang D. Klein, NAACL 2009 to appear

期待値を求める部分で、全てのサンプルを使って求めるのではなく
(1) Incremental EM (iEM), 今の期待値の計算部分の、各データに対応する部分を一個ずつ取り替えていく
(2) Stepwise EM (sEM), 期待値の計算部分を、確率的に近似し(確率的勾配法に似ている)、今持っている期待値を補間をとりながら更新していく

の二通りで近似する。どちらも実装は単純であり、データ数に比例する計算量で行える。
ただ、iEMは全てのデータを保存しておかないのが問題であり、sEMではそれは必要なく、必要な作業領域量はデータ数に依存せずスケーラブルである。ただ、sEMの場合確率的に近似する部分でのステップ幅の設定が必要である。だがそれも彼らの実験結果によれば、大部分でロバストに動くらしい。

実験では品詞推定、単語分割、機械翻訳における単語アライメントで試しており単語分割などにおいて2回ぐらいで収束し(バッチだと10回ぐらい)、品詞推定の場合では20回ぐらい(バッチだと100回ぐらい)で収束し、大体バッチEMと同程度の精度が達成可能である。

これによる簡単な拡張も考えられるだろう。
例えば、EMにおいて、期待値を求める部分で、そのままの値を使わずに定数を引いたものを使って期待値を求めれば、Dirichlet Processの平均場近似と同じ効果が得られることが知られている([pdf:percy ACL tutorial 06])。これと同じことをオンラインEMでやろうと思えば、ちょっとずつ値を引けばいいだろう。


以下EM法の解説。自分で書いていていろいろ思い出した。
--
はじめに多項分布の最尤推定からはじめる。箱の中に赤、青、黄の三色の球がそれぞれ何個か入っていて、これらの数を推定したい。実際に100回試しにとってみて(毎回球は戻すとする)、赤、青、黄の数がそれぞれ50個、20個、30個だったとすると、直感的には箱の中の球の割合はそれぞれ0.5、0.2、0.3であると推定される。この推定は最尤推定といって、実際に観測した事象が最も起きやすい割合がこうだから、そう推定しているのである。数式で書けばパラメータθ(今回これは球の割合)を推定するにあたり、観測した事象Xの確率p(X|θ)を最大にするようなθを求めるのが最尤推定である

(これに対し、p(θ|X)を最大にするようなθを求めるのは、ベイズの定理よりp(θ|X) = P(X|θ)P(θ)/P(X)であり、P(X)はθに依存しないので、結果としてP(X|θ)P(θ)を最大にするようなθを求めればよい。P(θ)はθに対する事前分布であり、大抵、「こういうθがいいなぁ」と人が決める。この推定を事後確率最大化推定、MAP推定と呼ぶ)

先ほどの場合の最尤推定を真面目にしてみると、箱中の赤、青、黄が入ってある割合をそれぞれp(r), p(b), p(y)とした時、観測された事象の確率(組み合わせ数は無視)p(r)^50 * p(b)^20 * p(y)^30を最大にするのはp(r)=0.5, p(b)=0.2, p(y)=0.3の時なので、そう推定しているのである(a^bはaのb乗)。きちんと求めようと思えば、これの対数をとってp(r)+p(b)+p(y)=1の拘束条件をとって、最大化をラグランジュの未定乗数法を使って求めればでてくる。

#多項分布の場合の最尤推定は人間は自然にできるので、当たり前だといわれるが、まじめにやると結構大変

次に、箱が二つある場合を考える。それぞれの箱には違う割合で赤、青、黄の球が入っている。
一回の試行では、ある確率pで片方の箱を選び、(1-p)で違う方の箱を選ぶ。そして次にその箱から球を10個取り出す。

この試行を6回繰り返した時の、それぞれの赤、青、黄の数は以下のようであったとする。しかし、どちらの箱を選んだのかは分からない。

1回目 (2, 5, 3)
2回目 (7, 2, 1)
3回目 (2, 7, 1)
4回目 (1, 6, 3)
5回目 (2, 4, 4)
6回目 (7, 1, 2)

このような中でパラメータ(どちらの箱を選ぶかの確率p、及びそれぞれの箱の中の球の割合)を求めることを考える。

例えば、ウェブページのセッション毎のアクセスログがあり、各ページへのアクセス回数を記録していたとする。データを見ているとどうやら男性と女性のグループでアクセス傾向が違うと思われるのだが、それぞれのセッションでは男性か女性か分からないような場合を想定してほしい。こんな中でもそれぞれのアクセスがどっちで、それらのがどんなウェブページかを推定したいのが今回のケースにあたる。

このデータをよくみてみると、片方の箱には青が多めに入っていそうであり(1,3,4,5回目)、片方の箱には赤が多めに入っていそうである(2, 6回目)。そして、箱を選ぶ確率pはなんとなく4:2なので4/6になりそうである。この直観を最尤推定できちんとみてみる。

i回目に選んだ箱の番号をzi={A,B}で表すことにする。箱ziを選ぶ確率をp(zi)とし、p(xi|zi)を箱ziを選んだ時にxi=(r,b,y)個の玉の組み合わせが出てくる確率とする。この時の尤度p(X,Z)(X=x1,x2,..x6, Z=z1,z2,...z6)は
p(z1) p(x1|z1)
* p(z2) p(x2|z2)
...
* p(z6) p(x6|z6)
である。これが最大になるようなパラメータを求めたい。

しかし問題は、毎回どちらの箱を選んだか、つまりziが分からないことである。そこで、(少々強引だが)両方確率的に選んだと考える。つまり、p(z1) p(x1|z1)の代わりにp(A) p(x1|A) + p(B) p(x1|B) を尤度と考え、最大化する。p(x1, A) + p(x1, B) = P(x1)であり、P(X)を最大化していると考えられる。

{ p(A) p(x1|A) + p(B) p(x1|B) }
* { p(A) p(x2|A) + p(B) p(x2|B) }
...
* { p(A) p(x6|A) + p(B) p(x6|B) }

この確率を最大化するパラメータp()を求めたい。尤度の対数を計算することにすると次のL(X)を最大化するパラメータを求める最適化問題を解くことになる。

L(X) = ∑_{i} log ∑_{z} p(xi,z)

厄介な点は、対数の中に和の計算があることである(対数をとらなければいいと思えるが、乗算系の最適化は更に困難)。一般に、観測できない事象がある中で最尤推定したい場合には、このような対数の中に和が入る形になる。

で、ようやくEM。EMはこのような場合の最適化法である。オリジナルの場合のEMはQ関数が出てきてややこしいのだが、今回の多項分布のような場合は非常に簡単な計算になる(一般的には確率分布が指数型分布族の場合に簡単になる)

EMは次の方法で、L(X)が大きくなるようなパラメータを求める。ただし、L(X)が最大になるとは限らない。

---
1)パラメータの初期値を適当に設定する
2)次のE, Mステップを収束するまで交互に繰り返す。
Eステップ 今のパラメータでの、各潜在変数ziの期待値p(zi|xi)を求める
Mステップ この期待値を重みとして使って対数尤度(logP(X,Z))を最大化する
---

直感的にはEステップで今のパラメータを使った場合に、i番目のサンプルでziが使われた可能性はどれくらいかを求め、Mステップで先ほど求めたziを信じて尤度を最大化している。

EMはL(X)の下限を最大化しているのに相当する。

先ほどの例で試してみる。
最初に適当にパラメータを設定する。p(.|A)は箱Aの場合の赤、青、黄の割合のパラメータである。

p(A)=0.5 p(B)=0.5
p(.|A) = 0.2 0.4 0.4
p(.|B) = 0.4 0.4 0.2

はじめにEステップ。p(zi|xi)を求める

p(zi|xi) = p(xi|zi)p(zi)/P(xi)を求める。
p(xi)=p(xi|A)p(A)+p(xi|B)p(B)なので、

q(A,xi)=p(A)p(xi|A)
q(B,xi)=p(B)p(xi|B)
を求めれば

p(A|xi)=q(A,xi)/(q(A,xi)+q(B,xi))
p(B|xi)=q(B,xi)/(q(A,xi)+q(B,xi))

として計算できる

例えば、q(A,xi)=p(A)*p(xi|A)は
0.2 * (0.4^2 + 0.4^5 + 0.2^3)
である。実際の計算では掛け算を繰り返すと非常に小さな値になって問題が生じるので、対数で計算して、あとで戻してあげる。

次にMステップ。今求まったp(zi|xi)をつかってL(x)を最大化するような
パラメータp(.|A), p(.|B)を求める。今回の多項分布の場合は、単純にそれぞれの出現回数を
p(zi|xi)で重みづけした結果を出現回数だと思って、計算すればよい

例えば箱Aの赤の仮想的な出現回数 r'は
p(A|x1) * 2 +
p(A|x2) * 7 +
p(A|x3) * 2 +
p(A|x4) * 1 +
...

であり、これで求まったそれぞれの色の仮想的な出現回数r',b',y'をつかって
p(r|A) = r'/(r'+b'+y')
とすればよい。

これを繰り返すとp(z)とp(x|z)が推定できる。

実際には、EMは最大になるのを求められるとは限らないので様々な初期値から始める

---

| | Comments (345) | TrackBack (0)

« March 2009 | Main | May 2009 »