2013-05-01

LDA CVB0 の C++ 実装

LDA

ちょっと C++ で実装したいものがあるのだけど、その前に練習ということで LDA の CVB0 推論を C++ で実装してみた。

https://github.com/shuyo/iir/tree/master/lda/ldacvb0_cpp

VC++2010 でしかビルドしてないが、一応後で gcc でもビルドしてみようかなとか思っているところ。
ビルドには @herumi 謹製 cybozulib が必要。

https://github.com/herumi/cybozulib

mmap ででかいファイルを効率的に読み込んだり、トピック単語分布の上位20語を抽出するのに便利なクラスが用意されていたり、とかいったあたりでお世話になっている。

cybozulib は std::string と互換かつ UTF-8 を透過的に扱える文字列クラスもあるんで、本来ならマルチバイト対応したかったんだけど入力テキストは ascii のみ。
今回の実装を始めたときはまだ正規表現まわりで互換性が不足していたので。が、その中で色々要望を出して対応してもらったので、後でマルチバイト対応するかも。

入力コーパスは Python 版実装と比較できるように nltk のものをそのまま読める、ようにしようとした。

http://nltk.googlecode.com/svn/trunk/nltk_data/index.xml

Brown corpus ならここからダウンロードして、展開、POS のアノテーションが付いているので -p オプションを付けて渡す。
読み込ませたいファイルをすべて指定しないといけないのだが、 cygwin ならこう書ける。

ldacvb0.exe -p brown/????

reuters corpus も同様。こちらは POS 付いてないので -p オプション不要。ただし、コマンドラインが長くなりすぎるので全てのファイルを渡すことはできない。
いろいろ手抜きですまん。

速度的には Python + numpy の20倍強といったところかな。
十分速いっちゃあ速いんだけど、ダブルバッファとか頑張ってしたのに、期待したほど速くはないというのが正直な感想。
うまく使えるアルゴリズムでの numpy は結構速いよってことなんだろうけど。

速度は20倍だけど、実装コストも20倍って感じ。
C/C++ で書くとなるとどうしても無駄なコピーをしてはならないという強迫観念に襲われてしまって、なんか素直な書き方ができなくなるんだよなあ。経験不足なんだろうけど。

2013-02-28

Labeled LDA (Ramage+ EMNLP2009) の perplexity 導出と Python 実装

LDA 機械学習

3年前に実装したものの github に転がして放ったらかしにしてた Labeled LDA (Ramage+ EMNLP2009) について、英語ブログの方に「試してみたいんだけど、どういうデータ食わせたらいいの？」という質問コメントが。
うーん、そうね、そういうところ書いてないから全然わかんないよね。手に入りやすいコーパスをだれでも食わせられるようにしてあると、やっぱり評価しやすいよね。
ということで、Labeled LDA に nltk のコーパスを食わせるスクリプトをさっくり書いてみた。

https://github.com/shuyo/iir/blob/master/lda/llda_nltk.py

同じディレクトリに llda.py がある状態で、そのまま起動すれば nltk.corpus.reuters から 100件のドキュメントをサンプリングして、適当なパラメータで Labeled LDA の推論して、ラベル-単語分布を出力する（ LDA と違ってトピック単語分布ではないところがミソ）。
Labeled LDA はマルチラベルに対応しており、nltk.corpus.reuters も各文書が複数カテゴリを持つので、データとしてはぴったり。ただし、カテゴリ毎の文書数が違いすぎるので、ちゃんと均したデータを作ったほうがいいかもしれない、と思いつつそこは手抜き。

自分の llda.py のコードを久しぶりに見たが、まあ約3年も前に書いたコードなものだから全般的に直したい衝動が湧き上がってくるのはおいとくとしても、perplexity のメソッドにただ一行 pass って書いてあるのはさすがにひどすぎる。この機会にこれくらいはさすがにやっつけておこう。
そこで Labeled LDA(Ramage+ EMNLP2009) を読みかえしてみたら、perplexity の式が論文内では導出されていなかった。
llda.py を実装したのは 2010年の 7月…… PRML 読書会ラス２が行われ、TokyoNLP の第１回が開催され、gihyo.jp の連載も始まったばかりの頃……。ははあ、さてはお前自力で perplexity の式が導出できなかったんだな（生暖かい目）。
というわけでさっくり導出＆実装を perplexity 対応に更新。

https://github.com/shuyo/iir/blob/master/lda/llda.py

