深入简出理解scipy.sparse.csr_matrix和scipy.sparse.csc_matrix

一、导入

  在用python进行科学运算时，常常需要把一个稀疏的np.array压缩，这时候就用到scipy库中的sparse.csr_matrix函数和sparse.csc_matric函数。

  其中csr_matrix：Compressed Sparse Row marix,而csc_matric：Compressed Sparse Column marix。

---

二、引例

这里先放上一个官网的经典例子：

from scipy.sparse import csr_matrix
import numpy as np
indptr = np.array([0,2,3,6])
indices = np.array([0,2,2,0,1,2])
data = np.array([1,2,3,4,5,6])
csr_matrix_0 = csr_matrix((data,indices,indptr),shape=(3,3))
print(csr_matrix_0.toarray())

输出：

这段代码什么意思？
我们不妨按照下面的思路来解释。

---

三、解释csr_matrix

1. 存储元素

  假设我们有一个稀疏矩阵X：
$X=\begin{bmatrix}1 & 0 & 2 \\0 & 0 & 3 \\4 & 5 & 6 \end{bmatrix}$X=⎣⎡​104​005​236​⎦⎤​

  我们现在要把它压缩处理。

  很容易想到，如果我们只存储非零的元素，就可以实现压缩。上面矩阵的非零元素一共有6个，分别是1、2、3、4、5、6（按行的顺序读），这里我们先用一个名为data的np.array来把这6个数据存储下来：data = np.array([1,2,3,4,5,6])。

  但是这个矩阵是3行3列的（共9个元素），我们怎么知道上面那6个元素该放在什么位置？于是我们需要一个索引来指引我们非零元素的位置。

2. 列索引

  我们现在创建一个列索引来解决这个问题，对于上面6个非零元素，我们依次找到它在稀疏矩阵中位于哪一列：

元素	所在列
1	0
2	2
3	2
4	0
5	1
6	2

  这里，我们用一个名为indices的np.array来把这6个元素所在列位置存储下来：indices = np.array([0,2,2,0,1,2])。

  可是光知道列，我们怎么知道某个元素到底在哪一行？难道要再用一个np.array来记录元素所在行位置？从而使用“每个元素分别存储其行列位置”来达到数据压缩。这效果能好吗？

  其实，这样操作可以压缩，但并不是什么好的压缩方法，因为它还有数据冗余。比如就拿上面那个例子来说，对于稀疏矩阵
$X=\begin{bmatrix}1 & 0 & 2 \\0 & 0 & 3 \\4 & 5 & 6 \end{bmatrix}$X=⎣⎡​104​005​236​⎦⎤​

  它的第0行有2个元素，分别是1和2，那么我们分别存储它们的行列信息[0,0]和[0,2]，这里就会发现，行的信息数据冗余了（第0行这个数据我们存储了两遍）。如果拿第2行来举例，效果会更直观，对于元素4，5，6，他们都是同一行，那么如果这样存储，那么“第2行”这个信息相会被重复存储3遍，何其浪费！

  所以这种想法是不可取的，我们需要更高效地压缩数据。

3. 游标指针

  据以上信息，我们可以发现，无论是data还是indices，它们的个数相同，而且一一对应，即：data[0]=1对应的列就是indices[0]=0，data[1]=2对应的列就是indices[1]=2等。

  那么我们就产生了想法，因为元素是按顺序排下来的，且两者一一对应，那么如果我们能把indices / data 分几份（即几行），这样就可以解决问题了。（因为indices和data是一一对应的，所以分割位置相同即可以）

  比如说，现在data = np.array([1,2,3,4,5,6])，我们这样分[1 2 | 3 | 4 5 6]。其中第0组[1 2]就是第0行的元素，第1组[3]就是第1行的元素，第2组[4 5 6]就是第2行的元素。
  与此同步，indices = np.array([0,2,2,0,1,2])，这样分[0 2 | 2 | 0 1 2]。其中第0组[0 2]就是第0行元素的列索引，第1组[2]就是第1行元素的列索引，第2组[0 1 2]就是第2行元素的列索引。

  那么如何对indices / data进行分组呢？肯定是用一个东西来存储indices的分组下标，这里很容易想到用游标指针来指示。用一个名为indptr的np.array来存储每组（每行）起始元素的位置即可（这是常用的套路）。

  如下：

