のんびり読書日記

日々の記録をつらつらと

久しぶりに開発環境を整理してみる

普段使ってるdotfileをgithubに上げるために、Linuxの開発環境を少し整理してみた。

vimrcを整理

vimrcはかなり前に整理したものをずっと使ってたんだけど、周りの人の環境を聞いたり、VIMテクニックバイブルを読んだり>して以下のように設定を追加してみた。

まだほとんど使いこなせてないけど、unite, neocomplcacheが便利。知らない間にvimもすごい便利になってたんだなあ…。

Vimテクニックバイブル ?作業効率をカイゼンする150の技

Vimテクニックバイブル ?作業効率をカイゼンする150の技

tcshからbashに移行

Linux上でのシェルは以前からtcshを使っていたんだけど、最近のディストリビューションだとデフォルトではtcshが入っていないことが多く、毎回tcshを入れるのが面倒になってきたのでbashに移行してみた。tcshでもあまり大したことはしてなかったので、移行といってもcshrcの内容をbashrcに移していくだけなんだけど。

zshに移行した方が面白そうだけど、それはまた気が向いたらで。

確率的勾配降下法による行列分解を試してみた

前々回のNMF(Non-negative Matrix Factorization)に続いて行列分解ネタです。言語処理学会全国大会のチュートリアル「推薦システム -機械学習の視点から-」で紹介されていた、確率的勾配降下法による行列分解を試してみました。チュートリアルの資料は公開されていないようでしたので、元論文の方のリンクを張っておきます。実際には同じ著者の別の論文を引用されてましたが、僕には下の論文の方が分かりやすかったのでこっちで。

作成したコードは以下に置いてあります。行列演算にEigenを使用しているので、実行する場合は事前にEigenをインストールしておいてください。また入力データにはMovieLensの100Kのデータを使用したので、こちらも合わせてダウンロードしてください。

比較と勉強を兼ねて、以下の3つの手法で行列を分解するクラスを入れています。手法の切り替え用オプションは用意してないので、適宜ソースコードを修正して下さい。

  • 勾配降下法(gradient descent)
    • データ全体から勾配を計算。大規模データには向かない
  • 確率的勾配降下法(stochastic gradient descent)
    • データ全体ではなくデータ1個ずつ行列の更新を実行。省メモリ
  • 確率的勾配降下法にバイアス項を追加(stochastic gradient descent)
    • 評価値のゆらぎをバイアスとして補正
      • ユーザAは高評価をつけやすい、アイテムBは低評価ばかりつけられているなど。
    • 今回は全体の平均、ユーザごとの評価の揺らぎ、アイテムごとの評価のゆらぎを使用

数式など詳細については上記の論文をご参照下さい。

実行する時は以下のようにします。MovieLensの100kのデータは交差検定用に初めからトレーニング・テストのペアが5セット用意されているので、それを使って5分割交差検定を実行します。またレートの予測の評価にはNetflix Prizeでも使用されている、二乗平均平方根誤差(RMSE)を表示しました。ちなみにNetflix Prizeの2009年度の優勝チームのRMSEは0.8558だったようです(元のデータが違うので単純に比較できないでしょうが)。

% g++ -Wall -O3 factorize_sgd.cc -o factorize_sgd
% ./factorize_sgd
Usage: ./factorize_sgd dir ncluster niter eta lambda

% ./factorize_sgd /path/movilens 10 100 0.1 0.1
Training data: /path/movielens/u1.base
Test data:     /path/movielens/u1.test
Factorizing input matrix ...
RMSE=0.975

Training data: /path/movielens/u2.base
Test data:     /path/movielens/u2.test
Factorizing input matrix ...
RMSE=0.968

Training data: /path/movielens/u3.base
Test data:     /path/movielens/u3.test
Factorizing input matrix ...
RMSE=0.967

Training data: /path/movielens/u4.base
Test data:     /path/movielens/u4.test
Factorizing input matrix ...
RMSE=0.964

Training data: /path/movielens/u5.base
Test data:     /path/movielens/u5.test
Factorizing input matrix ...
no data: row:863 col:1680
no data: row:896 col:1681
no data: row:916 col:1682
RMSE=0.960

