2010年8月5日木曜日

bigram feature について。

twitter で @takeda25 さんが指摘されていたのですが, CRF で f(y_{i-1},y_{i},x_{i}) という観測素性とラベル bigram の組みまで考慮した素性関数があまり使われないのは, 素性数が増えるからというよりは計算に時間が掛かるためだろうか? という事を私なりに考えてみました.

CRF の1事例に対するパラメータ推定に掛かるオーダーはラベル数L, 系列長Tの時に, O(L^2T) です.
forward-backward でラティス中の位置 i, ラベル j のノードの alpha を計算する際には以下の logsumexp の計算を行います.

for (k = 0; k < L; ++k)
alpha_{i,j} = logsumexp(alpha_{i,j}, alpha_{i-1,k}+cost)

この計算を全ての i, j について行います.

この時足し込んでいる cost が 観測素性+遷移素性のコストで, 観測素性 x に対してラベル bigram を考慮する場合, 遷移素性のコスト計算をする時に一緒に足す事になるため, 単純に観測された bigram 素性の数だけ加算回数が増えます.
ですが, logsumexp の呼び出し回数自体は変わりません. (最初この部分も定数倍になると勘違いしていました)

cost の加算部分だけが定数倍される程度なので, 計算時間にはそんなに影響が出ないのではないかと考えたのですが, 実は CNF だと結構影響が出ました...

CNF では cost の計算時に, logstic 関数による非線形な重み付けが必要になるため, exp の計算が入ります. exp の計算が素性数だけ増えるため, 試してみたら結構遅くなりました..

試してみたのは CoNLL2000 の学習データで, 使用したテンプレートは以下の通りです. 見事に全部 bigram テンプレートです.

# Bigram
B00:%x[-2,0,0]
B01:%x[-1,0,0]
B02:%x[0,0,0]
B03:%x[1,0,0]
B04:%x[2,0,0]
B05:%x[-1,0,0]/%x[0,0,0]
B06:%x[0,0,0]/%x[1,0,0]

B10:%x[-2,1,0]
B11:%x[-1,1,0]
B12:%x[0,1,0]
B13:%x[1,1,0]
B14:%x[2,1,0]
B15:%x[-2,1,0]/%x[-1,1,0]
B16:%x[-1,1,0]/%x[0,1,0]
B17:%x[0,1,0]/%x[1,1,0]
B18:%x[1,1,0]/%x[2,1,0]

B20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
B21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
B22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]

B


これで学習を行ってみたところ, 次のようになりました.

* iter = 1

time ./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s bigramtest.save --l2 --pool=1000000 --block=24 --lambda=2 --iter=1
labels: 22
features: 338552
bound: 3
ufeatures: 0
bfeatures: 76329
instance: 8936
gate functions: 20
parameters of gate functions: 31
uparameters: 0
bparameters: 40301712
model parameters: 40301743
epoch: 0 err:0.069070(14624/211727)
./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s 1323.90s user 41.84s system 96% cpu 23:40.72 total


ちなみに, テンプレートを元々の

U00:%x[-2,0,0]
U01:%x[-1,0,0]
U02:%x[0,0,0]
U03:%x[1,0,0]
U04:%x[2,0,0]
U05:%x[-1,0,0]/%x[0,0,0]
U06:%x[0,0,0]/%x[1,0,0]

U10:%x[-2,1,0]
U11:%x[-1,1,0]
U12:%x[0,1,0]
U13:%x[1,1,0]
U14:%x[2,1,0]
U15:%x[-2,1,0]/%x[-1,1,0]
U16:%x[-1,1,0]/%x[0,1,0]
U17:%x[0,1,0]/%x[1,1,0]
U18:%x[1,1,0]/%x[2,1,0]

U20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
U21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
U22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]

# Bigram
B


にすると, iter = 1 の学習時間は以下の通り.


time ./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s bigramtest.save --l2 --pool=1000000 --block=24 --lambda=2 --iter=1
labels: 22
features: 338552
bound: 3
ufeatures: 76328
bfeatures: 1
instance: 8936
gate functions: 20
parameters of gate functions: 31
uparameters: 1679216
bparameters: 528
model parameters: 1679775
epoch: 0 err:0.057314(12135/211727)
./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s 67.23s user 5.61s system 92% cpu 1:18.67 total


だいたい bigram テンプレート数倍だけ遅くなりました..
exp の計算が bigram テンプレート数倍必要になるので, そこが影響したようです.

CRF の場合は logistic 関数の計算無しで加算のみになるため, ここまでの影響は出ないと考えられるので, bigram 素性を使ったとしてもそこまで遅くなるという事は無いのではないかと思います.

ちなみにオール bigram テンプレートを使った時の学習に必要なメモリは約 450 MB でした.

bigram 素性を使うとタグ付け時にトークンの文脈依存性を上手く学習できたりしそうですし, 意味の無い素性ではないと思うのですが, 何であまり使われていないんでしょうね.