[pos]      0 1   2   3 4 5
data   =  [1 2 | 3 | 4 5 6]
indices = [0 2 | 2 | 0 1 2]

  第0组（第0行）在 indices / data 中的起始下标[0]；
  第1组（第1行）在 indices / data 中的起始下标[2]；
  第2组（第2行）在 indices / data 中的起始下标[3]。

  这里需要注意的是，我们还没有记录元素总个数，不然我们怎么知道有几个元素。所以我们把元素个数这个属性找一个地方安置，如下：把这些下标索引拼合起来再加上一个元素个数，得到一个名为indptr的np.array存储起来：indptr = np.array([0,2,3,6])。

4. 调用函数

indptr = np.array([0,2,3,6])
indices = np.array([0,2,2,0,1,2])
data = np.array([1,2,3,4,5,6])

  拥有了以上三个np.array信息的我们足以记录压缩矩阵了。

  这就是函数csr_matrix((data,indices,indptr),shape=(3,3))各个参数的含义，其中shape=(3,3)表示这是一个3行3列稀疏矩阵。

5. 总述

$X=\begin{bmatrix}1 & 0 & 2 \\0 & 0 & 3 \\4 & 5 & 6 \end{bmatrix}$X=⎣⎡​104​005​236​⎦⎤​

[pos]      0 1   2   3 4 5
data   =  [1 2 | 3 | 4 5 6]
indices = [0 2 | 2 | 0 1 2]

[pos]     0 1 2 3
indptr = [0,2,3,6]

第0行：

（1）起始元素游标indptr[0]=0，终止元素游标indptr[1]-1=2-1=1；
（2）对应indices中第0个到第1个元素，为[0,2]；
（3）对应data中第0个到第1个元素，为[1,2]；
（4）含义：第0行在第0、2列有非零元素存在，且元素分别为1、2。

第1行：

（1）起始元素游标indptr[1]=2，终止元素游标indptr[2]-1=3-1=2；
（2）对应indices中第2个到第2个元素，为[2]；
（3）对应data中第2个到第2个元素，为[3]；
（4）含义：第1行在第2列有非零元素存在，且元素为3。

第2行：

（1）起始元素游标indptr[2]=3，终止元素游标indptr[3]-1=6-1=5；
（2）对应indices中第3个到第5个元素，为[0,1,2]；
（3）对应data中第3个到第5个元素，为[4,5,6]；
（4）含义：第2行在第0、1、2列有非零元素存在，且元素分别为4、5、6。

第i行：

（1）起始元素游标indptr[i]，终止元素游标indptr[i+1]-1；
（2）对应indices中第indptr[i]个到第indptr[i+1]-1个元素，即indices[ indptr[i] : indptr[i+1] ]；（这里不用减1，Python基础知识）
（3）对应data中第indptr[i]个到第indptr[i+1]-1个元素，即data[ indptr[i] : indptr[i+1] ]；（同上不用减1）
（4）含义：第 i 行在第indices[ indptr[i] : indptr[i+1] ]列有非零元素存在，且元素分别为data[ indptr[i] : indptr[i+1] ]。

6. 官方解释

  第i行的列下标存储在indices[indptr[i]:indptr[i+1]]，并且它们对应的值存储在data[indptr[i]:indptr[i+1]]。如果没有提供shape参数，那么矩阵的维数就会根据 index arrays 自动推断。

  其实这个scipy.sparse.csr_matrix函数还有一些其他用法，这里就不再赘述，可以参考源码：

---

四、解释csc_matrix

  同样，函数csc_matrix也是大致相同的意思，只不过把csr_matrix函数的行列交换了一下。

  即对于同一稀疏矩阵：
