行列分解ライブラリredsvdを公開しました

大規模疎行列向けの行列分解ライブラリredsvdを公開しました．
redsvd

大規模疎行列向けの特異値分解や主成分分析，固有値分解を行うライブラリredsvdを公開しました．
修正BSDライセンスで公開しており，コマンドラインから使える他，C++ライブラリが用意されています．

例えば，行と列数がそれぞれ10万，非零の要素が1000万からなる疎行列に対する上位20位までの特異値分解を約2秒で処理します．

特異値分解とか，使っている技術の詳細とか応用事例を以下に簡単に紹介しましたので，興味のある方は参考にしてください．

特異値分解とは

まず行列を適当に復習します．行列Xの転置をX^tと表すことにします．またIを単位行列とし，Oを全ての成分が0である零行列とします．また，行列XX^t=IであるようなXを直交行列と呼びます．Xが直交行列の時，Xvはベクトルvを長さを変えずに回転させます．ここでは簡単のために，成分が全て実数である実行列のみを扱いますが，一般の場合では直交行列をエルミート行列と読み替えてください。

n行m列の行列Aが与えられた時，この行列はつぎのように必ず分解できます．この分解を特異値分解（SVD: Singular Value Decomposition)と呼びます．
A = USV^t
ただしU，Vはそれぞれn行n列，m行m列の直交行列であり，Sはn行m列の対角行列（対角成分以外は0）です．Sの対角成分は大きい順に並んでいるとします．これらを特異値と呼びます．UとVの列ベクトルをそれぞれ左特異ベクトル，右特異ベクトルと呼びます．Aを線形変換としてみた場合，A = USV^tという変換は任意の線形変換は，回転・軸毎のスケーリング・回転という三つの連続した操作に変換できることを意味します．

ちなみに
Ax = r xであるようなベクトルxをAの固有ベクトル，rを固有値と呼びます．それに対し特異値ベクトルの場合は
u_iA=s_i v_i
Av_i=s_i u_i
のような関係があります．但し，u_i, v_iはそれぞれU, Vのi列目の列ベクトルです．Aが対称行列の場合は，特異値分解A=USV^t=VSU^t=A^tとなり，V=Uとなって，固有値分解U^tAU = S（A=USU^t)に対応しますので，特異値分解は一般行列への特異値分解の一般化と見ることもできます．

この特異値分解は行列の"マスターアルゴリズム"であり，これさえできれば，線形方程式を効率的に解けたり，他の有用な行列分解（LU, QR）を効率的に求められたり，Aが観測行列の時、データ解析ができたりする非常に有用な分解です．例えば主成分分析は，主成分方向がVの各列であり，各サンプルの主成分スコアがUSの各列に対応するので特異値分解そのものです．レコメンデーションや画像索引，HMMの最尤解を求めるなどSVDの新しい使い方については最後の方で説明します．

特異値分解を理解するために別の表現をしてみます．先程の特異値分解で，Sが対角行列であることを利用して式を並び替えると，次のように表現できます
A = \sum_i s_i u_i v_i^t
但し，u_i, v_iはそれぞれU, Vのi列目の列ベクトルであり，s_iはi番目の特異値です．T_i = u_i v_i^tとした時T_iはAと同じサイズであり，また行と列の自由度は1であることから階数(rank)は1です．さらに，U, Vが直交行列であることからi≠jの時T_i^t T_j = Oとなり，{T_i}はお互い直交しています．よって，特異値分解はAをrankが1で直交している行列の重み付き和として表現していることになります．

この表現から導けることとして、Aの階数がrの時，Aはr個の行列の和A = \sum_{i=1...r} s_i u_i v_i^tで表せます。

これはさらに強いことがいえて，Sの要素のうち，i>rのときs_i=0とした行列をS~とした時，行列A~=US~Vは，階数がrである任意の行列の中でAとの二乗距離が最小であるような行列となっています．つまり階数rの中でAの最良の近似がA~となっていることになります．

S~の0となっているのに対応する部分をUとVから取り除き、U~,V~をそれぞれn行r列の直交行列とし，S~をr行r列で対角成分に特異値が並んだ対角行列とした時、A~=U~S~V~となります．

これをtruncated SVDと呼びます．