Result of cross validation: RMSE=0.967

各手法での結果を比較すると、こんな感じになりました。

手法 RMSE クラスタ 反復数 eta lambda
勾配降下法 1.012 10 100 0.001 0.001
確率的勾配降下法 0.967 10 100 0.1 0.1
確率的勾配降下法(バイアス項つき) 0.981 10 100 0.1 0.1

バイアス項を加えたときが一番精度がよくなるはずなのですが、普通の確率的勾配降下法よりも結果が悪くなってますね。。。うーん、理解が間違ってるのかバグってるのか。それと今回は全体の平均、ユーザごとの評価の揺らぎ、アイテムごとの評価のゆらぎをバイアス項として加えたのですが、論文にはこれ以外に、ユーザの行動履歴(映画なら鑑賞した映画のリストなど)を組み合わせる話や、評価時刻を考慮したモデルなども書かれていましたので、次回はその辺を試したいと思います。

あとeta, lambdaパラメータの設定がちょっと面倒で、適切に変更してやらないとすぐに値が膨れ上がって、正しい結果が得られなくなってしまいました。更新時にどうにか抑えることはできないのかなぁ。あとパラメータいじって最適なものを自分で探すってのはちょっとイヤですね。自動で最適なものを探すようにしたい。

行列演算部分はEigenを使ったのですが、今回ぐらいのだと大したことをしていないので、別に自前で書いても良さそう。無駄に他のライブラリへの依存しない方がいいかなぁ。いずれ修正してみます。

とりあえず次回はユーザの行動履歴を加えたバージョンに続きます。たぶん。

追記(2010/4/13)

コメント欄にてご指摘いただき、ソースコードを一部修正しました。ありがとうございます!>しましまさん

修正内容は以下の通りです。

  • ユーザ、アイテムのバイアス項の設定部分の修正
  • 元のコードでは各バイアス項(全体の平均、ユーザごとの評価のゆらぎ、アイテムごとの評価のゆらぎ)を初めに設定してその値をそのまま使っていたのを、行列の更新時に同時にバイアス項も更新させる手法を追加

変更後の各手法の結果は以下の通りです。

手法 RMSE クラスタ 反復数 eta lambda
勾配降下法 1.012 10 100 0.001 0.001
確率的勾配降下法 0.967 10 100 0.1 0.1
確率的勾配降下法(バイアス項つき・更新なし) 0.979 10 100 0.1 0.1
確率的勾配降下法(バイアス項つき・更新あり) 0.962 10 200 0.1 0.1

バイアス項を各ループで更新するバージョンが一番結果がよくなりました!論文どおりでほっとしました^^; もっとしっかりサーベイしないとダメですねぇ…。

Bayesian Setsの特許について

別にブログに書いてもしょうがないかなーと思っていたのですが、同じような目に遭う方がいるかもしれないのでちょろっとだけ書いておきます。

先日Stupaという関連文書検索システムを公開したのですが、その中で使用していたBayesian Setsというアルゴリズムが既に特許を取得されているため、公開を停止してほしいってメールが来ました。以前に公開したBayesian SetsのCPANモジュールAlgorithm::BayesianSetsも同様に下ろしてほしいとのことでした。特許の内容は以下のページに書いてあります。

特許の出願者がBayesian Setsの論文の著者と大学の機関のようなので、おそらく論文発表の前に出願したのではないかと思います。請求項の内容などをすべて詳細に読んだわけではないのですが、やはりBayesian Setsのアルゴリズムについての特許のようです。ちなみに今回メールを送ってきたのは出願者ではなく、出願者からライセンスを受けたという会社でした。

とりあえずソフトウェアのライセンス項目に特許に関する注意を加えたのですが、その後またメールが来て、特許が取得されたアルゴリズムを実装する事自体が特許取得者の許可がいるとのことで、やはり削除してほしいとのことでした。あまり揉めるのも嫌なので、結局StupaからBayesian Setsに関するコードをすべて除いて、リリースしなおしました。検索部分はIDF+内積、IDF+Cosine距離に変更しましたが、Bayesian Setsと比較してそこまで精度が落ちる分けでもないので、まあ問題はないかなと。

