SciPy稀疏矩阵核心指南：格式选型、性能优化与实战应用

假设你手上有一张1000万行乘1000万列的关联矩阵，里面真正有数据的格子可能只有一百万上下。直接用普通数组存？光按float64算就是 800TB 内存，服务器当场罢工。换成 SciPy 稀疏矩阵，同一份数据通常只要几十 MB。夸张吗？不夸张，这就是scipy.sparse存在的意义。

这篇指南我会把构造、格式选型、常用操作、实际落地场景和调试优化一次性讲清楚。适合刚接触稀疏矩阵的小白，也适合用了一段时间但总在性能上踩坑的同学。重点不是罗列 API，而是告诉你每一种选择背后的原因，以及哪些写法会让性能直接崩掉。

1. 为什么稀疏数据必须用稀疏矩阵，先算一笔账

1.1 内存账本：普通数组和稀疏矩阵差了六个数量级

先来一个最直观的计算。假设矩阵形状是(10_000_000, 10_000_000)，也就是 1000 万行乘 1000 万列，元素总数是10^7 * 10^7 = 10^14。

用float64存，每个元素 8 字节，总字节数是10^14 * 8 = 8 * 10^14字节，也就是大约 800TB。这个量级不是“服务器够不够大”的问题，而是根本不现实。绝大多数关联矩阵、评分矩阵、图邻接矩阵，真实数据密度往往在 1e-6 以下，也就是说百万个元素里才有一个非零值。

同样的矩阵，假设非零元素只有 100 万个。用 CSR 格式（Compressed Sparse Row，压缩行存储）大致是这样：

• data数组：100 万 * 8 字节 = 8MB
• indices数组：100 万 * 4 字节 = 4MB（按 int32 估算）
• indptr数组：1000 万 + 1 个元素 * 4 字节 ≈ 40MB

合计约 52MB。这里还没算 Python 对象本身的少量开销，但整体就是从 800TB 降到 52MB 的差距。如果你用的是float32，data还能再省一半。

我见过不少项目，明明数据是高度稀疏的，却因为初期图省事用了普通numpy数组，结果把集群内存打满、任务频繁 OOM。项目上线前换成scipy.sparse，单机就跑完了。这往往是优化成本最低、收益最大的一次改动。

1.2 省内存只是表面，计算复杂度才是关键

很多人以为稀疏矩阵只是“压缩存储”，这个理解不够。稀疏矩阵真正的价值是让算法复杂度从“矩阵维度”转移到“非零元素个数”。

以矩阵乘向量为例。普通稠密矩阵：y = A @ x，需要O(n^2)次乘加。如果是 1000 万维，那就是10^14次浮点运算。而稀疏矩阵乘稠密向量只需要扫描非零元素，复杂度是O(nnz)，也就是 100 万次级别。这两种差距大到直接影响算法能不能落地。

这个性质在迭代算法里尤其关键。比如链接分析里常用的幂迭代，每一轮都是A @ p；图卷积、谱聚类、PageRank 风格的算法，最重的计算部分几乎都是稀疏矩阵乘向量（SpMV）。只要 A 保持稀疏、p 保持稠密，一轮迭代就是O(nnz)。如果你不小心把 A 转成稠密矩阵，几轮迭代后计算量就爆炸了。

所以“稀疏”不只是一个存储策略，它是让大规模数值计算变得可行的核心手段。

1.3 这些高频场景，天天都在生产稀疏矩阵

什么样的数据天然是稀疏的？几乎可以列出一长串：

• 推荐系统的用户-物品评分矩阵，绝大多数用户只对极少数物品有过交互。
• 社交网络、知识图谱的邻接矩阵，一个节点只连接海量节点中的极小一部分。
• 自然语言处理里的词共现矩阵，词汇表动辄几十万，但每个词共同出现的词有限。
• 有限元方法中的刚度矩阵，单元只与相邻单元耦合，矩阵呈带状或块状稀疏。
• 图像处理中的差分算子、离散拉普拉斯算子，也几乎都是稀疏的。

这些场景有一个共同点：如果你用稠密矩阵去描述，绝大多数存储和计算都在“瞎忙活”。稀疏矩阵并不是什么小众工具，而是数据工程、科学计算领域真正的日常。

