SparseMatrixのデータ構造入門

★(後日記)この内容はちょっとわかりにくいし内容もあまり意味ないのでタイトルの内容について急いで知りたい人はこの後に書いた日記、こちらの改善された説明の方を参照してくれ。

行列の計算、とりわけ大規模なものともなると、世に出回っている定番のライブラリを使う方向にしか思考が行かない人が多いようだ。
ところが皆が皆そうでは、行列計算手法の進化なんて期待できない。本格的なビジネス目的であったり、パフォーマンスがとにかく最優先で余談を許さないような状況ならともかく、多少なりともオリジナリティを優先させられる余地があるのなら、その計算方法をなんとか自前で考えて見てはいかがだろうか。予備知識なしで自分で考えてみたものでも、結構使えるものになることがある。使えるとそれを改良するのが楽しくなってくる。するとさらにのめり込んで、もしかしたら新しいパラダイムが確立できるかも知れない。ここでは、c++stlのmapを使った、SparseMatrixのデータ構造について考えてみよう。

SparseMatrixのデータ構造については古く、CRSとかCCSに代表される多くの定番手法が知られているが、まずはそういった予備知識なしで考えてみよう。教科書ありき、先人の知恵ありきではない。まずは自分の頭でゼロから問題を考え、悩み、ある程度の高みに到達しながらもさらに先に進もうとして行き詰ったときにはじめて、先人にすがればよいのだ。では行ってみよう。
要は疎行列である。殆どの成分がゼロで、数%程度の要素が非ゼロの値をとるのである。ならば巨大行列全体(100000x100000とか!)のメモリを確保したりするのではなく、非ゼロの要素のみを覚えておけばよいではないか。
stl のpairやらmapを知っていればすぐに、

map< pair< int, int >, float >  m_sparse_map;

というデータ構造が思い浮かぶ。これの意味は、非ゼロ要素の「行id, 列idのペア」を、その「値」に関連付けて、map の中に保存しておくのである。この中に含まれない(i, j)のペアは、その要素がゼロであることを意味する。この方式なら、非ゼロ要素分しかメモリを必要とせず、A[i][j]の要素は何かと問い合わせられたときにも、即座にその値を返すことができる。
ところがこれ(だけ)では不十分だ。行列の基本的な演算では多くの場合、ある行に存在する値の全セットや、ある列に存在する値の全セットが必要となる。例えば行列にベクトルをかけたり、行列に行列をかけたりする場合がそうだ。
そのような要請にも対処できるように、ここではまた別のデータ構造を考えてみよう。手っ取り早くは以下のようなものになるだろう。

multimap< int, pair< int, float > > m_row_map;
multimap< int, pair< int, float > > m_col_map;

これは前者は、キーとする「行id」 に、その行の各非ゼロ要素の「列idと値のペア」を関連付け、後者では、キーとする「列id」に、その列の各非ゼロ要素の「行idと値のペア」を関連づけている。各行や各列には勿論複数の非ゼロ要素が存在するだろうから、multimapを使用した。
いかがだろうか。これで一番最初の m_sparse_mapでは、i行j列の値の問い合わせに即座に答えることができるし、続く m_row_map, m_col_map では、任意i行の全非ゼロ要素とその列id、および任意j列の全非ゼロ要素とその行idを、即座にセットで取り出すことが可能だ。(言うまでもないが、multimap::equal_rangeを使う)
ちなみに、筆者による上記方法の実装で、とある2000x2000程度のサイズの行列どうしの掛け合わせ演算速度は、通常のplainな配列ベースのデータ構造の場合に比べて、少なくとも20倍以上の高速化が達成された。これは行列のサイズが大きくなれば、100倍、1000倍ものパフォーマンスアップになることが予想されよう。
…ところで、m_row_map、 m_col_mapの実装において、multimapなんて持ち出す必要があるだろうか?単純な1次元配列の組み合わせでじゅーぶんなのではなかろうか?そうすれば一行で書くことはできなくなるが、パフォーマンスはあがるのではないだろうか?と、ここまでくれば、有名なCCS,CRSといったデータ構造の説明を紐解き、理解することは、比較的簡単なはずだ。
また、異なる行どうしを入れ替えたい場合はどうなるだろうか。などなど、いろいろと実用的な処理の実装を考えながらもっと優れたデータ構造について思いを巡らせるのはとても楽しい過程であるが、続きはまた後日、書き足すこととしよう。