ちなみに, CRFSuite では高速化の為にいろいろな工夫がされているそうで, その工夫の中にラベル bigram と観測素性の組み合わせは考慮しないという前提に基づいた手法があるのか, bigram 素性を扱うと高速化手法のどれかが使えなくなってしまうっぽい?
ソースコードを落としてきてちら見してみたのですが, 人のコードを読むのが苦手な私ではすぐには問題箇所が分からず, 「おー! コード超キレイ!」とそんな感想しか言えませんでした...

時間のある時にじっくり読んで勉強させてもらいます!

2010年7月9日金曜日

a hybrid markov/semi-markov cnf

名前だけ知りつつもどんな物なのか知らなかったので、5月に semi-markov crf の論文 を読んでみました。

semi-crf は入力系列に対する最適なセグメンテーションを学習する学習器です。
crf が P(Y|X, W) を求めていたのに対して、semi-crf では P(S|X, W) を推定します。

文章で書くと、2つの違いは

* crf では入力系列 X に対して生成されうるすべてのラベル系列 Y' に対して、正解系列とそれ以外の系列を弁別するよう学習する
* semi-crf では入力系列 X に対して生成されうるすべてのセグメンテーション結果 S' に対して、正解のセグメンテーションとそれ以外のセグメンテーションを弁別するよう学習する

ということになります。

semi-crf では入力系列が与えられた際に、t_j , ... , u_j までを1つのセグメント s_j として、1つのラベルを割り当てます。
この時、t_j 、 u_j は常に以下を満たす必要があります。

1 <= t_j <= u_j <= |S| , t_{j+1} = u_j + 1

これは、セグメントの終了位置と次のセグメントの開始位置は一致するよということです。
で、crf で使われていた素性関数 f の代わりに segment feature function g ( g^k(y_j,y_{j-1},X,t_j,u_j) )を使います。
といっても、あまり大きく変わっているという事は無くて、crf の素性関数がトークン x_i とラベル y_i のペアに対して 1 or 0 を返していたのが、セグメント s_i とラベル y_i のペアに変わった程度です。

大きな違いは以上なのですが、パラメータ推定をする際のラティスの形が違います。

crf ではラティスはきれいな格子でしたが、semi-crf では各位置 i について、i で始まるセグメントと i で終了するセグメントのみ接続を許すため、結構複雑な形をしています。
ただ、これは形態素解析の物と同じなので 自然言語処理-基礎と応用- の1章、形態素解析の 1.1.5 節あたりを読むと分かりやすいと思います。

パラメータ推定は crf と同じ方法が使えます。

とりあえずここまで理解したところで、cnf で作ってみよう! と考えて作り始めたのですが、いざ作ってみて失敗に気づきました。

semi-crf ではセグメント単位で同じラベルをふるため、ラベル付けは IOB とかじゃ無く、トークン列に対して AAA、BBB のようにつければ良いやと思っていたのですが、いざ動かす所で、正解の与え方ではまりました。
というのも、AA, AA, AA というラベル列を正解コーパスでつけたつもりでも、AAA, AAA、AAAAAA のどれとも区別が出来なくなってました。。

それで、方向を修正し、結局ラベル付けは IOB タグに限定する事にして(cnf ではそんな事はないのですが...セグメントの開始位置を区別しつつ、cnf と同じフォーマットのコーパスを入力に取れるようにする方法が他に思いつきませんでした。)、せっかく IOB を使うのだから、従来使えていた cnf と同じ素性も使えるようにしてみました。

IOB タグを使い、cnf で使っていたトークン素性を使えるようにするという事は、あるトークン x がセグメントの開始 B になりやすいだとか、また逆に B にはなりにくいという様な情報を学習できます。
さらに、 cnf では素性テンプレートを使い、カレントの位置 i の n 個前のトークンが何であったか等も素性として扱えるため、セグメントの外にあるトークンの情報も扱えるという利点があります。

まさにこれと同じ事をやっているっぽい論文

A Hybrid Markov/Semi-Markov Conditional Random Field for Sequence Segmentation

があったので、名前をお借りして a hybrid markov/semi-markov cnf と名付けてみました。
上の論文自体はまだちゃんと読んで理解したわけではないので、自信を持って同じとも違うとも言えないのですが、ちらっと読んだところ、私と同じくセグメント素性にトークン素性を追加したのではなく、セグメント素性をトークン素性の和で表しているように見えます。
そうするとコストの計算時に、トークン素性のコストを予め計算しておいて、各セグメントで共通する部分についてはキャッシュを使うという事が出来るようになります。
ちゃんと読んでからでないと断言はできませんが、あってたらそこが私の作った学習器との違いになりそうです。

とりあえず、コードは以前 cnf をアップした所にそのまま追加しました。

http://code.google.com/p/cnf/

semi* が追加部分です。

素性テンプレートは以下のようになっています。

# Unigram Segment
S:%s[15,0]

# Bigram Segment
T