2. 七种稀疏格式该怎么选，别只会 csr_matrix 和 coo_matrix

scipy.sparse提供了七种稀疏矩阵格式，很多人从头到尾只用过csr_matrix和coo_matrix，很大原因是教程只教这两种。但格式选错了，性能会差很多。

2.1 七种格式速查表

用一张表把每种格式的存储特点、擅长操作和典型用途列清楚：

大多数情况下，你只需要重点理解 COO、CSR、CSC 三种，然后记住一个原则：用 COO 或 DOK/LIL 去构造，用 CSR/CSC 去做计算。

2.2 构造阶段与计算阶段分开用格式，不是多此一举

为什么不能一步到位直接用 CSR 构建？因为每一行数据的插入成本不一样。

CSR 的存储结构是“行指针 + 列索引 + 数值”。它天生是为批量查询和计算设计的，但如果你要往中间插入一个新元素，就可能要把行内所有元素往后挪，最坏情况是每次插入都是O(nnz)。用循环往里填 10 万个元素，运行时间会非常难看。

COO 的思路完全相反。它只维护三个数组：row、col、data。新增一个元素就是往数组尾部追加，成本极低。它不适合直接做高效计算，但转成 CSR/CSC 只需要一次转换。

DOK 是基于字典散列的，适合“随机单点赋值”，比如A[3, 5] = 2.0。LIL 按行存储列表，适合逐行构造和行内切片。

我自己的习惯是：

• 数据能一次性凑齐三元组（行、列、值），优先 COO，然后.tocsr()。
• 必须按逻辑随机写入，比如模拟过程中动态填值，优先 DOK。
• 按行一轮一轮扩展，比如逐行读数据，优先 LIL。
• 无论哪种方式构造完，正式参与运算之前统一转成 CSR 或 CSC。

2.3 COO、CSR、CSC 之间怎么互相倒，成本高不高

A.tocsr()、A.tocsc()是频率最高的两个转换。它们的内部过程本质上是一次排序和合并：把三元组按照行号、列号排序，聚合同一下标。复杂度通常接近O(nnz log nnz)，取决于当前格式。对于百万级非零元素，一次转换耗时通常是几十毫秒到几百毫秒级别，完全可接受。

但要注意，这种转换不能放在热循环里。我见过有人在一个迭代算法里每一轮都做A.tocsc()，结果每轮都在做无用功。正确做法是：一次性转换，把转换后的矩阵保存下来反复用。

还有一个小技巧：A.T是视图操作，代价很小，不要把转置当成昂贵的复制操作。CSR 转成 CSC，本质上常常就是“转置视图再做一个内部转换”，比你想的便宜。

3. 构造与转换中的反直觉细节，这里最容易踩坑

3.1 最忌讳的写法：Python 循环逐点赋值

很多新手会写出类似这样的代码：

from scipy.sparse import dok_matrix import numpy as np A = dok_matrix((10000, 10000), dtype=np.float64) for i in range(100_000): A[i % 10000, (i * 7) % 10000] = i

这看起来没什么问题，但每次赋值都是一次 Python 方法调用，同时 DOK 内部散列表还要处理一次哈希插入。10 万次可能还能忍，如果是几百万次非零元素，这个循环就会成为整个程序的瓶颈。

同样的数据，用 COO 批量构造：

from scipy.sparse import coo_matrix import numpy as np n = 100_000 row = np.arange(n) % 10000 col = (np.arange(n) * 7) % 10000 data = np.arange(n).astype(np.float64) A = coo_matrix((data, (row, col)), shape=(10000, 10000)).tocsr()

所有下标计算都在向量化层完成，Python 循环完全消失。我用同等规模数据对比过，批量方式通常比逐点循环快几十倍以上。数值越大，差距越夸张。

记住一个铁律：能用数组构造就不要用循环，能批量就不要逐点。

3.2 三条可靠的构造路径

按照数据来源不同，我会选择不同的构造路线。

第一种：你已经有一张“行、列、值”的表，比如从数据库查出了三万行(user_id, item_id, rating)。

