2016 — 稀疏矩陣與 CSR#

案件摘要#

一個 \(10^9 \times 10^9\) 的矩陣若用稠密格式儲存需要約 \(8 \times 10^{18}\) 位元組——全世界的硬碟加起來也裝不下。但若它像真實世界的矩陣一樣，每行只有幾十個非零元素，只需儲存這些非零值加上「它們在哪裡」的索引即可。CSR（Compressed Sparse Row，壓縮稀疏行）格式從 1960 年代的手工稀疏技術出發，經 SciPy 2001–2016 年的標準化，成為今日 Web 規模計算、有限元與圖神經網路的通用指紋格式。這是一樁「不儲存零，只儲存結構」的結案報告。

前因 -- 為什麼會有這個案子#

• 1950–60 年代，有限元方法（FEM，Turner、Clough 1956 起源；Zienkiewicz 1967 教科書）與結構分析產生巨大而稀疏的矩陣：離散化後每個節點只與相鄰幾個節點耦合。
• 1960 年代，IBM 的稀疏矩陣技術（如 1969 年 Tewarson 的整理、IBM 的 AS/SS 程式庫）開始系統化發展「只存非零元素」的資料結構；1973 年 Tewarson《Sparse Matrices》出版。
• 1970 年代，George 與 Liu 的巢狀剖分（nested dissection，1973/1981）證明稀疏消去的填充（fill-in）可以透過好的排序大幅控制。
• 1998 年 PageRank 問世：Web 超連結矩陣 \(10^{10}\) 規模、每列平均約 10 個非零——PageRank 迭代 \(v_{k+1} = \alpha A v_k + \dots\) 必須靠稀疏矩陣-向量乘法才跑得動。
• 2001 年 SciPy 問世（Travis Oliphant 等），把稀疏矩陣格式標準化為 Python 生態的一部分；2016 年前後 SciPy.sparse 成為科學 Python 的標配，並成為圖神經網路（GNN）框架的底層。

線索與推理 -- 數學式、程式、理論#

線索一：三種格式 —— COO、CSR、CSC 的分工#

COO（coordinate format）：三個陣列 (row, col, data) 各存非零元素的行、列、值。最簡單、允許重複項（相加自動合併），是「從零建構」的格式。

CSR（Compressed Sparse Row）：三個陣列 - data[j]：第 \(j\) 個非零值（按行優先排列） - indices[j]：該值在哪一列 - indptr[i]：第 \(i\) 行的非零元素在 data 中的起始位置

行 \(i\) 的非零元素就是 data[indptr[i]:indptr[i+1]]。這是「行操作最快」的格式：行切片、矩陣-向量乘、迭代求解皆高效。

CSC（Compressed Sparse Column）：CSR 的轉置版（按列壓縮），是「列操作最快」的格式——SciPy 的求解器內部常用。

三者的取捨就是本案的偵訊邏輯：COO 建構、CSR 計算、CSC 解方程。

線索二：稀疏矩陣-向量乘法 \(y = Ax\) —— 只碰非零#

CSR 格式下的 SpMV（sparse matrix-vector multiplication）：

for i = 0 .. m-1:
    s = 0
    for j = indptr[i] .. indptr[i+1]-1:
        s += data[j] * x[indices[j]]
    y[i] = s

運算量正比於非零元素個數 nnz(A)，而非 \(mn\) 。儲存量為 \(8\,\text{nnz} + 4\,\text{nnz} + 4(m+1)\) 位元組（double + int32）。對 \(10^9 \times 10^9\) 、每行 10 個非零的矩陣：稠密需 \(8 \times 10^{18}\) B，CSR 只需約 \(1.2 \times 10^{11}\) B ≈ 120 GB——從「不可能」變成「一個叢集可承載」。

關鍵問題是快取與並行：CSR 的行內運算天然獨立（各行可並行），但 x[indices[j]] 的間接索引造成記憶體存取不規則。現代解法包括 CSR 的分塊變體（BCSR/Blocked CSR，把鄰近非零打包）與 GPU 上的 SpMV 核（cuSPARSE）。SpMV 是 PageRank、共軛梯度法、GNN 訊息傳遞的共同熱點。

線索三：填充與排序 —— 稀疏消去的隱藏戰場#

直接對稀疏矩陣做 LU 分解會產生填充（fill-in）：原本為零的位置在消去過程中變成非零。以 1973 年 George 的巢狀剖分為例，對 \(\sqrt{n} \times \sqrt{n}\) 網格的有限元矩陣：

