为什么稀疏矩阵能加速计算：从压缩存储到高效求解

当处理的数据中，绝大多数元素的值为零时，这个数据就是“稀疏的”。在科学计算、机器学习和网络分析中，稀疏矩阵无处不在。如果按常规方式存储和计算，会浪费大量内存和算力。核心加速策略是：只存储和计算非零元素。

1. 识别与定义稀疏矩阵

• 判断标准：一个 $m \times n$ 的矩阵，如果其中零元素的个数远多于非零元素的个数（例如非零元素占比低于1%），就可视为稀疏矩阵。
• 常规存储的问题：二维数组存储一个 $10000 \times 10000$ 的矩阵，即使99%是零，也需要分配 $10^8$ 个存储单元，其中99%被浪费。计算时，需要遍历所有元素，与零进行大量无效运算。
• 目标：设计一种压缩存储方案，跳过所有零值，将计算资源集中于有意义的非零数据上。

2. 实现压缩存储：三种经典格式

压缩存储的核心是，用三个数组（或两个数组加一个计数）来记录非零元素的值、行位置和列位置。

2.1. COO 坐标格式

这是最直观的格式，直接记录每个非零元素的“坐标”。

• 数据结构：创建三个等长的数组。
  1. values[]：存储所有非零元素的值。
  2. row_indices[]：存储每个非零元素所在的行号。
  3. col_indices[]：存储每个非零元素所在的列号。
• 构造步骤：
  1. 遍历原始矩阵，找出所有非零元素。
  2. 将每个非零元素的值添加到 values 数组末尾。
  3. 记录该元素的行号和列号，分别添加到 row_indices 和 col_indices 数组末尾。
• 示例：矩阵 [[0, 5, 0], [0, 0, 3]] 可表示为：
  • values = [5, 3]
  • row_indices = [0, 1]
  • col_indices = [1, 2]
• 适用场景：非常适合动态构建矩阵（即逐个添加非零元素）。

2.2. CSR 压缩稀疏行格式

这是最常用的格式，特别适合按行进行矩阵-向量乘法。

• 数据结构：用三个数组表示。
  1. values[]：存储所有非零元素的值，按行顺序连续存放。
  2. col_indices[]：存储每个非零元素所在的列号（与 values 中的元素一一对应）。
  3. row_ptr[]：存储每行第一个非零元素在 values 数组中的起始位置。其长度为矩阵行数加一。
• 构造步骤：
  1. 创建 row_ptr 数组，大小为 行数 + 1，并初始化 row_ptr[0] = 0。
  2. 遍历矩阵每一行，统计该行非零元素的个数 nnz。
  3. 计算下一行的起始位置：row_ptr[i+1] = row_ptr[i] + nnz。
  4. 再次遍历矩阵，按行顺序将非零元素的值填入 values 数组，将其列号填入 col_indices 数组。
• 查找第 i 行的非零元素：
  1. 定位起始位置：start = row_ptr[i]
  2. 定位结束位置：end = row_ptr[i+1]
  3. 该行所有非零元素为 values[start] 到 values[end-1]，它们的列号为 col_indices[start] 到 col_indices[end-1]。
• 适用场景：高效的行遍历和矩阵-向量乘法。

2.3. CSC 压缩稀疏列格式

与CSR对称，是按列压缩的格式。

• 数据结构：类似CSR，但以“列”为单位压缩。
  1. values[]：按列顺序存储非零元素。
  2. row_indices[]：存储每个非零元素所在的行号。
  3. col_ptr[]：存储每列第一个非零元素在 values 数组中的起始位置。
• 构造步骤：参照CSR格式，只需将所有“行”的操作替换为“列”。
• 适用场景：需要频繁按列访问矩阵时，如在某些优化算法中。

3. 核心加速：稀疏矩阵运算

压缩存储的意义在于让运算跳过零值。

3.1. 稀疏矩阵-向量乘法

这是最核心的运算，很多迭代算法（如求解线性方程组）依赖它。以CSR格式的矩阵 $A$ 乘以稠密向量 $x$ 得到向量 $y$ 为例。

• 算法步骤：
  1. 初始化结果向量 $y$ 的所有元素为 $0$。
  2. 遍历矩阵的每一行 i：
     a. 获取该行非零元素的索引范围：start = row_ptr[i], end = row_ptr[i+1]。
     b. 遍历 j 从 start 到 end-1：
     • 获取列索引：col = col_indices[j]
     • 执行累加：$y[i] = y[i] + values[j] \times x[col]$
  3. 输出向量 $y$。
• 加速原理：整个乘法只进行了 nnz（非零元素总数）次乘加操作，而非 $m \times n$ 次。如果矩阵稀疏度为99%，运算量可降至原来的1%。

3.2. 稀疏线性方程组求解

许多科学计算问题最终归结为求解 $Ax = b$，其中 $A$ 是大型稀疏矩阵。

• 直接法：如基于高斯消去的LU分解。专用算法（如SuperLU）会设计消去顺序，避免在零的位置产生非零元（填入），保持矩阵结构的稀疏性，从而减少存储和计算。
• 迭代法：如共轭梯度法。其核心迭代步骤就是稀疏矩阵-向量乘法。只要能高效完成 $Av$，就能高效求解。这正是压缩存储的用武之地。

4. 实践指南：在Python中使用scipy.sparse

SciPy库提供了完善的稀疏矩阵工具。以下是如何创建、操作和求解一个稀疏线性方程组。

• 创建CSR矩阵：

  import numpy as np
  from scipy import sparse
  
  # 直接使用数据列表创建
  row = np.array([0, 0, 1, 2, 2])
  col = np.array([0, 2, 1, 0, 2])
  data = np.array([1, 5, 3, 4, 6])
  A_sparse = sparse.csr_matrix((data, (row, col)), shape=(3, 3))
  
  # 或者从稠密矩阵转换
  A_dense = np.array([[1, 0, 5], [0, 3, 0], [4, 0, 6]])
  A_sparse = sparse.csr_matrix(A_dense)

• 执行稀疏矩阵-向量乘法：

  x = np.array([1, 2, 3])
  y = A_sparse.dot(x)  # 高效计算，只处理非零元

• 求解稀疏线性方程组 $Ax = b$：

  from scipy.sparse.linalg import spsolve
  
  b = np.array([1, 2, 3])
  x_solution = spsolve(A_sparse, b)  # 内部使用优化的稀疏求解器