from scipy.sparse import coo_matrix, csr_matrix rows = [0, 1, 2, 0] cols = [1, 0, 2, 2] vals = [1.0, 2.0, 3.0, 4.0] A = coo_matrix((vals, (rows, cols)), shape=(3, 3)).tocsr()

第二种：程序模拟过程中需要随机更新单个元素，用 DOK。

from scipy.sparse import dok_matrix A = dok_matrix((1000, 1000), dtype=np.float64) A[123, 456] = 2.5 A[0, 0] = 1.0 # 构造完成后转成 CSR A = A.tocsr()

第三种：按行追加数据，用 LIL 逐行填，再转 CSR。

from scipy.sparse import lil_matrix A = lil_matrix((1000, 1000)) A[100, :] = np.random.rand(1000) A = A.tocsr()

从已有稠密矩阵转稀疏也可以，只有一个小门槛：

sparse_A = csr_matrix(dense_A)

但我一般只在小矩阵时这么干。如果稠密矩阵本身大到快把内存撑爆，再一次性转稀疏，那只是又制造了一次内存峰值。正确思路是从源头就直接构建稀疏结构。

3.3 重复下标、未排序和显式零，三个不容易注意的细节

COO 构造时，如果同一个(row, col)出现了多次，SciPy 默认会把值累加。这个行为有时候是好事，比如你在装配有限元矩阵时，同一个位置会从多个单元贡献刚度值，累加正是你想要的。

但有时候你并不希望累加，那就先手动处理：

A.sum_duplicates()

另一个坑是索引未排序。CSR 格式并不强制要求行内索引有序，但很多运算在假设排序时才会更快。你可以在转换后主动调用：

A.sort_indices()

还有一个容易被忽略的坑：显式零。稀疏矩阵里的0跟“没有元素”是两回事。如果你因为某种原因写入了0，它会占用data里的一个位置，拖慢计算。清理方法：

A.eliminate_zeros()

这三个方法都不是花架子，处理大规模数据时能显著影响后续性能和存储体积。

4. 高效操作与计算实战，这些 API 行为要知道

4.1 切片和索引：行快列慢，转置比你想象中便宜

CSR 格式的对齐单位是“行”。所以行切片非常高效：

row_block = A[100:200, :]

这会返回一个行数更少的 CSR 矩阵，内部通过indptr直接截取区间，速度很快。

但列切片就不同了：

col_block = A[:, 100:200]

CSR 要遍历每一行去抽列，效率明显差。如果代码里频繁出现列切片，最好的办法是把矩阵转成 CSC：

B = A.tocsc() col_block = B[:, 100:200]

CSC 就是 CSR 的“转置视角”，列操作天然高效。记住这个原则就行：行操作用 CSR，列操作用 CSC，别让切片去硬扛不适合自己的格式。

还有一个容易混的点：对稀疏矩阵做A[0]，返回的仍然是一个稀疏矩阵，不是 numpy 数组。如果你真的要拿它参与某些只能处理数组的计算，才需要.toarray()。而toarray()的结果是稠密向量，长度等于列数，如果列数很大，内存会瞬间涨起来，这个行为一定要心里有数。

4.2 逐元素操作与统计：先问一句你的矩阵是什么格式

稀疏矩阵加、减、逐元素乘都遵循常规语义。但有个细节：当两个格式不同的矩阵做运算时，SciPy 可能会自动转换格式，这个转换可能带来额外开销。所以尽量让参与运算的两个矩阵使用同一种格式。

逐元素乘法用multiply，注意它和矩阵乘法@完全不同：

C = A.multiply(B) # 逐元素相乘，保持稀疏 D = A @ B # 矩阵乘法，结果可能是稀疏或稠密

统计类操作也要小心。A.sum(axis=0)返回的结果类型带有矩阵形态，不是普通一维数组。想拿到可用的一维数组，要这样处理：

col_sums = np.asarray(A.sum(axis=0)).ravel() row_sums = np.asarray(A.sum(axis=1)).ravel()

一个最常见的反例是np.sum(A)。如果你把稀疏矩阵直接传给 numpy 的 sum 函数，numpy 为了统一接口，会先把矩阵转为稠密，然后求和。小矩阵无所谓，大矩阵直接内存爆炸。稳妥写法永远是A.sum()。