また日本ではまだ公表のみで特許の取得には至っていないようですが、特許が権利化された場合は出願日まで遡って権利が発生するらしいので、後々のことを考えるとやはり使わない方がよさそうです。そもそもGoogle Code上で公開しているので、日本の特許だけで閉じて考えられるのか非常に怪しいところですが。

特許についての知識は全くないので、いきなり今回のメールが来たときはかなりビックリしました。公開してあまり経たずにこんなことが起きたので、モチベーションもガタ落ちしたり…。地雷踏んだなぁという感じなのですが、こういうことってよくあるものなんですかね?企業でがっつり取り組む気なら、あらかじめ知財担当にお願いして特許を調べてもらうこともできると思いますが、個人でオープンソースを作っている場合とかでは特許を調べ上げるのは結構大変そうな気がします。日本の特許だけじゃなくて海外の特許までが影響を受けるかもしれないとなると、調べるのにかなりの労力が必要になりそうですし。今回のような場合ならまずは、論文の著者に特許を取得するのかどうか聞くのがいいのかな?

結論は特にないのですが、すぐに頭を切り替えて別の対応策を取れるように打たれ強くなりたいな、と思う出来事でした。本当はもっと事前対策を打てることが一番なのですが^^;

NMF(Non-negative Matrix Factorization)を試してみた

先週は言語処理学会の全国大会に参加してきたのですが、チュートリアル「推薦システム -機械学習の視点から-」で紹介されていた行列分解に興味が湧いたので実装しようと奮闘中です。とりあえず行列いじりの練習に、手元の本に説明があるNMF(Non-negative matrix factorization)を実装してみました。参考にした書籍は「Rで学ぶクラスタ解析」と「集合知プログラミング」です。

NMFはn×mの行列をn×k, k×mの2行列の積に分解する手法で、分解後の行列はクラスタリング結果としても利用できます。例えば文書の特徴をbag-of-wordsで表した文書-単語行列を分解すると、文書の各クラスタへの重要度を表す行列と、特徴語の各クラスタへの重要度を表す行列が出力されます。より詳しくは前述の書籍をご参照下さい。

作成したコードは以下の場所に置いてあります(汎用性のないコードだなぁ…)。行列演算部分はEigenを使ってみました。

ちなみに入力データはSparseMatrixに入れているのですが、分解後の行列やtemporaryの行列をDenseのMatrixにしているせいか、大きめのデータ(5万文書くらい)を入れると落ちてしまいます。メモリ不足が原因かな…?

実際に実行する時は以下のように(Eigenを事前にインストールしておく必要があります)。入力データは1行に1文書の特徴を持つタブ区切りテキストです。入力データは6x9の文書-特徴語行列なのですが、NMFでは1文書のデータを列ベクトルで表現するのが普通らしいので、入力を読み込んだ後に転置して9x6の特徴語-文書行列として処理しています。

% g++ -Wall -O3 nmf.cc -o nmf

% cat input.tsv         # 6x9の行列(6文書, 9種類の特徴)
阿佐田   J-POP       10   J-R&B       6   ロック  4
小島     ジャズ       8   レゲエ      9
古川     クラシック   4   ワールド    4
田村     ジャズ       9   メタル      2   レゲエ  6
青柳     J-POP        4   ロック      3   HIPHOP  3
三輪     クラシック   8   ロック      1

% ./nmf input.tsv 3 50  # 9x3, 3x6 の行列に分解
Reading input data
Factorizing input matrix
 loop: 10
 loop: 20
 loop: 30
 loop: 40
 loop: 50
Input matrix was factorized. ( V = W * H )
=== W matrix ===
J-POP   0.0004  0.4423  0.0000
J-R&B   0.0000  0.2320  0.0000
ロック  0.0521  0.1962  0.0000
ジャズ  0.0000  0.0001  0.5020
レゲエ  0.0000  0.0000  0.4493
クラシック      0.5053  0.0000  0.0000
ワールド        0.0051  0.0017  0.0010
メタル  0.0036  0.0047  0.0013
HIPHOP  0.0059  0.0076  0.0001

