python的高级数组之稀疏矩阵

 

1. 稀疏矩阵的定义：

具有少量非零项的矩阵（在矩阵中，若数值0的元素数目远多于非0元素的数目，并且非0元素分布没有规律时，）则称该矩阵为稀疏矩阵；相反，为稠密矩阵。非零元素的总数比上矩阵所有元素的总数为矩阵的稠密度。

稀疏矩阵的两个动机：稀疏矩阵通常具有很大的维度，有时甚大到整个矩阵（零元素）与可用内存不想适应；另一个动机是避免零矩阵元素的运算具有更好的性能。

1. 稀疏矩阵的格式

存储矩阵的一般方法是采用二维数组，其优点是可以随机地访问每一个元素，因而能够容易实现矩阵的各种运算。对于稀疏矩阵，采用二维数组的存储方法既浪费大量的存储单元来存放零元素，又要在运算中浪费大量的时间来进行零元素的无效运算。因此必须考虑对稀疏矩阵进行压缩存储（只存储非零元素）。

scipy.sparse模块提供了许多来自于稀疏矩阵的不同存储格式。这里仅描述最为重要的格式csr、csc和lil。csr、csc是用于矩阵-矩阵和矩阵-向量运算的有效格式，lil格式用于生成和更改稀疏矩阵。python不能自动创建稀疏矩阵，所以要用scipy中特殊的命令来得到稀疏矩阵。

(1) 压缩稀疏行（csr，compressed sparse row）：或csr_matrix  按行对矩阵进行压缩的。

   csr使用了三个数组，分别为数值、行偏移（表示某一行的第一个元素在数值里面的起始偏移位置，在行偏移的最后补上矩阵总的元素个数）、列号。csr是一种编码的方式

一维数组data（数值）:有序地存储了所有的非零值，它具有与非零元素同样多数量的元素，通常由变量nnz表示。

一维数组indptr（行偏移量）：包含了证书使得indptr[i]是data中元素的索引，它是行i中的第一个非零元素。如果整个行i为零，则indptr[i]==indptr[i+1]

如初始矩阵有m行，则len(indptr)==m+1

一维数组indices（列号:）: 其使用如下方式包含列索引信息:indices[indptr[i]:indptr[i+1]]是一个具有行i中非零元素的列索引的整数数组。len(indice)==len(data)==nnz

 

备注：列索引表示数值所在的列号，从0开始。

      数组data：包含矩阵中的非零元素，以行优先的形式保存。

    行偏移：csr中行索引被压缩，没有行索引，这里用行偏移表示行索引。

实例：

             

如上图所示：data=(1,7,2,8,5,3,9,6,4)

            indices=(0,1,1,2,0,2,3,1,3)    #列索引

            indptr=(0,2,4,7,9)  #行偏移（表示某一行的第一个元素在数值里面的起始偏移位置，在行偏移的最后补上矩阵总的元素个数）

在python中使用：

import numpy as np

from scipy.sparse import csr_matrix

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])

a=csr_matrix((data, indices, indptr), shape=(3, 3)).toarray() #生成csr格式的矩阵

print(a)    #运行结果：

[[1 0 2]

 [0 0 3]

 [4 5 6]]

     解析：第i行的列索引存储在indices[indptr[i]:indptr[i+1]]中,对应的值为data[indptr[i]:indptr[i+1]]。即例如第0行的列索引为indices[0:2]=[0,2]（第i行中非零元素的列索引组成的整数数组）,值为data[0:2]=[1,2];第1行的列索引为indices[2:3]=[2],值为data[2:3]=[3]…

(2) 稀疏列矩阵csc（compressed sparse column）,用于csc格式的类型为：csc_matrix  按列对矩阵进行压缩的。

  与csr格式相比唯一的不同点是indptr和indices数组的定义，该定义与列有关。

csc格式的实例：

import numpy as np

import scipy.sparse as sp

a=np.array([[1,0,2,0],[0,0,0,0],[3,0,0,0],[1,0,0,4]])

as=sp.csc_matrix(a)