perplexity 計算のために導出した式もさらしておこう。間違いあればご指摘歓迎。
いつもの LDA っぽくやるには、ドキュメント-トピック分布とトピック-単語分布を推定する必要がある。Labeled LDA (Ramage+ EMNLP2009) とは異なるが、それらを慣れた記号θとφで表すと、

$\theta_{dz}=\mathbb{E}[p(z|d)]=\begin{cases}\frac{n_{dz}+\alpha}{n_{d\cdot}+M_d\alpha}\hspace{20px}&\text{if }z\in\lambda^{(d)}\\0&\text{otherwize}\end{cases}$

$\phi_{zw}=\mathbb{E}[p(w|z)]=\frac{n_{zw}+\eta}{n_{z\cdot}+V\eta}$

λ^(d) というのは文書 d がもつラベルに対応するトピックの集合、M_d はその集合のサイズ。
Labeled LDA は要はカテゴリラベルごとにトピックが割り当てられて、各文書のトピックが観測済みのカテゴリラベルによる制約を受けるモデルなので、ひもづけられていないθ_dz が消える必要があり、ちょっとややこしくなっている。
この実装 llda.py では common カテゴリを勝手に増やして、すべての文書のラベルに追加している。こうすることで「多くの文書で同様な出現分布を持つ単語」は common カテゴリに入る。つまり stop words がキレイに common カテゴリに入る。ここは Labeled LDA の見どころ。

ちなみに Labeled LDA では１つのカテゴリには１つのトピックが対応付けられるわけだが、もちろんそれが最適であるとは限らない。複数のラベルでトピックを共有したほうがよい場合も十分考えられる（coffee カテゴリと tea カテゴリは共通のトピックがかなりあるだろう）。そこらへんも同時に推論するのが以前紹介した Dirichlet Process with Mixed Random Measures (Kim+ ICML2012) だったりする。

DP-MRM も実装する気まんまんだったのだが、いやもうここだけの話、HDP-LDA よりはるかにめんどくさい。自分の HDP-LDA 実装を参考に眺めているうちに、久しぶりにコード見たけどこれはひどいな、あー実装し直したい（どこかで聞いたな……）という衝動に負けて hdplda.py をゼロから実装しなおしてお腹いっぱいになってしまい。まあそのうちごにょごにょ。

2012-12-10

Online VB inference for HDP (Wang+ 2011) を実装してみたけど

ノンパラベイズ LDA

[Wang+ AISTATS 2011] Online Variational Inference for the Hierarchical Dirichlet. Process

HDP-LDA を Online な VB で解くという論文。
Teh さんらの VB 推論の論文も読んだんだけど、実にアクロバティックな分解が出てきて、いやこれはどうなんだろう……という気分になってしまい、実装してみる気にはならなかった。
一方、Wang さんらのこれはとても素直な式展開で、そんな簡単でいいんだーという感じ。自力で全部の式導出してみたけど、引っかかるところも特になし。

とはいえやっぱり VB なので、きつい独立の仮定がどれくらい精度に影響があるのか気になるところ。
アルゴリズムはとても簡単なので、これは試してみるっきゃあないでしょう、と実装してみた。
論文の pseudo code では u, v, λ を保持しておくように書かれていたが、それだけでは likelihood が計算出来ないので ζと $\varphi$ も持ちまわっている*1。

https://github.com/shuyo/iir/blob/master/lda/hdp_online.py

HDP-LDA の CGS 推論はこの前2回目の実装してみた *2んだけど、中核部分だけで300行くらいあって、とにかく間違いやすくて、めんどくさくて、遅くて。
一方 Online HDP の方は中核部分は100行未満、 dish や table が増えたり減ったりしないので、なんと楽なことか。推論も速い。
と、ここまではいいことだらけだけど、一番肝心の結果が定量的にも定性的にも全然ダメ。

そうすると、どこか間違えちゃったんだろうな、と当然思うところ。
幸い、Wang さんら自身が公開されている Python 実装があるので、これを読んでみたのだが、論文とは似ても似つかない実装になっていて、どういうこと状態。
そんなわけで実装がまるで似ていないのだが、一応動かしてみたところ、結果の「定量的にも定性的にも全然ダメ」というところは同じだった。

パラメータを最適なものを選べば良くなったりするんだろうかとか、もうちょっと粘るつもりだけど、なにか有益なアドバイスがもらえたり、みたいなことをちょっぴり期待してとりあえず実装を晒しておく。

ああ、そうそう。[Wang+ 11] では likelihood を求めるための具体的な式が書かれていないので、自分で導出した。
ここで間違えているという可能性もあるわけなので、晒しておこう。

トピック-単語分布 : $\bar{\phi}_{kw}\propto\lambda_{kw}$

