Hatena::ブログ(Diary)

計算論的メンタルレキシコン このページをアンテナに追加 RSSフィード

2011-07-17

行列(画像)分解アルゴリズムGoDec (Zhou+, ICML2011)の実装を公開しました

つい2週間ほど前,機械学習のトップカンファレンスICMLが開催されました.その中のGoDecという行列分解アルゴリズムを実装したので公開します.このアルゴリズムは,簡単にいえば「外れ値抜き特異値分解」で,昨日のICML読み会で発表しました.論文はこれです.

厳密な版(Nai:ve GoDec)は遅いですが実装は非常に簡単です.遅い版でも,数百x数百ピクセルの小さな画像であれば,十分実用的な速度で動くので,実装して試してみた次第です.GoDecの論文では,この厳密な版(Nai:ve GoDec)が線形収束することを証明した上で,さらに,実用的に早くなるように(証明はないようですが)乱択化低ランク近似で高速化したFast GoDecも提案されています.僕のソースコードはNew BSDライセンスにしておきましたが,GoDec自体の特許関係は全く調べておりません.もしかしたら,特許が取られているかもしれませんので,商用にGoDecを用いることには注意が必要です.

ソースコードと使い方

コンパイルには,Eigenが必要です.C++0xを使っているので,gcc 4.6.0以上奨励です.

git clone git@github.com:niam/godec.git
cd godec
./waf configure
./waf build
build/src/godec -r 8 -k 762 hogehoge.png

これは,簡単にいえば,「762個までの外れ値を許容して,hogehoge.pngを行列だと思ってランク8で低ランク近似しろ」ということです.

例えば,以下の画像をhogehoge.pngとして入力したとします.

f:id:niam:20110717191416p:image

すると,X.png, L,png, S.png, LpS.pngがカレントディレクトリに出来ます.X.pngは,元の画像をグレイスケールに変換したもので,これが,分解対象の行列になります.

f:id:niam:20110717191415p:image

これをランク8の行列Lと,要素数762の行列Sに分解し,L+Sを実行したものがLpSです.

f:id:niam:20110717191413p:image

GoDecは,特異値分解に加えて外れ値の数kの自由度を持つので,特異値分解と比較するためには,自由度の数を合わせないとフェアではありません.200x180行列において,1ランク分の情報を保存するには200次元の特異ベクトルu+180次元の特異ベクトルv+特異値σの381個の自由度が必要になります.従って,2ランク分の情報がちょうど762自由度なので,単純な特異値分解で上位8+2=10で低ランク近似した場合と自由度がおなじになります.元の画像を特異値分解で上位10で低ランク近似したものが,下の画像です.

f:id:niam:20110717191414p:image

GoDecでは,頭の上の木漏れ日による光の点が外れ値として処理されハッキリ表示できているのに対し,特異値分解(SVD)では周りの色に紛れてしまっていることが分かります.

SVDと比較する場合,この自由度を合わせる作業は面倒くさいので,自動的に自由度を合わせてくれるモードも作りました.

build/src/godec -r 10 --rdiff=2 --comparesvdmode hogehoge.png

とやると,上の実験設定と同じ実験設定になります.つまり,特異値分解の方のランクが10で,それと2ランク差がつくようにGoDecのランクと,許容する外れ値の数kを自動的に設定してくれます.comparesvdmodeでは,10ランクで近似したsvd.pngも作られます.

アルゴリズムの解説

行列を,小さいランクの行列で近似することを低ランク近似と言います.低ランク近似の代表的な方法が特異値分解です.しかし,特異値分解は,行列内に,外れ値のような他の要素よりずっと大きい/小さい要素があったりすると,その要素を上手く再現出来なかったり,また,その要素に引きずられて,その要素の周りがにじむように近似誤差が大きくなってしまったりすることが知られています.

このアルゴリズムは,簡単に言えば「外れ値抜きの特異値分解」です.

  1. 外れ値を抜いて特異値分解でランクrで低ランク近似する
  2. もとの行列と低ランク近似した行列を見比べて,要素の差(の絶対値)が大きい「外れ値」要素を差の大きい方からk個見つけてきて新しい外れ値にする

を繰り返したものが,線形収束することを示しているのがこの論文です.

僕の発表資料はこちらです.解釈に間違いがある可能性がありますので,詳細は論文本体を参照されることをお勧めします.

2011-07-02

C++0xで10行で書けるマージンパーセプトロン

C++0xの練習として,マージンパーセプトロンを書きました.データ読み込みやコマンドライン処理を含めてもperceptron_binary.cppという90行にみたない1ファイルだけで済ませてあり,一応,コピペなどで気軽に使えることを意識しています.

ソースコード使い方

