【Numpy核心编程攻略：Python数据处理、分析详解与科学计算】2.23 稀疏矩阵：CSR格式的存储与运算

2.23 稀疏矩阵：CSR格式的存储与运算

目录

2.23.1 稀疏矩阵存储格式

2.23.1.1 稀疏矩阵的定义

稀疏矩阵是一个大多数元素为零的矩阵。在实际应用中，稀疏矩阵广泛存在于科学计算、机器学习、图论等领域。直接存储稀疏矩阵会浪费大量的内存空间，因此需要使用高效的存储格式。

2.23.1.2 常见的稀疏矩阵存储格式

• COO (Coordinate) 格式：通过行索引、列索引和值来存储非零元素。
• CSR (Compressed Sparse Row) 格式：通过行指针、列索引和值来存储非零元素。
• CSC (Compressed Sparse Column) 格式：通过列指针、行索引和值来存储非零元素。
• DOK (Dictionary of Keys) 格式：通过字典存储非零元素。
• LIL (List of Lists) 格式：通过链表存储非零元素。

2.23.2 CSR格式解析

2.23.2.1 CSR格式的基本概念

CSR 格式通过三个数组来存储稀疏矩阵：

1. data：存储所有非零元素的值。
2. indices：存储所有非零元素的列索引。
3. indptr：存储每一行第一个非零元素在 data 和 indices 中的位置。

2.23.2.2 CSR格式的表示

假设有一个稀疏矩阵：

[ A = \begin{bmatrix}
1 & 0 & 0 \
0 & 2 & 3 \
0 & 0 & 4
\end{bmatrix} ]

其 CSR 格式表示为：

• data = [1, 2, 3, 4]
• indices = [0, 1, 2, 2]
• indptr = [0, 1, 3, 4]

2.23.2.3 代码实现

import numpy as np
from scipy.sparse import csr_matrix

# 定义稀疏矩阵
data = [1, 2, 3, 4]  # 非零元素的值
indices = [0, 1, 2, 2]  # 非零元素的列索引
indptr = [0, 1, 3, 4]  # 每一行第一个非零元素的位置

# 创建 CSR 格式的稀疏矩阵
A = csr_matrix((data, indices, indptr), shape=(3, 3))

# 输出稀疏矩阵
print(A.toarray())  # 输出: [[1 0 0]
                   #        [0 2 3]
                   #        [0 0 4]]

2.23.3 矩阵向量乘法优化

2.23.3.1 矩阵向量乘法的基本原理

矩阵向量乘法是线性代数中的基本操作，对于稀疏矩阵，可以利用 CSR 格式进行优化，以减少计算量和内存占用。

2.23.3.2 CSR格式下的矩阵向量乘法

假设稀疏矩阵 ( A ) 为 CSR 格式，向量 ( x ) 为密集向量，矩阵向量乘法 ( y = A \cdot x ) 的步骤如下：

1. 遍历每一行的非零元素。
2. 计算每一行的贡献值。
3. 累加所有行的贡献值，得到最终结果 ( y )。

2.23.3.3 代码实现

import numpy as np
from scipy.sparse import csr_matrix

# 定义稀疏矩阵
data = [1, 2, 3, 4]  # 非零元素的值
indices = [0, 1, 2, 2]  # 非零元素的列索引
indptr = [0, 1, 3, 4]  # 每一行第一个非零元素的位置

# 创建 CSR 格式的稀疏矩阵
A = csr_matrix((data, indices, indptr), shape=(3, 3))

# 定义密集向量 x
x = np.array([1, 2, 3])

# 计算矩阵向量乘法 y = A * x
y = A.dot(x)  # 计算矩阵向量乘法

# 输出结果
print(y)  # 输出: [1 8 12]

2.23.3.4 优化效果

使用 CSR 格式进行矩阵向量乘法可以显著减少计算时间和内存占用，特别是在矩阵中非零元素占比很小的情况下。

import time

# 生成一个大型稀疏矩阵
n = 10000
m = 10000
density = 0.01  # 非零元素占比
data = np.random.rand(int(n * m * density))  # 非零元素的值
indices = np.random.randint(0, m, size=int(n * m * density))  # 非零元素的列索引
indptr = np.zeros(n + 1, dtype=int)  # 每一行第一个非零元素的位置
for i in range(n):
    start = indptr[i]
    end = start + np.random.binomial(m, density)
    indptr[i + 1] = end

A = csr_matrix((data, indices, indptr), shape=(n, m))

# 生成一个密集向量 x
x = np.random.rand(m)

# 测量密集矩阵向量乘法的时间
dense_A = A.toarray()
start_time = time.time()
y_dense = dense_A.dot(x)
dense_time = time.time() - start_time

# 测量稀疏矩阵向量乘法的时间
start_time = time.time()
y_sparse = A.dot(x)
sparse_time = time.time() - start_time

# 比较时间
print(f"密集矩阵向量乘法时间: {dense_time:.4f} 秒")
print(f"稀疏矩阵向量乘法时间: {sparse_time:.4f} 秒")

2.23.4 与SciPy的协同

2.23.4.1 SciPy 中的稀疏矩阵模块

SciPy 提供了 scipy.sparse 模块，用于处理稀疏矩阵。该模块支持多种稀疏矩阵存储格式，并提供了高效的稀疏矩阵运算方法。

2.23.4.2 代码实现

2.23.4.2.1 创建稀疏矩阵

import numpy as np
from scipy.sparse import csr_matrix