# Unigram
U00:%x[-2,0,0]
U01:%x[-1,0,0]
U02:%x[0,0,0]
U03:%x[1,0,0]
U04:%x[2,0,0]
U05:%x[-1,0,0]/%x[0,0,0]
U06:%x[0,0,0]/%x[1,0,0]

U10:%x[-2,1,0]
U11:%x[-1,1,0]
U12:%x[0,1,0]
U13:%x[1,1,0]
U14:%x[2,1,0]
U15:%x[-2,1,0]/%x[-1,1,0]
U16:%x[-1,1,0]/%x[0,1,0]
U17:%x[0,1,0]/%x[1,1,0]
U18:%x[1,1,0]/%x[2,1,0]

U20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
U21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
U22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]

# Bigram
B


トークン素性については以前と変わりありません。

# Unigram Segment
S:%s[15,0]

# Bigram Segment
T

の部分が追加したテンプレートの仕様で、

S:%s[ length , gate function に与える重みの初期値 ]

が unigram segment のテンプレート、

T:%s[ length , gate function に与える重みの初期値 ]

が bigram segment のテンプレートになります。
length は考慮するセグメントの最大長です。

S はセグメントの S ですが、 T は良い頭文字が思いつかなかったので、とりあえず T にしました。
トークン素性の bigram の B 同様、セグメントのラベル遷移だけを素性にする場合は T だけ書いておけば OK です。

ソースを落としてきたら、

# make
# ./src/semicnflearn src/semitemplate data/conll2000/train.txt semicnf.save
# ./src/semicnftagger src/semitemplate semicnf.save data/conll2000/test.txt

みたいにして使います。

コーパスは従来の系列ラベリングで使えた物は大抵使えるんじゃないかと思いますが、セグメンテーションに適した手法なのでどうせならちゃんとしたタスクのコーパスで評価したいところですね。

調べてみた所、公開されている semi-crf の実装もあまり無いようですし、素性テンプレートの使える semi-markov cnf はそこそこ貴重なのではないかと自分で思ってみたりしています。

今回から、学習器の学習率のスケジューリングとL1 正則化は
Stochastic Gradient Descent Training for L1-regularized Log-linear Models with Cumulative Penalty
の手法に変更しました。

@unnonouno さんのブログ
http://unnonouno.blogspot.com/2010/02/l1-sgd.html
で紹介されていた手法で、説明が秀逸ですw

[2010.7.21 追記]
template の仕様を変更しました。
以前の仕様では、unigram segment と bigram segment で各1つずつしかテンプレートを記述できなかったのですが、複数のテンプレートを記述できるようにしました。
それに伴って、

S:%s[ length , gate function に与える重みの初期値 ]

から、

S:%s[ length, col, 重みの初期値]

と、col を指定できるようにテンプレートのマクロも変更しています。

おまけで、新しい template の例を以下に。

# Unigram Segment
S01:%s[15,0,0]
S02:%s[15,1,0]

# Bigram Segment
T

# Unigram
U00:%x[-2,0,0]
U01:%x[-1,0,0]
U02:%x[0,0,0]
U03:%x[1,0,0]
U04:%x[2,0,0]
U05:%x[-1,0,0]/%x[0,0,0]
U06:%x[0,0,0]/%x[1,0,0]

U10:%x[-2,1,0]
U11:%x[-1,1,0]
U12:%x[0,1,0]
U13:%x[1,1,0]
U14:%x[2,1,0]
U15:%x[-2,1,0]/%x[-1,1,0]
U16:%x[-1,1,0]/%x[0,1,0]
U17:%x[0,1,0]/%x[1,1,0]
U18:%x[1,1,0]/%x[2,1,0]

U20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
U21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
U22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]

# Bigram
B

2010年6月13日日曜日

グラフラプラシアンで推薦

以前縁あって小町さんと一緒に仕事をさせてもらい論文に名前を載せてもらったのですが、会社だけでなく自宅でもちょっと使いたいなーということもあり、実装してみることにしました。

参考にしたのは以下の論文です。
ラプラシアンラベル伝播による検索クリックスルーログからの意味カテゴリ獲得

元論文と違うのは、インスタンス-パターン行列の要素を単純な頻度から別の尺度に変えている点です。
元々そのまんま実装してみたところ、非常にレアな場合なのですが、ジェネリックパターン1つのみと共起するようなインスタンスがあった場合に、これが上位に出やすくなるという問題が発生し、どうにかできないかなと模索していたところ、小町さんからアドバイスを頂き、それを基に手を加えています。

とりあえず動作検証のためにMovieLens Data Setsを使って実験してみました。

最初にデータのフォーマットをツールの入力形式へ変更。
perl swap.pl < ../ml-data/u.data | LC_ALL='C' sort > mlens.data


#!/opt/local/bin/perl

# user id | item id | rating | timestamp
use warnings;
use strict;

