LDA 実装の比較

(1年前に公開し忘れていた記事を下書き消化のために投稿しておきます)

LDA の実装を比較のためにいくつか触ってみました。だいたい以下のような選択肢かなというところですが、メモ的に書いた各言語の LDA 実行までの手順などをまとめておきます。

	利用言語	利用形態	負荷拡張性	アルゴリズム
Spark MLlib 1.6.0	Scala,Java, Python,R	ライブラリ	並列/分散	変分ベイズ, EM
gensim	Python	ライブラリ	並列	変分ベイズ
GibbsLDA++	シェル	コマンド	-	Gibbs Sampling
R	R	ライブラリ	-	Gibbs Sampling

実行コスト

それぞれ実行させて時間を計測していますが、下記に示すようにそれぞれ実行条件がかなり違っています。誤解のないように読み取ってください。

• 実行環境 (物理/仮想/Linux/Windows)
• 実行条件 (ファイルから読み込むか既にメモリ上にあるのか)
• 対象データの前処理具合 (BoW の生成からかなど)

形態素解析を行い不要な形態素を除去した $13,033$ 文書 $20,780$ 単語に対して、$k=100$, $iter=50$ の条件で実行。

Spark MLlib

アルゴリズム	Nodes	Cores	実行時間
EM	5	80	224.092
EM	1	8	81.854
EM	1	1	112.606
変分ベイズ	1	8	220.147
変分ベイズ	1	1	310.367

実行時間には RDD がローカルファイルを読み込んで BoW を作成する処理も含まれているため他より時間がかかっています。分散実行した結果が最も遅くなっているのは LDA の計算コストそのものよりもシャッフルやノード間のデータ転送コストが大きいためでしょうね。

gensim

アルゴリズム	Nodes	Cores	実行時間
変分ベイズ	1	4	15.396
変分ベイズ	1	1	20.576

Python の gensim は変分ベイズのみのようですが workers を増やして並列処理が可能で速度も速くパフォーマンス的なバランスが良さそうです。

GibbsLDA++

アルゴリズム	Nodes	Cores	実行時間
Gibbs Sampling	1	1	58.993

シングルスレッドに加えてローカルファイルの入出力処理が入りながら $k=10$ 付近では gensim よりも高速でしたが、$k$ が大きくなると遅くなるのでしょうか。スケールは厳しいかも知れませんがちょっとしたアドホックな実行向きですかね。

R言語

アルゴリズム	Nodes	Cores	実行時間
Collapsed Gibbs Sampling	1	1	24.247

計算自体は意外にも GibbsLDA++ より速かった R の lda.collapsed.gibbs.sampler {lda} (ただしローカルファイルを読み込んでライブラリに渡せる形式に変換する処理は遅かったのでトータルで見たら遅い)。R は実用速度よりも lda {lda} や lda.cvb0 {lda} など実装やアルゴリズムが多く試せる部分に期待すると良い印象です。

それぞれの実装を使ってみる

gensim

# -*- coding: utf-8 -*-
# Python gensim による LDA のサンプル
from gensim import corpora, models
import time

# 文書データの読み込み (形態素解析/名詞抽出済みの1行1文書で単語が空白で区切られたテキストファイル)
texts = [ ]
for line in open('docs.gibbs', 'r'):
    texts.append(line.split())

# 辞書の作成 (id:単語:出現数)
dictionary = corpora.Dictionary(texts)
dictionary.save_as_text('./docs.dic')
# dictionary = corpora.Dictionary.load_from_text("./docs.dic")

# コーパス作成
corpus = [dictionary.doc2bow(text) for text in texts]
corpora.MmCorpus.serialize('corpus.mm', corpus)
# corpus = corpora.MmCorpus('corpus.mm')

# LDA の計算
t0 = int(time.time() * 1000)
lda = models.ldamodel.LdaModel(corpus, num_topics=100, iterations=50)
t1 = int(time.time() * 1000)
print t1 - t0

# LDA の計算 (マルチプロセッサ対応版)
t0 = int(time.time() * 1000)
lda = models.ldamulticore.LdaMulticore(corpus, num_topics=10, iterations=50, workers=4)
t1 = int(time.time() * 1000)
print t1 - t0

GibbsLDA++

GibbsLDA++ 及びその Java 移植版の JGibbsLDA はコマンドライン型の LDA 実装です。Gibbs Sampling を実装していますが並列処理 (マルチコア対応) や分散処理には対応していないようです。