• 自然排序：填充 \(O(n^2)\)
• 巢狀剖分排序：填充 \(O(n \log n)\)

這是稀疏直接法的核心成果。間接法（共軛梯度、GMRES、多重網格）則完全避開填充，每步只需 SpMV——與 CSR 格式完美匹配。SciPy 中 scipy.sparse.linalg.spsolve（SuperLU/CXSparse 底層）與 cg/gmres 分別代表這兩條路線。

程式碼範例：CSR 格式手寫稀疏矩陣-向量乘法，對照 SciPy#

import numpy as np
import scipy.sparse as sp

def csr_matvec(indptr, indices, data, x):
    m = len(indptr) - 1
    y = np.zeros(m)
    for i in range(m):
        s = 0.0
        for j in range(indptr[i], indptr[i + 1]):
            s += data[j] * x[indices[j]]
        y[i] = s
    return y

# 建構一個稀疏矩陣：5x5 三對角 + 隨機非零
rng = np.random.default_rng(7)
rows, cols, vals = [], [], []
for i in range(5):
    for j in range(max(0, i - 1), min(5, i + 2)):   # 三對角帶
        rows.append(i); cols.append(j); vals.append(rng.normal())

A = sp.coo_matrix((vals, (rows, cols)), shape=(5, 5)).tocsr()
print("COO 建構 -> CSR:\n", A.toarray())
print("indptr =", A.indptr, "\nindices =", A.indices, "\ndata =", A.data.round(3))

x = rng.normal(size=5)
y1 = csr_matvec(A.indptr, A.indices, A.data, x)   # 手寫 SpMV
y2 = A @ x                                        # SciPy 內建
print("手寫 y =", y1.round(6))
print("SciPy y =", y2.round(6))
print("‖y手寫 - ySciPy‖ =", np.linalg.norm(y1 - y2))

# 規模對比：10000x10000、每行 10 個非零
n, nnz_per_row = 10_000, 10
nnz = n * nnz_per_row
rows = np.repeat(np.arange(n), nnz_per_row)
cols = rng.integers(0, n, size=nnz)
vals = rng.normal(size=nnz)
B = sp.coo_matrix((vals, (rows, cols)), shape=(n, n)).tocsr()
print(f"nnz={nnz}, CSR 記憶體 ≈ {(B.data.nbytes + B.indices.nbytes + B.indptr.nbytes)/1e6:.1f} MB")
print(f"稠密儲存需 = {n*n*8/1e6:.0f} MB  (倍數 ≈ {(n*n*8)/(B.data.nbytes+B.indices.nbytes+B.indptr.nbytes):.0f}×)")

輸出顯示手寫 SpMV 與 SciPy 內建完全一致；規模對比顯示 CSR 相比稠密節省約 800 倍記憶體——「只存非零」正是讓 \(10^9\) 規模計算成為可能的關鍵指紋。

結案 -- 後果與影響#

• SciPy.sparse 成為科學 Python 的標準稀疏格式：coo/csr/csc/dia/lil 等格式分工明確，spsolve、eigsh、cg 底層接上 SuperLU、ARPACK 等成熟 C/Fortran 函式庫。
• 有限元（FEniCS、deal.II 的 Python 介面）、結構力學、偏微分方程數值解全面以 CSR 為樞紐格式。
• 2000–2010 年代 PageRank 與 Web 圖計算（Hadoop/Spark 時代）把 SpMV 推到 \(10^{10}\) 規模。
• 2016 年起圖神經網路（Kipf & Welling 的 GCN, 2016）的訊息傳遞 \(\sigma(\tilde{A}XW)\) 本質上就是 SpMV——PyTorch Geometric、DGL 的底層即稀疏格式與稀疏乘法。
• 影響至今：推薦系統的隱式反饋矩陣、生物網路、交通網路、社會網路分析，凡有「大而稀」結構的地方就有 CSR 的指紋。

關鍵人物與文獻#

• R. P. Tewarson, Sparse Matrices, Academic Press (1973)。
• A. George, Nested dissection of a regular finite element mesh, SIAM J. Numer. Anal. 10, 345 (1973)。
• T. E. Oliphant, Python for Scientific Computing, Comput. Sci. Eng. 9(3), 10 (2007)；SciPy.sparse 文件 (2001–2016)。
• T. N. Kipf & M. Welling, Semi-Supervised Classification with Graph Convolutional Networks, ICLR (2017, arXiv 2016)。