3.2 稀疏矩阵存储与求解
文档摘要
3.2 稀疏矩阵存储与求解 本节摘要：scipy.sparse 用"只存非零元素"的思路处理非零占比极低的矩阵，COO、CSR、CSC 三种格式分别擅长构建、行操作与列操作。当矩阵上万阶、非零占比低于百分之几时，稀疏存储能省下几个数量级的内存，配合 sparse.linalg 的 cg、gmres 迭代求解器，才能把方程组求解推进到百万未知数。本节还会给出 diags 构造器、稀疏乘向量性能实测、预条件加速收敛与 eigs 挑算特征值的完整用法，并给出判断是否值得稀疏化的量化标准，让选型有据可依。 3.2 稀疏矩阵存储与求解 学习目标 阅读完本节，你应当能够： 用非零占比与矩阵规模两个数字判断一个问题是否值得用稀疏存储； 说清 COO、CSR、CSC 的存储结构、转换方法与各自擅长的操作；
3.2 稀疏矩阵存储与求解
本节摘要：scipy.sparse 用"只存非零元素"的思路处理非零占比极低的矩阵，COO、CSR、CSC 三种格式分别擅长构建、行操作与列操作。当矩阵上万阶、非零占比低于百分之几时，稀疏存储能省下几个数量级的内存，配合 sparse.linalg 的 cg、gmres 迭代求解器，才能把方程组求解推进到百万未知数。本节还会给出 diags 构造器、稀疏乘向量性能实测、预条件加速收敛与 eigs 挑算特征值的完整用法，并给出判断是否值得稀疏化的量化标准，让选型有据可依。
学习目标
阅读完本节，你应当能够：
- 用非零占比与矩阵规模两个数字判断一个问题是否值得用稀疏存储；
- 说清 COO、CSR、CSC 的存储结构、转换方法与各自擅长的操作；
- 用 diags、eye、random 等构造器快速搭建常见稀疏矩阵；
- 解释稀疏矩阵乘向量为什么快，以及哪些操作会悄悄把矩阵变稠密；
- 用 cg 与 gmres 求解大规模方程组，说明预条件如何改变收敛速度；
- 用 eigs 与 eigsh 只计算大规模矩阵的少数特征值。
一、问题与直觉
上一节结尾我们留了个问题：矩阵大到稠密存储装不下时怎么办。现在把数字摆出来。一个 1000 乘 1000 网格上的二维问题，未知数是 100 万。稠密存储要 100 万乘 100 万个浮点数，约 8 TB，任何一台普通机器都装不下。但这类问题的每个未知数只和上下左右四个邻居相关，矩阵每行只有五个非零元素——100 万个非零，约 80 MB。
这就是稀疏矩阵的全部动机：不存零，只存非零。先跑一段最小代码，看看稀疏矩阵长什么样。
import numpy as np from scipy.sparse import csr_matrix dense = np.array([[1, 0, 2], [0, 0, 3], [4, 0, 5]]) S = csr_matrix(dense) print(S)
输出是一份坐标清单：行号、列号、值，零元素一个都不出现。这种"打印出来只有非零项"的直观印象，正是稀疏存储的核心理念。
需要先说明的一点是：稀疏与否应该在建模阶段就决定，而不是等矩阵生成后再补救。从稠密矩阵转稀疏，内存早就花出去了；从数据源直接以坐标形式组织非零项，才是稀疏路线的正确打开方式。scipy.sparse 里所有格式共享同一套运算接口，加减乘、转置、切片、求和都写成一样的形式，格式切换只影响速度不影响正确性——这让我们可以先跑通再优化。
二、核心原理
2.1 三种主力格式
scipy.sparse 提供六种格式，主力是 COO、CSR、CSC，另外三种——LIL、DIA、BSR——在特定场景有专长。三种主力格式存储的是同一份数据，区别在组织方式。
COO（坐标格式）用三个等长数组分别存行号、列号、值，每个非零元素一行记录，构建最直观，适合从数据文件或计算过程一次性汇总坐标后创建。它不适合算术运算，做乘法、切片时都会先内部转换。
CSR（压缩行格式）把非零元素按行排序压缩：data 存值，indices 存列号，indptr 存每行的起点。行切片、行求和、矩阵乘向量都直接在这套结构上完成，是通用性最好的格式。
CSC（压缩列格式）是 CSR 的转置视角：按列压缩，列切片快，行切片慢。同一个矩阵，格式换一下，操作的快慢就换边。
三种格式可以随时互转：
from scipy.sparse import coo_matrix row = np.array([0, 0, 1, 2, 2, 2]) col = np.array([0, 2, 2, 0, 1, 2]) data = np.array([1, 2, 3, 4, 5, 6]) C = coo_matrix((data, (row, col)), shape=(3, 3)) Csr = C.tocsr() # 转成行压缩 Csc = C.tocsc() # 转成列压缩 print(Csr.indptr) # 行指针数组
转换本身有成本，一般是一次性的。构建阶段用 COO 或 LIL，计算阶段转成 CSR 或 CSC，这是最常见的节奏。
另外三种格式看一眼就好，用到时再回来翻。LIL 用每行的列表存非零，逐项赋值便宜，适合在循环里增量搭矩阵，但算术运算慢，搭完记得转走；DIA 按对角线存储，只对规则带状矩阵高效，有限差分里偶尔出现；BSR 把矩阵切成小块再按 CSR 的方式压缩，有限元里的刚度矩阵常带块结构，用它省索引开销。选型的口诀是：构建找 COO，增量找 LIL，运算找 CSR 或 CSC，特殊结构再谈 DIA 与 BSR。
图：稀疏矩阵三种格式对比
2.2 构造器：diags、eye、random
从稠密矩阵转稀疏只是入门，工程里更多是用构造器直接搭。diags 按偏移量一次构造多条对角线，是有限差分问题的标配：
from scipy.sparse import diags, eye, random n = 10 D = diags([-1, 2, -1], offsets=[-1, 0, 1], shape=(n, n)) I = eye(n, format='csr') R = random(n, n, density=0.05, random_state=42)
offsets 里的 0 是主对角，1 是上一条，负一是下一条。把数组换成 [-1, 0, 1, n, -n] 和 [-1, 4, -1, -1, -1]，就是 3.3 实战里五对角矩阵的原型。random 的 density 参数直接控制非零占比，测试格式性能时很好用。
2.3 性能：什么时候值得用稀疏
稀疏不是免费的午餐，格式有索引开销，运算有软件层成本。判断标准就两条：非零占比和矩阵规模。
| 非零占比 | 矩阵规模 | 建议 |
|---|---|---|
| 低于百分之五 | 上万阶 | 稀疏明显划算，内存差几个数量级 |
| 低于百分之五 | 千阶以下 | 稀疏可用，但优势不明显 |
| 高于百分之十 | 任意 | 稠密更快，稀疏的索引开销白花 |
| 带状且规则 | 任意 | DIA 或 CSR 都行，看操作类型 |
存储量可以量化对比。稠密要 n 的平方个浮点数；稀疏要非零个数个"值加索引"对。交叉点在非零占比约一半处，但工程上考虑到软件开销，非零占比百分之几以上就不值得折腾了。
性能上最典型的差距在矩阵乘向量：稠密是 n 的平方次乘加，稀疏只扫非零元素，复杂度与非零个数成正比。我们用三千阶三对角矩阵实测一把：
import time n = 3000 A_sp = diags([-1, 2, -1], offsets=[-1, 0, 1], shape=(n, n), format='csr') A_de = A_sp.toarray() v = np.random.default_rng(0).random(n) t0 = time.perf_counter(); y1 = A_de @ v; t_de = time.perf_counter() - t0 t0 = time.perf_counter(); y2 = A_sp @ v; t_sp = time.perf_counter() - t0 print(f"稠密 matvec: {t_de*1000:.2f} 毫秒, 稀疏 matvec: {t_sp*1000:.2f} 毫秒")
在多数机器上稀疏版本会快一个数量级左右。注意别把 n 放大——稠密那份 3000 阶矩阵已经占约 72 MB，再翻几倍内存先扛不住。
矩阵乘向量的优势同样适用于矩阵乘矩阵，但结果矩阵的稀疏度取决于两个矩阵的结构：稀疏乘稀疏可能还是稀疏，稀疏乘稠密则必然稠密，稠密乘稠密更是白费了稀疏的功夫。切片操作也吃格式：CSR 的行切片是取一段行指针范围，成本与切片行数相关；CSC 的行切片要跨列扫描，慢一个量级。列切片则正好反过来。所以写代码前先想清楚这步是行操作还是列操作，选对格式，比事后优化省心得多。
2.4 迭代求解器与预条件
稀疏矩阵的直接法也有：spsolve 用 SuperLU 做稀疏 LU 分解，对十万阶以内的问题又快又稳，不用管收敛参数，是首选。但规模再往上，直接法分解过程中会产生填充——零元素变成非零——内存可能失控，这时要换迭代法。
迭代求解器从初始猜测出发，反复做矩阵乘向量，直到残差降到阈值。每步成本与非零个数成正比，这是它能在百万阶规模存活的根本原因。scipy.sparse.linalg 里：
- cg 共轭梯度法，要求矩阵对称正定，收敛快，内存占用小；
- gmres 广义最小残差法，对一般矩阵适用，但可能需要重启；
- bicgstab 双共轭梯度稳定法，非对称矩阵的常见选择。
迭代法的软肋是收敛速度。矩阵条件数差时，几千步都压不到阈值。预条件就是给方程组做一次近似求逆的变换，把条件数拉下来。最常用的是不完全 LU 分解 spilu，把 A 近似分解成 L 乘 U，但砍掉填充，成本可控：
from scipy.sparse.linalg import cg, spilu M = spilu(A) Mx = lambda v: M.solve(v) x2, info2 = cg(A, b, M=Mx, tol=1e-10)
效果往往是迭代次数从上千降到几十。对角预条件取对角线倒数，更便宜，对付条件数温和的问题够用。
2.5 eigs：稀疏特征值
稠密特征值分解要算全部特征值，对大矩阵不现实。eigs 只算指定的少数几个：k 是数量，which 选头部还是尾部，sigma 提供位移反迭代——想算靠近某个值的特征值时用它。对称矩阵用 eigsh 更快更稳，svds 对应稀疏奇异值。用法上只要记住"稠密特征值全算，稀疏特征值挑着算"就够了。
典型场景是图论里的拉普拉斯矩阵：图很大，但每个节点的邻居有限，邻接矩阵天然稀疏，而谱聚类、网络稳定性分析只关心最小的几个特征值。这种问题用 eigs 是唯一可行的路线，全量分解的规模根本不允许。
from scipy.sparse.linalg import eigs, eigsh vals, vecs = eigs(A_sp, k=6, which='LM') # 模最大的六个 vals_s, vecs_s = eigsh(A_sp, k=6, sigma=0) # 对称矩阵，靠近零的六个
三、工程实践要点
3.1 格式选型速查
| 需求 | 格式 | 理由 |
|---|---|---|
| 从坐标数据构建 | COO | 三个数组直接给，无需排序 |
| 增量逐项修改 | LIL | 单元素赋值便宜 |
| 行切片、行求和、matvec | CSR | 行结构紧凑，扫一遍完事 |
| 列切片、列求和 | CSC | 与 CSR 对偶 |
| 对角与带状 | DIA | 对角存取天然高效 |
| 块结构（有限元） | BSR | 块级压缩，减少索引开销 |
3.2 常见问题
问：为什么我的稀疏矩阵运算变慢了？
答：多半是格式不对口，或者操作结果本身变稠密了。行操作请确认是 CSR，列操作请确认是 CSC；两个稀疏矩阵相乘的结果稀疏度取决于结构，别想当然。
问：迭代求解器不收敛怎么办？
答：先查条件数，再上预条件，最后才调 maxiter 和 tol。maxiter 给到十万步还收敛不了的，问题基本在矩阵本身，不在迭代器。
问：spsolve 和 cg 结果不一样，信哪个？
答：两边都是正规实现，差异通常来自容差——spsolve 是直接法，精度接近机器极限；cg 默认容差在 1e-5 量级。把 cg 的 tol 调到 1e-10 再比，差距就会消失。直接法与迭代法的取舍不是精度，而是规模与内存。
⚠️ 常见坑：对稀疏矩阵调用 toarray，或直接塞进 numpy 函数，矩阵被隐式转成稠密，几 GB 内存瞬间耗尽，程序往往不是报错而是被系统杀掉。需要检查小结果时，先切片再转稠密，或者用 data 数组直接看值。
💡 关键直觉：稀疏是存储问题也是算法问题——COO、CSR、CSC 解决"怎么存"，cg、gmres 解决"怎么算"。只有把迭代法配套上，稀疏化才真正完整，因为直接法在超大矩阵上的分解填充会偷偷吃掉你的内存预算。
一节小结
- 要点一：非零占比低于百分之几、规模上万阶时，稀疏存储才明显划算。
- 要点二：COO 适合构建，CSR 适合行操作与矩阵乘向量，CSC 适合列操作，转换一次到位。
- 要点三：diags 用偏移量一次搭出多对角矩阵，是偏微分方程离散的标准构造法。
- 要点四：稀疏矩阵乘向量复杂度与非零个数成正比，这是迭代法能跑百万阶的速度来源。
- 要点五：十万阶以内优先 spsolve 直接法；更大规模换迭代法，避免分解填充失控。
- 要点六：cg 要求对称正定，gmres 更通用，预条件能把迭代次数从上千压到几十。
- 要点七：eigs 只算指定的少数特征值，对称矩阵用 eigsh，别再全量分解。
下一节把这两节的知识装进一个真实问题：二维泊松方程的有限差分求解，稠密与稀疏两条路线同场竞技。