文書-トピック分布 : $\bar{\pi}_{jk}=\frac{1}{N_j}\sum_t \varphi_{jtk}\sum_n \zeta_{jnt}$

あと、Wang さんらの実装が出力する *.topics ファイルから、単語-トピック分布を出力するスクリプトも貼り付けておく。
*.topics ファイルは論文の中で λ にあたるパラメータ(実装の中の m_lambda は、論文の中のλ-ηに対応しているのでややこしい)をべたっと出力しているので、上の式が正しければかんたんにトピック単語分布が得られる。

import sys, codecs, numpy

# topic-word distribution
with open(sys.argv[1], 'rb') as f:
    b = numpy.array([[float(x) for x in s.strip().split()] for s in f])
phi = b / b.sum(axis=1)[:,None]
K, W = phi.shape

with codecs.open('ap/vocab.txt', 'rb', 'utf-8') as f:
    voca = [s.strip() for s in f]
assert W == len(voca)

for k in xrange(K):
    print "\n-- topic: %d" % k
    for w in numpy.argsort(-phi[k])[:20]:
        print "%s: %f" % (voca[w], phi[k,w])

*1:後に出てくる Wang さんらの実装では、イテレーションの中でその時のパラメータで書く文書の likelihood を計算してそれらの変数の値は捨てていた

*2:この実装は Tokyo.ML 用のネタなので今のところはスルー

2012-08-23

階層ディリクレ過程を実装してみる (5) HDP-LDA の更新式を導出 (パープレキシティ)

LDA HDP 機械学習ノンパラベイズ

HDP-LDA の更新式を実装できたら、それが正しく動いているのか、そして収束したかどうかを確認するために perplexity を求めたいところだが、こちらも例によって論文には詳細な数式は書き下されていないので、最後にこれをやっつけよう。

$\displaystyle\exp\left(-\frac 1I\;\log p(w_1,\cdots,w_I|\text{Training corpus})\right)$

論文ではパープレキシティはこのような式で定義されている。
LDA の眷属は bag-of-words モデルなので、単語はすべて i.i.d. であることを使って展開していくと、

$p(w_1,\cdots,w_I|\text{Training corpus})=\prod_{w\in\text{Test data}}p(w|\boldsymbol{X})$
$=\prod_{w\in\text{Test data}}\sum_z p(z|\boldsymbol{X})p(x=w|z,\boldsymbol{X})$

ここで事後分布 p(z|X) と p(x|z,X) はそれぞれ次のようになる。

$p(z_{ji}|\boldsymbol{X})\sim\text{DP}\left(\alpha+n_j,\;\;\frac{\alpha}{\alpha+n_j}H+\frac{1}{\alpha+n_j}\sum_k n_{jk}\delta_{\phi_k}\right)$
$p(x|z=k,\boldsymbol{X})\sim\text{Dir}\left(\beta+n_{k1},\;\cdots,\;\beta+n_{kV}\right)$

ディリクレ分布の事後分布は PRML 等でもおなじみのものだが、ディリクレ過程の事後分布はなかなか見たこと無い人も多いかと。かの Ferguson の原典にはもちろんちゃんと書かれているんだけど、普通それ読まないし読めない。
で、なんかいいのないのかなあと探してたら、まさに Teh さんのサイトで Springer の Encyclopedia of Machine Learning の Dirichlet Process の項の原稿なのかな？　が公開されていた。

(Yee Whye Teh) Dirichlet Process

百科事典の項目なのでさすがに導出過程などは全部省かれて入るものの、ディリクレ過程の性質が CRF やら stick breaking やら、そして事後分布までひと通り網羅されていて、もちろん今風の記法で書かれているので、とりあえず DP を扱う上で知っておかないといけないことを抑えておきたいという人にはおすすめかと。

