疎行列格納形式の話

久しぶりに何かまともな記事を書こうと思ったので書く。

自然現象をうまく計算可能な形に持っていくと大体連立一次方程式になる、らしい。

で、連立一次方程式は次元数的に規模は違えど大体が要素の殆どが零要素になっている疎行列(Sparse Matrix) になる、らしい。

で、規模が大きくなるにつれて零要素がどんどん増えていくので、格納効率は非常に悪い。

零要素は可能な限り無くして、メモリを効率的に扱いというあたりで疎行列格納方式を使う。

補足:疎行列を普通の行列の格納方法でメモリに入れようとした場合

例えば、5万次元の行列を普通にメモリに入れようとすることを考える。

float 型(4Byte) で、5万の行列を普通に確保しようとすると、

(sizeof(float) * 50000^2) / 1024 / 1024
= (4 * 50000^2) / 1024 / 1024 / 1024
= 9.536 (GByte)

ということで約9.5GBのメモリを消費する。

そのうち80%程度が零要素だと考えると、「ものすごいメモリの無駄遣い」というのが分かると思う。

本題に戻る

実際の行列を使うのも良いけど、結構めんどくさいので、適当な4次の正方行列を用意する。

零要素が40%ぐらい(7/16)。

これを、零要素がほぼ無い状態にしたい。

Compressed Sparse Row (Compressed Row Storage) 方式

研究室内だと CRS という呼び方で括弧内の呼び方が多いのだが、英語の論文(NVIDIA が出してるのとか)を読むと CSR の方がよく呼ばれているので、一応 CSR として話す。

CSR は、その通り、行成分を圧縮する。

これしか言うことないので、実際に上の行列を CSR に変換した後の状態を示す。

val が行列の成分、col が列の番号、 row が行成分となる。

row が次元数 + 1 のサイズを持つのは、C++イテレータ的に言うと、begin と end の関係が欲しいので。

これを C 言語で計算しようとすると、こんな感じになる。

int n = 4; // 次元数
float* val = ... ;
float* row = ... ;
float* col = ... ;
float* b = ... ;
float* x = ... ;

// Ax = b
for(int i = 0 ; i < n ; ++i) {
  b[i] = 0;
  for(int j = row[i] ; j < row[i + 1] ; ++j) {
    b[i] += val[j] * x[ col[j] ];
  }
}

Coordinate 形式

この形式は非零要素だけファイルに出力するとき、自然と似たような操作をしていると思う。

非零要素の行、列番号と値だけを保持すれば良い、ということ。

僕自身は計算に使うというよりは、ファイルへの保存時に使うので計算は省略。

ファイルへ出力する時は、ファイルの1行に行、列番号、値を出力すれば良い。

0 0 1
0 2 2
1 1 3
...

他にも CSR の逆で列を圧縮した CCR や、ELL, DIA, などがある。

その辺の話は NVIDIA からレポートが出ていて(http://www.nvidia.com/object/nvidia_research_pub_001.html)、そのレポートを原点とした CUDA で使えるライブラリで cusp-library - Generic Parallel Algorithms for Sparse Matrix and Graph Computations - Google Project Hosting というのもある。

このはてダでも ELL ぐらいは紹介しようと思う。