実世界から得られる多くの行列は，行数，列数に比べ実際の自由度は小さい場合が多いです．先程の文書・単語行列の例でいえば，同じような単語集合が同じ文書で出現するなどです．そこで上位の特異値のみからなるSVD，つまりtruncated SVDを行うことで、実際のAの情報をコンパクトに表現することができます。

また、多くの成分が0であるような行列を疎行列と呼びます．例えば，行が文書，列が単語に対応し，i行j列の成分が文書iに単語jが出現した時1，そうでなければ0であるとして，作られた文書-単語行列Aは非常に疎になります（英語Wikipediaの例だと非零の割合は0.1%）．レコメンデーションで使われるような行が購入者，列が購入したアイテムに対応する行列も疎行列の代表例です．従来の疎行列向けの固有値分解では，クリノフ部分空間を利用した手法が多く利用されていました。

redsvdが利用しているアルゴリズム

redsvdは次の論文で紹介されている乱択化アルゴリズムを元にしています．

"Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions", N. Halko, P.G. Martinsson, J. Tropp, arXiv 0909.4061 [pdf] [nips tutorial slide]

オリジナルのアルゴリズムが一方向のみにサンプリングしているのに対し，redsvdでは行と列の両方向でサンプリングをしていますのでそのまま論文で述べられている解析は適用できませんが、多くの場合うまくいっているようです。