話を元に戻して。
これらの事後分布の分布を perplexity に放り込むには厳密には更に周辺化とかうんぬんしなくといけなくて、そんな計算出来そうにない積分しだしたら頭が痛くなってくるので、それぞれの期待値で代用していいことにしよう(つまり点推定したものを使うわけだが、だいたいこれは p(z|X) と p(w|z,X) が独立であるという仮定を入れることに相当する）。*1

そこで期待値 θ_jk, φ_kv を求めるとそれぞれ

$\theta_{jk}=\mathbb{E}[p(z_{ji}=k|X)]=\begin{cases}\frac{1}{n_j+\alpha}\left(n_{jk}+\alpha\cdot\frac{m_k}{m_{\cdot}+\gamma}\right)\hspace{20px}&\text{if }k\text{ is used}\\ \frac{1}{n_j+\alpha}\left(\alpha\cdot\frac{\gamma}{m_{\cdot}+\gamma}\right)&\text{if }k\text{ is new,}\\\end{cases}$

$\phi_{kv}=\mathbb{E}[p(x_{ji}=v|z_{ji}=k,\;X)]=\begin{cases}\frac{n_{kv}+\beta}{n_k+V\beta}\hspace{20px}&\text{if }k\text{ is used}\\\frac 1V&\text{if }k\text{ is new}\\\end{cases}$

となる。
あとはこのθとφを使って $\prod_{x_{ji}\in\text{Test data}}\sum_k\theta_{jk}\phi_{kx_{ji}}$ を求めれば良い。

ちなみに本来なら perplexity を求めるのに用いるテストデータは訓練データと関係ない未知のドキュメントであることが望ましいわけだが、p(z|X) を未知のドキュメントに対して求めるのは骨が折れる。
そこでコーパスの各ドキュメントの単語の 90% を訓練に使って、推定されたθ_jkと残りの 10% の単語で perplexity を求めるということがよくやられるようだ。

今回の記事まとめ、そして DP-MRM の実装のために１年前の HDP-LDA の自分の実装を見なおしていたわけだけど、当時は全く理解せずにかろうじてわかったつもりの部分から実装しているから、いやあ本当にひどいコードでよく動いているなと自分でも関心（苦笑）。
今だったらだいぶすっきり書けると思うので、コードの見通し優先で実装し直したいところ。まあそうやって実装しなおしても、１年後くらいに「うわ、これはひどい……」とか言ってそうな気もするけどｗ

*1:という流れだと解釈しているんだけど、あってるかな？

2012-08-20

DP-MRM (Kim+ ICML12) の数式を確認中

ICML LDA HDP ノンパラベイズ機械学習

Dirichlet Process with Mixed Random Measures (Kim+ ICML2012, 以下 DP-MRM) を実装してみようかと論文の計算式をひと通り確かめているのだが、いくつか気になることが。

k, l のサンプリング

論文の (9), (10) 式にあげられている k, l の全条件付き事後分布。

$\displaystyle p(k_{jt}=k, l_{jt}=l|\boldsymbol{k}_{-jt}, \boldsymbol{l}_{-jt},\text{rest})\propto\frac{m_{jk}+\eta}{m_{j}+K\eta}\times\frac{m_{kl}}{m_k+\gamma_k}f_{kl}(x_{ji})$

末尾の x_ji は x_jt の間違い、というのは前回の記事ですでに書いたとおりだが、さらに r_j を考慮にいれないといけないはず。
導出は HDP の時の導出を参考にしてもらうと簡単だが、ちょっと手抜き気味に書くと、

$\displaystyle \begin{eqnarray} p(k_{jt}=k, l_{jt}=l|\boldsymbol{k}_{-jt}, \boldsymbol{l}_{-jt},\text{rest}) &=&\frac{p(\boldsymbol{x}_{jt}|k_{jt}=k,l_{jt}=l,\cdots)p(k_{jt}=k,l_{jt}=l|\cdots)}{p(\boldsymbol{x}_{jt}|\cdots)} \\ &\propto& f_{kl}(\boldsymbol{x}_{jt})\cdot p(l_{jt}=l|k_{jt}=k,\cdots)\cdot p(k_{jt}=k|\cdots) \end{eqnarray}$

ここで $p(l_{jt}=l|k_{jt}=k,\cdots)$ の扱いは論文の通りで OK だが、 $p(k_{jt}=k|\cdots)$ はr_j、つまり label(j) の影響を受けることが触れられていない。
したがって

$p(k_{jt}=k|\boldsymbol{k}_{-jt}, \boldsymbol{l}_{-jt}, \boldsymbol{t})\propto\begin{cases}m_{jk}^{-jt}+\eta\hspace{10px}&\text{if }r_{jk}=1,\\ 0&\text{if }r_{jk}=0\end{cases}$

が正しいと思われる。

perplexity

論文の 4.3 章で perplexity を求めている。

$p(\boldsymbol{x}^{\text{test}})=\prod_{ji\in\boldsymbol{x}^{\text{test}}}\frac 1S\sum_{s=1}^S\sum_{k=1}^K\sum_l \theta_{jkl}^{(s)}\psi_{jkl}^{(s)},$
$\displaystyle \theta_{jkl}^{(s)}=\frac{n_{jkl}+\alpha\{m_{kl}/(m_{k}+\gamma_k)\}}{n_{jk}+\alpha}$

ψの方は問題ないので省略。
当初、ここで求めているθは事後分布 $p(\theta_{j\cdot}=\phi_l^k|\boldsymbol{t},\boldsymbol{k},\boldsymbol{l})$ (の期待値)かと思ったのだが（もちろん文脈的にもそうあるべき）、右辺を見ると l で summation を取ったらほぼ 1 になる(未知の l を選ぶ確率の分小さい)ので、あれー違うのかなあとしばらく悩んでしまった。でも perplexity の式では k でも周辺化しているし……。
というわけでこのθの式は正しくはこうなるはず。

$\displaystyle \theta_{jkl}^{(s)}=\frac{n_{jkl}+\alpha\{m_{kl}/(m_{k}+\gamma_k)\}}{n_{j}+|\text{label}(j)|\alpha}$

$|\text{label}(j)|$ は document j に振られたラベルの個数。 $\sum_k r_jk$ と書いてもいいかも。
論文では、DP-MRM の perplexity は Labeled LDA との比較にのみ用い、またその際のデータセットはシングルラベルに限定しているようであり、そのときには論文の式とこの新しい式は一致するので、結果自体には問題はない気はする。
が、やはりθの式はマルチラベルのことも考慮したものであるほうがいいだろう。

とっとと英語で書けって？　そーねー。ま、実装済んだら。

【以下 2012/8/22 追記】
上に書いたのは間違い。まじめに DP の事後分布を評価したら、こうなった。

$\displaystyle \theta_{jkl}^{(s)}=\frac{1}{n_{j}+\alpha}\left(n_{jkl}+\alpha\cdot\frac{m_{jk}+\eta}{(m_{j}+|\text{label}(j)|\eta)}\cdot\frac{m_{kl}}{m_{k}+\gamma_k}\right)\hspace{20px} \text{if}\; r_{jk}=1$

そうだよな、eta が出てこないなんて変だと思ったんだ。
ちなみにシングルラベルの場合、j に対して k が1つしか対応しないので、その k に対して n_jk=n_j, m_jk=m_j, |label(j)|=1 となり、[Kim+ 2012] のθと上の式は一致する。

【以下 2012/8/23 追記】
都合良く勘違いしていたが、よくよく確認したら、[Kim+ 2012] で perplexity を求めているのはマルチラベルの dataset に対してだった。
そうすると、論文に書いてあるθの式では確率が不当に高くなり、DP-MRM の perplexity として本来よりも良い値が出てしまう懸念があるのだが、大丈夫なのだろうか。

あと実装してて気づいたが、l についても周辺化するのだから、l=l_new に対するθ, ψも必要だ。
これをもし忘れてたら、今度は perplexity が本来の値より悪くなる。まあ上の「間違いθ」よりは影響はだいぶ小さい気はするが。

$\displaystyle \theta_{jkl_{\text{new}}}=\frac{\alpha\gamma_k}{(n_{jk}+\alpha)(m_k+\gamma_k)}$
$\displaystyle \psi_{kl_{\text{new}}v}=\frac 1W$

しっかし「実装してみて、やっと気づく」ことが多すぎる……。全部が全部実装なんてしてられないので、もうちょっと実装する前に気づければいいんだけど（苦笑）。

2012-08-15

DP-MRM (Kim+ ICML12) の更新式も導出

ICML LDA HDP ノンパラベイズ機械学習

昨日の HDP-LDA の更新式導出のついでに、ICML 読み会で紹介した教師ありノンパラ LDA である (Kim+ ICML2012) Dirichlet Process with Mixed Random Measures (DP-MRM) の更新式も導出しておこう。
DP-MRM のモデルについては ICML 読み会での発表資料参考。

[Kim+ ICML2012] Dirichlet Process with Mixed Random Measures : A Nonparametric Topic Model for Labeled Data from Shuyo Nakatani

www.slideshare.net

Teh さんの HDP の論文では、実装に必要になる f_k の一般式までしか書かれていなくて、HDP-LDA の場合の書き下しが行われていなかったのだが、DP-MRM でも同じように f_kl という式に帰着して、でもその f_kl は具体的には書き下されていないというスタイルになっている。いやそんなところは真似しなくても。
なので、その f_kl を展開してみよう。

$\displaystyle f_{kl}(x_{ji})=\frac{\int f(x_{ji}|\phi_l^k)\prod_{x_{j'i'}\in\boldsymbol{x}_{kl}}f(x_{j'i'}|\phi_l^k)h(\phi_l^k)d\phi_l^k}{\int\prod_{x_{j'i'}\in\boldsymbol{x}_{kl}}f(x_{j'i'}|\phi_l^k)h(\phi_l^k)d\phi_l^k}$
where $\boldsymbol{x}_{kl}=\{x_{ji}|k_{jt_{ji}}=k,\;l_{jt_{ji}}=l\}$

こうやって天下りで与えられると意味不明だが、、例によって $p(x_{ji}|\boldsymbol{t}, \boldsymbol{k}, \boldsymbol{l}, \boldsymbol{x}^{-ji})$ を f_kl と書いているだけである。

DP-MRM でも HDP-LDA と同じく f(x|φ) は単語の多項分布 Multi(φ)、h(φ) は単語分布を引くディリクレ分布 Dir(β) という設定なので、これを上の式に放り込む。
v=x_ji, t=t_ji, k=k_jt, l=l_jt ということにしつつ、 $n_lkv = \sharp\{x_ji|v=x_ji, t=t_ji, k=k_jt, l=l_jt\}$ 、-ji が付いた場合はそこから x_ji を除いたもの、といういつもの記法も使って計算していく。
HDP-LDA とほとんど同じ流れなので、細かいところは省略。

$\displaystyle f_{kl}(x_{ji})=\frac{\int\phi_{lkv}\prod_{w}\phi_{lkw}^{n_{lkw}^{-ji}}\prod_{w}\phi_{lkw}^{\beta-1}d\phi_l^k}{\int\prod_{w}\phi_{lkw}^{n_{lkw}^{-ji}}\prod_{w}\phi_{lkw}^{\beta-1}d\phi_l^k}$
$\displaystyle =\frac{\Gamma(\beta+n_{lkv}^{-ji}+1)\prod_{w\neq v}\Gamma(\beta+n_{lkw}^{-ji})}{\Gamma(V\beta+n_{lk}^{-ji}+1)}\cdot\frac{\Gamma(V\beta+n_{lk}^{-ji})}{\prod_w\Gamma(\beta+n_{lkw}^{-ji})}$
$\displaystyle =\frac{\beta+n_{lkv}^{-ji}}{V\beta+n_{lk}^{-ji}}$

$\displaystyle f_{kl_{\text{new}}}(x_{ji})=\int f(x_{ji}|\phi)h(\phi)d\phi=\frac 1V$

$\displaystyle f_{kl}(\boldsymbol{x}_{jt})=\frac{\prod_w\Gamma(\beta+n_{lkw}^{-jt}+n_{jtw})}{\Gamma(V\beta+n_{lk\cdot}^{-jt}+n_{jt})} \cdot \frac{\Gamma(V\beta+n_{lk\cdot}^{-jt})}{\prod_w\Gamma(\beta+n_{lkw}^{-jt})}$ 【2012/08/22 修正】

$\displaystyle f_{kl_{\text{new}}}(\boldsymbol{x}_{jt})=\frac{\Gamma(V\beta)\prod_v\Gamma(\beta+n_{\cdot\cdot v}^{jt})}{\Gamma(V\beta+n_{\cdot\cdot\cdot}^{jt})\prod_v\Gamma(\beta)}$

ああ、そう言えば、[Kim+ ICML2012] の式 (9), (10) に出てくる $f_{kl}(x_{ji})$ は $f_{kl}(x_{jt})$ の間違い。

2012-08-14

階層ディリクレ過程を実装してみる (4) k の全条件付き分布と HDP-LDA の更新式を導出

LDA HDP 機械学習ノンパラベイズ

来週、サイボウズ・ラボユース合宿なんてのがあったりする(昨年の様子ちら見)。
ちょっと缶詰っぽい感じでコードをガジガジ書く機会になるので、この前の ICML 読み会で紹介した DP-MRM でも実装してみるかー、でも HDP-LDA の実装の細かいところ忘れちゃった（てへ）、TokyoML で発表する約束もしてるし、ちょっと計算やり直してみよっかな、という感じで Teh さんの HDP の論文をまたぞろひっぱりだして読んでたりする。
しかし、初めて読んだときは全くわからなくて、計算も何度も間違えて、なんかもうひどく苦労した印象しかなかったのだが、今読みかえすと、あれれこんなに簡単だったっけ……？　ディリクレ分布のお化けみたいな計算を10回くらいこなしたおかげか、なんか迷わずさくさく計算できちゃってあの苦労はなんだったんだ状態。人間の適応力って大したもんですね。

で、計算をまとめた自分のブログを見返したら、t のサンプリングまで導出しておいて、「 k のサンプリングはまた次回」で終わっている。これはひどい。

というわけで、１年前から放置されていた記事の続き「 k のサンプリング式 (34) の導出」をようやく書こう。

$p(k_{jt}=k|\boldsymbol{t}, \boldsymbol{k}^{-jt}) \propto \begin{cases}m_{\cdot k}^{-jt}f_k^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})\hspace{10px}&\text{if } k \text{ previously used,}\\\gamma f_{k^{\text{new}}}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})&\text{if }k=k^{\text{new}}\end{cases}$ 　　(34)

