今日，製本用の修士論文を提出しました． ６年つづいた学生生活も終わりが近づいてきました．
この６年間にはいろんなことがありました．
前半の４年間は駒場へ通っていました． 2007 年 10 月からは数学科にいました． そこにはガチでやばい数学erがいて，そうでない人も段々頭を書き換えられて数学語ばかり喋るようになっていくやばい場所でした． 特別ガチでやばい二人と輪読した結び目理論の本は飛び抜けて難しかったですが，あれはワクワクしました． 僕は数学科生として結局は中途半端な感じになり，コンピューターを前からやってみたかったのもあって，情報理工に進むことにしました．
少し戻って，2007 年の駒場祭の帰り，友人と２人でゲームをつくる学内サークルを作りました(現ノンリニア)． はじめの頃に集まった十数人のメンバーに全力で支えられた１年間のことを，今でもよく覚えています． 思いつきと勢いで始めたサークルがこんなに大きくなって今でも回っていて，不思議な気分です． 先日の追い出しコンパ，とても楽しかったです！
後半の２年間は本郷へ通っていました． 情報理工に進み，今月まで中川研に所属しています． 最初の一年は機械学習の勉強でした． パーセプトロンや SVM すら聞いたことのない状態からはじめて，いつも頭の中がフワフワしている気分でした． 系列学習をやっていましたが，ちょっと理解がフワフワすぎで，うまくいきませんでした． 去年の２月ごろからは近傍探索のハッシュを勉強し始めました． こちらは割と性に合っていました． 機械学習に対する認識も少しずつこなれてきて，気分的には今やっとスタートラインに立てたかなってくらいの感じです．
インターンやバイトもしました (僕にしては珍しく積極的です)． おととしの夏には Google に通いました． Mozc チームは技術レベルが半端なく高くて，そこでの仕事は貴重な体験でした． 作ったものが OSS で公開されたのが嬉しかったです． そこで知り合った iwiwi くんに (しつこくお願いして) 紹介してもらって，直後の 11 月からは PFI でアルバイトをしています． こちらは密で濃くて尖っていて，とても楽しく刺激的です． 去年の３月には pixiv のインターンに参加しました． 僕が会ってこなかったタイプの人がたくさんいて，新しい刺激をたくさんもらいました．
今年に入って，初めての論文 (修論) を書き，初めて会議論文を投稿しました． 結果がどうなるかドキドキですが，ともかく大学生活はほぼ終わりました． 今やっと羽をのばしはじめたところ，あと一ヶ月も残っていません．小旅行的サムシングの予定が２つほどありますが，他の時間何しようかなあとぼんやり考えています．
４月からは社会人です． PFI で研究開発メンバーとしてがんばっていきます． よろしくお願いいたしますっ！
2011年12月24日土曜日
Eigenの実対称行列の固有値問題をARPACKで解く
ハイパーお久しぶりです． 修論に追われていて，先輩方に言われた「正月はもう一週間前に終わった」という言葉をかみしめています (とてもにがい)．
今日は，ARPACK を使って Eigen の実対称行列の固有値・固有ベクトルを求めるために薄いラッパーを書いたので，ご紹介いたします．
基本的な使い方は README に書いてあります． 求めたい固有値の数と種類 (実数 (or 絶対値) の意味で大きい (or 小さい) 方からいくつか，あるいは次数の意味で両側から合わせていくつか) を指定して，行列を渡すだけです． 行列は掛け算にだけ用いられるので，Eigen の Expression Template を直接渡すこともできます． また README の 2 番目の例のように，対象となる行列をベクトルにかける演算を自分で書いて渡せば，その固有値を求めることができます． これは複数の大規模疎行列の和や積の形をした行列をそのまま扱う場合などに有効です．
中身の話をすると，基本的には ARPACK の関数 dsaupd, dseupd, ssaupd, sseupd のラッパーです． これらの使い方についてはソースコード (dsaupd.f, ..., sseupd.f) の先頭に詳しく書かれていますが，数値計算マニュアル：線型代数 (3) にこれの翻訳と解説が書かれていて，わかりやすいです (と言っても ARPACK 自体は Fortran で書かれた古のコードなので，いろいろはまって苦労しました)． ARPACK の C++ ラッパーとしては公式に ARPACK++ というものがありますが，これがいかんせん古く，今のコンパイラーで通るように直すのが大変そうだったので，今回欲しかった実対称行列の固有値問題に絞ってラッパーを書いたものが arpaca です．
最近は，大規模固有値問題を解くためのライブラリーは ARPACK 以外にもあるようです． 特に有力そうなのは SLEPc ですが，これは MPI を用いた並列計算が前提のようで，本当に巨大な疎行列の固有値問題を解く場合には使えそうですが，お手軽さはなさそうです． 他にも anasazi は Krylov-Schur アルゴリズムや Davidson アルゴリズムのブロック化版などを実装しています． これは Trilinos という大きなパッケージの一部らしく，今回はあまり見ていません．
アルゴリズムの話をします． ARPACK (ARnoldi PACKage の略) が実装しているのは Implicitly Restarted Arnoldi (IRA) method という手法です． Arnoldi method は一般の行列をヘッセンベルグ行列に変換する手法ですが，実対称行列やエルミート行列に対しては対称三重対角行列となり効率が良く，特に Lanczos method の名前で知られています． Lanczos method は対象となる行列 A の Krylov 部分空間 {q, Aq, A2q, ...} を求めることで三重対角化を行いますが，implicit restart を行うことでこのうち必要となる固有ベクトルを含む低次元の部分空間を効率よく求めることができて，出てきた小さな三重対角行列の固有値分解は高速に求めることができる，というのが IRA の方法です． ここらへんの話は Matrix Computations や Fundamentals of Matrix Computations が詳しいです． まだ詳しく読んでいませんが，二つ目の本には IRA のさらに変形版である Krylov-Schur アルゴリズムや，中途半端な位置にある固有値も求めることができる Jacobi-Davidson method なども載っています．
時間があったらここらへん自分でも実装してみたいです． 実装の際には，Eigen::SelfAdjointEigenSolver の中に三重対角行列の固有値ソルバーが書かれていて，参考になりそうです． これは internal なので直接呼び出すのはアレかもしれませんが…… (外に露出させろというバグも提出されていますが，対応されるかどうか)． Eigen の疎行列モジュールは，3.1.0-alpha1 で今まで unsupported/SparseExtra にあったソルバーがたくさんメイン入りするなどしているので，将来的に大規模固有値問題もサポートされないかなあと期待はしています．
というわけで，メリクリさまでした！
2011年8月19日金曜日
Suffix Array を作る - SA-IS の実装
Suffix Array は今若者の間で人気のデータ構造です． マイ suffix array を実装することで，オシャレ度がアップしてモテ系になり，女子力も上がると言われています． その中でも今特に，手軽でクールな SA-IS (アルファベットサイズ固定の下で線形時間で省メモリで suffix array が作れる今最強のアルゴリズム) の実装がブームです． 僕もブームに便乗して，実装してみました．
ところで，SA-IS は流行っているので，日本語でもすでに様々なところで記事が書かれています (日付順)．
仕掛けが気になる人はまずここら辺を眺めて，え，そこどうなってんの，ってなったら下の論文を読むのが良いです．
SA-IS の元論文は 2009 年ですが，上の論文 (まだドラフト，PDF が G. Nong さんのページにある) に C のコードがついてるので，それを見ながらがわかりやすいです (というかソース読まないとわかんない)．
この論文についている実装ですが，'\0' を番兵に使うのでバイナリとかだと使えません． ちゃんとした実装は sais が良さそうで，上の G. Nong さんのページからも refer されています． esaxx - C++ enhanced suffix array template library - Google Project Hosting でもこれが使われています．
僕が書いた実装は
beam2d/sara - GitHub
にあります． 本体はヘッダー sara.h 一個です． 使い方は簡単で，文字列
上の論文にある C のコード (100 行弱) よりは長い (160 行強) ですが，番兵なしで，入力と出力がテンプレート化されてます (sais もそうですが)． 文字の型は組み込み整数型ならなんでもいいですが，値の範囲は 0 から指定した値 ch_max までに収まっている必要があります．
検証と評価用のサンプルが sara_main.cc (sara) です． 今回僕が唐突に SA-IS を実装し始めたのはレポートのためなんですが (冒頭の件？ 何の話ですか？)，その中では The Canterbury Corpus にあるデータを使って検証と実行時間の計測を行い，
アルゴリズムをひと通り説明します． が，正直がんばってわかりやすく説明するげんきが残っていないので，適当です！！！ わからなくなったら上のブログとか論文とかついてるソースコードとかを見るのが良いです．
各 suffix に対して，右隣の suffix より (辞書式順序で) 大きいか小さいかで L-type と S-type というラベルを振ります． L/S を並べると例えば mmiissiissiippii$ という文字列 ($ は終端文字で最小の文字で必ず S-type とする) なら LLSSLLSSLLSSLLLLS となります． 最終的な suffix array はまず先頭文字についてソートされているので，先頭文字ごとのバケツを作ります． suffix 同士の先頭文字が等しい場合，L-type の方が S-type より必ず先に来ることに注意します． つまり各バケツはさらに前半・後半で L-type / S-type にわかれます．
SA-IS ではまず L-type だけをソートします． 先頭文字と type が等しい場合，右隣の suffix の大小がそのままもとの suffix の大小になることを利用して，「まだ見ていない L-type の右隣の suffix の中で最小のものを含むソートされたリスト」を持っておき，最小のものを取り出し，左隣の suffix が L-type なら対応するバケツに入れて，リストに加える，という操作をリストが空になるまで繰り返します． このリストは，後述する初期値を除いて構築中の suffix array そのものが使えます (まだ埋まっていない部分は飛ばしながら見ます)． 次に S-type では逆向きに，今できたソートされた L-type suffix を大きい方から順に見て，左隣の suffix が S-type ならバケツの後端から順に加える，という操作を繰り返します．
これをやるためには，L-type をソートするときに「まだ見ていない L-type の右隣の suffix の中で最小のものを含むソートされたリスト」を常に持っている必要があります． これは，L の右隣にいる S (LMS: Left Most S-type という) 全部のソートされたリストをはじめに作り，上の方法で最小のものを取り出しては左隣の L-type を入れる，という操作を繰り返せばできます．
理由ですが，L-type が並んでいるところ (...LLLLLS... みたいな． この S が LMS) では右の suffix ほど小さくなります． 特にこれらの L-type の直後に来る LMS の suffix はこの L-type suffix のどれよりも小さいです． 同じように，S-type が並んでいるところでは右に行くほど大きくなります． つまり LMS は suffix の大小という意味で谷になっています． なので谷を全部持っておけば，必ず最小のものが含まれるし，取り出した最小の suffix の一個前 (山を後ろに一歩登ったところ) を加えていけば「まだ見ていない L-type の右隣の suffix の中で最小のものを含むリスト」が常に保たれることになります． S-type のソートでも同じように，今度は山のてっぺんから左に下っていく感じで進めていけばよく，これは単純にソート済みの L-type suffix とすでにソートされた S-type suffix を後ろから順に見ていけばよいです．
そのためにはじめにソートされた LMS-suffix のリストを作るのですが，ここが SA-IS 最大のポイントになります． 各 LMS 文字に対して，そこから次の LMS 文字までの substring を LMS-substring と呼びます (右端の LMS 文字を含みます)． LMS-substring は各文字と L/S-type とのペアに関する辞書式順序で比較します． LMS-substring をソートして，各 LMS-substring に対して LMS-substring 同士の大小と一致するような整数を振り，それを出現順に並べた文字列 S1 を考えます． すると，各 LMS-suffix の順序と対応する S1 の suffix の順序が一致します． なので S1 の suffix array を作れば LMS-suffix のソートされたリストが得られます． S1 の suffix array は，重複する文字がない場合にはただ文字をソートすればよいです (これは実は LMS-substring を整数に変換するところですでにやっています)． S1 に重複する文字がある場合には SA-IS を再帰的に適用します． LMS 文字は全体の半分以下しかないので，問題サイズは半分以下になっており，再帰させても線形時間で終了します．
まだ続きます． 上でサラっと LMS-substring のソートと書いた部分です． ここでも四段落上で書いたのとほとんど同じやり方が使えます． 今度は LMS-substring ごとの suffix を考え，それを全部集めたものをソートします． L/S-type はもとのと同じものが使えます． この場合，最初の「左隣が L-type な suffix のうち最小のものを含むソート済みリスト」の初期値として，LMS 文字一文字からなる suffix をバケツに入れたものが使えます． あとは上と同じように L-type，S-type の順でソートすれば，全 LMS-substring の suffix のソートが完了します． LMS-substring に対応する部分だけを取り出せば，LMS-substring のソートが得られます．
まとめると，
アルゴリズムの概要はこんな感じです． 何言ってるのかわからない部分もあると思いますが，上のリンクとか論文とかを参考にしてください． 本当は山と谷の図をたくさん書いて，具体的な文字列で手でアルゴリズムを回してみるとわかりやすいです． 論文は，ソースと本文を行ったり来たりしてればわかると思います． やってることはシンプルで，ソースも短くて (100 行弱)，コメントもたくさん入っています．
ところで，SA-IS は流行っているので，日本語でもすでに様々なところで記事が書かれています (日付順)．
- SAIS(Suffix Array - Induced Sorting) - EchizenBlog-Zwei
- SA-IS: SuffixArray線形構築 - sileの日記
- SA-IS - (iwi) ｛ 反省します - TopCoder部
- 接尾辞配列(Suffix Array)の効率的な実装手法について - xxxxxeeeeeの日記
仕掛けが気になる人はまずここら辺を眺めて，え，そこどうなってんの，ってなったら下の論文を読むのが良いです．
- G. Nong, S. Zhang and W. H. Chan, Two Efficient Algorithms for Linear Time Suffix Array Construction, IEEE Transactions on Computers, To Appear
SA-IS の元論文は 2009 年ですが，上の論文 (まだドラフト，PDF が G. Nong さんのページにある) に C のコードがついてるので，それを見ながらがわかりやすいです (というかソース読まないとわかんない)．
実装と実験
この論文についている実装ですが，'\0' を番兵に使うのでバイナリとかだと使えません． ちゃんとした実装は sais が良さそうで，上の G. Nong さんのページからも refer されています． esaxx - C++ enhanced suffix array template library - Google Project Hosting でもこれが使われています．
僕が書いた実装は
beam2d/sara - GitHub
にあります． 本体はヘッダー sara.h 一個です． 使い方は簡単で，文字列
str に対して sara::make(str, ch_max) で suffix array が返ってきます．上の論文にある C のコード (100 行弱) よりは長い (160 行強) ですが，番兵なしで，入力と出力がテンプレート化されてます (sais もそうですが)． 文字の型は組み込み整数型ならなんでもいいですが，値の範囲は 0 から指定した値 ch_max までに収まっている必要があります．
検証と評価用のサンプルが sara_main.cc (sara) です． 今回僕が唐突に SA-IS を実装し始めたのはレポートのためなんですが (冒頭の件？ 何の話ですか？)，その中では The Canterbury Corpus にあるデータを使って検証と実行時間の計測を行い，
std::lexicographical_compare と std::sort を使うだけの単純な suffix array 構築と比較しました． 詳細は省きますが，普通の文書だと 4 倍くらいの差しかなくて，割とクイックソート速いのねという感じです． 一方，似たパターンが連続するバイナリデータだと愚直な文字列比較ソートは激遅で，例えば第一世代 Core i5 2.67GHz で， ptt5 (FAX, 501KB) は std::sort だと 28.1 秒， sara::make だと 0.040 秒でした． 繰り返しを多く含む人工データではさらに差が広がります．適当な解説
アルゴリズムをひと通り説明します． が，正直がんばってわかりやすく説明するげんきが残っていないので，適当です！！！ わからなくなったら上のブログとか論文とかついてるソースコードとかを見るのが良いです．
各 suffix に対して，右隣の suffix より (辞書式順序で) 大きいか小さいかで L-type と S-type というラベルを振ります． L/S を並べると例えば mmiissiissiippii$ という文字列 ($ は終端文字で最小の文字で必ず S-type とする) なら LLSSLLSSLLSSLLLLS となります． 最終的な suffix array はまず先頭文字についてソートされているので，先頭文字ごとのバケツを作ります． suffix 同士の先頭文字が等しい場合，L-type の方が S-type より必ず先に来ることに注意します． つまり各バケツはさらに前半・後半で L-type / S-type にわかれます．
SA-IS ではまず L-type だけをソートします． 先頭文字と type が等しい場合，右隣の suffix の大小がそのままもとの suffix の大小になることを利用して，「まだ見ていない L-type の右隣の suffix の中で最小のものを含むソートされたリスト」を持っておき，最小のものを取り出し，左隣の suffix が L-type なら対応するバケツに入れて，リストに加える，という操作をリストが空になるまで繰り返します． このリストは，後述する初期値を除いて構築中の suffix array そのものが使えます (まだ埋まっていない部分は飛ばしながら見ます)． 次に S-type では逆向きに，今できたソートされた L-type suffix を大きい方から順に見て，左隣の suffix が S-type ならバケツの後端から順に加える，という操作を繰り返します．
これをやるためには，L-type をソートするときに「まだ見ていない L-type の右隣の suffix の中で最小のものを含むソートされたリスト」を常に持っている必要があります． これは，L の右隣にいる S (LMS: Left Most S-type という) 全部のソートされたリストをはじめに作り，上の方法で最小のものを取り出しては左隣の L-type を入れる，という操作を繰り返せばできます．
理由ですが，L-type が並んでいるところ (...LLLLLS... みたいな． この S が LMS) では右の suffix ほど小さくなります． 特にこれらの L-type の直後に来る LMS の suffix はこの L-type suffix のどれよりも小さいです． 同じように，S-type が並んでいるところでは右に行くほど大きくなります． つまり LMS は suffix の大小という意味で谷になっています． なので谷を全部持っておけば，必ず最小のものが含まれるし，取り出した最小の suffix の一個前 (山を後ろに一歩登ったところ) を加えていけば「まだ見ていない L-type の右隣の suffix の中で最小のものを含むリスト」が常に保たれることになります． S-type のソートでも同じように，今度は山のてっぺんから左に下っていく感じで進めていけばよく，これは単純にソート済みの L-type suffix とすでにソートされた S-type suffix を後ろから順に見ていけばよいです．
そのためにはじめにソートされた LMS-suffix のリストを作るのですが，ここが SA-IS 最大のポイントになります． 各 LMS 文字に対して，そこから次の LMS 文字までの substring を LMS-substring と呼びます (右端の LMS 文字を含みます)． LMS-substring は各文字と L/S-type とのペアに関する辞書式順序で比較します． LMS-substring をソートして，各 LMS-substring に対して LMS-substring 同士の大小と一致するような整数を振り，それを出現順に並べた文字列 S1 を考えます． すると，各 LMS-suffix の順序と対応する S1 の suffix の順序が一致します． なので S1 の suffix array を作れば LMS-suffix のソートされたリストが得られます． S1 の suffix array は，重複する文字がない場合にはただ文字をソートすればよいです (これは実は LMS-substring を整数に変換するところですでにやっています)． S1 に重複する文字がある場合には SA-IS を再帰的に適用します． LMS 文字は全体の半分以下しかないので，問題サイズは半分以下になっており，再帰させても線形時間で終了します．
まだ続きます． 上でサラっと LMS-substring のソートと書いた部分です． ここでも四段落上で書いたのとほとんど同じやり方が使えます． 今度は LMS-substring ごとの suffix を考え，それを全部集めたものをソートします． L/S-type はもとのと同じものが使えます． この場合，最初の「左隣が L-type な suffix のうち最小のものを含むソート済みリスト」の初期値として，LMS 文字一文字からなる suffix をバケツに入れたものが使えます． あとは上と同じように L-type，S-type の順でソートすれば，全 LMS-substring の suffix のソートが完了します． LMS-substring に対応する部分だけを取り出せば，LMS-substring のソートが得られます．
まとめると，
- 各 suffix に対して L/S-type のラベルを振る (ここは文字列を後ろから前へ一回走査するだけでできます)．
- LMS-substring をソートする． これはまず LMS 文字一文字からなる LMS-substring の suffix をバケツに入れ，そこから上記の方法で L-type，S-type の順にソートする．
- ソートされた LMS-substring に整数値を振る． 等しい LMS-substring 同士は同じ値になるように，小さい方から順に 0, 1, 2, ... と振る． これをもとの文字列での LMS の出現順と同じ順序で並べて S1 を作る．
- 整数値が被らなければ (= 重複する LMS-substring がなければ) すでに作った LMS-substring のソートが LMS-suffix のソートを与える． そうでなければ S1 に対して SA-IS を再帰的に適用して S1 の suffix array を作り，対応する LMS-suffix のソートを得る．
- LMS-suffix のソートされたものを対応する各バケツの後端に入れる． これを初期状態として，上記の方法で L-type，S-type の順にソートする．
アルゴリズムの概要はこんな感じです． 何言ってるのかわからない部分もあると思いますが，上のリンクとか論文とかを参考にしてください． 本当は山と谷の図をたくさん書いて，具体的な文字列で手でアルゴリズムを回してみるとわかりやすいです． 論文は，ソースと本文を行ったり来たりしてればわかると思います． やってることはシンプルで，ソースも短くて (100 行弱)，コメントもたくさん入っています．
2011年7月31日日曜日
正則な確率行列の逆行列
確率行列とは，行和が 1 な非負の正方行列のことです．確率行列 X に対して，Xi j を状態 i から状態 j への遷移確率と見なせば，X はマルコフ連鎖の遷移確率と見なすこともできます．
今正則な確率行列 X が与えられたとします．X の逆行列 X-1 は行和が 1 です．証明は簡単で，全成分が 1 である列ベクトルを 1 と書くことにすると，X の行和が 1 であることから 1 = X1 で，この式の両辺に左から X-1 をかければ X-11 = 1 となります．
これだけでは X-1 が確率行列とは言えません．一般に成分が負になる場合があるからです．逆行列が確率行列になるための条件は何でしょうか？ とりあえず考えつく方法として 3 つの解法があるようです．
一番簡単なのは，X と X-1 を遷移確率行列だと思うことです．X が非決定的な遷移確率行列 (i に対して複数の j で正な値を取ることがある) の場合，X で遷移したあと，X-1 で遷移すると，確率 1 でもとの場所に戻ってくる，などという奇妙なことは絶対に起こりません．つまり，逆行列が確率行列になるのは，X が決定的な遷移確率行列の場合だけです．行列の言葉で言い換えると X が置換行列ということになります．
正面から証明する方法もあります．X-1 が非負行列で，X のある行 i に対して 2 つの列 j, j' で正の値をとっていると仮定します．XX-1 は単位行列なので，i 以外のすべてのインデックス i' に対して，X の i 行目と X-1 の i' 列目は内積が 0 です．これらは共に非負ベクトルなので，X-1 の (j, i'), (j', i') 成分はともに 0 となります．つまり X-1 の j, j' 行目は第 i 成分を除いて 0 となり，この 2 行は一次従属になります．これは X-1 が正則であることに反します．よって X はどの行も 1 つの成分を除いてゼロであり，置換行列になります．
最後に，固有値からも証明できるようです．確率行列はスペクトル半径が 1 です．X-1 は X の固有値の逆数を固有値に持つので，これらが共に確率行列となるためには，両方ともすべての固有値が長さ 1 である必要があります．つまり X は等長変換を表す行列 (直交行列) ということです．第 i 成分だけが 1 の行ベクトル v を考えると，vX のユークリッドノルムが 1 である必要があります．つまり X の第 i 行はユークリッドノルムが 1 となりますが，行和が 1 な非負ベクトルでユークリッドノルムも 1 となるのは，ある一つの成分だけが 1 というベクトル，しかありません．よって X は置換行列となります．
とまあ，きっかけがあってそんなことをちょっと考えていました．
今正則な確率行列 X が与えられたとします．X の逆行列 X-1 は行和が 1 です．証明は簡単で，全成分が 1 である列ベクトルを 1 と書くことにすると，X の行和が 1 であることから 1 = X1 で，この式の両辺に左から X-1 をかければ X-11 = 1 となります．
これだけでは X-1 が確率行列とは言えません．一般に成分が負になる場合があるからです．逆行列が確率行列になるための条件は何でしょうか？ とりあえず考えつく方法として 3 つの解法があるようです．
一番簡単なのは，X と X-1 を遷移確率行列だと思うことです．X が非決定的な遷移確率行列 (i に対して複数の j で正な値を取ることがある) の場合，X で遷移したあと，X-1 で遷移すると，確率 1 でもとの場所に戻ってくる，などという奇妙なことは絶対に起こりません．つまり，逆行列が確率行列になるのは，X が決定的な遷移確率行列の場合だけです．行列の言葉で言い換えると X が置換行列ということになります．
正面から証明する方法もあります．X-1 が非負行列で，X のある行 i に対して 2 つの列 j, j' で正の値をとっていると仮定します．XX-1 は単位行列なので，i 以外のすべてのインデックス i' に対して，X の i 行目と X-1 の i' 列目は内積が 0 です．これらは共に非負ベクトルなので，X-1 の (j, i'), (j', i') 成分はともに 0 となります．つまり X-1 の j, j' 行目は第 i 成分を除いて 0 となり，この 2 行は一次従属になります．これは X-1 が正則であることに反します．よって X はどの行も 1 つの成分を除いてゼロであり，置換行列になります．
最後に，固有値からも証明できるようです．確率行列はスペクトル半径が 1 です．X-1 は X の固有値の逆数を固有値に持つので，これらが共に確率行列となるためには，両方ともすべての固有値が長さ 1 である必要があります．つまり X は等長変換を表す行列 (直交行列) ということです．第 i 成分だけが 1 の行ベクトル v を考えると，vX のユークリッドノルムが 1 である必要があります．つまり X の第 i 行はユークリッドノルムが 1 となりますが，行和が 1 な非負ベクトルでユークリッドノルムも 1 となるのは，ある一つの成分だけが 1 というベクトル，しかありません．よって X は置換行列となります．
とまあ，きっかけがあってそんなことをちょっと考えていました．
2011年7月2日土曜日
eLog - C++ で簡易ログと簡単時間計測
C++ で使える簡単なロガーを書きました．ロガーと書いてしまいましたが，そんな大げさなものではなくて，ほとんど std::cerr の代わりみたいなものです．見た目はほとんど google-glog です．
beam2d/elog - GitHub
ヘッダーだけなので，すぐ使えます．/usr/local/include へのインストールは ./waf install でできます．もちろん，自分のコードツリーにコピーしてきてもいいです．MIT License です．
使い方．
ファイル名と行数も出力します．そのうち日付時刻とかも入れたいです．
あと簡単なベンチマーク機能があります．アルゴリズムをパパっと実装して速度を見るのに便利です．こっちは iwiwi くんの benchmark を参考にしています．
ログの構文は，周りの邪魔にならない簡素なものがいいと思っていて，google-glog はその意味で好きですが，マクロです．マクロは名前空間代わりのプレフィックスを付けた方が親切ですが，ログの構文は簡素な方がいいので，上のようになっています．その代わり，そういうグローバルスコープの単純な名前は極力減らす方針で書いています．
あと google-glog はリンクがいるし，最初に Init なんたらを呼ばなきゃいけないし，gflags との関係がややこしいので，何も考えずに使えるものが欲しかったということもあります．
型ごとに verbosity を設定するのはほぼ思いつきなので，便利なのかどうかわからないです．
というわけで，大したものではないですが，ぜひ．
beam2d/elog - GitHub
ヘッダーだけなので，すぐ使えます．/usr/local/include へのインストールは ./waf install でできます．もちろん，自分のコードツリーにコピーしてきてもいいです．MIT License です．
使い方．
#include <elog/elog.h>
class ModuleA {};
class ModuleB {};
class MyLogger : public LOG::Logger { ... };
int main() {
LOG(INFO) << "基本は google-glog と同じ";
LOG(ERROR) << "LogLevel (severity) は INFO, WARN, ERROR, FATAL";
try {
LOG(FATAL) << "FATAL は";
} catch (const LOG::FatalLogError&) {
LOG(ERROR) << "例外を投げる";
}
LOG() << "INFO は省略できる";
// ロガーのレベルを変更できる
// LOG は括弧を付けなければ名前空間
LOG::SetDefaultLoggerLevel(LOG::ERROR);
LOG(INFO) << "出力されない";
LOG(WARN) << "出力されない";
LOG::SetDefaultLoggerLevel(LOG::INFO);
LOG(ModuleA, 0) << "型ごとに verbosity をつけてログを吐ける";
// 型ごとに異なる verbosity を設定できる (初期値は 0)
LOG::SetDefaultLoggerVerbosity<ModuleA>(2); // verbosity 2 以下を有効に
LOG(ModuleA, 2) << "出力される";
LOG(ModuleA, 3) << "出力されない";
LOG(ModuleB, 1) << "出力されない";
unsigned char uc = 65;
LOG() << "(un)signed char は整数として出力される: " << uc;
LOG() << "char は普通のストリームと同様: " << 'A';
// アサートもある
CHECK(true) << "出力されない";
try {
CHECK(false) << "出力され，そして";
} catch (const LOG::CheckError&) {
LOG(ERROR) << "例外を投げる";
}
// デフォルトのロガーは std::clog に吐くが，他のストリームに変えるには:
std::ostream& ostream = std::cout;
LOG::StreamLogger cout_logger(ostream);
LOG::SetLogger(cout_logger);
// この場合 SetDefaultLoggerLevel 等は使わない
cout_logger.set_level(LOG::WARN);
cout_logger.SetTypeVerbosity<ModuleA>(1);
// ロガーの実装は自分で作って差し替えられる (Logger の定義は logger.h)
MyLogger my_logger;
LOG::SetLogger(my_logger);
}ファイル名と行数も出力します．そのうち日付時刻とかも入れたいです．
あと簡単なベンチマーク機能があります．アルゴリズムをパパっと実装して速度を見るのに便利です．こっちは iwiwi くんの benchmark を参考にしています．
#include <elog/benchmark.h>
// BENCHMARK のログ出力はすべて INFO レベルです
int main() {
// タイトルを出してから経過時間を計測して，最後に出力
BENCHMARK(my_bench) {
for (int i = 0; i < 100000; ++i);
LOG() << "途中でも出力できる: " << my_bench.GetTime() << " sec";
for (int i = 0; i < 100000; ++i);
}
// 複数の計測結果をまとめる
LOG::BenchmarkSuite suite("title");
BENCHMARK(suite, case_1) {
for (int i = 0; i < 100000; ++i);
LOG() << "こっちも同じく: " << case_1.GetTime() << " sec";
for (int i = 0; i < 100000; ++i);
}
BENCHMARK(suite, case_2, "開始時にメッセージを出力できる") {
for (int i = 0; i < 200000; ++i);
}
// 計測結果を表にして出力
suite.LogChart();
}
とすると出力は
my_bench: ... [INFO] bench.cc(8): 途中でも出力できる: 0.000623941 sec my_bench: 0.00120592 sec case_1: ... [INFO] bench.cc(17): こっちも同じく: 0.000519037 sec case_1: 0.00106406 sec case_2: 開始時にメッセージを出力できる... case_2: 0.00103211 sec title | time (sec) --------+----------- case_1 | 0.00106 case_2 | 0.00103 --------+----------- (other) | 0.00006 total | 0.00216
ログの構文は，周りの邪魔にならない簡素なものがいいと思っていて，google-glog はその意味で好きですが，マクロです．マクロは名前空間代わりのプレフィックスを付けた方が親切ですが，ログの構文は簡素な方がいいので，上のようになっています．その代わり，そういうグローバルスコープの単純な名前は極力減らす方針で書いています．
あと google-glog はリンクがいるし，最初に Init なんたらを呼ばなきゃいけないし，gflags との関係がややこしいので，何も考えずに使えるものが欲しかったということもあります．
型ごとに verbosity を設定するのはほぼ思いつきなので，便利なのかどうかわからないです．
というわけで，大したものではないですが，ぜひ．
2011年6月24日金曜日
数理輪講で Anchor Graph Hashing の紹介をしてきました
今日、大学の数理情報学輪講で発表してきました。
今回は来週開催される ICML 2011 の論文 Hashing with Graphs の紹介です。Anchor Graph Hashing (AGH) という近傍探索のための最新のハッシング手法を提案しています (なんでアルゴリズム名をタイトルにしなかったんだろう……)。一応スライドを読めばなんとなくわかるんじゃないかと思います。
AGH は Spectral Hashing (SH) と同じ最適化問題を考えていて、それを解いてから out-of-sample extension の話に持っていくという流れも同じです。違うのは最適化問題を低ランク近似して解くところと、そのおかげで Nyström method という手法が使えて、SH みたいに分布に変な仮定を置いて固有関数問題をガリガリ解かなくて済むところ。データの次元、データ数のそれぞれについて線形時間のアルゴリズムになっています。
低ランク近似は、第一著者の Liu が ICML2010 提案している Anchor Graph に基づいています。この Anchor Graph が良い性質をたくさん持っていて、アンカーを用いたマルコフ連鎖的な解釈もあって、ハッシング以外にもいろいろ使えるんじゃないかと思います。
スライドでは論文の後半で述べられている階層的手法 (2-AGH) を紹介できませんでしたが、これは半分のビットは AGH で作って、残りは前半の各ビットに対して、うまく分けられなかった・分けるべきでなかった部分を補正するようなビットを入れる、という手法です。AGH は SH と同じく PCA みたいな手法なので、ビット数を増やすと固有値の低い (ノイズの多い) 主成分まで使うことになり、結果として MAP 精度が低下してしまう問題があります。階層的ハッシングでは半分のビットだけ AGH にして、後半はアドホックに作ることで、この問題をある程度回避しようとしています。が、半分は AGH で作るので、実験では 2-AGH でも下がったりしてます。
SH と同じ問題なので、多様体学習的な要素が入ってきています。なので与えられた距離関数についての最近傍探索ではないです。実際、実験結果では単純に全比較をした場合よりも精度が高いです。面白いのは、Spectral Embedding したあとに単純な全比較をする方法 (SE l_2 Scan) よりも階層的ハッシング (2-AGH) の方が精度が良いということ。近傍探索については、遠い点同士の類似度はざっくりゼロに落としちゃう方が精度が上がるということなんでしょうか。
Nyström method のあたりはちょっと難しいですが、SH に出てくる Laplace-Beltrami 作用素の固有関数問題のあたりに比べればだいぶわかりやすいと思います。[Bengio+, 04] をパラパラ読めば結構分かるんではないかと。固有ベクトルから固有関数を近似する方法として、結構使われている？みたいです。
ともかく Anchor Graph Hashing は面白いので、気になる方はぜひ論文も読むといいと思います。
rinko2011-agh
View more presentations from beam2d
今回は来週開催される ICML 2011 の論文 Hashing with Graphs の紹介です。Anchor Graph Hashing (AGH) という近傍探索のための最新のハッシング手法を提案しています (なんでアルゴリズム名をタイトルにしなかったんだろう……)。一応スライドを読めばなんとなくわかるんじゃないかと思います。
AGH は Spectral Hashing (SH) と同じ最適化問題を考えていて、それを解いてから out-of-sample extension の話に持っていくという流れも同じです。違うのは最適化問題を低ランク近似して解くところと、そのおかげで Nyström method という手法が使えて、SH みたいに分布に変な仮定を置いて固有関数問題をガリガリ解かなくて済むところ。データの次元、データ数のそれぞれについて線形時間のアルゴリズムになっています。
低ランク近似は、第一著者の Liu が ICML2010 提案している Anchor Graph に基づいています。この Anchor Graph が良い性質をたくさん持っていて、アンカーを用いたマルコフ連鎖的な解釈もあって、ハッシング以外にもいろいろ使えるんじゃないかと思います。
スライドでは論文の後半で述べられている階層的手法 (2-AGH) を紹介できませんでしたが、これは半分のビットは AGH で作って、残りは前半の各ビットに対して、うまく分けられなかった・分けるべきでなかった部分を補正するようなビットを入れる、という手法です。AGH は SH と同じく PCA みたいな手法なので、ビット数を増やすと固有値の低い (ノイズの多い) 主成分まで使うことになり、結果として MAP 精度が低下してしまう問題があります。階層的ハッシングでは半分のビットだけ AGH にして、後半はアドホックに作ることで、この問題をある程度回避しようとしています。が、半分は AGH で作るので、実験では 2-AGH でも下がったりしてます。
SH と同じ問題なので、多様体学習的な要素が入ってきています。なので与えられた距離関数についての最近傍探索ではないです。実際、実験結果では単純に全比較をした場合よりも精度が高いです。面白いのは、Spectral Embedding したあとに単純な全比較をする方法 (SE l_2 Scan) よりも階層的ハッシング (2-AGH) の方が精度が良いということ。近傍探索については、遠い点同士の類似度はざっくりゼロに落としちゃう方が精度が上がるということなんでしょうか。
Nyström method のあたりはちょっと難しいですが、SH に出てくる Laplace-Beltrami 作用素の固有関数問題のあたりに比べればだいぶわかりやすいと思います。[Bengio+, 04] をパラパラ読めば結構分かるんではないかと。固有ベクトルから固有関数を近似する方法として、結構使われている？みたいです。
ともかく Anchor Graph Hashing は面白いので、気になる方はぜひ論文も読むといいと思います。
2011年5月28日土曜日
PCA の練習
練習に http://archive.ics.uci.edu/ml/datasets/Iris のデータをさくっと PCA してみました。
g++4.6 -std=c++0x でコンパイルできます。RedSVD を使う方はコメントアウトしてあります (svd, pca_dat の定義を消してコメントを戻せば RedSVD 版になる)。RedSVD の方は上位 2 次元だけ出力しています。Eigen は魔法ですね。チュートリアルも親切。でもリファレンス読むのは難しい。
#include <iostream>
#include <string>
#include <utility>
#include <vector>
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <eigen3/Eigen/Dense>
// #include <redsvd/redsvd.hpp>
using namespace std;
using namespace boost;
int main() {
vector<pair<string, vector<float>>> data;
for (string line; getline(cin, line); ) {
vector<string> row;
split(row, line, is_any_of(","));
if (row.size() < 2) continue;
auto p = make_pair(move(row.back()), vector<float>());
row.pop_back();
for (const auto& val : row) p.second.push_back(lexical_cast<float>(val));
data.push_back(move(p));
}
Eigen::MatrixXf data_matrix(data.size(), data[0].second.size());
for (size_t i = 0; i < data.size(); ++i) {
for (size_t j = 0; j < data[i].second.size(); ++j) {
data_matrix(i, j) = data[i].second[j];
}
}
Eigen::JacobiSVD<Eigen::MatrixXf> svd(
data_matrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
const auto& pca_dat = svd.matrixU() * svd.singularValues().asDiagonal();
// const auto& pca_dat = REDSVD::RedPCA(data_matrix, 2).scores();
for (int i = 0; i < pca_dat.rows(); ++i) {
cout << data[i].first;
for (int j = 0; j < pca_dat.cols(); ++j) cout << ',' << pca_dat(i, j);
cout << endl;
}
}g++4.6 -std=c++0x でコンパイルできます。RedSVD を使う方はコメントアウトしてあります (svd, pca_dat の定義を消してコメントを戻せば RedSVD 版になる)。RedSVD の方は上位 2 次元だけ出力しています。Eigen は魔法ですね。チュートリアルも親切。でもリファレンス読むのは難しい。
登録:
投稿 (Atom)