print(as)

print(as.data)

print(as.indptr)

print(as.indices)

print(as.nnz)    #运行结果：

[1 3 1 2 4]

[0 3 3 4 5]     #注意此处，同一矩阵csr格式的indptr为[0 2 2 3 5]

[0 2 3 0 3]

5

(3) 基于行的链表格式：lil（row-based linked list format）

 1. 链表稀疏格式在列表数据中以行方式存储非零元素，

列表data: data[k]是行k中的非零元素的列表。如果该行中的所有元素都为0，则它包含一个空列表。

列表rows: 是在位置k包含了在行k中的非零元素列索引列表。

lil格式的同一示例：

import numpy as np

import scipy.sparse as sp

a=np.array([[1,0,2,0],[0,0,0,0],[3,0,0,0],[1,0,0,4]])

as=sp.lil_matrix(a)

print(as.data)

print(as.rows)

print(as.nnz)  #运行结果：

[list([1, 2]) list([]) list([3]) list([1, 4])]

[list([0, 2]) list([]) list([0]) list([0, 3])]

5

 2. 用lil格式更改和切割矩阵：

lil格式最适合切片的方法，即以lil格式提取子矩阵，并通过插入非零元素来改变稀疏模式。

例如：提取

import numpy as np

import scipy.sparse as sp

a=np.array([[1,0,2,0],[0,0,0,0],[3,0,0,0],[1,0,0,4]])

as=sp.lil_matrix(a)

print(as)

bs=as[0:2,0:3]      #切片提取0，1行，0，1，2列组成的子矩阵

print(bs)

print(bs.data)

print(bs.rows)   

#运行结果：

  (0, 0)        1

  (0, 2)        2

[list([1, 2]) list([])]

[list([0, 2]) list([])]

更改：插入新的非零元素会自动更新属性

as[0,1]=17

print(as.data)

print(as.rows)

print(as.nnz)  

#结果： [list([1, 17, 2]) list([]) list([3]) list([1, 4])]

[list([0, 1, 2]) list([]) list([0]) list([0, 3])]

6

1. 生成稀疏矩阵：

numpy包的命令eye、identity、diag和rand都有其对应的稀疏矩阵，这些命令需要额外的参数来指定所得矩阵的稀疏矩阵格式。

import numpy as np

import scipy.sparse as sp

print(sp.eye(20,20,format = 'lil'))

print(sp.spdiags(np.ones((20,)),0,20,20,format = 'csr'))

print(sp.identity(20,format = 'csc'))

print(sp.rand(20,200,density=0.1,format='csr'))    #sp.rand命令需要额外的参数来描述生成随机矩阵的密度。

1. 稀疏矩阵方法

将稀疏矩阵类型转换为另一种类型和数据或数组的方法：

as.toarray  #转换稀疏矩阵类型为数组

as.tocsr

as.tocsc

as.tolil

 

#通过issparse、isspmatrix_lil、isspmatrix_csc、isspmatrix_csr等方法检查稀疏矩阵的类型。

import numpy as np

import scipy.sparse as sp

a=np.array([[1,0,2,0],[0,0,0,0],[3,0,0,0],[1,0,0,4]])

def sparse_sin(a):

    if not (sp.isspmatrix_csr(a) or sp.isspmatrix_csc(a)):

        a=a.tocsr()

        a.data=sin(a.data)

        return(a)

b=sparse_sin(a)

print(b)

 

#稀疏矩阵方法的dot，用于矩阵-矩阵或者矩阵-向量乘法运算，返回csr_matrix或numpy array

例如：import numpy as np

import scipy.sparse as sp

a=np.array([[1,0,2,0],[0,0,0,0],[3,0,0,0],[1,0,0,4]])

as=sp.csr_matrix(a)

b=np.array([1,2,3,4])

c=as.dot(b)     #结果为：[ 7  0 3  17]

print(c)

c=as.dot(as)        #结果仍为稀疏矩阵

print(c)

d=np.dot(as,b)

print(d)     #不能返回期望的结果