=== H matrix (transposed) ===
阿佐田  0.0000  0.5881  0.0000
小島    0.0000  0.0000  0.4753
古川    0.2006  0.0000  0.0002
田村    0.0006  0.0008  0.4256
青柳    0.0088  0.2118  0.0000
三輪    0.4024  0.0000  0.0000

ここではクラス多数を3にして行列の分解を行いました。分解後の1つ目の行列は特徴語の各クラスタでの重要度を表し、2つ目の行列は文書の各クラスタでの重要度を表します。一番重要度の高いクラスタがその文書が属しているクラスタと考えられるので、各クラスタは「阿佐田、青柳」「小島、田村」「古川、三輪」のように分かれていることが分かります。

ちなみにbayonでクラスタリングするとこんな感じで、クラスタの構成は同じになりますね。

% bayon -n 3 input.tsv
1       田村    小島
2       阿佐田  青柳
3       三輪    古川

まだEigenの使い方がよく分かってないとか、そもそも線形代数が苦手でどうしようという感じなのですが、とりあえずウォーミングアップは済んだので、次はもうちょい突っ込んだもの作りたいですね。

Twitter Streaming APIでデータ収集

Twitterからデータを引っ張ってきたいと前から思ってたので、TwitterStreaming APIを試し中。とりあえず1日分(2010/02/10 12:00 〜 2010/02/11 12:00)のデータを引っ張ってきてみました。ドキュメントはほとんど読んでないままやってるので、いろいろ間違ってるかも。

実際に引っ張ってくるコードはこんな感じ。ユーザ名、TweetのID、日付、Tweetの文面をタブ区切りで出力します。Config::Pitについてはここを参照。

#!/usr/bin/perl

use strict;
use warnings;
use AnyEvent::Twitter::Stream;
use Config::Pit;
use Data::Dumper;
use Encode qw(encode);

my $config = pit_get('twitter.com');
my $method = 'sample';

my $cv = AnyEvent->condvar;
my $listener = AnyEvent::Twitter::Stream->new(
    username => $config->{username},
    password => $config->{password},
    method   => $method,
    #track    => 'google',
    on_tweet => sub {
        my $tweet = shift;
        eval {
            if ($tweet->{id} && $tweet->{user}{screen_name}
                && $tweet->{text} && $tweet->{created_at}) {
                $tweet->{text} = encode 'utf-8', $tweet->{text}
                    if utf8::is_utf8($tweet->{text});
                $tweet->{text} =~ s/\n/ /g;
                printf "%s\t%d\t%s\t%s",
                    $tweet->{user}{screen_name}, $tweet->{id},
                    $tweet->{created_at}, $tweet->{text};
                print "\n";
            }
        };
        if ($@) {
            warn Dumper $@;
#            exit 1;
        }
    },
    on_keepalive => sub {
        warn "ping\n";
    },
    on_eof => $cv,
);
$cv->recv;

1日分のTweetをカウントしてみると下の結果になった。あと各Tweetから何かしら特徴を求めたいので、とりあえず日本語だけに限定して、カタカナ・平仮名・漢字のいずれを含む文章のうち、MeCabで一般名詞か固有名詞がとれるものだけを抽出してみた。中国語も入っちゃってそうだけど、数は少ないと思うので無視で。

- Tweet ファイルサイズ
Tweet 1718239件 225MB
日本語を含むTweet 53537件 9.7MB
日本語で名詞を含むTweet 53195件 9.7MB

Twitterの投稿量が1日で170万件だけってわけないと思うので、かなりの量が省かれちゃってるのかな?実際はどれくらいあるんだろう。

あとオマケで1時間おきのTweetの数をカウントしてみた。

時(日本時間) Tweet
0 88150
1 88423
2 84763
3 82074
4 81279
5 79663
6 81182
7 82106
8 81632
9 83521
10 86534
11 85347
12 79614
13 68669
14 59306
15 49191
16 46703
17 45748
18 46131
19 48526
20 53187
21 62295
22 72542
23 81651
  • 日本語で名詞を含むもの