ビルド

GibbsLDA++: A C/C++ Gibbs Sampling LDA から GibbsLDA++-0.2.tar.gz をダウンロードして解凍しビルドを行います。

$ tar zxvf GibbsLDA++-0.2.tar.gz
$ cd GibbsLDA++-0.2/
$ make all

g++ が利用可能であれば make all でビルドが完了し src/lda コマンドができあがります。

error: ‘atof’ was not declared in this scope, error: ‘printf’ was not declared in this scope が発生する場合はここ を参考に 2 ファイルに #include を追加。

$ vim src/util.cpp
...
#include <stdio.h>
#include <stdlib.h>  // 追加
#include <string>
...
$ vim src/lda.cpp
...
*/

#include <stdio.h>  // 追加
#include "model.h"
...
$ make clean; make all
make -C src/ -f Makefile clean
make[1]: Entering directory `/export/home/t_takami/git/gibbslda++/GibbsLDA++-0.2/src'
...
make[1]: Leaving directory `/export/home/t_takami/git/gibbslda++/GibbsLDA++-0.2/src'

$ src/lda --help
Please specify the task you would like to perform (-est/-estc/-inf)!
Command line usage:
        lda -est -alpha <double> -beta <double> -ntopics <int> -niters <int> -savestep <int> -twords <int> -dfile <string>
        lda -estc -dir <string> -model <string> -niters <int> -savestep <int> -twords <int>
        lda -inf -dir <string> -model <string> -niters <int> -twords <int> -dfile <string>

データの準備

先頭行に文書数、以降「1行=1文書」となるように、文書に含まれる単語を空白区切りで羅列したファイルを作成します。なお空行 (特徴を表す単語数が 0 の文書) が存在すると Invalid (empty) document! のエラーが発生します。

$ head -n 5 docs.gibbs
13033
ストレージ お願い 選択 カード カード 場合 場合 回答 方法 方法 可能 内部 機種 容量 録画 録画 録画 録画 設定 時間 保存 保存 保存 保存
世間 世間 世間 大事 気 ズレ 常識 常識 常識 常識 常識 常識 自分 自分 自分
キャラクター 服 操作 ゲーム 金髪 キャラ 機種 銀髪
最高 最高 妄想 デート デート
電話 電話 電話 別れ 本当 プレゼント プレゼント 好き 好き お願い お願い 涙 涙 涙 気分 成田 会話 陰性 目 関係 関係 部屋 恋愛 ウルグアイ ウルグアイ ウルグアイ ウル グアイ フクロウ 最終 感情 裏側   …
…

実行結果

コマンドライン上でトピック数 $k$ や繰り返し数を指定して実行します。

$ src/lda -est -niters 50 -ntopics 10 -twords 5 -dfile docs.gibbs
Sampling 50 iterations!
Iteration 1 ...
...
Iteration 50 ...
Gibbs sampling completed!
Saving the final model!

実行が終わるといくつかのファイルができあがっているはずです。

model-final.others はこのモデルが作成された時の実行パラメータ、wordmap.txt には単語に振られた ID が格納されています。

$ cat model-final.others
alpha=5.000000
beta=0.100000
ntopics=10
ndocs=13033
nwords=20779
liter=50

$ cat wordmap.txt
20779
Tシャツ 4601
あい 1829
あいこ 19897
あいさつ 2125
...

model-final.twords はトピックごとに最も特徴を表している単語を -twords オプションで指定した個数羅列します (-twords オプションを省略したするとこのファイルは作成されません)。

$ cat model-final.twords
Topic 0th:
        人   0.047465
        感じ   0.019363
        気   0.018178
        他   0.016968
        普通   0.016004
Topic 1th:
        仕事   0.043820
        時間   0.024824
        家   0.019440
        今   0.017962
        母   0.016881
Topic 2th:
        場合   0.033522
        月   0.018820
        会社   0.018083
        保険   0.015252
        お願い   0.012468
…

model-final.phi はトピックごとの単語出現確率。1行=1トピック に対応し 1 行内のカラム数はコーパス上の単語数と一致します。つまり「$k$行(全単語数)列」の行列になります。model-final.twords の元となる数値データですね。

$ head -c 100 model-final.phi
0.000002 0.000002 0.000002 0.000002 0.000799 0.000002 0.000002 0.002415 0.000002 0.000002 0.000002 0

model-final.theta は文書ごとにどのトピックに属するかを表す確率。1行=1文書 に対応し 1カラム=1トピックに相当します。つまり「(文書数)行$k$列」の行列データです。