$X=\begin{bmatrix}1 & 0 & 2 \\0 & 0 & 3 \\4 & 5 & 6 \end{bmatrix}$X=⎣⎡​104​005​236​⎦⎤​

  csc_matrix函数：
  存储元素（按列的顺序读元素）：data = np.array([1,4,5,2,3,6])
  列索引：indices = np.array([0,2,2,0,1,2])
  游标指针：indptr = np.array([0,2,3,6])

  套路一样，不再详细解释。

---

五. 例题

例1：
使用csr_matrix函数压缩稀疏矩阵
$X=\begin{bmatrix}1 & 0 & 2 \\0 & 0 & 3 \\4 & 5 & 6 \end{bmatrix}$X=⎣⎡​104​005​236​⎦⎤​

from scipy.sparse import csr_matrix,csc_matrix
import numpy as np
indptr = np.array([0,2,3,6])
indices = np.array([0,2,2,0,1,2])
data = np.array([1,2,3,4,5,6])
csr_matrix_0 = csr_matrix((data,indices,indptr),shape=(3,3))
print(csr_matrix_0.toarray())

例2：
使用csc_matrix函数压缩稀疏矩阵
$X=\begin{bmatrix}1 & 0 & 2 \\0 & 0 & 3 \\4 & 5 & 6 \end{bmatrix}$X=⎣⎡​104​005​236​⎦⎤​

from scipy.sparse import csr_matrix,csc_matrix
import numpy as np
indptr = np.array([0,2,3,6])
indices = np.array([0,2,2,0,1,2])
data = np.array([1,4,5,2,3,6])
csr_matrix_0 = csc_matrix((data,indices,indptr),shape=(3,3))
print(csr_matrix_0.toarray())

例3：
使用csc_matrix函数压缩稀疏矩阵
$X=\begin{bmatrix}1 & 0 & 4 \\0 & 0 & 5 \\2 & 3 & 6 \end{bmatrix}$X=⎣⎡​102​003​456​⎦⎤​

from scipy.sparse import csr_matrix,csc_matrix
import numpy as np
indptr = np.array([0,2,3,6])
indices = np.array([0,2,2,0,1,2])
data = np.array([1,2,3,4,5,6])
csr_matrix_0 = csc_matrix((data,indices,indptr),shape=(3,3))
print(csr_matrix_0.toarray())

注意：对比例1、例2和例3，发现当矩阵是方阵的时候，两个函数创建矩阵的转置关系。

例4：
使用csr_matrix函数压缩稀疏矩阵
$X=\begin{bmatrix}0 & 0 & 6 & 0 & 2 & 0 \\1 & 0 & 0 & 8 & 0 & 4 \\0 & 5 & 3 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 1 & 1 \end{bmatrix}$X=⎣⎢⎢⎡​0100​0050​6030​0801​2001​0401​⎦⎥⎥⎤​

from scipy.sparse import csr_matrix,csc_matrix
import numpy as np
indptr = np.array([0,2,5,7,10])
indices = np.array([2,4,0,3,5,1,2,3,4,5])
data = np.array([6,2,1,8,4,5,3,1,1,1])
csr_matrix_0 = csr_matrix((data,indices,indptr),shape=(4,6))
print(csr_matrix_0.toarray())

例5：
使用csc_matrix函数压缩稀疏矩阵
$X=\begin{bmatrix}0 & 0 & 6 & 0 & 2 & 0 \\1 & 0 & 0 & 8 & 0 & 4 \\0 & 5 & 3 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 1 & 1 \end{bmatrix}$X=⎣⎢⎢⎡​0100​0050​6030​0801​2001​0401​⎦⎥⎥⎤​

from scipy.sparse import csr_matrix,csc_matrix
import numpy as np
indptr = np.array([0,1,2,4,6,8,10])
indices = np.array([1,2,0,2,1,3,0,3,1,3])
data = np.array([1,5,6,3,8,1,2,1,4,1])
csr_matrix_0 = csc_matrix((data,indices,indptr),shape=(4,6))
print(csr_matrix_0.toarray())