例によって省略されている x を左辺の given に補って、k_jt に絡んでいる {x_ji} たちを軸に展開していくことになる。
そこで $\boldsymbol{x}_{jt}=\{x_{ji}|t_{ji}=t\}$ , $-\boldsymbol{x}_{jt}=\{x_{ji}|t_{ji}\neq t\}$ と書くことにすると、

$p(k_{jt}=k|\boldsymbol{t}, \boldsymbol{k}^{-jt},\boldsymbol{x})=\frac{p(\boldsymbol{x}_{jt}|k_{jt}=k,\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt})p(k_{jt}=k|\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt})}{p(\boldsymbol{x}_{jt}|\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt})}$

分子第1項は、t のサンプリングのところで使った $f_k^{-x_{ji}}(x_{ji})=p(x_{ji}|\boldsymbol{t},\boldsymbol{k},\boldsymbol{x}^{-ji})$ の x_ji がちょうど $\boldsymbol{x}_{jt}$ に置き換わったものなので、 $f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})$ と書くことにする。つまり t_ji のサンプリングはそれと結びつく x_ji 一つで済むが、k_jt のサンプリングはテーブル t に座っているすべての x_ji たちが関わってきてしまう、ということ。
実装にあたってはこの $f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})$ の評価が一番めんどくさい。