my %dic;
my $file = "../ml-data/u.item";
open (F, $file) or die;
while ()
{
chomp;
my @t = split(/\|/);

unless (defined $dic{$t[0]})
{
$dic{$t[0]} = $t[1];
}
}
close (F);
while (<>)
{
chomp;
my @tokens = split(/\t/);
print $dic{$tokens[1]}, "\t$tokens[0]\t$tokens[2]\n";
}


入力ファイルのフォーマットは次のようにTSV形式で<インスタンス-パターン-評価値>として、インスタンスでソートしておきます。


'Til There Was You (1997) 152 4
'Til There Was You (1997) 178 3
'Til There Was You (1997) 223 1
'Til There Was You (1997) 299 2
'Til There Was You (1997) 342 1
'Til There Was You (1997) 416 3
'Til There Was You (1997) 530 2
'Til There Was You (1997) 532 3
'Til There Was You (1997) 782 2
...


で、実行。

./quetchup -f mlens.data -o test.save

デフォルトでの実行は alpha = 0.0001, iteration = 10, エッジのカットに使う閾値 0.1 としています。

./quetchup -c sample.txt -r test.save | head -50

sample.txt

Snow White and the Seven Dwarfs (1937)


白雪姫がお気に入りとすると、以下のようにディズニー系の映画を上位にレコメンドすることができます。


[Instance]
Snow White and the Seven Dwarfs (1937):0.999801
Cinderella (1950):5.19753e-07
Fantasia (1940):5.06323e-07
Dumbo (1941):4.17243e-07
Beauty and the Beast (1991):3.94127e-07
Aladdin (1992):3.88318e-07
Mary Poppins (1964):3.71035e-07
Jungle Book, The (1994):3.67852e-07
Sword in the Stone, The (1963):3.67118e-07
Lion King, The (1994):3.64219e-07
Pinocchio (1940):3.60772e-07
Robin Hood: Prince of Thieves (1991):3.5665e-07
Aristocats, The (1970):3.51889e-07
Fox and the Hound, The (1981):3.34974e-07
Star Trek VI: The Undiscovered Country (1991):3.24104e-07
Alice in Wonderland (1951):3.19468e-07
Homeward Bound: The Incredible Journey (1993):2.97177e-07
Sound of Music, The (1965):2.97109e-07
Walk in the Sun, A (1945):2.96338e-07
Very Natural Thing, A (1974):2.96338e-07
Bedknobs and Broomsticks (1971):2.93205e-07
Star Trek IV: The Voyage Home (1986):2.90837e-07
Gone with the Wind (1939):2.89412e-07
Nightmare Before Christmas, The (1993):2.85079e-07
Star Trek: Generations (1994):2.8431e-07
Winnie the Pooh and the Blustery Day (1968):2.82242e-07
101 Dalmatians (1996):2.81954e-07
Swiss Family Robinson (1960):2.6257e-07
Aladdin and the King of Thieves (1996):2.54459e-07
Grease (1978):2.53118e-07
African Queen, The (1951):2.52688e-07
Ghost (1990):2.52399e-07
Hunt for Red October, The (1990):2.4501e-07
Miracle on 34th Street (1994):2.42372e-07
Home Alone (1990):2.4004e-07
Top Gun (1986):2.35881e-07
Abyss, The (1989):2.35302e-07
Little Women (1994):2.33329e-07
Old Yeller (1957):2.32874e-07
Pocahontas (1995):2.31783e-07
Mrs. Doubtfire (1993):2.31762e-07
Jumanji (1995):2.31231e-07
Pete's Dragon (1977):2.29198e-07
Christmas Carol, A (1938):2.28481e-07
New York Cop (1996):2.28269e-07
Maltese Falcon, The (1941):2.26072e-07
Goofy Movie, A (1995):2.25177e-07
James and the Giant Peach (1996):2.23962e-07
Star Trek: The Motion Picture (1979):2.23878e-07


ですが、よく見るとスタートレックとかが入り込んでいます。
そこで、次のように嫌いなものリストを用意します。

sample2.txt

Walk in the Sun, A (1945)
Star Trek: Generations (1994)


./quetchup -c sample.txt -n sample2.txt -r test.save | head -50


