第11章 稀疏矩阵计算 - 科学计算的核心 (Sparse Matrix Computation - The Core of Scientific Computing)
一、导读
11.1 本章学习目标
通过本章的学习，读者将能够：
理解稀疏矩阵的定义、特性及其在科学计算和工程应用中的普遍性
掌握常见的稀疏矩阵存储格式（CSR、COO、ELL、HYB等）及其优缺点
深入理解稀疏矩阵向量乘法（SpMV）的并行化策略和优化技术
学会根据不同的稀疏模式选择合适的存储格式和kernel实现
掌握共享内存在SpMV优化中的应用，减少全局内存访问
理解稀疏矩阵的其他核心操作：转置、加法、乘法、预处理
能够分析SpMV的性能瓶颈，识别内存带宽和不规则访问的影响
了解现代稀疏线性代数库（如cuSPARSE）的功能和最佳实践
11.2 为什么本章如此重要
稀疏矩阵是指大部分元素为零的矩阵。在许多科学和工程应用中，稀疏矩阵无处不在：有限元分析中的刚度矩阵、图算法中的邻接矩阵、机器学习中的特征矩阵、电路模拟中的导纳矩阵、自然语言处理中的词袋模型等。这些矩阵的规模可能达到数百万甚至数十亿行，但非零元素的比例通常低于1%。
对于稀疏矩阵，使用稠密存储和计算方法会浪费大量的存储空间和计算资源。例如，一个100万×100万的矩阵，如果只有0.1%的非零元素，使用双精度稠密存储需要8TB内存，而稀疏存储可能只需要几百MB。更重要的是，稀疏矩阵向量乘法（SpMV）可以跳过零元素，显著减少计算量。
然而，稀疏计算也带来了独特的挑战。非零元素的分布可能是不规则的，导致内存访问模式不可预测。不同的存储格式在空间效率、访问模式和并行性方面有不同的权衡。选择合适的数据结构和算法对于性能至关重要。
在GPU上实现高效的稀疏计算尤其具有挑战性。GPU擅长规则的、大规模的并行计算，而稀疏数据的不规则性可能导致负载不均衡、内存访问分散和分支分歧。因此，稀疏计算是测试GPU架构灵活性和编程技巧的试金石。
在实际项目中，稀疏计算的性能直接影响整体应用的效率。在迭代求解器中，SpMV通常是每步迭代的主要开销；在图神经网络中，稀疏矩阵操作占据了大部分计算时间；在推荐系统中，稀疏矩阵分解是核心算法。因此，掌握稀疏计算优化技术对于构建高性能应用至关重要。
本章将从稀疏矩阵的基本概念出发，详细介绍各种存储格式的特点，然后深入探讨SpMV的并行化策略和优化技术，最后讨论其他稀疏操作和库函数的使用。
11.3 核心问题引导
什么是稀疏矩阵？如何量化稀疏度？稀疏计算的挑战是什么？
CSR、COO、ELL、HYB等存储格式各有什么特点？何时应该选择哪种格式？
SpMV的基本计算模式是什么？为什么它在GPU上难以高效实现？
如何使用共享内存优化SpMV？行分块和列分块策略各有什么优劣？
如何处理不规则的稀疏模式导致的负载不均衡问题？
稀疏矩阵转置的作用是什么？如何高效实现？
cuSPARSE库提供了哪些稀疏操作？如何与自定义kernel结合使用？
新兴的稀疏计算硬件（如Tensor Core对稀疏的支持）如何改变优化策略？
二、核心概念详解
11.1 稀疏矩阵的定义与挑战
稀疏矩阵的定义：
稀疏矩阵是指大部分元素为零的矩阵。虽然没有严格的数学定义，但在实践中，如果非零元素的比例低于10%，通常认为是稀疏的。在某些应用中（如有限元分析），非零比例可能低至0.01%或更低。
稀疏度的量化：
$$sparsity = \frac{\text{number of zero elements}}{\text{total elements}}$$
$$density = 1 - sparsity = \frac{\text{number of non-zero elements}}{\text{total elements}}$$
稀疏计算的核心挑战：
不规则的内存访问：非零元素的位置可能是任意的，导致随机访问模式，缓存命中率低。
负载不均衡：不同行的非零元素数量可能差异很大，导致线程间工作量不均。
索引开销：除了存储非零值，还需要存储位置信息（行索引、列索引），增加了存储和访问开销。
并行化困难：传统的稠密矩阵算法依赖于规则的行列结构，难以直接应用于稀疏矩阵。
格式转换成本：不同的算法可能需要不同的存储格式，格式转换本身也有开销。
11.2 稀疏矩阵存储格式
选择合适的存储格式是稀疏计算的第一步。不同的格式在空间效率、访问模式和适用算法方面有不同的权衡。
COO（Coordinate Format）：
COO是最简单的稀疏格式，使用三个数组存储：
values[]：非零元素的值row_indices[]：每个非零元素的行索引col_indices[]：每个非零元素的列索引
Example matrix:
[1 0 0 0]
[0 0 2 0]
[0 3 0 0]
[0 0 0 4]
COO representation:
values: [1, 2, 3, 4]
row_indices: [0, 1, 2, 3]
col_indices: [0, 2, 1, 3]优点：简单直观，易于构建和扩展。
缺点：冗余存储（每个元素都需要行索引），不适合快速行访问。
适用场景：矩阵构建、格式转换。
CSR（Compressed Sparse Row）：
CSR是最常用的稀疏格式，使用三个数组：
values[]：非零元素的值（按行优先顺序）col_indices[]：每个非零元素的列索引row_offsets[]：每行的起始位置在values数组中的索引（长度为行数+1）
Example matrix (same as above):
values: [1, 2, 3, 4]
col_indices: [0, 2, 1, 3]
row_offsets: [0, 1, 2, 3, 4] // row i has values from row_offsets[i] to row_offsets[i+1]-1优点：紧凑的行存储，适合行访问操作（如SpMV）。
缺点：插入新元素困难，列访问效率低。
适用场景：SpMV、行操作为主的算法。
CSC（Compressed Sparse Column）：
CSC是CSR的列版本，将行和列的角色互换：
values[]：非零元素的值（按列优先顺序）row_indices[]：每个非零元素的行索引col_offsets[]：每列的起始位置
优点：适合列访问操作。
缺点：行访问效率低。
适用场景：列操作、稀疏矩阵乘法的中间步骤。
ELL（ELLPACK）：
ELL将每行的非零元素填充到固定长度，存储为二维数组：
values[num_rows][max_nnz_per_row]：非零值col_indices[num_rows][max_nnz_per_row]：列索引
如果某行的非零元素少于最大值，用零填充。
优点：规则的访问模式，适合GPU并行化。
缺点：如果非零元素分布不均匀，填充浪费严重。
适用场景：非零元素分布均匀的稀疏矩阵。
HYB（Hybrid）：
HYB结合了ELL和COO的优点：
- ELL部分存储每行的大部分非零元素（固定数量）
- COO部分存储超出ELL容量的额外元素
优点：平衡了规则访问和灵活性。
缺点：实现复杂，需要调优ELL宽度参数。
适用场景：非零元素分布不均匀但大部分行有相似数量的矩阵。
DIA（Diagonal Storage）：
DIA存储矩阵的对角线：
values[num_diagonals][num_cols]：每条对角线的值offsets[]：每条对角线的偏移量
优点：对于带状矩阵非常高效。
缺点：只适用于特定结构的矩阵。
适用场景：有限差分法产生的带状矩阵。
格式对比表：
| 格式 | 空间效率 | SpMV效率 | 插入效率 | GPU友好性 |
|---|---|---|---|---|
| COO | 中 | 低 | 高 | 中 |
| CSR | 高 | 高 | 低 | 中 |
| CSC | 高 | 中 | 低 | 中 |
| ELL | 低* | 高 | 低 | 高 |
| HYB | 中 | 高 | 低 | 高 |
| DIA | 高** | 高*** | 低 | 高 |
*取决于填充程度
**对于带状矩阵
***仅适用于带状矩阵
11.3 稀疏矩阵向量乘法（SpMV）
SpMV是稀疏计算中最核心的操作，定义为：
$$y = A \times x$$
其中A是稀疏矩阵，x和y是稠密向量。
CSR格式的SpMV：
// Serial SpMV using CSR
void spmv_csr(const float* values, const int* col_indices,
const int* row_offsets, const float* x, float* y,
int num_rows) {
for (int i = 0; i < num_rows; i++) {
float sum = 0.0f;
for (int j = row_offsets[i]; j < row_offsets[i + 1]; j++) {
sum += values[j] * x[col_indices[j]];
}
y[i] = sum;
}
}GPU并行化策略：
策略1：一行一线程
每个线程负责计算一行的结果。这是最直观的映射。
__global__ void spmv_csr_one_thread(const float* values,
const int* col_indices,
const int* row_offsets,
const float* x, float* y,
int num_rows) {
int row = blockIdx.x * blockDim.x + threadIdx.x;
if (row < num_rows) {
float sum = 0.0f;
for (int j = row_offsets[row]; j < row_offsets[row + 1]; j++) {
sum += values[j] * x[col_indices[j]];
}
y[row] = sum;
}
}优点：实现简单，无竞争。
缺点：如果行的非零元素数量差异大，负载不均衡严重。短行的线程很快完成，长行的线程成为瓶颈。
策略2：多行warp协作
一个warp（32线程）协作计算一行。使用warp内的并行归约来累加结果。
__global__ void spmv_csr_one_warp(const float* values,
const int* col_indices,
const int* row_offsets,
const float* x, float* y,
int num_rows) {
int warp_id = (blockIdx.x * blockDim.x + threadIdx.x) / 32;
int lane_id = threadIdx.x % 32;
if (warp_id < num_rows) {
int start = row_offsets[warp_id];
int end = row_offsets[warp_id + 1];
int nnz = end - start;
float sum = 0.0f;
// Each thread in warp processes a subset of the non-zeros
for (int j = start + lane_id; j < end; j += 32) {
sum += values[j] * x[col_indices[j]];
}
// Warp-level reduction
#pragma unroll
for (int offset = 16; offset > 0; offset >>= 1) {
sum += __shfl_down_sync(0xFFFFFFFF, sum, offset);
}
// First thread in warp writes result
if (lane_id == 0) {
y[warp_id] = sum;
}
}
}优点：更好地利用warp并行性，适合中等长度的行。
缺点：对于非常长的行，warp可能不够；对于非常短的行，仍有浪费。
策略3：基于ELL的SpMV
ELL格式的规则结构使其非常适合GPU并行化。
#define MAX_NNZ_PER_ROW 32
__global__ void spmv_ell(const float* values,
const int* col_indices,
const float* x, float* y,
int num_rows) {
int row = blockIdx.x * blockDim.x + threadIdx.x;
if (row < num_rows) {
float sum = 0.0f;
for (int j = 0; j < MAX_NNZ_PER_ROW; j++) {
int idx = row * MAX_NNZ_PER_ROW + j;
float val = values[idx];
if (val != 0.0f) { // Skip padding zeros
int col = col_indices[idx];
sum += val * x[col];
}
}
y[row] = sum;
}
}优点：完全合并的内存访问，无分支分歧（如果所有行长度相同）。
缺点：填充浪费，不适合非均匀分布的矩阵。
11.4 共享内存优化
行分块策略：
将连续的多行加载到共享内存中，重用向量x的元素。
__global__ void spmv_shared_mem(const float* values,
const int* col_indices,
const int* row_offsets,
const float* x, float* y,
int num_rows) {
extern __shared__ float shared_x[];
int row = blockIdx.x * blockDim.x + threadIdx.x;
if (row < num_rows) {
// Load relevant portion of x into shared memory
// (Simplified: in practice, need to determine which x elements are needed)
float sum = 0.0f;
for (int j = row_offsets[row]; j < row_offsets[row + 1]; j++) {
sum += values[j] * x[col_indices[j]];
}
y[row] = sum;
}
}这个优化的效果取决于向量x的重用程度。如果相邻行访问相似的列，共享内存可以显著减少全局内存访问。
预取优化：
在计算当前元素的同时，预取下一个元素的列索引和值。
// Software prefetching in SpMV
float sum = 0.0f;
int next_j = row_offsets[row];
float next_val = values[next_j];
int next_col = col_indices[next_j];
for (int j = next_j + 1; j < row_offsets[row + 1]; j++) {
// Use prefetched values
sum += next_val * x[next_col];
// Prefetch next
next_val = values[j];
next_col = col_indices[j];
}
sum += next_val * x[next_col];11.5 负载不均衡处理
动态调度：
使用CUDA动态并行或工作队列来动态分配行给空闲线程。
行重排序：
将行数相近的行分组在一起，使得同一Block或warp内的行长度相似。这可以通过对行按非零元素数量排序来实现。
混合格式：
使用HYB格式，将规则部分用ELL处理，不规则部分用COO处理。
11.6 其他稀疏操作
稀疏矩阵转置：
将CSR转换为CSC（或反之）。这在某些算法中是必要的预处理步骤。
// CSR to CSC conversion (simplified)
void csr_to_csc(const float* csr_values, const int* csr_col_indices,
const int* csr_row_offsets,
float* csc_values, int* csc_row_indices,
int* csc_col_offsets,
int num_rows, int num_cols, int num_nnz) {
// Step 1: Count non-zeros per column
memset(csc_col_offsets, 0, (num_cols + 1) * sizeof(int));
for (int i = 0; i < num_nnz; i++) {
csc_col_offsets[csr_col_indices[i]]++;
}
// Step 2: Exclusive scan to get column offsets
// (Use prefix sum algorithm)
// Step 3: Scatter elements to CSC format
// (Process each row of CSR and place in correct column position)
}稀疏矩阵加法：
两个稀疏矩阵相加，需要合并它们的非零元素。如果使用CSR格式，需要同时遍历两个矩阵的行。
稀疏矩阵乘法：
两个稀疏矩阵相乘是最复杂的稀疏操作之一。结果的非零模式难以预先确定，通常需要临时数据结构。
预处理（Preconditioning）：
在迭代求解器中，预处理用于加速收敛。不完全LU分解（ILU）是常用的预处理方法。
11.7 cuSPARSE库
NVIDIA的cuSPARSE库提供了全面的稀疏线性代数操作：
- 格式转换：CSR、COO、CSC、ELL、HYB之间的转换
- SpMV：针对各种格式的优化实现
- SpMM：稀疏矩阵稠密矩阵乘法
- SpSV：稀疏三角求解
- SpMM：稀疏矩阵乘法
- 预处理：ILU、Jacobi等
使用cuSPARSE可以大幅简化开发，并获得经过充分优化的性能。
#include <cusparse.h>
void cusparse_spmv(cusparseHandle_t handle,
const float* csr_values, const int* csr_col_indices,
const int* csr_row_offsets,
const float* x, float* y,
int num_rows, int num_nnz) {
cusparseMatDescr_t mat_descr;
cusparseCreateMatDescr(&mat_descr);
cusparseSetMatType(mat_descr, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(mat_descr, CUSPARSE_INDEX_BASE_ZERO);
float alpha = 1.0f;
float beta = 0.0f;
cusparseScsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
num_rows, num_rows, num_nnz,
&alpha, mat_descr,
csr_values, csr_row_offsets, csr_col_indices,
x, &beta, y);
cusparseDestroyMatDescr(mat_descr);
}三、程序员的视角
3.1 在实际项目中的应用
有限元分析：结构力学、热传导、流体动力学等仿真使用稀疏矩阵表示离散化的偏微分方程。SpMV是迭代求解器的核心操作。
图算法：图的邻接矩阵是稀疏的。PageRank、社区检测、最短路径等算法可以表达为稀疏矩阵操作。
推荐系统：用户-物品交互矩阵是高度稀疏的。矩阵分解和协同过滤依赖稀疏计算。
自然语言处理：词袋模型和TF-IDF矩阵是稀疏的。文档聚类和分类使用稀疏向量操作。
3.2 常见陷阱与误区
格式选择不当：对于非均匀分布的矩阵使用ELL格式，导致大量填充浪费。应根据数据特点选择格式。
忽视索引开销：只关注值的存储，忘记索引也占用空间和带宽。对于非常稀疏的矩阵，索引开销可能超过值的开销。
负载不均衡：使用一行一线程策略处理非均匀矩阵，导致严重的线程闲置。应考虑warp协作或混合格式。
重复格式转换：在不同算法步骤之间频繁转换格式，转换开销超过计算收益。应尽量减少转换次数。
3.3 与其他层级的联系
与内存层次的联系：SpMV是典型的内存带宽受限操作。共享内存和缓存的利用对性能至关重要。
与前缀和的联系：CSR到CSC的转换需要使用前缀和来计算列偏移。
与归约的联系：Warp级SpMV使用归约来累加行结果。
四、代码实战
4.1 CSR格式的SpMV完整实现
#include <cuda_runtime.h>
#include <stdio.h>
// CSR SpMV with one thread per row
__global__ void spmv_csr_kernel(const float* values,
const int* col_indices,
const int* row_offsets,
const float* x, float* y,
int num_rows) {
int row = blockIdx.x * blockDim.x + threadIdx.x;
if (row < num_rows) {
float sum = 0.0f;
for (int j = row_offsets[row]; j < row_offsets[row + 1]; j++) {
sum += values[j] * x[col_indices[j]];
}
y[row] = sum;
}
}
void gpu_spmv_csr(const float* d_values, const int* d_col_indices,
const int* d_row_offsets, const float* d_x, float* d_y,
int num_rows) {
int blockSize = 256;
int gridSize = (num_rows + blockSize - 1) / blockSize;
spmv_csr_kernel<<<gridSize, blockSize>>>(d_values, d_col_indices,
d_row_offsets, d_x, d_y,
num_rows);
}
int main() {
// Example: 4x4 matrix with 4 non-zeros
int num_rows = 4;
int num_nnz = 4;
float h_values[] = {1.0f, 2.0f, 3.0f, 4.0f};
int h_col_indices[] = {0, 2, 1, 3};
int h_row_offsets[] = {0, 1, 2, 3, 4};
float h_x[] = {1.0f, 1.0f, 1.0f, 1.0f};
float h_y[4];
float *d_values, *d_x, *d_y;
int *d_col_indices, *d_row_offsets;
cudaMalloc((void**)&d_values, num_nnz * sizeof(float));
cudaMalloc((void**)&d_col_indices, num_nnz * sizeof(int));
cudaMalloc((void**)&d_row_offsets, (num_rows + 1) * sizeof(int));
cudaMalloc((void**)&d_x, num_rows * sizeof(float));
cudaMalloc((void**)&d_y, num_rows * sizeof(float));
cudaMemcpy(d_values, h_values, num_nnz * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_col_indices, h_col_indices, num_nnz * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(d_row_offsets, h_row_offsets, (num_rows + 1) * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(d_x, h_x, num_rows * sizeof(float), cudaMemcpyHostToDevice);
gpu_spmv_csr(d_values, d_col_indices, d_row_offsets, d_x, d_y, num_rows);
cudaMemcpy(h_y, d_y, num_rows * sizeof(float), cudaMemcpyDeviceToHost);
printf("SpMV result:\n");
for (int i = 0; i < num_rows; i++) {
printf("y[%d] = %.0f\n", i, h_y[i]);
}
// Expected: [1, 2, 3, 4]
cudaFree(d_values);
cudaFree(d_col_indices);
cudaFree(d_row_offsets);
cudaFree(d_x);
cudaFree(d_y);
return 0;
}4.2 使用cuSPARSE
#include <cusparse.h>
#include <stdio.h>
int main() {
cusparseHandle_t handle;
cusparseCreate(&handle);
// Create sparse matrix in CSR format
// ... (setup code)
// Perform SpMV
// ... (using cusparseScsrmv or newer API)
cusparseDestroy(handle);
printf("cuSPARSE SpMV demonstrated.\n");
return 0;
}五、深入思考
5.1 历史背景
稀疏计算的历史可以追溯到20世纪50年代，早期的大型科学计算就遇到了稀疏矩阵问题。CSR格式在1970年代成为标准，至今仍是使用最广泛的格式。
GPU加速稀疏计算的研究始于2000年代末。早期的工作集中在优化CSR格式的SpMV，后来发展出ELL、HYB等更适合GPU的格式。cuSPARSE库的发布标志着稀疏计算在GPU上的成熟。
近年来，Tensor Core开始支持稀疏矩阵乘法，利用结构化稀疏性获得额外加速。这代表了稀疏计算的新方向。
5.2 前沿方向
结构化稀疏性：利用特定的稀疏模式（如块对角、带状）获得更高的性能。NVIDIA的AMPERE架构支持2:4结构化稀疏，可获得2倍加速。
图神经网络：GNN中的稀疏矩阵操作是研究热点。如何高效处理动态变化的稀疏模式是一个挑战。
压缩稀疏格式：使用压缩技术进一步减少存储和带宽需求。
量子稀疏计算：量子计算机在处理某些稀疏问题上可能有指数级优势。
5.3 延伸阅读
NVIDIA cuSPARSE Library Documentation
Bell et al., "Implementing Sparse Matrix-Vector Multiplication on Throughput-Oriented Processors", SC 2009
Saad, "Iterative Methods for Sparse Linear Systems"
Vuduc et al., "Fast Sparse Matrix-Vector Multiplication on GPUs"
六、实验笔记
实验1：不同存储格式的SpMV性能对比
目的：比较CSR、ELL、HYB格式在SpMV中的性能。
步骤：
将同一矩阵转换为三种格式
测量每种格式的SpMV执行时间
分析非零分布对性能的影响
预期：对于均匀分布的矩阵，ELL最快；对于非均匀分布，HYB或CSR更好。
实验2：负载不均衡的影响
目的：量化行长度变化对SpMV性能的影响。
步骤：
生成具有不同非零分布的测试矩阵
测量一行一线程策略的执行时间
比较warp协作策略的效果
预期：行长度变化越大，一行一线程的性能越差。warp协作可以缓解但不完全解决问题。
七、本章小结
稀疏矩阵计算是科学计算和工程应用的核心。本章介绍了稀疏矩阵的定义、存储格式和核心操作。
CSR是最常用的稀疏格式，适合行访问操作。ELL格式具有规则的访问模式，适合GPU并行化。HYB格式结合了两者的优点，适合非均匀分布的矩阵。
SpMV是稀疏计算中最核心的操作。一行一线程策略简单但可能负载不均衡；warp协作策略更好地利用并行性；ELL格式的SpMV具有最佳的内存访问模式。
共享内存可以用于重用向量元素，减少全局内存访问。负载不均衡是稀疏计算的主要挑战，可以通过行重排序、动态调度或混合格式来缓解。
cuSPARSE库提供了全面的稀疏操作实现，在生产代码中应优先使用。
八、习题
格式转换：给定一个4×4矩阵，写出其CSR、COO和ELL表示。
SpMV复杂度：对于有M个非零元素的稀疏矩阵，SpMV的时间复杂度是多少？与稠密矩阵向量乘法相比如何？
ELL填充分析：如果一个1000×1000矩阵的平均每行非零数为5，但最大为100，ELL格式的填充率是多少？
负载不均衡：设计一个测试用例，使一行一线程的SpMV性能最差。解释原因。
共享内存优化：在什么情况下共享内存对SpMV的帮助最大？给出一个具体的稀疏模式示例。
稀疏矩阵转置：实现CSR到CSC的完整转换函数。
HYB参数调优：对于给定的矩阵，如何选择HYB格式中ELL部分的宽度？
结构化稀疏：解释NVIDIA Ampere的2:4结构化稀疏特性。它如何加速SpMV？
性能建模：估算CSR格式SpMV的全局内存访问次数。假设矩阵有M个非零元素。
综合实现：实现一个完整的HYB格式SpMV kernel，包括ELL部分和COO部分的处理。