時(日本時間) Tweet
0 3856
1 2876
2 1856
3 1100
4 705
5 610
6 624
7 930
8 1153
9 1619
10 1972
11 2269
12 2506
13 2286
14 2195
15 2069
16 2310
17 2508
18 2666
19 2666
20 3051
21 3419
22 3830
23 4119

さーて次は実際にTwitterのデータを使って何かつくろう。関連Twitter検索作ろうかな。

bayonを使って画像からbag-of-keypointsを求める

クラスタリングツールbayonOpenCVを使って、画像からbag-of-keypointsを特徴量として抽出する手順について書きたいと思います。bag-of-keypointsは自然言語処理でよく使用されるbag-of-words(文章を単語の集合で表現したもの)と同じようなもので、画像中の局所的な特徴量(keypoint)の集合で画像の特徴を表します。bag-of-wordsと同じ形式ですので言語処理と同じように、bag-of-keypointsデータを使ってクラスタリングツールに適用したり、転置インデックスに載せたりといったことが可能になります。

今回は画像からbag-of-keypointsを取り出し、そのデータを使ってbayonで画像集合をクラスタリングするところまでやってみます。ちなみに画像処理は完全に素人で、この記事もニワカ知識で書いているので、間違っている箇所やもっと効率のいい手法があるかもしれません。その際は記事にコメントいただければ幸いですm(_ _)m

bag-of-keypoint抽出手順の概要

  1. 各画像から局所特徴量を抽出する
  2. 局所特徴量をクラスタリングし、クラスタの中心をvisual words(代表的なパターン)とする
  3. 各局所特徴量に対し一番近いvisual wordを求め、画像全体の特徴量をvisual wordsのヒストグラム(bag-of-keypoints)で表現する

画像の局所特徴量はSIFTを使います。SIFTでは画像中の特徴的な点(keypoint)をいくつか抽出します。各keypointは128次元の特徴ベクトルとなり、抽出されるkeypointの数は画像によって異なります。このようにkeypointごとに特徴ベクトルを持つため、画像全体の特徴ベクトルを表すにはそのまま使用するのは困難です。

そこで局所特徴量をクラスタリングして各クラスタの中心ベクトルを特徴的なパターン(ここではVisual Wordsと呼ぶ)とし、各局所特徴量は最も近いVisual Wordで置きかえることで、画像全体の特徴をVisual Wordの頻度(ヒストグラム)で表現することができます。このVisual wordのヒストグラムをbag-of-keypointと呼びます。

SIFTやbag-of-keypointsについては、以下の資料が非常に分かりやすかったので、より詳しく知りたい方はそちらをご参照下さい。

必要なもの

OpenCVのバージョンが2.0以降の場合は、SIFT Feature Detectorをコンパイルする時は以下のページを参考にしてソースコードを修正しておいてください。

画像データはとりあえずcaltech101を使用しますが、別で用意できるのであれば何でも構いません。

1.局所特徴量を抽出してbayonの入力を作成

まずはSIFT Feature Detectorを使って各画像からSIFTを特定し、bayonに適用する入力ファイルを作成します。以下のように実行してください。"/path/sift/bin/siftfeat"はSIFT Feature Detectorをビルドした時に生成される実行プログラムのパスです。

% mkdir data
% find /path/caltech101 -name "*.jpg" | ./save_sift.pl /path/sift/bin/siftfeat data/siftmap.tsv > data/input.tsv

SIFTの各座標を1つのドキュメントとして出力します。このとき各座標に固有のIDを割り振るので、各画像中のSIFT特徴量がどのIDに採番されたかを data/siftmap.tsv に保存しておきます。

2. 局所特徴量をクラスタリングしてVisual Wordsを特定

次に上記のSIFT特徴量をクラスタリングします。得られたクラスタの中心ベクトルがVisual wordsになります。

ただし実行環境によっては、すべてのSIFT特徴量をクラスタリングするにはメモリが足りない場合があります。そこでランダムにいくつかの座標を抜き出して、その座標のみでクラスタリングを行うようにします。下記では入力ファイルからランダムに10分の1行選択してます。本来なら各カテゴリから均等に選択するようにした方がよさそうですが、今回は全体からランダムに抽出します。

