GPU求解稀疏线性系统速度远慢于CPU的原因排查

问题背景 #

使用SciPy（CPU）和CuPy（GPU）求解规模约为100000×100000的稀疏矩阵A，以及100000×8的稀疏矩阵B构成的线性系统，结果CPU耗时0.689秒，GPU耗时37.774秒，GPU性能远不如CPU。

核心代码 #

import cupyx
import cupyx.scipy.sparse.linalg
import time
import scipy
import scipy.sparse.linalg
import pathlib

file_dir = str(pathlib.Path(__file__).parent.resolve())

# ------------ cpu (scipy) ------------------------------- 
A_host = scipy.sparse.load_npz(file_dir+'/A.npz')
B_host = scipy.sparse.load_npz(file_dir+'/B.npz')

t00 = time.time()
P_host = scipy.sparse.linalg.spsolve(A_host, B_host)
t01 = time.time()

# ------------ gpu (cupy) -------------------------------  
A_device = cupyx.scipy.sparse.csr_matrix(A_host)
B_device = cupyx.scipy.sparse.csr_matrix(B_host)

t10 = time.time()
# 2nd parameter of the below function must be cupy.ndarray
P_seg = cupyx.scipy.sparse.linalg.spsolve(A_device, B_device.toarray())
t11 = time.time()

# ------------ timing results -------------------------------
print('cpu:', t01-t00)
print('gpu:', t11-t10)

关键原因分析 #

1. 输入格式强制转换导致稀疏性丢失 #

SciPy的spsolve支持直接传入稀疏矩阵B，而CuPy的spsolve要求右侧必须是稠密数组，代码中通过B_device.toarray()将稀疏B转为稠密矩阵，带来两个核心问题：

• 额外内存开销：稀疏B转为稠密后，内存占用大幅提升，增加GPU内存拷贝与管理的耗时；
• 算法优化失效：稀疏求解器针对稀疏右侧项的优化完全无法利用，只能执行通用的稠密矩阵运算，效率骤降。

2. 求解器底层实现的成熟度差异 #

SciPy的spsolve依赖SuperLU等经过数十年优化的CPU稀疏代数库，这些库针对各类稀疏矩阵结构（对称性、带宽、非零元分布）做了深度调优，适配CPU的缓存与多线程调度。而CuPy的spsolve底层实现相对较新，针对当前矩阵结构的优化不足，无法发挥GPU的并行优势。

3. 并行度与问题规模不匹配 #

GPU的优势在于大规模并行计算，但当前问题中B仅包含8个独立线性系统，且矩阵A的稀疏结构可能导致单步运算的并行度不足——GPU的大量核心无法被有效利用，反而因线程调度、内存访问延迟等问题拖慢速度。相比之下，CPU的多线程求解器能高效适配这类中等规模的稀疏矩阵运算。

4. 数据转换的隐性开销 #

代码中B_device.toarray()在GPU内部完成稀疏到稠密的转换，这一操作本身就会消耗额外时间，而CPU端直接使用稀疏矩阵求解，完全避免了这类转换开销。

优化建议 #

• 规避稀疏转稠密操作：若CuPy的spsolve不支持稀疏B，可将B拆分为单个列向量逐个求解，或直接调用CuSPARSE底层接口实现稀疏右侧项的求解；
• 选择适配的GPU求解器：针对矩阵特性（如是否对称正定），选用迭代求解器（如CG、GMRES）替代直接法spsolve——迭代法在GPU上通常比直接法更适合大规模稀疏矩阵；
• 提前完成数据部署：确保所有数据提前加载到GPU，避免计时范围包含不必要的跨设备数据拷贝；
• 分析矩阵结构：检查A的非零元分布、带宽等特性，若矩阵适合GPU并行，可尝试使用AmgX等专门的GPU稀疏求解库。

内容的提问来源于stack exchange，提问作者Andy