分子第2項は CRF なり Polya なりの式に出ている通りである。

$p(k_{jt}=k|\boldsymbol{t},\boldsymbol{k}^{-jt},-\boldsymbol{x}_{jt}) \propto \begin{cases}m_{\cdot k}^{-jt}\hspace{10px}&\text{if } k \text{ previously used,}\\\gamma&\text{if }k=k^{\text{new}}\end{cases}$

分母は正規化項に押し込んでしまえば、(34) 式が得られる。

ここまで base measure H の密度関数を h(φ)、 emission F(φ) の密度関数を f として計算してきた。
つまり HDP は任意の分布 H, F(φ) に対してモデリング可能な一般的な枠組みだったわけだ（実際には $f_k^{-x_{ji}}(x_{ji})$ の積分が計算可能であるために H は F(φ) の共役事前分布であることが望ましい)。
この HDP で特に base measure H を単語分布 H=Dir(β), F(φ) をトピック-単語間の多項分布 F(φ)=Multi(φ) と置いた場合が HDP-LDA である。
というわけでいよいよ HDP-LDA の実装に必要な更新式を出してみる。

$f_{k}^{-x_{ji}}(x_{ji})$ にそれぞれの密度関数を放り込もう。

$f(x_{ji}=v|\phi_k)=\phi_{kv}$
$h(\phi)=\frac 1Z \prod_v \phi_v^{\beta-1}$ , where $Z=\frac{\prod_v \Gamma(\beta)}{\Gamma(\sum_v \beta)}$