% ./rand_line.pl data/input.tsv 10 > data/input_10.tsv

では次にクラスタリングを行います。bayonを使って以下のように実行してください。

% bayon -l 1.5 -c data/centroid.tsv --clvector-size 128 data/input_10.tsv > /dev/null

クラスタの中心ベクトル(これが visual words になります)を data/centroid.tsv に保存しておきます。その際デフォルトでは50次元までしかベクトルの要素を保存しないので、中心ベクトルのサイズをSIFTと同じ128次元に指定しておく必要があります。またクラスタリング結果そのものは特に使わないので、/dev/nullに流し込んで廃棄します。

各SIFT特徴量に最も近いvisual wordを特定しbag-of-keypointsを作成

次に各SIFT特徴量に対し、一番類似しているvisual wordを特定します。これはbayonの-Cオプションでさくっとできます。一番近いものだけ分かればいいので、--classify-size=1 を指定して類似するベクトルの出力を1つだけにしておきます。結構時間がかかるので気長に待ちます。

% bayon -C data/centroid.tsv --classify-size 1 data/input.tsv > data/classify.tsv

最後に各画像の特徴量をvisual wordsのヒストグラム(bag-of-keypoints)で表現します。

% ./assign_vwords.pl data/siftmap.tsv data/classify.tsv > data/bok.tsv

得られたbag-of-keypointsは以下のようになります。

% cat data/bok.tsv
/path/caltech101/101_ObjectCategories/accordion/image_0001.jpg 20  1   39  1   59  1   77  1   79  1   120 6   122 1   144 2   167 1   200 1
...

先頭が画像ファイル名で、それ以降は "visualword_id1 \t count1 \t visualword_id2 \t count2 ..." となります。visualword_id部分が例えば「りんご」「大根」のような単語に変えてみれば、これはbag-of-wordsとまったく同じ表現であることが分かります。

ここでちょっと思ったのですが、bag-of-keypointsは別に一番似ているものの頻度じゃなくて、似ているもの上位5個ぐらいの類似度(0〜1)を合わせたベクトルでもいいんじゃないかな...?まあとりあえず今回は文献のやり方に従います。

画像集合をクラスタリング

作成したbag-of-keypointsを使って、画像集合をbayonでクラスタリングしてみます。caltechは総数が9143個ですが、とりあえず500個のクラスタに分割します。またbayonを実行する際、--idfオプションをつけて特徴ベクトルに補正をかけることにします。

% bayon -n 500 --idf data/bok.tsv > data/cluster.tsv

少しクラスタリング結果を眺めてみた印象だと、人の顔はよく解析できているけど、ボロボロのクラスタが結構な割合であるので、もう少し改善が必要な感じでした。bayonの精度がいまいちなのか、bag-of-keypointsの結果がいまいちなのかをちゃんと検証できてないのがダメなのですが、まあとりあえずある程度は使えるかなと。

なかなかうまくできたクラスタをいくつか載せておきます。


まとめ

bayonを使って画像からbag-of-keypointsを抽出しました。また抽出結果を使用してbayonで画像集合をクラスタリングし、ある程度の精度が得られていることを確認しました。bag-of-wordsと同じ形式ですので、他にも既存の検索エンジンや分類器に入れたりといったこともできると思います。ゆくゆくは類似画像検索が作れればいいなーと妄想してます。

また今回は局所特徴量にSIFTを使用したのですが、SIFTは特許が取られているので商用で使う場合はちょっと厳しそうです。そこでSIFTと同じような局所特徴量としてSURFというアルゴリズムがあるようですので、現在これを使ってbag-of-keypointsを求めるよう実装中です。SURFはOpenCVにも入っているので比較的簡単に使えます。

あと今回はコマンドラインから一つ一つ呼び出すようにしてましたが、実際に使うときは面倒なので、もう少しまとめた形で実装し直して公開したいですね。

追記:SURFも特許がとられてたいたため、商用利用は許可がいるようです。

結局のところ商用で局所特徴量は何を使うのがいいんですかね?知識が足りなさ過ぎる…。