1. 引言为什么需要稀疏矩阵在科学计算和工程应用中我们常常会遇到维度极高但绝大多数元素为零的矩阵这类矩阵被称为稀疏矩阵。例如在有限元分析、社交网络图、推荐系统、自然语言处理的词袋模型以及电路仿真中矩阵的非零元素占比可能不足1%。如果使用标准的NumPy二维数组numpy.ndarray存储这样的矩阵将造成巨大的内存浪费和计算开销。SciPy的scipy.sparse模块专门为解决这一问题而设计它提供了多种高效的稀疏矩阵存储格式和丰富的线性代数运算功能。本文目标带你深入理解SciPy稀疏矩阵的核心概念、七种存储格式的适用场景并通过丰富的Python代码实例掌握其创建、转换、运算及实际应用。2. SciPy稀疏矩阵的核心存储格式SciPy主要提供了七种稀疏矩阵格式可分为两大类COO, CSR, CSC通用格式适用于大多数代数运算。LIL, DOK, BSR, DIA专用格式针对特定场景优化如增量构建、块状结构、对角线密集。下表对比了这七种格式的关键特性格式全称优点缺点典型用途COOCoordinate Format构建简单转换灵活不支持直接算术运算和切片从坐标数据快速构建矩阵CSRCompressed Sparse Row高效的行切片、矩阵向量乘法列切片慢修改结构代价高线性代数运算、迭代法求解CSCCompressed Sparse Column高效的列切片、矩阵向量乘法行切片慢修改结构代价高与CSR类似但列操作多时使用LILList of Lists支持灵活的增量修改内存开销大运算效率较低逐步构建或修改稀疏矩阵DOKDictionary of KeysO(1)时间访问单个元素内存开销大不适用于大规模矩阵测试、小规模矩阵的随机访问BSRBlock Sparse Row针对稠密子块优化缓存友好块大小需已知通用性较差有限元分析中的块状结构DIADiagonal Storage存储对角线元素极其高效仅适用于对角线附近非零的矩阵常微分方程数值解中的三对角矩阵3. 创建稀疏矩阵五种常用方法下面通过代码演示如何创建稀疏矩阵。3.1 从坐标列表创建COO格式import numpy as np from scipy import sparse 定义矩阵维度 (5x5) rows np.array([0, 1, 2, 3, 4, 0, 2, 4]) cols np.array([0, 1, 2, 3, 4, 2, 0, 2]) data np.array([1., 2., 3., 4., 5., 0.5, 1.5, 2.5]) 创建COO格式稀疏矩阵 coo_matrix sparse.coo_matrix((data, (rows, cols)), shape(5, 5)) print(COO矩阵:) print(coo_matrix.toarray()) # 转换为稠密数组查看 print(f非零元素数量: {coo_matrix.nnz}) print(f内存占用字节: {coo_matrix.data.nbytes coo_matrix.row.nbytes coo_matrix.col.nbytes})3.2 从稠密数组转换# 创建一个稠密的5x5矩阵大部分元素为0 dense_matrix np.array([ [1, 0, 0, 0, 2], [0, 0, 3, 0, 0], [0, 4, 0, 0, 0], [0, 0, 0, 5, 0], [6, 0, 0, 0, 7] ]) 转换为CSR格式最常用的运算格式 csr_from_dense sparse.csr_matrix(dense_matrix) print(\n从稠密数组转换的CSR矩阵:) print(csr_from_dense) print(行指针数组 (indptr):, csr_from_dense.indptr) print(列索引数组 (indices):, csr_from_dense.indices) print(值数组 (data):, csr_from_dense.data)3.3 创建特殊稀疏矩阵# 创建单位矩阵 identity sparse.eye(5, formatcsr) print(\n5x5单位矩阵 (CSR格式):) print(identity.toarray()) 创建对角线矩阵 diag_matrix sparse.diags([1, 2, 3], offsets[-1, 0, 1], shape(5, 5), formatcsr) print(\n三对角矩阵:) print(diag_matrix.toarray()) 创建随机稀疏矩阵10%密度 random_sparse sparse.random(5, 5, density0.1, formatcsr, random_state42) print(\n随机稀疏矩阵 (密度0.1):) print(random_sparse.toarray())4. 格式转换与属性操作不同格式间可以高效转换并可以获取矩阵的各种属性。# 从COO转换为CSR和CSC csr_matrix coo_matrix.tocsr() csc_matrix coo_matrix.tocsc() print(fCOO - CSR 转换成功: {csr_matrix.format}) print(fCOO - CSC 转换成功: {csc_matrix.format}) 获取矩阵基本信息 print(f\n矩阵形状: {csr_matrix.shape}) print(f矩阵数据类型: {csr_matrix.dtype}) print(f非零元素数量: {csr_matrix.nnz}) print(f稀疏度: {1 - csr_matrix.nnz / (csr_matrix.shape[0] * csr_matrix.shape[1]):.2%}) 转换为其他数据结构 dense_array csr_matrix.toarray() # 转为NumPy数组 coo_back csr_matrix.tocoo() # 转回COO格式 dok_matrix csr_matrix.todok() # 转为DOK格式便于随机访问 检查矩阵属性 print(f\n是否为CSR格式: {sparse.isspmatrix_csr(csr_matrix)}) print(f是否为稀疏矩阵: {sparse.issparse(csr_matrix)})5. 线性代数运算实战稀疏矩阵支持大多数NumPy风格的线性代数运算。5.1 矩阵乘法# 创建两个稀疏矩阵 A sparse.random(5, 3, density0.3, formatcsr, random_state1) B sparse.random(3, 4, density0.4, formatcsr, random_state2) 稀疏矩阵乘法结果仍为稀疏矩阵 C A.dot(B) # 等价于 A B print(矩阵A (5x3):) print(A.toarray()) print(\n矩阵B (3x4):) print(B.toarray()) print(\n矩阵乘法结果 C A B (5x4):) print(C.toarray()) print(fC的格式: {C.format}) 稀疏矩阵与稠密向量乘法 x np.random.randn(4) b C.dot(x) print(f\n稠密向量 x: {x}) print(f矩阵向量乘积 b C x: {b})5.2 求解线性系统from scipy.sparse.linalg import spsolve import time 创建一个对称正定稀疏矩阵模拟泊松方程离散化 n 100 main_diag 2 * np.ones(n) off_diag -1 * np.ones(n-1) A_poisson sparse.diags([main_diag, off_diag, off_diag], offsets[0, -1, 1], shape(n, n), formatcsr) 创建右侧向量b b_vec np.ones(n) 使用稀疏求解器 start time.time() x_solution spsolve(A_poisson, b_vec) sparse_time time.time() - start print(f稀疏求解 {n}x{n} 系统耗时: {sparse_time:.4f} 秒) 对比稠密求解小规模演示大规模会内存不足 if n 1000: A_dense A_poisson.toarray() start time.time() x_dense np.linalg.solve(A_dense, b_vec) dense_time time.time() - start print(f稠密求解 {n}x{n} 系统耗时: {dense_time:.4f} 秒) print(f稀疏求解加速比: {dense_time/sparse_time:.1f}x) print(f解的最大绝对误差: {np.max(np.abs(x_solution - x_dense)):.2e})5.3 特征值计算from scipy.sparse.linalg import eigs 计算最大特征值适用于大规模稀疏矩阵 这里用小矩阵演示 A_small sparse.random(50, 50, density0.1, formatcsr, random_state42) A_small A_small A_small.T # 使其对称 计算前3个最大特征值 eigenvalues, eigenvectors eigs(A_small, k3, whichLM) # LM: Largest Magnitude print(前3个最大特征值:) for i, val in enumerate(eigenvalues): print(f λ{i1} {val.real:.4f} {val.imag:.4f}j) 验证特征值分解: A * v ≈ λ * v for i in range(3): v eigenvectors[:, i] Av A_small.dot(v) λv eigenvalues[i] * v error np.linalg.norm(Av - λv) / np.linalg.norm(λv) print(f特征值 {i1} 的相对误差: {error:.2e})6. 实战应用图与社交网络分析稀疏矩阵是图邻接矩阵的自然表示。以下示例分析一个小型社交网络。import matplotlib.pyplot as plt import networkx as nx 创建有向图的邻接矩阵5个节点 行源节点列目标节点1表示有连接 adjacency_data np.array([1, 1, 1, 1, 1, 1]) adjacency_rows np.array([0, 1, 2, 1, 3, 4]) adjacency_cols np.array([1, 2, 3, 4, 2, 0]) adj_matrix sparse.csr_matrix( (adjacency_data, (adjacency_rows, adjacency_cols)), shape(5, 5) ) print(社交网络邻接矩阵:) print(adj_matrix.toarray()) 计算每个节点的出度和入度 out_degree adj_matrix.sum(axis1).A1 # .A1将矩阵展平为1D数组 in_degree adj_matrix.sum(axis0).A1 print(f\n节点出度: {out_degree}) print(f节点入度: {in_degree}) 计算二阶邻居朋友的朋友 two_step_neighbors adj_matrix.dot(adj_matrix) print(\n二阶连接矩阵朋友的朋友:) print(two_step_neighbors.toarray()) 可视化图 G nx.from_scipy_sparse_array(adj_matrix, create_usingnx.DiGraph) pos nx.spring_layout(G) plt.figure(figsize(6, 4)) nx.draw(G, pos, with_labelsTrue, node_colorlightblue, node_size800, arrowsize20, font_weightbold) plt.title(社交网络图可视化) plt.show()7. 性能对比与最佳实践7.1 内存占用对比import sys 创建一个1000x1000密度0.01的矩阵 n 1000 density 0.01 nnz int(n * n * density) 稀疏矩阵CSR格式 rows np.random.randint(0, n, sizennz) cols np.random.randint(0, n, sizennz) data np.random.randn(nnz) sparse_mat sparse.csr_matrix((data, (rows, cols)), shape(n, n)) 稠密矩阵仅用于对比实际可能内存不足 dense_mat sparse_mat.toarray() 计算内存占用 sparse_memory (sparse_mat.data.nbytes sparse_mat.indices.nbytes sparse_mat.indptr.nbytes) dense_memory dense_mat.nbytes print(f矩阵大小: {n}x{n}) print(f非零元素数量: {nnz} (密度: {density:.1%})) print(fCSR稀疏矩阵内存: {sparse_memory / 1024:.1f} KB) print(f稠密矩阵内存: {dense_memory / 1024 / 1024:.1f} MB) print(f内存节省比例: {dense_memory/sparse_memory:.0f}x)7.2 格式选择指南构建阶段使用COO或LIL格式它们支持灵活的增量添加。计算阶段转换为CSR行操作多或CSC列操作多格式以获得最佳性能。特殊结构对角线密集用DIA块状结构用BSR需要随机访问用DOK。通用建议CSR是最常用的通用格式大多数SciPy稀疏线性代数函数都针对CSR/CSC优化。7.3 常见陷阱与解决方案# 陷阱1对CSR矩阵进行频繁的元素修改 bad_csr sparse.random(5, 5, density0.2, formatcsr) # 错误做法直接赋值效率极低 # bad_csr[2, 3] 10 # 不推荐 正确做法先转换为LIL或DOK修改再转回CSR good_lil bad_csr.tolil() good_lil[2, 3] 10 good_csr good_lil.tocsr() 陷阱2忽略矩阵格式对切片性能的影响 csr_mat sparse.random(100, 100, density0.1, formatcsr) csc_mat csr_mat.tocsc() CSR行切片快 row_slice_csr csr_mat[10:20, :] # 快 row_slice_csc csc_mat[10:20, :] # 慢 CSC列切片快 col_slice_csc csc_mat[:, 10:20] # 快 col_slice_csr csr_mat[:, 10:20] # 慢 print(记住CSR适合行操作CSC适合列操作)8. 总结SciPy稀疏矩阵是处理大规模稀疏数据的利器。通过本文你应该掌握七种存储格式的特点与适用场景特别是COO、CSR、CSC的核心区别。五种创建方法能够从各种数据源构建稀疏矩阵。线性代数运算包括矩阵乘法、线性系统求解和特征值计算。