本文最后更新于 258 天前,其中的信息可能已经过时,如有错误请发送邮件到wuxianglongblog@163.com
稀疏矩阵
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 格式在乘法计算时比较高效
构建稀疏矩阵
from scipy.sparse import *
import numpy as np
创建一个空的稀疏矩阵:
coo_matrix((2,3))
<2x3 sparse matrix of type ''
with 0 stored elements in COOrdinate format>
也可以使用一个已有的矩阵或数组或列表中创建新矩阵:
A = coo_matrix([[1,2,0],[0,0,3],[4,0,5]])
print A
(0, 0) 1
(0, 1) 2
(1, 2) 3
(2, 0) 4
(2, 2) 5
不同格式的稀疏矩阵可以相互转化:
type(A)
scipy.sparse.coo.coo_matrix
B = A.tocsr()
type(B)
scipy.sparse.csr.csr_matrix
可以转化为普通矩阵:
C = A.todense()
C
matrix([[1, 2, 0],
[0, 0, 3],
[4, 0, 5]])
与向量的乘法:
v = np.array([1,0,-1])
A.dot(v)
array([ 1, -3, -1])
还可以传入一个 (data, (row, col))
的元组来构建稀疏矩阵:
I = np.array([0,3,1,0])
J = np.array([0,3,1,2])
V = np.array([4,5,7,9])
A = coo_matrix((V,(I,J)),shape=(4,4))
print A
(0, 0) 4
(3, 3) 5
(1, 1) 7
(0, 2) 9
COO 格式的稀疏矩阵在构建的时候只是简单的将坐标和值加到后面,对于重复的坐标不进行处理:
I = np.array([0,0,1,3,1,0,0])
J = np.array([0,2,1,3,1,0,0])
V = np.array([1,1,1,1,1,1,1])
B = coo_matrix((V,(I,J)),shape=(4,4))
print B
(0, 0) 1
(0, 2) 1
(1, 1) 1
(3, 3) 1
(1, 1) 1
(0, 0) 1
(0, 0) 1
转换成 CSR 格式会自动将相同坐标的值合并:
C = B.tocsr()
print C
(0, 0) 3
(0, 2) 1
(1, 1) 2
(3, 3) 1
求解微分方程
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
from numpy.linalg import solve, norm
from numpy.random import rand
构建 1000 x 1000
的稀疏矩阵:
A = lil_matrix((1000, 1000))
A[0, :100] = rand(100)
A[1, 100:200] = A[0, :100]
A.setdiag(rand(1000))
转化为 CSR 之后,用 spsolve
求解 $Ax=b$:
A = A.tocsr()
b = rand(1000)
x = spsolve(A, b)
转化成正常数组之后求解:
x_ = solve(A.toarray(), b)
查看误差:
err = norm(x-x_)
err
6.4310987107687431e-13
sparse.find 函数
返回一个三元组,表示稀疏矩阵中非零元素的 (row, col, value)
:
from scipy import sparse
row, col, val = sparse.find(C)
print row, col, val
[0 0 1 3] [0 2 1 3] [3 1 2 1]
sparse.issparse 函数
查看一个对象是否为稀疏矩阵:
sparse.issparse(B)
True
或者
sparse.isspmatrix(B.todense())
False
还可以查询是否为指定格式的稀疏矩阵:
sparse.isspmatrix_coo(B)
True
sparse.isspmatrix_csr(B)
False