[Instance]
Snow White and the Seven Dwarfs (1937):0.9998
Cinderella (1950):2.96839e-07
Fantasia (1940):2.81548e-07
Dumbo (1941):2.54658e-07
Aladdin (1992):2.4008e-07
Lion King, The (1994):2.30137e-07
Beauty and the Beast (1991):2.10186e-07
I Don't Want to Talk About It (De eso no se habla) (1993):2.01313e-07
Jurassic Park (1993):1.96396e-07
Coldblooded (1995):1.95535e-07
Robin Hood: Prince of Thieves (1991):1.74025e-07
Fox and the Hound, The (1981):1.71627e-07
African Queen, The (1951):1.69855e-07
Sound of Music, The (1965):1.69781e-07
Mary Poppins (1964):1.64425e-07
Wizard of Oz, The (1939):1.6242e-07
Swan Princess, The (1994):1.56013e-07
Homeward Bound: The Incredible Journey (1993):1.53814e-07
Mamma Roma (1962):1.52174e-07
Good Man in Africa, A (1994):1.50378e-07
Some Like It Hot (1959):1.46221e-07
Citizen Kane (1941):1.42743e-07
Anne Frank Remembered (1995):1.40035e-07
Dances with Wolves (1990):1.36046e-07
It's a Wonderful Life (1946):1.3514e-07
Hunt for Red October, The (1990):1.3501e-07
Young Frankenstein (1974):1.34608e-07
Rebecca (1940):1.32806e-07
Winnie the Pooh and the Blustery Day (1968):1.30989e-07
It Happened One Night (1934):1.30592e-07
M (1931):1.30189e-07
Maltese Falcon, The (1941):1.29532e-07
Lassie (1994):1.28929e-07
Foreign Correspondent (1940):1.28287e-07
Grease (1978):1.2729e-07
Alice in Wonderland (1951):1.23873e-07
Quiz Show (1994):1.22688e-07
Singin' in the Rain (1952):1.22222e-07
Sword in the Stone, The (1963):1.18912e-07
Lawrence of Arabia (1962):1.15236e-07
To Kill a Mockingbird (1962):1.14732e-07
Swiss Family Robinson (1960):1.14698e-07
Drop Dead Fred (1991):1.13925e-07
Vertigo (1958):1.13817e-07
Top Gun (1986):1.12622e-07
Bonnie and Clyde (1967):1.12262e-07
Blue Angel, The (Blaue Engel, Der) (1930):1.10351e-07
Annie Hall (1977):1.09473e-07
Turning, The (1992):1.07052e-07


と、無事スタートレックとか消えてくれました(スタートレック好きな人ごめんなさい!)。
ちなみに、head で出力を絞っているので見えていませんが、出力にはスコア付きのパターンリストも含まれます。
これで、どんなパターンで各インスタンスがつながっているのかも分かりやすいです。

お気に入りリストが1つしかなくても十分な精度が出そうだし、複数あっても無問題。
嫌いなものあればそれに近いものを推薦候補から外しやすい。

この手法で推薦とか悪くないなと思う今日この頃な訳です。

ソースコードを公開したいのですが、Quetchup アルゴリズムは一応特許が取られているはずなので、それがどんな内容か分からないと不安だったりします。。

2010年3月29日月曜日

gibbs sampling

3月に入って netwalker を購入しました。

電車通勤の時間を利用してコードを書きたいなと思っていたのですが、なかなかよろしい感じです。
といってもキーボードは使いやすいとは言い難く、長いコードは書きたくないですが...

ということで、netwalker で作成したプログラム第一号を公開します。
第一号は gibbs sampling を使ったモチーフ抽出アルゴリズムの実装です。

前々から gibbs sampling について調べてはいたのですが、適度な練習問題が無く実装はしていませんでした。
今月になって購入したバイオインフォマティクスの数理とアルゴリズムにちょうど良い例があったので、それを実装してみました。

そもそもモチーフって何という話ですが、一言で言うとタンパク質の配列パターンの事を指すそうです。
私は専門ではないので詳しい事は分かりませんが、ここでは与えられた複数のタンパク質配列から、特徴的な配列パターンを見つけてくる事をモチーフ抽出と言ってよさそうです。

具体的にどうモチーフを抽出するかですが、載っている例では、各タンパク質A、C、T、Gの出現確率をP(A)、P(C)、P(T)、P(G)、部分配列中の位置 j に各タンパク質が現れる確率を f_j(A)、f_j(C)、f_j(T)、f_j(G)とした時に、相対エントロピーを最大化するように部分配列を選んでくる事でモチーフを抽出します。

これをローカルマルチプルアラインメントと呼ぶそうですが、この問題はまともにやると NP 困難なので、gibbs sampling を使って効率的に解くよ、というお話しになっています。

という事で、アルゴリズムは以下の通り。

1. 各入力配列 s_k から, 長さ L の部分配列 t_k を一様ランダムに選択 (k = 1,...,n)
2. 入力配列から, s_i を一様ランダムに選択
3. f'_j(a) を t_i 意外の t_k から計算したモチーフ領域の j 番目の列に置ける文字 a の頻度とする ( f'_j(a) = (|{k|t_k[j] = a, k != i}|) / (n-1))
4. s_i から長さ L の部分配列 t'_i を, Πj f'_j(t'_i[j])/p(t'_i[j]) に比例する確率でランダムに選択
5. t_i を t'_i で置き換える
6. 2-5 を十分な回数繰り返す

ステップ 3. の頻度は確率の間違えかな?と思うのですが、一応そのまま書いておきます。
ステップ 4. では求めた事後分布に従って新しくモチーフ候補を選択し直すわけですが、ここはランダムで選択しないと初期値から動かなくなるので、まず計算したスコアの配列を確率分布として使えるように総和が1となるように正規化してから、PRMLの下巻に紹介されている棄却サンプリングを使いました。
提案分布には一様分布を使っています。

