SciPy中稀疏矩阵的处理

今天药忘吃喽~ 2022-05-09 08:58 316阅读 0赞

在矩阵中，若数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律时，则称该矩阵为稀疏矩阵。    ——来自百度百科。

为什么会用到稀疏矩阵，最近在做协同过滤算法时，调用评分图和信任图，数据的稀疏程度达到99.9%，这样的数据存储到内存中，0会占据大量的内存，本想无所谓，但奈何内存放不下这样的数据量，无奈进行稀疏矩阵的存储与计算。记录下学习笔记。

知识来源：[Scipy Lecture Nodes 中文版][Scipy Lecture Nodes]

先举一个小栗子，展示双精度数据大小和内存之间的关系图

import numpy as np
    import matplotlib.pyplot as plt
    
    x = np.linspace(0, 1e6, 10)
    plt.plot(x, 8.0 * (x**2) / 1e6, lw=5)
    plt.xlabel('size n')
    plt.ylabel('memory [MB]')
    plt.show()

显示结果

![70][]

## 稀疏矩阵的存储： ##

稀疏矩阵作为一个矩阵，绝大多数都是0，为空。存储所有的0是浪费。

所以想办法进行压缩。节省大量内存空间。

稀疏矩阵的可视化图（样例）：

![graph.png][]![graph_g.png][]![graph_rcm.png][]

### scipy.sparse中有七类稀疏矩阵: ###

**csc\_matrix: 压缩列格式  
csr\_matrix: 压缩行格式  
bsr\_matrix: 块压缩行格式  
lil\_matrix: 列表的列表格式  
dok\_matrix: 值的字典格式  
coo\_matrix: 座标格式 (即 IJV, 三维格式)  
dia\_matrix: 对角线格式**

每一个类型适用于一些任务  
许多都利用了由Nathan Bell提供的稀疏工具 C ++ 模块

*   *  算术操作的默认实现
        
         *  通常转化为CSR
         *  为了效率而子类覆盖
     *  形状、数据类型设置/获取
     *  非0索引
     *  格式转化、与Numpy交互(toarray(), todense())
     *  ...
 *  属性:
    
     *  mtx.A - 与mtx.toarray()相同
     *  mtx.T - 转置 (与mtx.transpose()相同)
     *  mtx.H - Hermitian (列举) 转置
     *  mtx.real - 复矩阵的真部
     *  mtx.imag - 复矩阵的虚部
     *  mtx.size - 非零数 (与self.getnnz()相同)
     *  mtx.shape - 行数和列数 (元组)
 *  数据通常储存在Numpy数组中

### 对角线格式 (DIA) ###

形状 (n\_diag, length) 的密集Numpy数组的对角线。固定长度距离当离主对角线比较远时会浪费空间。\_data\_matrix的子类 (带数据属性的稀疏矩阵类)  
**每个对角线的偏移：**0 是主对角线，负偏移 = 下面，正偏移 = 上面  
**快速矩阵 \* 向量 (sparsetools)**（\*代表矩阵相乘）  
**构建器接受 :**1、密集矩阵 (数组)；2、稀疏矩阵；3、形状元组 (创建空矩阵)；4、(数据, 偏移) 元组  
**没有切片、没有单个项目访问**  
**用法 : **通过有限微分解偏微分方程；有一个迭代求解器

### 举例DIA矩阵 ###

【1】

import numpy as np
    import scipy.sparse as sparse
    import matplotlib.pyplot as plt
    
    data = np.array([[1, 2, 3, 4]]).repeat(3, axis=0)
    offsets = np.array([0, -1, 2])
    mtx = sparse.dia_matrix((data, offsets), shape=(4, 4))
    print(mtx)

![70 1][]     ![70 2][]

由此可以看出，mtx是一个4\*4的矩阵，***mtx的值，我是真不知道为什么，跪求高人指点，网友说看下面的图能看懂了，可能我脑回路有问题，我是真没看懂***。

![è¿éåå¾çæè¿°][1360851187_6918.jpg]

【2】右图为大神给的解释图，**但是我还是有点问题**

![70 3][]![70 4][]

dia\_matrix源码：

![70 5][]

【3】矩阵相乘

![70 6][]

### 列表格式(LIL) ###

**基于行的连接列表**。1、每一行是一个Python列表（排序的）非零元素的列索引；2、行存储在Numpy数组中 (dtype=np.object)；3、非零值也近似存储  
**高效增量构建稀疏矩阵**  
**构建器接受 :**1、密集矩阵 (数组)；2、稀疏矩阵；3、形状元组 (创建一个空矩阵)

### 举例LIL矩阵 ###

【1】

![70 7][]

【2】，更多的切片和索引