4.3 矩阵乘法里的“隐形转换”与“稠密化”

矩阵乘法是稀疏矩阵最常用的操作，但它的结果有时是稀疏的，有时是稠密的，取决于你乘的是什么。

• A @ B，如果 A、B 都是稀疏矩阵，结果是稀疏矩阵。
• A @ x，如果 x 是普通一维 numpy 数组，结果是一个稠密的一维数组。
• A @ X，如果 X 是二维稠密矩阵，结果是一个二维稠密数组。

这个行为本身合理，但要注意一个常见陷阱：如果你想让稀疏乘以稀疏保持稀疏，必须保证两边都是稀疏格式。不要把数组随手传入@，否则结果会被稠密化。

举个例子，做对称归一化时，我会这样写：

from scipy.sparse import diags import numpy as np row_sum = np.asarray(A.sum(axis=1)).ravel() inv_sqrt = diags(1.0 / np.sqrt(row_sum + 1e-10)) L = inv_sqrt @ A @ inv_sqrt

这里inv_sqrt是稀疏对角矩阵，A是 CSR，整个乘法链保持稀疏。如果你把它写成inv_sqrt.toarray() @ A.toarray() @ inv_sqrt.toarray()，内存和计算量会同步爆炸。

还有一个需要留意的场景：A.T @ A在推荐系统里非常常见，结果是物品-物品共现矩阵。它即使保持稀疏，非零元素也可能比你预期多得多，因为高热度物品集中会让结果变密。建议做完后看一下.nnz，检查稀疏度是否还在可控范围。

5. 三个落地场景里的常用套路，直接照抄

5.1 图算法：CSR 的 indptr 和 indices 就是天然邻接表

很多做图算法的人容易忽视这一点：CSR 格式的三件套indptr、indices、data，本质上就是压缩后的邻接表。

比如要遍历某个节点的所有邻居：

def neighbors(A, node): start, end = A.indptr[node], A.indptr[node + 1] for j in range(start, end): yield A.indices[j], A.data[j]

这一操作不需要任何额外存储，也避开了 Python 层面对动态结构的依赖，速度非常可观。

做图的归一化拉普拉斯时，也可以直接用稀疏操作：

from scipy.sparse import eye, diags A = csr_matrix(...) D = eye(A.shape[0]) # 或直接用度 D_inv_sqrt = diags(1.0 / np.sqrt(np.asarray(A.sum(axis=1)).ravel())) L_norm = eye(A.shape[0]) - D_inv_sqrt @ A @ D_inv_sqrt

只要每一步都保持稀疏格式，图规模在百万节点级别也能跑得动。

经典的链接分析迭代也可以这么写：

p = np.ones(n) / n alpha = 0.85 for _ in range(50): p = alpha * (A @ p) + (1 - alpha) / n

每一轮都是O(nnz)的矩阵乘向量，几千万元素的图在单机上也能完成。

5.2 推荐系统评分矩阵：行统计、列统计和相似度计算

用户-物品评分矩阵是稀疏矩阵最典型的应用。假设你有三张表：用户 id、物品 id、评分值，直接用 COO 构造：

rows = user_ids cols = item_ids vals = ratings R = coo_matrix((vals, (rows, cols)), shape=(n_users, n_items)).tocsr()

计算用户活跃度、物品热度：

user_activity = np.asarray(R.sum(axis=1)).ravel() item_popularity = np.asarray(R.sum(axis=0)).ravel()

计算物品相似度，最直接的方式是：

item_sim = R.T @ R

这里R.T是 CSC 视图，R.T @ R内部优化得当，比手动写两层循环不知道快多少倍。

如果你在做矩阵分解类的迭代算法，需要拿到某个用户的交互物品列表，利用 CSR 的indices就能快速拿到：

start, end = R.indptr[user_id], R.indptr[user_id + 1] items = R.indices[start:end] ratings = R.data[start:end]

这种做法比循环里每次R[user_id].toarray()高效得多，尤其在用户总数和物品总数都很大的时候。

5.3 数值求解：组装稀疏矩阵，交给 spsolve

