疎行列の格納形式と演算

なぜ疎行列を専用形式で持つのか

有限要素法・有限差分法で偏微分方程式を離散化すると、係数行列はほぼ必ず疎行列（sparse matrix）になります。ある格子点・ある自由度が相互作用するのは近傍の一部だけであり、行列の大多数の要素はゼロだからです。次数 n の行列を密行列（dense matrix）としてそのまま格納すると、メモリは要素数に比例して増え、n が百万規模になると密表現は現実的な容量を超えます。

疎行列格納形式の目的は単純です。ゼロを一切保存せず、非ゼロ要素の値とその位置（行・列インデックス）だけを保持する。非ゼロ数を nnz（number of non-zeros）と書くと、必要メモリは nnz にほぼ比例し、行列サイズが大きくなっても nnz が緩やかにしか増えない限り扱えます。有限差分法の5点ステンシルなら1行あたりの非ゼロはたかだか5個程度で、行数が増えても1行あたりの非ゼロ数はほぼ一定です。

COO ── 最も単純な三つ組表現

COO（Coordinate format）は非ゼロ要素を (行番号, 列番号, 値) の三つ組の羅列として持つ形式です。構造がシンプルで、要素の追加や行列の組み立て（アセンブリ）が容易なため、有限要素法で要素ごとの局所行列を積み上げてグローバル行列を作る段階でよく使われます。

COO の内部表現（3x3の例、非ゼロが4個）

  行列:
    [10  0  0]
    [ 0  0 20]
    [30 40  0]

  row  = [0, 1, 2, 2]
  col  = [0, 2, 0, 1]
  val  = [10, 20, 30, 40]

欠点は、行や列を軸にした演算（行ごとの内積など）をしようとすると、目的の行番号を持つ要素を毎回線形探索するか、事前にソートしておく必要があることです。COOは組み立て用の中間形式であり、実際の反復法ソルバではCSRやCSCに変換してから使うのが通例です。

CSR ── 行方向アクセスに最適化された形式

CSR（Compressed Sparse Row）は、行ごとに非ゼロをまとめ、各行の開始位置だけをポインタ配列で持つことで索引の重複を圧縮します。3本の配列で構成されます。

CSR の内部表現（同じ3x3行列）

  values      = [10, 20, 30, 40]        非ゼロの値（行優先順）
  col_index   = [0, 2, 0, 1]            各値の列番号
  row_pointer = [0, 1, 2, 4]            行iの値はvalues[row_pointer[i] : row_pointer[i+1]]

  例: 行2の非ゼロは values[2:4] = [30, 40]、列は col_index[2:4] = [0, 1]

row_pointer は行数+1個の要素しか持たず、COOのように行番号を非ゼロの数だけ重複して持つ必要がありません。これによりCOO比でメモリがさらに削減され、かつ「i行目の非ゼロを取り出す」操作が範囲アクセス1回で完結します。この性質のおかげで、行単位に内積を取る疎行列ベクトル積（SpMV）や、反復法ソルバ（共役勾配法、GMRESなど）で最も多用される形式になっています。

一方でCSRは列方向の抽出、つまり「j列目の非ゼロを全部取り出す」処理には向きません。全行を走査して該当する列番号を探す必要があり、コストが行数に比例してしまいます。

CSC ── 列方向アクセスに最適化された形式

CSC（Compressed Sparse Column）はCSRの行と列の役割を入れ替えた形式で、列ごとに非ゼロをまとめ、各列の開始位置をポインタ配列に持ちます。

CSC の内部表現（同じ3x3行列）

  values      = [10, 30, 40, 20]        非ゼロの値（列優先順）
  row_index   = [0, 2, 2, 1]            各値の行番号
  col_pointer = [0, 2, 3, 4]            列jの値はvalues[col_pointer[j] : col_pointer[j+1]]

CSCは列方向の抽出が高速なため、列指向のアルゴリズムに向きます。代表例がLU分解やコレスキー分解での列単位のピボット処理・列消去です。疎行列を扱う直接法ソルバ（直接法、direct solver）の多くはCSCを内部表現として採用しています。逆にCSCで行方向の抽出をすると、CSRで列方向を抽出するのと同じ非効率が生じます。

疎行列ベクトル積（SpMV）の演算特性

反復法ソルバの主要コストはほぼすべてSpMV（y = A·x の計算）に集約されます。CSR形式でのSpMVは次のような単純な二重ループになります。

CSR形式でのSpMV疑似コード

  for i in 0..n-1:
      sum = 0
      for k in row_pointer[i] .. row_pointer[i+1]-1:
          sum += values[k] * x[col_index[k]]
      y[i] = sum

この計算の性能特性を理解する鍵は「演算密度（arithmetic intensity）」、つまり読み込んだデータ1バイトあたり何回の浮動小数点演算をするかです。SpMVでは非ゼロ1個につき乗算1回・加算1回の計2FLOPを行う一方、読み込むデータは値（8バイト）と列索引（4バイト）に加え、x[col_index[k]] という間接参照（インダイレクトアクセス、indirect access）が発生します。

さらに厄介なのが x[col_index[k]] の間接参照です。col_index の値は行列の非ゼロパターンに依存してばらばらであり、アクセスする x の要素はキャッシュライン上で局所性を持ちません。行列の非ゼロパターンが帯（バンド）状に近ければ col_index は近い値が連続し、xへのアクセスも局所的でキャッシュヒット率が上がります。逆に非ゼロが行列全体に散らばっているとxへのアクセスはランダムに近くなり、キャッシュミスが多発します。

この理由から、疎行列ソルバの前処理としてしばしば行・列の並べ替え（リオーダリング、reordering。Cuthill-McKee法などが代表例）を行い、非ゼロパターンをできるだけ対角近くに集めてバンド幅を縮小します。これによりxへのアクセスの局所性が改善し、キャッシュヒット率とSpMV全体のスループットが向上します。ブロック化された非ゼロ構造を持つ行列では、複数の非ゼロをまとめて格納するブロックCSR（BCSR）を使うと、インデックスの再利用が増え間接参照の回数自体を減らせます。

まとめ

• 疎行列格納形式は非ゼロ要素だけを保持し、メモリ使用量を行列サイズのO(n²)からO(nnz)へ落とすための表現である。
• COOは組み立てに向くが行/列アクセスが遅く、CSRは行方向、CSCは列方向のアクセスに最適化されており、反復法はCSR、直接法はCSCを好んで使う。
• SpMVは非ゼロあたりの演算量が少なく再利用も乏しいため演算密度が低く、性能はほぼメモリ帯域とキャッシュ効率で決まるメモリバウンドな処理である。
• 間接参照によるベクトルxへのランダムアクセスがキャッシュミスの主因であり、リオーダリングやブロック化（BCSR）による非ゼロパターンの整理が実効性能を左右する。