![70 8][]

![70 9][]

### 字典格式 (DOK) ###

**Python字典的子类：**1、键是 (行, 列) 索引元组 (不允许重复的条目)；2、值是对应的非零值  
**高效增量构建稀疏矩阵**  
**构建器支持：**1、密集矩阵 (数组)；2、稀疏矩阵；3、形状元组 (创建空矩阵)  
**高效 O(1) 对单个元素的访问**  
**一旦创建完成后可以被高效转换为coo\_matrix**  
**算术很慢** **(循环用dict.iteritems())**  
**用法：**当稀疏模式是未知的假设或改变时

![70 10][]

### 座标格式 (COO) ###

**称为 ‘ijv’ 或 ‘triplet’ 格式：**1、三个NumPy数组: row, col, data；2、data\[i\]是在 (row\[i\], col\[i\]) 位置的值；3、允许重复值  
**构建稀疏矩阵的高速模式**  
**构建器接受：**1、密集矩阵 (数组)；2、稀疏矩阵；3、形状元组 (创建空数组)；4、(data, ij)元组  
**可以与CSR/CSC格式非常快的互相转换**  
**快速的矩阵 \* 向量 (sparsetools)**  
**快速而简便的逐项操作：**直接操作数据数组 (快速NumPy机制)  
**没有切片，没有算术** (直接)  
**使用：1、**在各种稀疏格式间的灵活转换；2、当转化到其他形式 (通常是 CSR 或 CSC), 重复的条目被加总到一起

case【1】

![70 11][]

### 行格式 (CSR) ###

**面向行**  
**三个Numpy数组: indices, indptr, data**

1.  indices是列索引的数组
2.  data是对应的非零值数组
3.  indptr指向indices和data开始的行
4.  长度是n\_row + 1, 最后一个项目 = 值数量 = indices和data的长度
5.  第i行的非零值是列索引为indices\[indptr\[i\]:indptr\[i+1\]\]的data\[indptr\[i\]:indptr\[i+1\]\]
6.  项目 (i, j) 可以通过data\[indptr\[i\]+k\], k是j在indices\[indptr\[i\]:indptr\[i+1\]\]的位置来访问

**快速矩阵向量相乘和其他算术 (sparsetools)**  
**构建器接受：**1、密集矩阵 (数组)；2、稀疏矩阵；3、形状元组 (创建空矩阵)；4、(data, ij) 元组；5、(data, indices, indptr) 元组  
**高效行切片，面向行的操作  
较慢的列切片，改变稀疏结构代价昂贵**  
**用途：**

1.  实际计算 (大多数线性求解器都支持这个格式)

### case ###

**【1】用(data, indices, indptr)元组创建:**

![70 12][]

解释下mtx.indptr，mtx.indices

\# 对于第i行，非0数据列是indices\[indptr\[i\]:indptr\[i+1\]\] 数据是data\[indptr\[i\]:indptr\[i+1\]\]  
\# 第0行，有非0的数据列是indices\[indptr\[0\]:indptr\[1\]\] = indices\[0:2\] = \[0,2\]  
\# 数据是data\[indptr\[0\]:indptr\[1\]\] = data\[0:2\] = \[1,2\],所以在第0行第0列是1，第0行第2列是2  
\# 第1行，有非0的数据列是indices\[indptr\[1\]:indptr\[2\]\] = indices\[2:3\] = \[2\]  
\# 数据是data\[indptr\[1\]:indptr\[2\] = data\[2:3\] = \[3\],所以在第1行第2列是3  
\# 第2行，有非0的数据列是indices\[indptr\[2\]:indptr\[3\]\] = indices\[3:6\] = \[0,1,2\]  
\# 数据是data\[indptr\[2\]:indptr\[3\]\] = data\[3:6\] = \[4,5,6\],所以在第2行第0列是4，第2行第1行是5，第2行第2列是6

**【2】用`(data, ij)`元组创建:**

![70 13][]

### 列格式 (CSC) (与CSR类似) ###

**三个Numpy数组: indices、indptr、data**

1.  indices是行
2.  索引的数组
3.  data是对应的非零值
4.  indptr指向indices和data开始的列
5.  长度是n\_col + 1, 最后一个条目 = 值数量 = indices和data的长度
6.  第i列的非零值是行索引为indices\[indptr\[i\]:indptr\[i+1\]\]的data\[indptr\[i\]:indptr\[i+1\]\]
7.  项目 (i, j) 可以作为data\[indptr\[j\]+k\]访问, k是i在indices\[indptr\[j\]:indptr\[j+1\]\]的位置