上の「ソースコード」や「使い方」のリンク先にもかいてありますが,gcc4.6さえ入っていれば,下記のように打てば使えます.基本的には,perceptron_binary.cppをコンパイルするだけです.gcc4.6が必要なのは,range-based forを使っているからです.

% g++ -std=c++0x -O2 perceptron_binary.cpp -o perceptron_binary
% g++ -O2 mknews20.cpp -o mknews20
% wget http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/news20.binary.bz2
% bzip2 -dc news20.binary.bz2|./mknews20
% ./perceptron_binary news20.train news20.test 10 0.6
Accuracy: 4819/4996=0.964572
Zerofication rate: 334121/1355194=0.246548

10行,というのはこちらの記事にインスパイアされました.

直接の参考にしたのは,kisa12012さんのこの記事や輪講発表資料です.

最初,パーセプトロンを実装していたところ,id:sleepy_yoshi さんのTokyoNLP #5で発表された資料でマージンパーセプトロンも簡単に出来ることを知り,ついでにやってみたという感じです.


プログラムの解説

主要なマージンパーセプトロンの部分は,以下の10行です.

struct perceptron_t{
  weight_t weight;  val_t    margin;
  perceptron_t(int size, val_t margin_=0.0): weight(size),  margin(margin_)   {}
  inline clabel_t classify(const datum_t& datum) const{ return weight*datum>0.0?1:-1; }
  inline void     train   (const datum_t& datum){
    if(((val_t)datum.first)*(weight*datum)<=margin)
      for(const auto& p : datum.second)
        weight[p.first] += ((val_t)datum.first)*p.second;
  }
};

