稀疏矩阵

Scipy 提供了稀疏矩阵的支持（scipy.sparse）。

稀疏矩阵主要使用 位置 + 值 的方法来存储矩阵的非零元素，根据存储和使用方式的不同，有如下几种类型的稀疏矩阵：

类型 描述
bsr_matrix(arg1[, shape, dtype, copy, blocksize]) Block Sparse Row matrix
coo_matrix(arg1[, shape, dtype, copy]) A sparse matrix in COOrdinate format.
csc_matrix(arg1[, shape, dtype, copy]) Compressed Sparse Column matrix
csr_matrix(arg1[, shape, dtype, copy]) Compressed Sparse Row matrix
dia_matrix(arg1[, shape, dtype, copy]) Sparse matrix with DIAgonal storage
dok_matrix(arg1[, shape, dtype, copy]) Dictionary Of Keys based sparse matrix.
lil_matrix(arg1[, shape, dtype, copy]) Row-based linked list sparse matrix

在这些存储格式中：

  • COO 格式在构建矩阵时比较高效
  • CSC 和 CSR 格式在乘法计算时比较高效

构建稀疏矩阵

In [1]:

  1. from scipy.sparse import *
  2. import numpy as np

创建一个空的稀疏矩阵：

In [2]:

  1. coo_matrix((2,3))

Out[2]:

  1. <2x3 sparse matrix of type '<type 'numpy.float64'>'
  2.     with 0 stored elements in COOrdinate format>

也可以使用一个已有的矩阵或数组或列表中创建新矩阵：

In [3]:

  1. A = coo_matrix([[1,2,0],[0,0,3],[4,0,5]])
  2. print A
  1.   (0, 0)    1
  2.   (0, 1)    2
  3.   (1, 2)    3
  4.   (2, 0)    4
  5.   (2, 2)    5

不同格式的稀疏矩阵可以相互转化：

In [4]:

  1. type(A)

Out[4]:

  1. scipy.sparse.coo.coo_matrix

In [5]:

  1. B = A.tocsr()
  2. type(B)

Out[5]:

  1. scipy.sparse.csr.csr_matrix

可以转化为普通矩阵：

In [6]:

  1. C = A.todense()
  2. C

Out[6]:

  1. matrix([[1, 2, 0],
  2.         [0, 0, 3],
  3.         [4, 0, 5]])

与向量的乘法：

In [7]:

  1. v = np.array([1,0,-1])
  2. A.dot(v)

Out[7]:

  1. array([ 1, -3, -1])

还可以传入一个 (data, (row, col)) 的元组来构建稀疏矩阵：

In [8]:

  1. I = np.array([0,3,1,0])
  2. J = np.array([0,3,1,2])
  3. V = np.array([4,5,7,9])
  4. A = coo_matrix((V,(I,J)),shape=(4,4))

In [9]:

  1. print A
  1.   (0, 0)    4
  2.   (3, 3)    5
  3.   (1, 1)    7
  4.   (0, 2)    9

COO 格式的稀疏矩阵在构建的时候只是简单的将坐标和值加到后面，对于重复的坐标不进行处理：

In [10]:

  1. I = np.array([0,0,1,3,1,0,0])
  2. J = np.array([0,2,1,3,1,0,0])
  3. V = np.array([1,1,1,1,1,1,1])
  4. B = coo_matrix((V,(I,J)),shape=(4,4))
  5. print B
  1.   (0, 0)    1
  2.   (0, 2)    1
  3.   (1, 1)    1
  4.   (3, 3)    1
  5.   (1, 1)    1
  6.   (0, 0)    1
  7.   (0, 0)    1

转换成 CSR 格式会自动将相同坐标的值合并：

In [11]:

  1. C = B.tocsr()
  2. print C
  1.   (0, 0)    3
  2.   (0, 2)    1
  3.   (1, 1)    2
  4.   (3, 3)    1

求解微分方程

In [12]:

  1. from scipy.sparse import lil_matrix
  2. from scipy.sparse.linalg import spsolve
  3. from numpy.linalg import solve, norm
  4. from numpy.random import rand

构建 1000 x 1000 的稀疏矩阵：

In [13]:

  1. A = lil_matrix((1000, 1000))
  2. A[0, :100] = rand(100)
  3. A[1, 100:200] = A[0, :100]
  4. A.setdiag(rand(1000))

转化为 CSR 之后，用 spsolve 求解 $Ax=b$：

In [14]:

  1. A = A.tocsr()
  2. b = rand(1000)
  3. x = spsolve(A, b)

转化成正常数组之后求解：

In [15]:

  1. x_ = solve(A.toarray(), b)

查看误差：

In [16]:

  1. err = norm(x-x_)
  2. err

Out[16]:

  1. 6.4310987107687431e-13

sparse.find 函数

返回一个三元组，表示稀疏矩阵中非零元素的 (row, col, value)：

In [17]:

  1. from scipy import sparse
  2.  
  3. row, col, val = sparse.find(C)
  4. print row, col, val
  1. [0 0 1 3] [0 2 1 3] [3 1 2 1]

sparse.issparse 函数

查看一个对象是否为稀疏矩阵：

In [18]:

  1. sparse.issparse(B)

Out[18]:

  1. True

或者

In [19]:

  1. sparse.isspmatrix(B.todense())

Out[19]:

  1. False

还可以查询是否为指定格式的稀疏矩阵：

In [20]:

  1. sparse.isspmatrix_coo(B)

Out[20]:

  1. True

In [21]:

  1. sparse.isspmatrix_csr(B)

Out[21]:

  1. False

原文: https://nbviewer.jupyter.org/github/lijin-THU/notes-python/blob/master/04-scipy/04.08-sparse-matrix.ipynb