$ head -n 5 model-final.theta
0.081081 0.216216 0.067568 0.067568 0.162162 0.081081 0.067568 0.108108 0.081081 0.067568
0.076923 0.076923 0.123077 0.092308 0.076923 0.076923 0.107692 0.076923 0.076923 0.215385
0.086207 0.103448 0.103448 0.086207 0.137931 0.086207 0.103448 0.086207 0.086207 0.120690
0.090909 0.090909 0.109091 0.127273 0.090909 0.090909 0.090909 0.090909 0.090909 0.127273
0.035971 0.028777 0.111511 0.323741 0.050360 0.086331 0.248201 0.039568 0.028777 0.046763

model-final.tassign は元データの 1行=1文書 ごとの [単語ID:トピック] 形式の羅列。

$ head -n 5 model-final.tassign
0:1 1:4 2:1 3:7 3:7 4:8 4:7 5:5 6:1 6:1 7:4 8:1 9:4 10:1 11:4 11:4 11:4 11:4 12:1 13:0 14:1 14:1 14:1 14:1
15:9 15:9 15:9 16:3 17:6 18:6 19:2 19:9 19:2 19:2 19:9 19:9 20:9 20:9 20:9
21:2 22:9 23:1 24:4 25:6 26:9 9:4 27:4
28:9 28:2 29:9 30:3 30:3
31:4 31:6 31:6 32:3 33:2 34:3 34:3 35:3 35:3 1:5 1:2 36:3 36:3 36:3 37:6 38:2 39:6 40:6 41:6 42:3 42:6 43:6 44:3 …

GibbsLDA++ はほぼ数値モデルを作成するところまでですので twords 以外の結果は自分で文書や単語とマッピングする必要があります。

R 言語

個人的に言語機能が取っつきにくくて敬遠気味の R 言語ですが lda.collapsed.gibbs.sampler {lda} を見ると Collapsed Gibbs Sampling のいくつかのモデルが実装されています。

並列処理 (マルチコア対応) や分散処理には対応していないようです。

インストール

まず CentOS 7.2 へ R 言語をインストールします (以下の rpm コマンドで 404 となる場合は URL のファイル部分を削除して該当するバージョンを探してください)。

$ sudo rpm -ihv http://ftp.riken.jp/Linux/fedora/epel/7/x86_64/e/epel-release-7-5.noarch.rpm
$ sudo yum install R

R のインストールが完了したら対話インターフェースを起動して LDA パッケージをインストールします。途中で LDA パッケージを /usr/lib64/R/library にインストールしようとしますが root ユーザではないので個人ライブラリを指定しています。またダウンロード元を聞かれますので Japan (Tokyo) CRAN mirror を選択します。

以下ではデモを実行するために reshape2 と ggplot2 を追加でインストールしています。

$ R
R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-redhat-linux-gnu (64-bit)
...
> install.packages("lda")
...
Would you like to use a personal library instead?  (y/n) y
...
Selection: 13
...
> install.packages("reshape2")
> install.packages("ggplot2")
> require("lda")
> demo(lda)

まだ ggplot2 の仕様変更に追従していないようで stat_count() must not be used with a y aesthetic. が発生してグラフは表示されません。qplot() の部分を正しく手動で起動すれば上手く行くと思いますが、LDA の計算処理自体は終わっていますので無視して進めます。

> top.topic.words(result$topics, 5)
     [,1]         [,2]        [,3]        [,4]       [,5]       [,6]
[1,] "algorithm"  "learning"  "model"     "learning" "neural"   "learning"
[2,] "model"      "algorithm" "report"    "neural"   "networks" "network"
[3,] "algorithms" "paper"     "models"    "networks" "network"  "feature"
[4,] "data"       "examples"  "bayesian"  "paper"    "learning" "model"
[5,] "results"    "number"    "technical" "network"  "research" "features"
     [,7]        [,8]          [,9]        [,10]
[1,] "knowledge" "problem"     "paper"     "learning"
[2,] "system"    "genetic"     "decision"  "reinforcement"
[3,] "reasoning" "control"     "algorithm" "paper"
[4,] "design"    "performance" "results"   "problem"
[5,] "case"      "search"      "method"    "method"

結果は出ていますよね。

データの準備

lda.collapsed.gibbs.sampler {lda} を使用するために vocab, documents に相当するデータファイルを準備します。