そこで v=x_ji, $n_{kw}^{-ji}$ を「トピック k, 単語 w をもつ x_mn たちの個数(ただし x_ji を除く)」とおくと、

$\displaystyle\begin{eqnarray} f_k^{-x_{ji}}(x_{ji}) &=& \frac{\int f(x_{ji}|\phi_k)\prod_{j'i'\neq ji,z_{j'i'}=k}f(x_{j'i'}|\phi_k)h(\phi_k)d\phi_k}{\int \prod_{j'i'\neq ji,z_{j'i'}=k}f(x_{j'i'}|\phi_k)h(\phi_k)d\phi_k} \\ &=&\frac{\int\phi_{kv}\prod_{w}\phi_{kw}^{n_{kw}^{-ji}}\prod_{w}\phi_{kw}^{\beta-1}d\phi_{k}}{\int\prod_{w}\phi_{kw}^{n_{kw}^{-ji}}\prod_{w}\phi_{kw}^{\beta-1}d\phi_{k}} \\ &=&\frac{\int\prod_{w}\phi_{kw}^{I(w=v)+n_{kw}^{-ji}+\beta-1}d\phi_{k}}{\int\prod_{w}\phi_{kw}^{n_{kw}^{-ji}+\beta-1}d\phi_{k}} \end{eqnarray}$