**快速矩阵向量相乘和其他算术 (sparsetools)**  
**构建器接受：**1、密集矩阵 (数组)；2、稀疏矩阵；3、形状元组 (创建空矩阵)；4、(data, ij)元组；5、(data, indices, indptr)元组  
**高效列切片、面向列的操作  
较慢的行切片、改变稀疏结构代价昂贵  
用途:**

1.  实际计算 (巨大多数线性求解器支持这个格式)

### case ###

**【1】用(data, indices, indptr)元组创建:**

![70 14][]

解释下mtx.indptr，mtx.indices

\# 对于第i列，非0数据行是indices\[indptr\[i\]:indptr\[i+1\]\] 数据是data\[indptr\[i\]:indptr\[i+1\]\]  
\# 第0列，有非0的数据行是indices\[indptr\[0\]:indptr\[1\]\] = indices\[0:2\] = \[0,2\]  
\# 数据是data\[indptr\[0\]:indptr\[1\]\] = data\[0:2\] = \[1,2\],所以在第0列第0行是1，第0列第2行是2  
\# 第1列，有非0的数据行是indices\[indptr\[1\]:indptr\[2\]\] = indices\[2:3\] = \[2\]  
\# 数据是data\[indptr\[1\]:indptr\[2\] = data\[2:3\] = \[3\],所以在第1列第2行是3  
\# 第2列，有非0的数据行是indices\[indptr\[2\]:indptr\[3\]\] = indices\[3:6\] = \[0,1,2\]  
\# 数据是data\[indptr\[2\]:indptr\[3\]\] = data\[3:6\] = \[4,5,6\],所以在第2列第0行是4，第2列第1行是5，第2列第2行是6

**【2】用`(data, ij)`元组创建:**

![70 15][]

### 块压缩行格式 (BSR) ###

**CSR带有密集的固定形状的子矩阵而不是纯量的项目**  
**块大小(R, C)必须可以整除矩阵的形状(M, N)**  
**三个Numpy数组: indices、indptr、data**

1.  indices是每个块列索引的数组
2.  data是形状为(nnz, R, C)对应的非零值

**快速矩阵向量相乘和其他的算术 (sparsetools)**  
**构建器接受：**1、密集矩阵 (数组)；2、稀疏矩阵；3、形状元组 (创建空的矩阵)；4、(data, ij)元组；5、(data, indices, indptr)元组  
**许多对于带有密集子矩阵的稀疏矩阵算术操作比CSR更高效很多**  
**用途:**

1.  类似CSR

### case ###

【1】**用`(data, ij)`元组创建:**

![70 16][]

**【2】用块大小`(data, indices, indptr)`的元组创建:**

![70 17][]

[Scipy Lecture Nodes]: https://wizardforcel.gitbooks.io/scipy-lecture-notes/content/10.html
[70]: /images/20220509/0a58fe2ede784044ae28646513823a3f.png
[graph.png]: http://www.scipy-lectures.org/_images/graph.png
[graph_g.png]: http://scipy-lectures.github.io/_images/graph_g.png
[graph_rcm.png]: /images/20220509/249a6d89a5f14e12ba26e4de523447ec.png
[70 1]: https://img-blog.csdn.net/20181005145858428?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hhbmc5MTY=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 2]: /images/20220509/ae8735c3878c4098b40990a1e6884b74.png
[1360851187_6918.jpg]: https://img-my.csdn.net/uploads/201302/14/1360851187_6918.jpg
[70 3]: https://img-blog.csdn.net/20181005212946722?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hhbmc5MTY=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 4]: /images/20220509/4a7388763f464d859f1a0f58cb77932c.png
[70 5]: /images/20220509/9e84815c57a84352af5311bae80b3d2e.png
[70 6]: /images/20220509/119872151bcd49158318c79f143c75b7.png
[70 7]: /images/20220509/ddb72eca6bfd4c55b2727d3c1d290191.png
[70 8]: /images/20220509/38d91c48d75448d6950e4f7f9f6fcce3.png
[70 9]: /images/20220509/60454d79fc104408aaa42cc934a4a760.png
[70 10]: /images/20220509/8d7845c8556d40d1b0eb2d6077eccdd1.png
[70 11]: /images/20220509/49468a67018747bca73c4103bb110feb.png
[70 12]: /images/20220509/426b652bb52f4872ab975e1654702e5c.png
[70 13]: /images/20220509/f20adb09ee99422ba434698d5e65a38d.png
[70 14]: /images/20220509/e3e324cf76824f8585936b917174be66.png
[70 15]: /images/20220509/6a90ef5f320c4a0d9a81122ec91ca120.png
[70 16]: /images/20220509/ea526ba1b39a44beb8e9ec58f35fc64b.png
[70 17]: /images/20220509/b69df928ff1641d4b74ff7b5927caef7.png