有限元、电路仿真、流体计算这类场景里，线性方程组的系数矩阵天然稀疏。常见做法是先把每个单元的贡献按三元组收集起来，装配成 COO，再转成 CSC 或 CSR 送去求解。

from scipy.sparse import coo_matrix, csc_matrix from scipy.sparse.linalg import spsolve # 装配阶段，反复追加三元组 rows = [] cols = [] vals = [] # ... 循环里不断往这三个列表追加 ... K = coo_matrix((vals, (rows, cols)), shape=(n_dof, n_dof)).tocsc() x = spsolve(K, b)

spsolve适合中小规模问题，背后会选择合适的稀疏 LU 或 Cholesky 分解。对于大规模矩阵，我更喜欢迭代求解器：

from scipy.sparse.linalg import cg, gmres x, info = cg(K, b, rtol=1e-8)

迭代求解器配合好的预条件，往往能处理百万甚至千万自由度的问题。另外还有一个宝藏工具LinearOperator，它允许你定义一个隐式的线性算子，不必显式构造矩阵，适合矩阵结构复杂但向量乘法简单的场景。

6. 调试、内存观测与最终优化清单

6.1 五分钟摸清一个稀疏矩阵的家底

拿到一个陌生稀疏矩阵，第一件事就是看清楚它是什么格式、多大、占多少内存。我会用一个小函数：

import numpy as np from scipy.sparse import issparse def inspect_sparse(A): if not issparse(A): print("not sparse matrix") return print("format:", A.getformat()) print("dtype:", A.dtype) print("shape:", A.shape) print("nnz:", A.nnz) print("density:", A.nnz / np.prod(A.shape)) total = 0 for name in ("data", "indices", "indptr", "row", "col"): if hasattr(A, name): total += getattr(A, name).nbytes print("approx memory: {:.2f} MB".format(total / 1024**2))

看nnz和density能快速判断数据稀疏度是否合理。如果某次操作后nnz暴涨，说明可能发生了意外的稠密化，这是排查性能问题时最敏锐的指标。

6.2 性能优化时的判断顺序

遇到稀疏矩阵操作很慢，我会按下面这个顺序排查：

1. 格式和操作是否匹配。行切片用 CSR，列切片用 CSC，别让指令做它不擅长的事。
2. 是不是在循环里反复做转换。tocsr()、tocsc()应该只做一次。
3. 是不是不小心把稀疏矩阵转成了稠密。搜索代码里的.toarray()、np.array(A)、np.sum(A)。
4. 是不是在 Python 层用循环处理非零元素。能向量化就向量化，能批量就批量。
5. 结果矩阵有没有过度膨胀。操作后看一眼.nnz，评估是否需要剪枝或截断。

很多时候性能问题不是“SciPy 不够快”，而是数据在格式和稠密状态之间反复横跳，把成本都花在了无意义的转换和搬运上。

6.3 我在实践里踩过的几个坑

第一个坑是np.sum(A)。处理大矩阵时，这一步会尝试把稀疏矩阵转为稠密，运气好是警告，运气差就是内存溢出。后来我强制自己所有矩阵统计都用A.sum()，没有例外。

第二个坑是A[0]的返回值类型。初学者很容易把它当成数组去跟np.array直接做运算，结果发现类型对不上，又顺手补一个.toarray()，然后又制造一次不必要的稠密化。除非真的需要展示或做小规模调试，否则尽量让数据保持稀疏格式。

第三个坑是循环里逐点赋值。早年写模拟代码，图省事写了个dok_matrix，结果一亿次非零赋值变成了灾难。改成 COO 批量构造之后，运行时间从“下班前跑完”变成“泡杯咖啡就出结果”。

第四个坑是坐标索引数据类型。矩阵维度特别大时，indptr和indices的类型会影响内存，也会影响某些操作的性能。如果发现索引数组占用的内存比data还多，就要检查是否值得用更短的整数类型，或者反过来，是否需要升级到更宽的类型来防止溢出。小细节，但大数据量下不能忽视。

稀疏矩阵的学习曲线不算陡，但每一个关键选择背后都有明确的性能动机。把这些细节吃透，后面处理大规模数据时会轻松得多。