2. 稀疏矩阵基础与存储优化

有限差分、有限元、电路网络等问题离散后得到的线性系统，绝大多数元素都是零。以二维泊松方程五点差分格式为例，若网格为 100×100，则未知量约 1 万个，而每行非零元只有 5 个左右，非零率不足 0.05%。如果仍用稠密二维数组存储，需要 10⁸ 个双精度数，约 800 MB；而按稀疏方式只存非零元，则只需约 5×10⁴ 个数值加索引，内存降到 MB 级。稀疏矩阵技术的第一价值，就是把存储和运算都压缩到"非零元"这一真实数据规模上。

数学表示：从三元组到压缩格式

稀疏矩阵在数学上仍是 A ∈ ℝ^m×n，只是非零元集合 nnz(A) 远小于 mn。最直观的存储是坐标格式 COO：用三个等长数组分别记录每个非零元的行号、列号和数值。COO 构建简单、便于从任意顺序的条目装配，但遍历时无法直接定位某一行的起点，随机访问效率低，通常作为"中间格式"使用。

真正用于计算的是压缩稀疏行 CSR。它把矩阵按行切分，用三个数组描述：values 按行优先顺序存放所有非零元数值；col_idx 存放对应的列号；row_ptr 长度为 m+1，row_ptr[i] 表示第 i 行在 values 中的起始位置，row_ptr[i+1]-row_ptr[i] 即该行非零元个数。第 i 行的非零元就是 values[row_ptr[i] … row_ptr[i+1]-1]。CSC 是 CSR 的转置视角：按列压缩，col_ptr 指向每列起点，适合按列访问或求解需要列遍历的算法。

CSR 的关键收益在于：行内元素连续存放，稀疏矩阵向量乘 y=Ax 可以写成对每行的一次连续扫描，缓存友好；同时 row_ptr 让"取第 i 行"成为 O(1) 定位。代价是插入新非零元需要移动数组，因此 CSR 适合"装配完成后一次性构建、之后反复使用"的场景，这正是有限差分矩阵的典型生命周期。

稀疏度与带宽：两个刻画结构的指标

稀疏度常定义为零元占比 1−nnz/(mn)，或直接用非零率 nnz/(mn) 描述。它决定存储压缩比，但不决定求解效率。带宽则刻画非零元离主对角线的最大距离：下带宽为 max(i−j)，上带宽为 max(j−i)，两者之和加一为总带宽。对一维三点差分，矩阵是三对角，带宽为 3；对 n×n 网格的五点差分，按行优先编号时，同列相邻节点的编号差为 n，因此半带宽为 n，总带宽约 2n+1。带宽越大，LU 分解产生的填充（fill-in）越多，直接法的计算量和内存都会显著上升。这也解释了为什么稀疏求解常配合 RCM、AMD 等重排序算法先降低带宽与填充。

存储效率对比与构建要点

以双精度数值（8 字节）和 32 位索引（4 字节）估算：稠密存储每个元素固定 8 字节；CSR 每个非零元占 8+4=12 字节，另加 row_ptr 的 4(m+1) 字节。对 nnz=5×10⁴、m=10⁴ 的二维五点差分，CSR 约 0.6 MB，压缩比超过 1000 倍。COO 则每个非零元需 4+4+8=16 字节，且无行定位能力，因此只建议在装配阶段临时使用。若索引超过 2³¹，需改用 64 位索引，内存会相应增加，这也是大规模问题中要留意的边界。

构建 CSR 的实用做法是"两遍扫描"：第一遍统计每行非零元个数，累加得到 row_ptr；第二遍按行填入 col_idx 与 values。下面以 n×n 网格五点差分（Dirichlet 边界，内部节点方程 4u_i,j−u_i−1,j−u_i+1,j−u_i,j−1−u_i,j+1=f_i,j）为例，用 Python 的 scipy 构建并核对内存：

import numpy as np
from scipy.sparse import csr_matrix

n = 100                      # 每边内部节点数，未知量 N = n*n
N = n * n
rows, cols, vals = [], [], []

def idx(i, j):               # 行优先编号
    return i * n + j

for i in range(n):
    for j in range(n):
        k = idx(i, j)
        rows.append(k); cols.append(k); vals.append(4.0)
        for di, dj in ((-1,0),(1,0),(0,-1),(0,1)):
            ii, jj = i + di, j + dj
            if 0 <= ii < n and 0 <= jj < n:
                rows.append(k); cols.append(idx(ii, jj)); vals.append(-1.0)

A = csr_matrix((vals, (rows, cols)), shape=(N, N))
print("nnz =", A.nnz)                       # 预期 5*N - 4*n（边界节点少项）
print("CSR 字节 ≈", A.data.nbytes + A.indices.nbytes + A.indptr.nbytes)
print("稠密字节 ≈", N * N * 8)
print("第 0 行非零列：", A.indices[A.indptr[0]:A.indptr[1]])

运行后可观察到：nnz 等于 5N−4n，因为四条边界上的节点缺少朝外的邻居；CSR 占用约 (12·nnz+4(N+1)) 字节，而稠密数组需 8N² 字节，两者相差三个数量级以上。遍历验证时，用 indptr 切片取出第 i 行的列号与数值，与按 (i,j) 手工推出的邻居集合比对，即可确认编号与装配方向一致。需要注意 scipy 的 indptr 即前文的 row_ptr，indices 即 col_idx，命名不同但语义相同。

最后要区分"存储压缩"与"计算高效"：CSR 解决了内存与行遍历问题，但带宽和填充仍由矩阵结构决定。实际工程中通常先由网格生成 COO 三元组，转为 CSR 后施加 RCM 等重排序降低带宽，再交给求解器。理解 row_ptr/col_idx/values 三者的对应关系，是后续阅读 SpMV 内核、剖析求解器内存行为的基础。