疎行列のデータ構造

最近、仕事で疎行列ばっか弄ってるのでまとめ。
COO形式、CSR形式、ELLPACK形式についてまとめます。

例えば以下の配列があるとします。

1	7	0	0
0	2	8	0
5	0	3	9
0	6	0	4

単純に行列として記録すると、上記を1次元に並べたデータ構造を使います。

values: 1 7 0 0 0 2 8 0 5 0 3 9 0 6 0 4

これとベクトルXを掛けてベクトルYを計算するプログラムはこうなります。

for (i = 0; i < M; ++i)
    for (j = 0; j < N; ++j)
        Y[i] += values[i*N+j] * X[j];

上記のデータだと0の部分が多いので、これを圧縮したのがCOO形式です。Coordinate形式の名の通り、座標を記録します。valuesに0が無くなるのが特徴です。

row:    0 0 1 1 2 2 2 3 3
column: 0 1 1 2 0 2 3 1 3
values: 1 7 2 8 5 3 9 6 4

これとベクトルXを掛けてベクトルYを計算する(Sparse Matrix-Vector multiplication、SpMVとよく呼びます)プログラムはこうなります。

for (i = 0; i < NZ; ++i)
    Y[row[i]] += values[i] * X[column[i]];

COO形式でもデータのrowの部分に同じ値が何度も出て来るので、これを圧縮したのがCSR形式です。Compressed Sparse Row形式の名前の通り、行を圧縮します。row_startの中に同じ数字が出てこなくなります。

row_start: 0 2 4 7 9
column:    0 1 1 2 0 2 3 1 3
values:    1 7 2 8 5 3 9 6 4

row_startは各行のデータが何番目から始まるかを示します。例えば、row_start[2]が4というのは、3行目(=2+1)のデータは、values[4]から始まるという意味です。

これのSpMVするプログラムはこうなります。

for (i = 0; i < M; ++i)
    for (j = row_start[i]; j < row_start[i+1]; ++j)
        Y[i] += values[j] * X[column[j]];

データは圧縮できていいんですが、計算に関節参照が増えます。ベクトル化しようとすると1行あたりのデータ数が一定じゃないのも辛い。

次は全然別のELLPACK形式。これはELLPACKというというelliptic problemsを解くためのシステムで利用されていたデータ構造ということで、この名前らしい。1行あたりのデータ数を決めて、パディングします。

column					values
0	1	*			1	7	*
1	2	*			2	8	*
0	2	3			5	3	9
1	3	*			6	4	*

*がパディングです。データは本来は縦方向に並べるのですが、今回は横方向に並べます。valuesの*は0を埋めておきます。1行あたりのデータ数が決まっているので、rowのデータは必要ありません。

column: 0 1 * 1 2 * 0 2 3 1 3 *
values: 1 7 0 2 8 0 5 3 9 6 4 0

これのSpMVするプログラムはこうなります。NPRが1行辺りのvalueの個数です。

for (i = 0; i < M; ++i)
    for (j = i*NPR; j < (i+1)*NPR; ++j)
        if (values[j] != 0)
            Y[i] += values[j] * X[column[j]];

if文で分岐せずに全て計算してしまっても構いません。そうするとこうなります。

for (i = 0; i < M; ++i)
    for (j = i*NPR; j < (i+1)*NPR; ++j)
        Y[i] += values[j] * X[column[j]];

valueの値の0を掛けるので、Xは何でもいいのです。ただし、領域外を参照しないようにcolumnの*には適当な値を設定する必要があります。0でもいいですが、キャッシュに無駄が出ないように一つ前と同じ値を埋めるのもいい方法です。同じ値を埋めると、データは以下のようになります。*の所にそれぞれ一つ前と同じ、1や2や3が入ります。

column: 0 1 1 1 2 2 0 2 3 1 3 3
values: 1 7 0 2 8 0 5 3 9 6 4 0

この方式の利点は規則的ってことです。またiとjのループを入れ替えると縦方向に処理ができます。大抵M >> NPRなので、この方がベクトル長が伸ばせて効率が良くなります。

欠点として、パディングの分データ量が増えます。1行あたりの0以外の値の数がほぼ均一なら問題はないのですが、長い行が混ざっていると無駄なデータが必要になってしまいます。

それを改善するのが、ELLPACKとCOOのハイブリッド方式です。1行につき2つのデータをELLPACK形式で保持し、残りのデータはすべてCOO形式にまとめます。COO形式が一つ一つのデータに依存関係がないため、どこからでも開始できるので、こうして他のデータを補完する形で使えるわけです。

ELLPACK部分
column: 0 1 1 2 0 2 1 3
values: 1 7 2 8 5 3 6 4

COO部分
row:    2
column: 3
values: 9

これのSpMVするプログラムはELLPACKとCOOを繋げた形になります。ここでは省略。

このハイブリッド形式は、ELLPACK部分が最適に作れて良さそうなデータに見えるのですが、1行あたりの0以外の値の数が均一である必要があるとか、要素の長い行は少なくないといけないとか、色々制限が多いので、実際の例ではあまり使われません。

続きはまた今度