コードと実行例は以下。


#!/usr/local/bin/perl

use warnings;
use strict;

my $L = 5;
my @dataset;
my %wordset;
my $words = 0;

while(<>)
{
chomp;
my $line = $_;
my @sequence = split(//,$line);
push @dataset,\@sequence;
}

# init words and probs
$words = &setwords(\%wordset, \@dataset);
&initwordprobs($words,\%wordset);
my $datasize = @dataset;

print "dataset\n";
foreach my $ar (@dataset)
{
print "@$ar\n";
}

# init motif set
my @pos;
foreach my $ar (@dataset)
{
my $p = int rand(@$ar-$L);
push @pos, $p;
#print $p,"\n";
}

# output first
print "init\n";
for (my $i = 0; $i < @dataset; $i++)
{
my $ar = $dataset[$i];
my $p = $pos[$i];
for (my $j = $p; $j < $p+$L; $j++)
{
print "$ar->[$j] ";
}
print "\n";
}

# iterration
for (my $iter = 0; $iter < 50; $iter++)
{
#print "epoch: $iter\n";
my $id = int rand($datasize);
my %wordsfreq;
# wordset * position
# +---+-------------+
# | |a|b|c|d|e|f|g|
# +---+-------------+
# | 0 | |
# | 1 | |
# | 2 | frequency |
# | . | |
# | L | |
# +---+-------------+
&initwordsfreq(\%wordset,$L,\%wordsfreq);
# count word frequency
&countwordfreq($id,\@dataset,\%wordsfreq,\@pos);
# pick up new motif i
my $ar = $dataset[$id];
$pos[$id] = &pickup($id,$ar,$L,\%wordsfreq, \%wordset);
}

# output
print "result\n";
for (my $i = 0; $i < @dataset; $i++)
{
my $ar = $dataset[$i];
my $p = $pos[$i];
for (my $j = $p; $j < $p+$L; $j++)
{
print "$ar->[$j] ";
}
print "\n";
}

# subroutine
sub pickup
{
my $id = shift;
my $ar = shift;
my $L = shift;
my $wordsfreqhr = shift;
my $wordsethr = shift;

my @scores;
my $total = 0;
for (my $p = 0; $p <= @$ar-$L; $p++)
{
my $score = 1;
my $index = 0;
for (my $j = $p; $j < $p+$L; $j++)
{
my $w = $ar->[$j];
$score *= $wordsfreqhr->{$index}->{$w}/$wordsethr->{$w};
$index++;
}
push @scores, $score;
$total += $score;
}
return int rand(@scores) if $total == 0;
for(my $i = 0; $i < @scores; $i++)
{
$scores[$i] /= $total;
}
#print "@scores\n";
my $flg = 1;
while($flg)
{
my $npos = int rand(@scores);
my $s = rand(1);
if ($s <= $scores[$npos])
{
$flg = 0;
return $npos;
}
}
}
sub countwordfreq
{
my $id = shift;
my $datasetar = shift;
my $wordsfreqhr = shift;
my $posar = shift;

for (my $i = 0; $i < @$datasetar; $i++)
{
next if ($i == $id);
my $ar = $datasetar->[$i];
for (my $j = $posar->[$i]; $j < $L+$posar->[$i]; $j++)
{
my $p = $j-$posar->[$i];
my $w = $ar->[$j];
$wordsfreqhr->{$p}->{$w}++;
}
}
for (my $l = 0; $l < $L; $l++)
{
foreach my $k (keys %{$wordsfreqhr->{$l}})
{
$wordsfreqhr->{$l}->{$k} /= (@$datasetar-1);
}
}
}

sub initwordsfreq
{
my $wordsethr = shift;
my $L = shift;
my $wordfhr = shift;
for (my $i = 0; $i < $L; $i++)
{
my %freqonpos;
foreach my $k (keys %$wordsethr)
{
$freqonpos{$k} = 0;
}
$wordfhr->{$i} = \%freqonpos;
}
}
sub initwordprobs
{
my $words = shift;
my $wordsethr = shift;

foreach my $k (keys %$wordsethr)
{
$wordsethr->{$k} = 1/$words;
}
}

sub setwords
{
my $wordshr = shift;
my $dataar = shift;
my $words = 0;

foreach my $ar (@$dataar)
{
foreach my $w (@$ar)
{
unless (defined $wordshr->{$w})
{
$wordshr->{$w} = defined;
$words++;
}
}
}
return $words;
}


サンプルデータ1

c a t c g a a a a a a
a a a a c a t c g a a
a a c a t c g a a a a
a a a a a a c a t c g


サンプルデータ2

t c a g a g a t t c g t a g
t t c a t t c g g g c g t
c g a t t c g c g a c t c
t g a a t t c g g a a


実行例1

# perl gibbs.pl < sample1.txt
dataset
c a t c g a a a a a a
a a a a c a t c g a a
a a c a t c g a a a a
a a a a a a c a t c g
init
g a a a a
a a a c a
a t c g a
a a a a a
result
c a t c g
c a t c g
c a t c g
c a t c g


実行例2

# perl gibbs.pl < sample2.txt
dataset
t c a g a g a t t c g t a g
t t c a t t c g g g c g t
c g a t t c g c g a c t c
t g a a t t c g g a a
init
t c a g a
c a t t c
t c g c g
a a t t c
result
a t t c g
a t t c g
a t t c g
a t t c g


上手い事特徴的な文字列を抽出できる事が確認できました。

2010年2月28日日曜日

Conditional Neural Fields on Google Code

CNF の著者の Jian Peng 氏に特許について質問をしてみたところ、問題ないということでしたので Google Code にプロジェクトを作成してコードを公開しました。

http://code.google.com/p/cnf/

あまりちゃんとした実装ではないので、使用は自己責任でお願いします。
間違ってるかもしれないので、間違いがあれば教えてくれると嬉しいです。

mercurial で管理しているので、以下のコマンドで落としてきて使用できます。

$ hg clone https://cnf.googlecode.com/hg/ cnf
$ cd cnf
$ make
$ ./src/cnflearn src/template data/conll2000/train.txt test.save
$ ./src/cnftagger src/template test.save data/conll2000/test.txt

cnflearn_main.cpp と cnftagger_main.cpp を見てもらえば分かる通り、実行用のプログラムはむちゃくちゃ手抜です。

logsumexp の計算部分には、NAIST の浅原先生の資料にあるコードを使わせて頂きました。

2010年2月16日火曜日

Conditional Neural Fields

年越ししてから既に2ヶ月が過ぎ、2月も終わりが見えかけてきた今日この頃です。
生存報告をかねて、少しだけ最近やっていた事を書いておきます。

去年の年末くらいに、面白そうな論文を見つけたのでそれを読みつつ、実装していました。
NIPS2009 で発表された論文です。

その名も Conditional Neural Fields 。
http://books.nips.cc/papers/files/nips22/NIPS2009_0935.pdf

名前から何か感じるところの有る人もいそうですが、これはCRFに隠れ層を加えて、非線形にした物になります。

自分のメモ用に、先に簡単に CRF についておきます。
--
CRF の説明はNLP2006のチュートリアル資料が割と分かりやすいです。
http://nlp.dse.ibaraki.ac.jp/~shinnou/lecture/nl/NLP2006.pdf

私はオンライン学習器でしかCRFを作った事が無いので、ここではそのつもりで書きます。
CRF では素性関数に対応するパラメータを学習するわけですが、その式は大雑把に書くと次のようになります。

Λ_t+1 = Λ_t + η * (正解の素性関数ベクトル - 期待素性関数ベクトル)

η は学習率、Λ は素性関数に対応するパラメータベクトルです。
正解の素性関数ベクトルは入力系列が与えられれば素性と正解のラベルのペアを取り出すだけで作れます。
要はどの素性とラベルのペアが何個有ったか、どのラベルとラベルの遷移が何回あったか、をカウントしてベクトルにしています。

期待素性関数ベクトルはちょっと計算が面倒ですが、現在のパラメータΛのもとで、入力系列から生起しうる全てのラベル系列を考えます。
各ラベル系列が生起した際の、各素性関数の期待値を計算し足し合わせた物が期待素性関数ベクトルです。
まともに全てのラベル系列を生成して試すともの凄い計算量になるので、これは Forward-Backward アルゴリズムを用いて計算するのが一般的です。
--

前置きが長くなりましたが、 CNF について書きます。
ちなみに間違ってるかもしれません。あまり期待しては行けません。
この記事は後でこっそり修正されてる可能性もあります。

CRF との違いは、観測素性の扱いになります。
CRF では、素性関数 f(x,y) は素性xとラベルy のペアを観測した際に1になるような関数です。

これが CNF では、f(h(X,t),y) になります。

数式は元論文のままだと比較しにくいので変えてます。
関数h はロジスティック関数で、X に対して非線形な値を取ります。
t は系列のカレントの位置を表しています。

f(h(X,t),y) は組み合わせ素性xとラベルy のペアを観測した際に、h(X,t)を返す関数と思ってください。
X が大文字になっている理由は、CNF では入力系列に対して隠れ層で組み合わせ素性を取り出すためです。
素性関数 f(h(X),y) は入力系列 X に対して、位置tに置ける組み合わせ素性 x とラベル y のペアに対して、h(X,t) を返します。

h(X,t) は中で何をしているのかと言うと、入力系列 X とパラメータθベクトルの内積を計算して、その値をロジスティック関数に乗せて非線形な値にしています。

こうする事で、単純に組み合わせ素性が観測されたというだけではなく、その素性に対して非線形な重みを与えています。

次に、パラメータの最適化ですが、実は素性関数が非線形になっているという事を無視すると、観測素性についてのパラメータと、遷移素性についてのパラメータの更新の式は一緒になります。

ただし、正解の素性関数ベクトルはどの素性とラベルのペアが何個有ったか、どのラベルとラベルの遷移が何回あったか、ではなく、素性に対して非線形な重みを計算し足し合わせたものになります。ラベルとラベルの遷移は相変わらずただのカウントです。
要は、カウントの仕方が1ずつ足していたのが、実数値ずつ足すようになった物です。

期待素性関数ベクトルも同じ要領です。 入力系列から生成されうる全てのラベル系列を考えて、各ラベル系列が生起した際の各素性関数の期待値を計算して足し合わせます。 ただし期待値を計算する際には素性に対して、これまで有った/無かったで済ませていたところにh関数の値が入ります。

パラメータ推定での1つ大きな違いは、これに加えて隠れ層のパラメータθについても更新することです。

θの更新は、次の式で行われます。
θt+1 = θt + Σt w_yt,g * ∂h/∂θ - E (Σt w_yt^,g * ∂h/∂θ)

wは素性関数に対応した重みパラメータです。
重みパラメータΛ = (W, U, θ) と思ってください。
Wは観測素性に対応する重みパラメータ
Uは遷移素性に対応する重みパラメータ
θは隠れ層のパラメータ
E は期待値です。

かなりやっつけな式なので詳しくは元論文を読んだ方が良いです。
要は、正解系列中の観測素性の重みベクトルと、入力系列から生成されうる全ラベル系列を考えた時の各観測素性の期待値ベクトルの差をなくすようにθを更新しています。

という事で、実装してみました。
実装は sgd + FOLOS です。

コード公開したいんだけど、権利関係やら何やら良く分からない。
論文の著者にメールして特許についてとか聞いといた方が良いと言われているのですが、そこで停止中。

--
ちなみに、CNF は組み合わせ素性、CRF は組み合わせ無し、のような印象を受けますが、CRF++やCrfSgdを使った事が有る方は分かると思いますが、CRF でも素性テンプレートを使って組み合わせ素性を扱えます。

素性テンプレートを使った組み合わせ素性の抽出はまさに CNF で言うところの中間層だと思っています。
なので、CNF と CRF の違いは実際のところロジスティック関数による非線形な値を観測素性に与えているところと、θの学習が入っているところだというのが私の理解です。

私の実装もそのまんまで、「ロジスティック関数による非線形な値を観測素性に与えているところと、θの学習を入れた」素性テンプレートを使える CRF です。

--
conll2000 のチャンキングで比較をしてみたので一応載せておきます。
一応素性テンプレートは揃えてありますが、内部での処理が微妙に違うため素性数は厳密に一緒ではありません。

まあ、見ての通りだと思いますが、恐らくロジスティック関数による非線形な値云々よりも、組み合わせ素性を使えるかどうかの方が重要で、CRF でもそれが出来れば問題ないんだろうなという感じでした。

元々Kernelなどを用いた非線形な分類の実体は、素性の組み合わせを考慮してより高次元空間へ写像して、その空間で線形分離しているわけで、素性テンプレートを使って組み合わせ素性を扱っている以上、あんま変わんない結果が出たとしても別段驚く事じゃない気もしますね。

* CrfSgd
[Epoch 50] -- wnorm: 9273.81 total time: 1664.69 seconds
Training perf: sentences: 8936 loss: 4473.67 objective*n: 9110.58
misclassifications: 794(0.375011%)
accuracy: 99.62%; precision: 99.32%; recall: 99.18%; FB1: 99.25
Testing perf: sentences: 2012 loss: 4429.27 objective*n: 5473.3
misclassifications: 1893(3.99561%)
accuracy: 96.00%; precision: 93.84%; recall: 93.62%; FB1: 93.73
Saving model file model.gz.
Done! 1664.69 seconds.

* CNF
epoch: 49 err:0.002985(632/211727)
./cnflearn template data/conll2000/train.txt testbb.save 3728.58s user 200.56s system 69% cpu 1:34:54.39 total

./cnftagger template testbb.save data/conll2000/test.txt | ./conlleval
processed 47377 tokens with 23852 phrases; found: 23718 phrases; correct: 22293.
accuracy: 96.02%; precision: 93.99%; recall: 93.46%; FB1: 93.73
ADJP: precision: 81.08%; recall: 75.34%; FB1: 78.11 407
ADVP: precision: 83.71%; recall: 80.72%; FB1: 82.19 835
CONJP: precision: 55.56%; recall: 55.56%; FB1: 55.56 9
INTJ: precision: 100.00%; recall: 50.00%; FB1: 66.67 1
LST: precision: 0.00%; recall: 0.00%; FB1: 0.00 0
NP: precision: 94.39%; recall: 93.88%; FB1: 94.14 12355
PP: precision: 96.81%; recall: 97.78%; FB1: 97.29 4859
PRT: precision: 83.15%; recall: 69.81%; FB1: 75.90 89
SBAR: precision: 87.72%; recall: 85.42%; FB1: 86.55 521
VP: precision: 93.95%; recall: 93.62%; FB1: 93.78 4642