# 生成一个大型稀疏矩阵
n = 10000
m = 10000
density = 0.01  # 非零元素占比
data = np.random.rand(int(n * m * density))  # 非零元素的值
indices = np.random.randint(0, m, size=int(n * m * density))  # 非零元素的列索引
indptr = np.zeros(n + 1, dtype=int)  # 每一行第一个非零元素的位置
for i in range(n):
    start = indptr[i]
    end = start + np.random.binomial(m, density)
    indptr[i + 1] = end

# 创建 CSR 格式的稀疏矩阵
A = csr_matrix((data, indices, indptr), shape=(n, m))

2.23.4.2.2 稀疏矩阵的基本操作

# 稀疏矩阵的形状
print(f"稀疏矩阵的形状: {A.shape}")

# 稀疏矩阵的非零元素数量
print(f"稀疏矩阵的非零元素数量: {A.nnz}")

# 稀疏矩阵的密度
print(f"稀疏矩阵的密度: {A.nnz / (A.shape[0] * A.shape[1])}")

# 稀疏矩阵的行指针
print(f"稀疏矩阵的行指针: {A.indptr}")

# 稀疏矩阵的列索引
print(f"稀疏矩阵的列索引: {A.indices}")

# 稀疏矩阵的非零元素值
print(f"稀疏矩阵的非零元素值: {A.data}")

2.23.4.3 稀疏矩阵的运算

2.23.4.3.1 矩阵加法

# 创建另一个稀疏矩阵
B = csr_matrix((data, indices, indptr), shape=(n, m))

# 矩阵加法
C = A + B

# 输出结果
print(f"稀疏矩阵 C 的形状: {C.shape}")
print(f"稀疏矩阵 C 的非零元素数量: {C.nnz}")

2.23.4.3.2 矩阵乘法

# 矩阵乘法
D = A.dot(B)

# 输出结果
print(f"稀疏矩阵 D 的形状: {D.shape}")
print(f"稀疏矩阵 D 的非零元素数量: {D.nnz}")

2.23.5 社交网络分析案例

2.23.5.1 案例背景

社交网络中，用户之间的关系可以用图来表示，图的邻接矩阵通常是稀疏矩阵。通过使用 CSR 格式的稀疏矩阵，可以更高效地进行图的操作和分析。

2.23.5.2 生成随机社交网络

import numpy as np
from scipy.sparse import csr_matrix

# 生成一个随机社交网络
n = 1000  # 用户数量
density = 0.01  # 关系密度

# 生成稀疏矩阵
data = np.ones(int(n * n * density))  # 非零元素的值
indices = np.random.randint(0, n, size=int(n * n * density))  # 非零元素的列索引
indptr = np.zeros(n + 1, dtype=int)  # 每一行第一个非零元素的位置
for i in range(n):
    start = indptr[i]
    end = start + np.random.binomial(n, density)
    indptr[i + 1] = end

# 创建 CSR 格式的稀疏矩阵
A = csr_matrix((data, indices, indptr), shape=(n, n))

2.23.5.3 图的操作

2.23.5.3.1 邻接矩阵转图

import networkx as nx

# 将稀疏矩阵转换为图
G = nx.from_scipy_sparse_matrix(A)

# 绘制图
plt.figure(figsize=(10, 10))
nx.draw(G, node_size=50, with_labels=False)
plt.title('Social Network Graph')
plt.show()

2.23.5.3.2 中心性分析

# 计算度中心性
degree_centrality = nx.degree_centrality(G)

# 计算介数中心性
betweenness_centrality = nx.betweenness_centrality(G)

# 输出中心性分析结果
print(f"度中心性: {degree_centrality}")
print(f"介数中心性: {betweenness_centrality}")

2.23.6 内存占用对比

2.23.6.1 测量内存占用

import sys

# 生成一个大型稀疏矩阵
n = 10000
m = 10000
density = 0.01  # 非零元素占比
data = np.random.rand(int(n * m * density))  # 非零元素的值
indices = np.random.randint(0, m, size=int(n * m * density))  # 非零元素的列索引
indptr = np.zeros(n + 1, dtype=int)  # 每一行第一个非零元素的位置
for i in range(n):
    start = indptr[i]
    end = start + np.random.binomial(m, density)
    indptr[i + 1] = end

# 创建 CSR 格式的稀疏矩阵
A = csr_matrix((data, indices, indptr), shape=(n, m))

# 生成相同大小的密集矩阵
dense_A = A.toarray()

# 测量内存占用
sparse_memory = sys.getsizeof(A.data) + sys.getsizeof(A.indices) + sys.getsizeof(A.indptr)
dense_memory = sys.getsizeof(dense_A)

print(f"稀疏矩阵的内存占用: {sparse_memory} 字节")
print(f"密集矩阵的内存占用: {dense_memory} 字节")
print(f"内存占用比例: {sparse_memory / dense_memory:.2%}")

2.23.6.2 结果分析

在上述代码中，稀疏矩阵的内存占用远小于密集矩阵的内存占用。特别是在矩阵的非零元素占比很小的情况下，使用稀疏矩阵可以显著节省内存。

2.23.7 总结与参考文献

2.23.7.1 总结

本文详细介绍了稀疏矩阵的存储格式，特别是 CSR 格式，以及如何在 NumPy 和 SciPy 中进行稀疏矩阵的运算。通过社交网络分析案例和内存占用对比，展示了 CSR 格式在实际应用中的高效性和优势。希望本文对您理解稀疏矩阵及其应用有所帮助。

2.23.7.2 参考文献

这篇文章包含了详细的原理介绍、代码示例、源码注释以及案例等。希望这对您有帮助。如果有任何问题请随私信或评论告诉我。