📚 课程目录第4章稀疏矩阵存储与运算

4. 稀疏矩阵存储与运算

预条件技术的性能上限，往往不在算法公式本身，而在稀疏矩阵的存储方式与运算实现。一个百万阶的二维五点差分矩阵，非零元约五百万，若按稠密方式存储需要 8×10¹² 字节量级，完全不可行；而按稀疏格式存储只需约 60 MB。本节围绕三种核心能力展开：如何用 CSR/CSC 紧凑地存、如何高效地算 y = Ax、以及为什么分解过程会产生填充元并改变复杂度。

一、CSR 与 CSC：用三个数组描述一个矩阵

压缩稀疏行（CSR）用三个一维数组描述矩阵：values 按行依次存放非零元数值，col_idx 存放对应列号，row_ptr 长度为 n+1，row_ptr[i] 到 row_ptr[i+1]-1 是第 i 行的数据区间。压缩稀疏列（CSC）完全对称，只是把行、列角色互换，col_ptr 指向每一列的非零元区间。

以一个 4×4 示例矩阵说明（示例数据，用于演示结构）：

A = [ 4  -1   0   0 ]
    [ -1  4  -1   0 ]
    [ 0  -1   4  -1 ]
    [ 0   0  -1   3 ]

CSR: values  = [4, -1, -1, 4, -1, -1, 4, -1, -1, 3]
     col_idx = [0, 1, 0, 1, 2, 1, 2, 3, 2, 3]
     row_ptr = [0, 2, 5, 8, 10]

CSC: values  = [4, -1, -1, 4, -1, -1, 4, -1, -1, 3]
     col_ptr = [0, 2, 5, 8, 10]   （本例对称，故与 row_ptr 相同）

存储量约为 2·nnz + n + 1 个整数/浮点数，与 nnz 成正比而非 n²。选择原则很直接：按行遍历选 CSR，按列遍历选 CSC。稀疏矩阵向量乘 y = Ax 天然按行累加，CSR 是首选；而稀疏三角回代、按列访问的预条件构造（如 ILU 的列更新）更适合 CSC。转置只需交换 row_ptr 与 col_ptr 的角色并重排索引，代价为 O(nnz)。

二、稀疏矩阵向量乘：CSR 上的行式累加

CSR 格式下 SpMV 的核心循环非常简洁：对每一行 i，把该行所有非零元 values[k] 与对应 x[col_idx[k]] 相乘并累加到 y[i]。计算量为 2·nnz 次浮点运算，访存量约为 nnz 个值加 nnz 个列索引，因此算术强度极低，性能通常受内存带宽而非浮点峰值限制。

实际实现时有几个关键点。第一，y 必须先清零或显式覆盖写，避免残留值污染结果。第二，行循环内对 x 的访问是间接寻址（gather），列号不连续会导致缓存命中率下降；可按行分块、对列号做局部排序以改善局部性。第三，多线程并行时按行划分天然无写冲突，因为每行只写 y[i]，这是 CSR 相对 CSC 在并行 SpMV 上的优势。若需按列并行，则应改用 CSC 并采用 y 的原子累加或私有副本归约。

验证 SpMV 正确性的低成本方法：构造已知解析解的向量，例如令 x 全为 1，则 y[i] 应等于第 i 行非零元之和；再与稠密矩阵乘法结果逐元素比对，误差应在舍入量级（如 1e-12 相对误差）。

三、填充元：稀疏分解的隐性代价

对稀疏矩阵做 LU 或 Cholesky 分解时，原本为零的位置可能在消元过程中变为非零，这就是填充元（fill-in）。其根源是消元产生的更新项 a_ij ← a_ij − a_ika_kj/a_kk：只要第 i 行与第 j 行同时在第 k 列有非零元，位置 (i,j) 就可能被填充。填充元直接决定分解的存储量与运算量，是预条件技术必须正视的问题。

控制填充的主要手段是重排：通过对称近似最小度（AMD）、嵌套剖分（ND）等排序改变消元顺序，使因子尽可能稀疏。另一个手段是丢弃策略：不完全 LU（ILU）只保留满足阈值条件的填充元，用可控的精度损失换取稀疏性与速度。判断填充规模时，可先统计符号分解得到的因子非零元数 nnz(L+U)，它与数值无关，只由结构和排序决定，适合在正式计算前做代价预估。

四、复杂度分析：从 O(nnz) 到 O(nnz(L+U))

把上述结构串起来，可以得到一条清晰的复杂度链条。SpMV 的时间复杂度为 O(nnz)，空间为 O(nnz)；一次预条件应用（如三角求解）的复杂度为 O(nnz(L+U))；而分解本身的复杂度取决于填充规模，对二维网格问题，采用嵌套剖分排序后 Cholesky 分解约为 O(n^1.5)，三维问题则升至 O(n²)。这意味着：预条件越强，分解代价越高，需要在迭代次数减少与单次迭代成本之间权衡。

工程上常用的判断流程是：先统计 nnz(A) 与符号分解得到的 nnz(L+U)，计算填充率 nnz(L+U)/nnz(A)；若填充率过高（例如超过 10 倍），说明排序不理想或问题本身稠密化严重，应改用 ILU(k) 或 ILU-T 等带丢弃的策略，或直接退回 Jacobi/SSOR 这类无填充预条件。对同一问题，可对比不同排序下的填充率与迭代次数，选择总耗时最小的组合，而不是盲目追求最强预条件。

常见误区：把 CSR 的 row_ptr 长度误写成 n。它必须是 n+1，末位元素等于 nnz，否则最后一行的区间无法确定。另一个误区是在并行 SpMV 中让多个线程写同一 y[i]，这会产生数据竞争，应按行划分任务。

小结本节的关键判断：存储格式服务于访问模式，CSR 适合行式 SpMV 与并行，CSC 适合列式操作与三角求解；SpMV 是带宽受限的 O(nnz) 运算；填充元是稀疏分解的隐性成本，必须通过重排与丢弃策略主动控制，并用符号分解提前预估代价。