今回利用している乱択化アルゴリズムが従来の乱択化を用いた手法とと違うのは，結果が非常に高い確率で真の値に一致するということです．オーバーサンプリングを利用することでさらに精度をあげることができます。redsvdを利用する場合は実際に欲しいrの値よりちょっと多め(+5ぐらい）に求めておき、上位だけを取り出すことで実現できます。詳細については上記論文を参照してください．

では具体的なアルゴリズムを紹介します．

n行m列の行列Aに対し上位r個の特異値に関する特異値分解をすることを考えてみます． 初めに各成分が平均0，分散1のGaussianからサンプルされたn行r列のGaussian行列Oを作ります． 次に，Y = A^t Oを求め，Yの各列を正規直交化します． この時AYY^t ≒ Aが成り立ちます． n行r列のB = AYを求めます． AはYによって列方向に圧縮していると考えられ，Y^tを書けることで元のAに復元できます．そして扱うのは圧縮した状態であるBです． 次に同様の方法で行方向に情報を圧縮します。 まず，各成分が平均0，分散1のGaussianからサンプルされたr行r列の行列Pをサンプリングし，Z= BPを求めます．そしてZの各列を正規直交化します． B = Z Z^t Bが成り立ちます． 最後にC = Z^t Bを求めます． Cはr行r列の行列であり，Z^tによってBを列方向に圧縮し，Zで復元できると考えられます．最後にCに対してSVDを従来手法で求め，C=USV^tを得ます． Cは非常に小さい行列なので，この部分は高速に行えます．

A = AYY^t = BY^t = ZZ^tBY^t = ZCY^t =ZUSV^tY^tとなります． ZU，YVは直交行列であり，それぞれAの左特異ベクトル，右ベクトルに対応します． またSはAの特異値が対角にならんだ対角行列です．

Eigenについて

以前C++の便利ツールとしてEigenを紹介しましたがredsvdはそれをフルに利用しています．例えばコア部分のコード(redsvd.hpp)は15行ぐらいです．疎行列の利用の仕方とかはredsvd.cppのconvertFV2Matとかを参考にしてみてください

特異値分解の応用

特異値分解は，解が一意に求まり（ベクトルが逆方向とか特異値が同じ場合の自由度はあるものの），しかもそれを効率的に求める手法が存在するという大きな利点があります．これは他のクラスタリング手法や隠れ層があるような手法には無い大きな特徴です．

SVD自体は古くからある技術であり，一度は様々な分野に応用されていましたが，最近再度，その重要性が注目されてきているようです

- 教師無学習による品詞推定
-- "SVD and Clustering for Unsupervised POS Tagging", M. Lamar, Y. Maron, M. Johnson and E. Bienenstock, ACL 2010
-- 各単語の前方コンテキストのBigram行列と，後方コンテキストのBigram行列をそれぞれSVDして，(USだけを利用して）隠れクラスタを求めた後にそれを特徴量としてクラスタリング．現在教師無学習の品詞推定で最高精度

- 索引
-- "Spectral Hashing", Y. Weiss, A. Torralba, R. Fergus, NIPS 2008
-- 以前記事にかいたり，WEB+DB Pressで記事を書きましたが，要は制約ありのグラフ分割問題に帰着して主成分分析で解く．

- 隠れマルコフモデル
-- "A Spectral Algorithm for Learning Hidden Markov Models", D. Hsu, S. M. Kakade and T. Zhang. COLT 2009
-- 隠れマルコフモデルの学習を従来のEMのような更新式で求めるのではなく，bigram,trigramの統計量からのSVD一発で求める
-- 高橋さんによる解説スライド [pdf]

もちろんスペクトラルクラスタリングなどにも利用できるでしょう。

今回のredsvdを利用して従来SVDや固有値分解が適用できなかったいろいろな手法に試していただければ幸いです．

redsvdの名前は普段飲んでいる飲み物由来でつけました。
なお，redsvdはPFIの20%プロジェクトで作りました．今後も何かできたら公開していきます．

C++の便利ツール・ライブラリ

フルタイムで働きはじめて4ヶ月。
いろんなことがありました。

今日はインターンが来ているということもあり日頃のC++コーディングライフの中で大変重用しているツールを紹介します。といってもどれも有名なツールでググれば解説がでてくるとは思いますので、一言ずつだけ紹介してみます。みなさんも何かよさげなライブラリ・ツールがありましたら教えてください。

- valgrind/callgrind/cachegrind

プログラムの実行結果を解析するツール群。まぁ、王道であえて紹介する必要はないかもしいませんが.。valgrindはプログラムのどこかでメモリが漏れているかどうかのチェックに使います．コードのどの部分で確保した領域がどこで漏れているかまで追跡することができます
　 valgrind --leak-check=full command

プログラムのどのが計算量的にボトルネックになっているかを調べるにはcallgrindを使います。ソースのどこで何クロック使っているかのレベルでわかります
　 valgrind --tool=callgrind command
　 callgrind_annotate -I=ソースの場所 --auto=yes callgrind.out.xxxx

cachegrindはcallgrindと同様にプログラムのどこで何回キャッシュミスしているかどうかを調べることができます。

- gtest

かゆいところに手が届くGoogle製C++テストフレームワークです。他のテストフレームワークを使ったことがないのでなんとも言えないですが、これでかなり満足しています。テストをちゃんと書かないとやっぱりだめですね。

- glog

Google製ロギングライブラリでサーバーなどを書いてログを吐く場合に使う。ずっと使い渋りをしていたが、使ってみるとこんなに便利なことはないと思う。面倒なログローテーションとかもやってくれるし落ちた時はどこで落ちたかも教えてくれる。gflagsを介した初期化をしたくないので、初期化の際、間に一つ挟んで使っています。

- re2

Google便利ツール三つ目、正規表現ライブラリ。これも直感的に使える。従来の多くの正規表現ライブラリは最悪時の計算量は爆発するが、re2は正規表現長に比例する計算量で必ず動作する保障があります。

- cmdline

tanakhさんが作ってくれたコマンドラインパーサー。fujimapでの使用例．gflagsは大きすぎるので、もっとコンパクトでヘッダだけのやつ作ってと頼んだらさらっと作ってくれた。

- waf

Pythonベースのビルドシステム。autotoolsがいくらたってもよくわからないところに現れた．といってもこれもよくわからないが最近良く使います。tanakhさんの解説が重宝します。

- openMP

共有メモリ型計算機環境下においてマルチスレッドなプログラム書くのに比べ圧倒的に簡単に並列化できる。ディレクティブをいれる形で並列部分を指示するのでopenMPが利用できない環境では勝手に無視されるのもよい。依存関係がないforではとりあえず使ってみてます。

- Eigen

最近特に使っている行列ライブラリ。C++上でmatlabのように直感的に書けて、しかもかなり速い(Eigen3βとかはGOTO BLAS並らしい）。SSE使えるときはつかってくれる。ヘッダファイルをコピーするだけでよい。解説記事がみあたらないので今度解説記事かいてみようかな。