vocab はコーパスを表していますので、まず全ての文書に含まれている (文字列としての) 単語 $N$ 個のリストを作成します。これにより各単語は vocab 上のインデックス $i$ で表されます。

{\tt vocab} = \{w_0, w_1, ..., w_{N-1}\} \\
{\tt vocab[}i{\tt]} = w_i

データファイルとしてはコーパスとして全文書内の単語を重複なく集めた 1行=1単語 の corpus_r.txt を作成します。

$ head -n 5 corpus_r.txt
あっせん
資本
人物
ドリップ
ランタノイド
$ wc -l corpus_r.txt
20778 corpus_r.txt

documents は文書 $d_k$ のリストになります。文書 $d_k$ は文書に含まれる単語のインデックス $i_{k,j}$ と、その単語のの文書内での出現数 $c_j$ の $m \times 2$ 行列で表されます。

{\tt documents} = \{ d_0, d_1, ..., d_k, ..., d_{m-1} \} \\
d_k = \begin{pmatrix}
i_{k,0} & i_{k,1} & ... & i_{k,n-1} \\
c_{k,0} & c_{k,1} & ... & c_{k,n-1}
\end{pmatrix}

データファイルとしては 1行=1文書 に相当し、文書中に含まれる単語の (0から始まる) インデックスと出現数をコロンで区切った bow_r.txt を用意します。

$ head -n 5 bow_r.txt
74:1    1109:1  1788:1  7000:2  10308:2 10552:1 12332:2 13489:1 14996:1 15448:1 15947:1 16354:4 17577:1 18262:1 19831:4
3256:3  5278:1  9039:1  12247:1 14529:6 17026:3
2181:1  4062:1  6270:1  6508:1  7405:1  8662:1  15448:1 18905:1
8045:2  9323:1  5934:2
288:3   624:1   691:1   820:2   1078:2  1109:2  1148:3  1251:1  2025:1  2050:1  2072:1  2090:1  2543:2  2626:1  2759:1…
$ wc -l bow_r.txt
13017 bow_r.txt

実行結果

以下のコードで実行します。LDA の処理よりデータファイルを解析して想定するデータ構造を作る方が大きいですが (R が得意な人ならもっと簡単に書けるかも?)。

require("lda")

vocab <- as.vector(read.table("corpus_r.txt", header=FALSE)$V1)

file <- file("bow_r.txt", "r")
docs <- list()
repeat {
  line <- readLines(con=file, 1)
  if(length(line) == 0) break
  doc <- NULL
  for(tc in strsplit(line, "\t")[[1]]){
    col <- c()
    for(c in strsplit(tc, ":")[[1]]){
      col <- c(col, as.integer(c))
    }
    if(is.null(doc))  doc <- cbind(col)
    else              doc <- cbind(doc, col)
  }
  docs <- append(docs, list(doc))
}
close(file)

result = lda.collapsed.gibbs.sampler(docs, 100, vocab, 50, 0.1, 0.1)

top.topic.words(result$topics, 5)

top.topic.words() を実行すると $k=100$ トピックそれぞれの特徴を表す単語が羅列されます。

     [,1]   [,2]     [,3]     [,4]     [,5]     [,6]   [,7]       [,8]
[1,] "問題" "期間"   "一番"   "目"     "部屋"   "発行" "アイドル" "箱"
[2,] "議員" "お願い" "天気"   "ロック" "家"     "今後" "ファン"   "部分"
[3,] "社会" "試用"   "中国"   "回転"   "トイレ" "お寺" "変"       "シール"
[4,] "宗教" "皆さん" "間違い" "不要"   "水"     "無理" "興味"     "ニス"
[5,] "主義" "正社員" "全部"   "中国"   "洗濯"   "内容" "素直"     "情報"
     [,9]   [,10]      [,11]    [,12]      [,13]    [,14]      [,15]  [,16]
[1,] "女性" "銀"       "画像"   "サイト"   "曲"     "お願い"   "回路" "仕事"
[2,] "人"   "近視"     "添付"   "方法"     "お願い" "マイ"     "連隊" "会社"
[3,] "男性" "大人"     "複数"   "登録"     "ピアノ" "理由"     "歩兵" "上司"
[4,] "男"   "スナック" "写真"   "ドメイン" "音楽"   "営業"     "電流" "人"
[5,] "女"   "絵本"     "アース" "ページ"   "感じ"   "ナンバー" "写真" "職場"
... (100トピックまで続く)