10行で書ける…というのは,class perceptron_tが10行ということです.すみません,さすがにc++でデータ読み込みまで含めて10行というのは無理です(汗 仕掛けはweight*datumにあります.weight_t型のoperator*をオーバーロードすることによって,重みベクトルと特徴量ベクトル内積をweight*datumで書けるようにしました.id:sleepy_yoshi さんや id:kisa12012 さんの資料に出てくるパーセプトロンのifをほとんどそのままの形で表せています.

数式:

  • if {y_t {¥bf w}_^¥top {¥bf x}_t ¥le ¥gamma}
    • {{¥bf w}_{i+1} ¥leftarrow {¥bf w}_i + y_t {¥bf x}_t}

プログラム

    if(((val_t)datum.first)*(weight*datum)<=margin)
      for(const auto& p : datum.second)
        weight[p.first] += ((val_t)datum.first)*p.second;

一方,更新式はプログラムっぽい形になってしまいました.これも,vectorのoperator +とoperator =をオーバーロードすれば数式っぽくかけると思うのですが,ナイーブに書くと途中でvectorのコピーが発生したり,データ中に現れる特徴だけではなくて,全ての特徴に対してforがかかったりしそうで,パフォーマンスが低下しそうです.プログラムで書けば2行ですむものにそこまで労力をかける気にもなりませんでした.


データフォーマットは,libsvmやliblinear,ollその他,みんなが使っている,1行1データで各行が"ラベル 特徴量ID:特徴量 特徴量ID:特徴量 ... 特徴量ID:特徴量"という形式です.

datum_tはpair<ラベル,vector<pair<特徴量ID, 値>>>で,trainのp.firstには特徴量IDが,p.secondには特徴量が入っています.datum.firstがラベルのy_iで,+1か-1を取ります.

weight_tはvector<val_t>の派生クラスで,operator*はweight_tの中で定義されています.

C++0xのRange-based forをfor(const auto& p : datum.second)で使ってますが,これ,この部分だけなら普通のC++でも書けますね.for(data_t::iterator it = datum.second.begin(); it!=datum.second.end();++it)みたいに.なので,分かりにくくなりますが,一応C++でも10行で書けます.C++0xのもうひとつの本領はラムダ式で,制度を求めるときに使ってます(後述).



C++0xの機能

Range-based for

今までC++でループを回すためには,for(int i=0;i<n++i)とintで回すか,for(c_t::iterator it=c.begin();it!=c.end;++it)とイテレータで回すかしかありませんでした.一方,今は他の多くの言語では,例えばPythonのfor(x in l)といった形のループが用意されています.C++0xでも,とうとうこの形でも書けるようになりました.

auto

先程のfor(c_t::iteratort it...のところのc::iteratorを書くのは,cがvector<pair<int, vector<pair<int,double> > >のように長いと泣けてきたりします.今まで,これを防ぐ代表的な方法は,typedefで短い型名c_tを付け直してやるぐらいしかありませんでした.これが,for(auto& it=c.begin()... )というように簡潔に書けます.

ラムダ式

今までも,operator()をオーバーロードすることによる関数オブジェクトは定義すれば作ることが出来ましたが,いちいち定義しなければならず,その場で作ることができませんでした.一方,Pythonでは,lambda x:x**2のように,その場で関数を作成する事ができます.C++でも,その場で関数オブジェクトを作成することができるようになりました.

今回は,集計のところで本領を発揮してくれています.

  printrate("Accuracy: ",  count_if(test_data.begin(), test_data.end(), 
    [&p](const datum_t& t){return p.classify(t)==t.first;}),  test_data.size());
  printrate("Zerofication rate: ", count_if(p.weight.begin(), p.weight.end(), 
    [](const val_t& w_i ){return (bool)(w_i==0.0);}) , (int)maxfid+1);

pがperceptron_t型で,パーセプトロンそのものです.printrate(文字列,int a, int b)は,文字列を表示した後,a,bをfloatに直して率を表示するだけの関数です.

ラムダ式はtest_dataに対するcount_ifの中で使っています.count_ifは,指定されたコンテナ(test_data)中の各要素について,渡されたpredicateを実行し,trueが返って時だけカウントするというものです.

[&p](const datum_t& t){return p.classify(t)==t.first;}という形をしています.t.firstには,正解が入っています.[&p]は現在のスコープからpという変数を使用できるという指定です.


マージンパーセプトロンについて

マージンパーセプトロンは,ものすごく乱暴にいってしまうとこういうものです:

  • 普通のパーセプトロンは「予測を間違えた時だけ更新する」という方針です.上のコードで,margin=0.0の場合がこれに当たります.
  • マージンパーセプトロンは,「予測があっていた時もスレスレだったら更新する」という方針です.そのために,スレスレ度のような「マージン」という量を定義します.そして,このマージンが一定以下だったら(=もう少しで間違えそう!,という時の「少し」がその量以下だったら),更新を行います.上のコードでは,マージンはそのままmarginです.

実験設定

実験は,news20.binaryで行いました.これは,既にBag-of-words形式で特徴量ベクトルになっている新聞記事を2つに分類する問題です.news20.binaryは,@hillbig さん作のollというオンライン学習ライブラリマニュアルに書かれている実験で使われています.

この実験データの作り方は,id:n_shuyo さんも以前,書かれています.この実験でも,この記事と同じデータの作り方をしています.

が,困ったことに,行の順番をランダムに並び替えるsort -Rが自分の環境にはなかったので,この際news20.binaryから上記事のnews20.trainとnews20.testを作成するmknews20という簡単なプログラムを作成しました.使い方は簡単でコンパイルして,bzip2 -dc news20.binary.bz2|./mknews20で,news20.trainとnews20.testが作成できます.上のgithubのREADME.mdに書いてあります.

実験結果

まず,普通のパーセプトロンです.10はイテレーション回数で,0.0がマージンです.biasはデフォルトでは入りません.毎イテレーションでデータをrandom_shuffleしています.

$ ./perceptron_binary news20.train news20.test 10 0.0                                                                                                                                                                           Accuracy: 4768/4996=0.954363
Zerofication rate: 1025984/1355194=0.757075

Zerofication rate(ゼロ化率)は,重みベクトル中の0の要素の数/重みベクトルの次元(=素性の次元)です.精度はそこそこですが,高いゼロ化率を達成できています.

マージンをあげると,ゼロ化率は急激に下がりますが,精度は1%ほど上がります.

$ ./perceptron_binary news20.train news20.test 10 1.0                                                                                                                                                                           
Accuracy: 4820/4996=0.964772
Zerofication rate: 680830/1355194=0.502386

マージンを上げることによって,しばらく精度が上がるが,3.0ぐらいをピークにどんどんゆっくりと精度が下がっていく.一方,ゼロ化率はマージンを上げることで急激に下がっていく.

$ ./perceptron_binary news20.train news20.test 10 2.0                                                                                                                                                                           
Accuracy: 4825/4996=0.965773
Zerofication rate: 557335/1355194=0.411258
$ ./perceptron_binary news20.train news20.test 10 3.0                                                                                                                                                                           
Accuracy: 4827/4996=0.966173
Zerofication rate: 483274/1355194=0.356609
$ ./perceptron_binary news20.train news20.test 10 4.0                                                                                                                                                                           
Accuracy: 4825/4996=0.965773
Zerofication rate: 423173/1355194=0.31226
$ ./perceptron_binary news20.train news20.test 10 5.0                                                                                                                                                                           
Accuracy: 4824/4996=0.965572
Zerofication rate: 379253/1355194=0.279851
$ ./perceptron_binary news20.train news20.test 10 6.0                                                                                                                                                                           
Accuracy: 4819/4996=0.964572
Zerofication rate: 334121/1355194=0.246548
$ ./perceptron_binary news20.train news20.test 10 7.0                                                                                                                                                                           
Accuracy: 4818/4996=0.964372
Zerofication rate: 312279/1355194=0.230431
$ ./perceptron_binary news20.train news20.test 10 8.0                                                                                                                                                                           
Accuracy: 4817/4996=0.964171
Zerofication rate: 286297/1355194=0.211259

マージンパーセプトロンは,マージンがとても大きい場合は,ほとんど常に更新するので,マージンが十分に大きいところでは確率的勾配降下法(SGD)と同じ更新をする.ただし,今回は学習率は1.0で固定であるが,SGDは段々学習率をイテレーション回数に反比例する形で減衰させることが普通である.しかも,SGDでは重みベクトルを0に引き戻す作用をする正則化項を入れることも通常行われるが,今回はいれていない.

$ ./perceptron_binary news20.train news20.test 10 100.0                                                                                                                                                                         
Accuracy: 4472/4996=0.895116
Zerofication rate: 67953/1355194=0.0501426
$ ./perceptron_binary news20.train news20.test 10 1000.0                                                                                                                                                                        
Accuracy: 3958/4996=0.792234
Zerofication rate: 68017/1355194=0.0501899

失敗する場合だけでなく,正解した場合も常に重みベクトルを更新すると,このように,精度ががた落ちになり,しかも,ゼロ化率も低くなってしまうことが分かる.

結論

正則化項のないところでは,マージンを大きくすると精度は上がるが,ゼロ化率は下がる.正則化項のあるところでは,もう少しましな挙動をすると思われるので,この場合は将来の課題としたい.

2011-06-26

TokyoNLP #6で発表しました〜言語アフリカ起源説〜

6/25にTokyoNLP #6で発表しました。id:nokuno さんがまとめてくださっています。

発表の中身は、世界の音素の多様性が言語がアフリカから広がる時の連続創始者効果を表しているという、Scienceの論文 を読むというものです。

僕自身、「この言語とこの言語には関連があるんじゃなかろうか」とか思ったときは、確かにまず音素を見るので、「音素を考える」という方針は悪くないと思っています。というわけで、音素について話しだしたら、音声学入門みたいな話に…(汗

全く自分の専門ではない話を、しかも言語処理の勉強会でするということで不安だったのですが、興味をもっていただけた方も多いようで、ほっとしています。

時間がオーバー気味だったので、最後の線形回帰とBICのところは端折ってしまいました。まぁ、最初から説明不足な点も多いスライドだったので、ここは、後で補足するかも知れません。

追記:@takeda25さんがスライドの異音のページに間違いを見つけてくださいました。スライドをreplaceしましたが、6/26 23:00現在、まだr反映されていないようです。6/27 5:00に見たら反映されていました。

僕以外の発表の感想はこちら。

@uchumikさん CRFと素性テンプレート

自分の発表の直前なので、あんまり聞けなかった…(汗

CRFと素性テンプレートのお話。

素性テンプレートは、去年、自分もちょっと使ったことがあるが、ちゃんと説明している文献が少ないので、貴重であると思う…と、思ったが、高村本には、そんな話まで載っていたんだっけ(汗

高村本の5章には、素性テンプレートっぽい話が載っているように見えますが、素性テンプレートに関してだけ言えば、やはり、こちらのスライドの方がよく説明されていると思います。

@niam 音素と発音、線形回帰とBICの話〜言語アフリカ起源説(Atkinson, Science 2011)を読む〜

僕の発表。

@shuyoさん はじめての生成文法(後編)

生成文法の色々ある理論が凄く分かりやすく説明されていて、大変おもしろかったです。

結局、やっていることは、「いくつかの言語で文と非文を分けられる文法を構築する」こと。

そうして得られた文法処理能力を、

1)人間が持っていること

2)持っていることを仮定したとしても、さらに生得的に持っていること

については、生成文法は説明が乏しいな、と感じました。

一定の文法処理能力を人間が生得的に持っていることは確かだと思います。

しかし、人間の文法処理が、本当にそのように行われているのかは分からない→真のモデルはよく分からない。さらに、その真のモデルを人間が生得的に持っていることの証明に乏しい。

僕は、実際には脳科学とかで頭の中の処理を見てみないと分からないんじゃね?と思っています。

ただ、文法理論の一部は、計算機で扱い易い構造を持っていたり、語学教育の場で文法を教えるのにもよい性質を持っているので、有用性は高いと思います。

@y_shindohさん 音声認識向け単語N-gram言語モデリング入門 ?今日からあなたも言語モデラー?

音声認識系の話はあまり知らなかったので、おもしろかった…chasenとか、普通に使われているんですね。

@nokuno 統計機械翻訳入門 その2 ?フレーズベース編?

SMTの話。高次のIBMモデルが何をやっているかのまとめスライドが素晴らしい。