Gamma Process Non-negative Matrix Factorization (GaP-NMF) in Julia

最近 Julia で遊んでいて、その過程で非負値行列因子分解(NMF)のノンパラ版の一つであるGamma Process Non-negative Matrix Factorization (GaP-NMF) を書いてみました。(まぁmatlabコードの写経なんですが)

https://github.com/r9y9/BNMF.jl

元論文: Bayesian Nonparametric Matrix Factorization for Recorded Music by Matthew D. Hoffman et al. in ICML 2010.

デモ

http://nbviewer.ipython.org/github/r9y9/BNMF.jl/blob/master/notebook/GaP-NMF.ipynb

適当な音声(音楽じゃなくてごめんなさい)に対して、GaP-NMFをfittingしてみた結果のメモです。$K=100$ で始めて、100回ほどイテレーションを回すと適度な数(12くらい)にtruncateしているのがわかると思います。予めモデルの複雑度を指定しなくても、データから適当な数を自動決定してくれる、ノンパラベイズの良いところですね。

ハマったところ

  • GIGの期待値を求めるのに必要な第二種変形ベッセル関数は、exponentially scaled versionを使いましょう。じゃないとInf地獄を見ることになると思います(つらい)。Juliaで言うなら besselkx で、scipyで言うなら scipy.special.kve です。

雑感

  • MatlabのコードをJuliaに書き直すのは簡単。ところどころ作法が違うけど(例えば配列の要素へのアクセスはmatlabはA(i,j)でJuliaはA[i,j])、だいたい一緒
  • というかJuliaがMatlabに似すぎ?
  • Gamma分布に従う乱数は、Distributions,jl を使えばめっちゃ簡単に生成できた。素晴らしすぎる
  • 行列演算がシンプルにかけてホント楽。pythonでもmatlabでもそうだけど(Goだとこれができないんですよ…)
  • 第二種変形ベッセル関数とか、scipy.special にあるような特殊関数が標準である。素晴らしい。

Python版と速度比較

bp_nmf/code/gap_nmf と比較します。matlabはもってないので比較対象からはずします、ごめんなさい

Gistにベンチマークに使ったスクリプトと実行結果のメモを置いときました https://gist.github.com/r9y9/3d0c6a90dd155801c4c1

結果だけ書いておくと、あらゆる現実を(ry の音声にGaP-NMFをepochs=100でfittingするのにかかった時間は、

Julia: Mean elapsed time: 21.92968243 [sec]
Python: Mean elapsed time: 18.3550617 [sec]

という結果になりました。つまりJuliaのほうが1.2倍くらい遅かった(僕の実装が悪い可能性は十分ありますが)。どこがボトルネックになっているのか調べていないので、気が向いたら調べます。Juliaの方が速くなったらいいなー

おわりに

GaP-NMFの実装チャレンジは二回目でした。(たぶん)一昨年、年末に実家に帰るときに、何を思ったのか急に実装したくなって、電車の中で論文を読んで家に着くなり実装するというエクストリームわけわからんことをしていましたが、その時はNaN and Inf地獄に負けてしまいました。Pythonで書いていましたが、今見るとそのコードマジクソでした。

そして二回目である今回、最初はmatlabコードを見ずに自力で書いていたんですが、またもやInf地獄に合いもうだめだと思って、matlabコードを写経しました。あんま成長していないようです(つらい)

Julia歴二週間くらいですが、良い感じなので使い続けて見ようと思います。