ここで多項ベータ積分 $\int\prod_i t_i^{\alpha_i-1}d\boldsymbol{t}=\frac{\prod_i\Gamma(\alpha_i)}{\Gamma(\sum_i\alpha_i)}$ (ただし t は t_i>0, Σt_i<1 なる単体上を動く) を思い出すと、

$\displaystyle=\frac{\Gamma(\beta+n_{kv}^{-ji}+1)\prod_{w\neq v}\Gamma(\beta+n_{kw}^{-ji})}{\Gamma(V\beta+n_{k\cdot}^{-ji}+1)}\cdot\frac{\Gamma(V\beta+n_{k\cdot}^{-ji})}{\prod_{w}\Gamma(\beta+n_{kw}^{-ji})}$

w=v 以外はキャンセルして、

$\displaystyle=\frac{\Gamma(\beta+n_{kv}^{-ji}+1)\Gamma(V\beta+n_{k\cdot}^{-ji})}{\Gamma(\beta+n_{kv}^{-ji})\Gamma(V\beta+n_{k\cdot}^{-ji}+1)}$

よってガンマ関数の特徴 $\Gamma(n+1)=n\Gamma(n)$ より

$\displaystyle f_k^{-x_{ji}}(x_{ji}) = \frac{\beta+n_{kv}^{-ji}}{V\beta+n_{k\cdot}^{-ji}}$ (ただし v=x_ji)

となる。あら簡単。
同様に $f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})$ についても計算することができる。

$n_{kw}=\#\{x_{mn}|k_{jt_{mn}}=k, x_{mn}=w \}$ 　(トピック k 単語 w である x_mn の個数)
$n_{kw}^{-jt}=\#\{x_{mn}|k_{jt_{mn}}=k, x_{mn}=w, (m,t_{mn})\neq(j,t) \}$ 　(トピック k 単語 w である x_mn で、テーブル (j,t) に座っていない個数)

とおくと、

$\displaystyle f_{k}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt}) = \frac{\prod_w\Gamma(\beta+n_{kw}^{-jt}+n_{jtw})}{\Gamma(V\beta+n_{k\cdot}^{-jt}+n_{jt})} \cdot \frac{\Gamma(V\beta+n_{k\cdot}^{-jt})}{\prod_w\Gamma(\beta+n_{kw}^{-jt})}$
$\displaystyle = \frac{\prod_w(\beta+n_{kw}^{-jt}+n_{jtw}-1)\cdots(\beta+n_{kw}^{-jt})}{(V\beta+n_{k\cdot}^{-jt}+n_{jt}-1)\cdots(V\beta+n_{k\cdot}^{-jt})}$
【2012/08/22 修正】

この２つを (32), (33), (34) 式に放り込めば、HDP-LDA の実装に必要な数式が揃う。

最後頑張ってガンマ関数を消しているが、実装で手抜きをしたいならガンマ関数を展開せず log-gamma 関数を使うのが楽かもしれない(対数を取らないと、浮動小数点の範囲をさくっとオーバーフローする)。
この HDP-LDA の実装では、ガンマ関数を展開した方の式で結果をキャッシュして使いまわしている(cur_log_base 〜 log_f_k_new_x_jt メソッドのあたり)が、Python 実装なら log-gamma 関数をばしっと呼んだほうが速かったかもしれない……。

【追記】
あ、f_{k^new} 忘れてた。

$f_{k^{\text{new}}}^{-x_{ji}}(x_{ji})=\int p(x_{ji}|\phi)p(\phi)d\boldsymbol{\phi}=\int \phi_v\frac{\Gamma(V\beta)}{\prod_w\Gamma(\beta)}\prod_w \phi_w^{\beta-1}d\boldsymbol{\phi}$
$=\frac{\Gamma(V\beta)}{\prod_w\Gamma(\beta)}\frac{\Gamma(\beta+1)\prod_{w\neq v}\Gamma(\beta)}{\Gamma(V\beta+1)}=\frac{1}{V}$

同様に

$f_{k^{\text{new}}}^{-\boldsymbol{x}_{jt}}(\boldsymbol{x}_{jt})=\frac{\Gamma(V\beta)\prod_w\Gamma(\beta+n_{jtw})}{\Gamma(V\beta+n_{jt})\Gamma(\beta)^V}$ 【2012/9/26 修正】
【／追記】

階層ディリクレ過程を実装してみる (5) HDP-LDA の更新式を導出 (パープレキシティ) - 木曜不足