1. 前言

數(shù)組和矩陣是數(shù)值計(jì)算的基礎(chǔ)元素。目前為止，我們都是使用NumPy的ndarray數(shù)據(jù)結(jié)構(gòu)來表示數(shù)組，這是一種同構(gòu)的容器，用于存儲(chǔ)數(shù)組的所有元素。

有一種特殊情況，矩陣的大部分元素都為零，這種矩陣被稱為稀疏矩陣。對(duì)于稀疏矩陣，將所有零保存在計(jì)算機(jī)內(nèi)存中的效率很低，更合適的方法是只保存非零元素以及位置信息。于是SciPy應(yīng)運(yùn)而生，為稀疏矩陣的表示及其線性代數(shù)運(yùn)算提供了豐富易用的接口。

2. 導(dǎo)入包

SciPy中提供了稀疏矩陣模塊scipy.sparse，為稀疏矩陣的表示及其線性代數(shù)運(yùn)算提供了豐富易用的接口。

import scipy.sparse as sp
import scipy.sparse.linalg
import scipy.linalg as la

3. 稀疏矩陣總覽

There are seven available sparse matrix types:
    1. csc_matrix: Compressed Sparse Column format
    2. csr_matrix: Compressed Sparse Row format
    3. bsr_matrix: Block Sparse Row format
    4. lil_matrix: List of Lists format
    5. dok_matrix: Dictionary of Keys format
    6. coo_matrix: COOrdinate format (aka IJV, triplet format)
    7. dia_matrix: DIAgonal format

• sp.coo_matrix（坐標(biāo)的列表）描述：將非零值及其行列信息保存在一個(gè)列表。構(gòu)造簡(jiǎn)單，添加元素方便，但訪問元素效率低下；
• sp.lil_matrix（列表的列表）：將每行的非零元素列索引保存在一個(gè)列表，將對(duì)應(yīng)值保存在另一個(gè)列表。支持切片操作，但不方便進(jìn)行數(shù)學(xué)運(yùn)算；
• sp.dok_matrix（值的字典）：將非零值保存在字典，非零值的坐標(biāo)元組作為字典的鍵。構(gòu)造簡(jiǎn)單，可快速添加刪除元素，但不方便進(jìn)行數(shù)學(xué)運(yùn)算；
• sp.dia_matrix（對(duì)角矩陣）：矩陣的對(duì)角線列表 。對(duì)于對(duì)角矩陣非常有效 ，但不適用非對(duì)角矩陣；
• sp.csc_matrix 和 sp.csr_matrix（壓縮列格式和壓縮行格式）：將值與行列索引的數(shù)組一起存儲(chǔ)。對(duì)于矩陣的向量乘法很高效，但構(gòu)造相對(duì)復(fù)雜；
• sp.bsr_matrix（塊稀疏矩陣）：與CSR類似，用于具有稠密子矩陣的稀疏矩陣。對(duì)于此類特殊矩陣很高效，但不適用一般矩陣。

注意：

• 為了有效地構(gòu)建矩陣，使用dok_matrix或者lil_matrix，lil_matrix類支持基本的切片和索引操作，語(yǔ)法與NumPy的arrays相似。
• COO格式也能有效率地構(gòu)建矩陣。盡管與NumPy有許多相似性，但是強(qiáng)烈不建議使用NumPy的函數(shù)直接對(duì)稀疏矩陣格式進(jìn)行操作，因?yàn)榭赡軐?dǎo)致不正確的結(jié)果。如果將NumPy的函數(shù)用在這些矩陣上，首先檢查SciPy在對(duì)應(yīng)的稀疏矩陣類上有沒有已經(jīng)實(shí)現(xiàn)的操作，或者使用toarray()方法將稀疏矩陣對(duì)象轉(zhuǎn)為NumPy的array。
• 實(shí)現(xiàn)乘法與轉(zhuǎn)置操作，則轉(zhuǎn)為CSC或CSR格式，lil_matrix格式是基于行的，所以轉(zhuǎn)為為CSR比CSC更有效率。所有的轉(zhuǎn)換在CSR,CSC和COO格式之間都是有效的，線性時(shí)間操作。

4. 稀疏矩陣詳細(xì)介紹

4.1 coo_matrix

coo_matrix是最簡(jiǎn)單的存儲(chǔ)方式。采用三個(gè)數(shù)組row、col和data保存非零元素的行下標(biāo)，列下標(biāo)與值。這三個(gè)數(shù)組的長(zhǎng)度相同一般來說，coo_matrix主要用來創(chuàng)建矩陣，因?yàn)閏oo_matrix無法對(duì)矩陣的元素進(jìn)行增刪改等操作，一旦創(chuàng)建之后，除了將之轉(zhuǎn)換成其它格式的矩陣，幾乎無法對(duì)其做任何操作和矩陣運(yùn)算。

為了創(chuàng)建sp.coo_matrix對(duì)象，需要?jiǎng)?chuàng)建非零值、行索引以及列索引的列表或數(shù)組，并將其傳遞給生成函數(shù)sp.coo_matrix。

values = [1, 2, 3, 4]
rows = [0, 1, 2, 3]
cols = [1, 3, 2, 0]
A = sp.coo_matrix((values, (rows, cols)), shape=[4, 4])
A

>>> A.toarray()
array([[1, 0, 0, 0],
       [0, 0, 0, 2],
       [0, 0, 3, 0],
       [4, 0, 0, 0]])
>>> type(A)
<class 'scipy.sparse.coo.coo_matrix'>
>>> type(A.toarray())
<class 'numpy.ndarray'>

SciPy的sparse模塊中稀疏矩陣的屬性大部分派生自NumPy的ndarray對(duì)象，同時(shí)也包括nnz（非零元素?cái)?shù)目）和data（非零值）等屬性。

A.shape, A.size, A.dtype, A.ndim

A.nnz, A.data

對(duì)于sp.coo_matrix對(duì)象，還可以使用row和col屬性來訪問底層的行列坐標(biāo)數(shù)組。

A.row, A.col

4.2 dok_matrix

dok_matrix適用的場(chǎng)景是逐漸添加矩陣的元素。dok_matrix的策略是采用字典來記錄矩陣中不為0的元素。所以字典的key存的是記錄元素的位置信息的元祖，value是記錄元素的具體值。

>>> S = sparse.dok_matrix((5, 5), dtype=np.float32)
>>> for i in range(5):
    	for j in range(5):
        	S[i,j] = i+j    # 更新元素
>>> S.toarray()
[[0. 1. 2. 3. 4.]
 [1. 2. 3. 4. 5.]
 [2. 3. 4. 5. 6.]
 [3. 4. 5. 6. 7.]
 [4. 5. 6. 7. 8.]]

優(yōu)點(diǎn)：對(duì)于遞增的構(gòu)建稀疏矩陣很高效，比如定義該矩陣后，想進(jìn)行每行每列更新值，可用該矩陣。當(dāng)訪問某個(gè)單元，只需要O(1)

缺點(diǎn)：不允許重復(fù)索引（coo中適用），但可以很高效的轉(zhuǎn)換成coo后進(jìn)行重復(fù)索引。

4.3 lil_matrix

lil_matrix適用的場(chǎng)景也是逐漸添加矩陣的元素。與dok不同，lil_matrix則是使用兩個(gè)列表存儲(chǔ)非0元素。data保存每行中的非零元素,rows保存非零元素所在的列。這種格式也很適合逐個(gè)添加元素，并且能快速獲取行相關(guān)的數(shù)據(jù)。

>>> l = sparse.lil_matrix((4, 4))
>>> l[1, 1] = 1
>>> l[1, 3] =2
>>> l[2, 3] = 3
>>> l.toarray()
array([[0., 0., 0., 0.],
       [0., 1., 0., 2.],
       [0., 0., 0., 3.],
       [0., 0., 0., 0.]])
>>> l.data
array([list([]), list([1.0, 2.0]), list([3.0]), list([])], dtype=object)
>>> l.rows
array([list([]), list([1, 3]), list([3]), list([])], dtype=object)

優(yōu)點(diǎn)：適合遞增的構(gòu)建成矩陣、轉(zhuǎn)換成其它存儲(chǔ)方式很高效、支持靈活的切片。

缺點(diǎn)：當(dāng)矩陣很大時(shí)，考慮用coo、算術(shù)操作，列切片，矩陣向量?jī)?nèi)積操作慢。

4.4 dia_matrix

如果稀疏矩陣僅包含非0元素的對(duì)角線，則對(duì)角存儲(chǔ)格式(DIA)可以減少非0元素定位的信息量。這種存儲(chǔ)格式對(duì)有限元素或者有限差分離散化的矩陣尤其有效。dia_matrix通過兩個(gè)數(shù)組確定： data和offsets。其中data對(duì)角線元素的值；offsets：第i個(gè)offsets是當(dāng)前第i個(gè)對(duì)角線和主對(duì)角線的距離。data[k:]存儲(chǔ)了offsets[k]對(duì)應(yīng)的對(duì)角線的全部元素。例子如下：

>>> data = np.array([[1, 2, 3, 4], [5, 6, 0, 0], [0, 7, 8, 9]])
>>> offsets = np.array([0, -2, 1])
>>> sparse.dia_matrix((data, offsets), shape=(4, 4)).toarray()
array([[1, 7, 0, 0],
       [0, 2, 8, 0],
       [5, 0, 3, 9],
       [0, 6, 0, 4]])

注意：offsets[0]=0表示第0個(gè)對(duì)角線與主對(duì)角線的距離為0，表示第0個(gè)對(duì)角線就是主對(duì)角線，data[0]就是第0個(gè)對(duì)角線的值。offsets[1]=-2表示第1個(gè)對(duì)角線與主對(duì)角線距離為-2，此時(shí)該對(duì)角線在主對(duì)角線的左下方，對(duì)角線上數(shù)值的數(shù)量為4-2=2，對(duì)應(yīng)的值為data[1, :2+1]，此時(shí)data[1, 3:]為無效的值，在構(gòu)造對(duì)角稀疏矩陣時(shí)不起作用。offsets[2]=1表示第2個(gè)對(duì)角線與主對(duì)角線距離為1，此時(shí)該對(duì)角線在主對(duì)角線的右上方，對(duì)角線上數(shù)值的數(shù)量為4-1=3，對(duì)應(yīng)的值為data[2, 1:]，此時(shí)data[2, :1]為無效的值，在構(gòu)造對(duì)角稀疏矩陣時(shí)不起作用。

4.5 csc_matrix & csr_matrix

csr_matrix是按行對(duì)矩陣進(jìn)行壓縮的，csc_matrix則是按列對(duì)矩陣進(jìn)行壓縮的。通過row_offsets,column_indices，data來確定矩陣。column_indices，data與coo格式的列索引與數(shù)值的含義完全相同，row_offsets表示元素的行偏移量。例子如下，

>>> 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((data, indices, indptr), shape=(3, 3)).toarray()
array([[1, 0, 2],
       [0, 0, 3],
       [4, 5, 6]])

注意：indices和data分別表示列索引和數(shù)據(jù)，第 i 行的列索引存儲(chǔ)在indices[indptr[i]:indptr[i+1]] 中，對(duì)應(yīng)的值為data[indptr[i]:indptr[i+1]]。即第0行的列索引為indices[0:2]=[0,2]，值為data[0:2]=[1,2]；第1行的列索引為indices[2:3]=[2]，值為data[2:3]=[3]…

CSR格式常用于讀入數(shù)據(jù)后進(jìn)行稀疏矩陣計(jì)算。

兩者的優(yōu)缺點(diǎn)互反：

• CSR優(yōu)點(diǎn)：高效的稀疏矩陣算術(shù)操作、高效的行切片、快速地矩陣向量?jī)?nèi)積操作；
• CSR缺點(diǎn)：緩慢地列切片操作（可以考慮csc）、轉(zhuǎn)換到稀疏結(jié)構(gòu)代價(jià)較高（可以考慮lil，dok）。
• CSC優(yōu)點(diǎn)：高效的稀疏矩陣算術(shù)操作、高效的列切片、快速地矩陣向量?jī)?nèi)積操作（不如csr，bsr塊）；
• CSC缺點(diǎn)：緩慢地行切片操作（可以考慮csr）、 轉(zhuǎn)換到稀疏結(jié)構(gòu)代價(jià)較高（可以考慮lil，dok）。

4.6 bsr_matrix

基于行的塊壓縮,通過row_offsets，column_indices，data來確定矩陣。與csr相比，只是data中的元數(shù)據(jù)由0維的數(shù)變?yōu)榱艘粋€(gè)矩陣（塊），其余完全相同。

>>> 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]).repeat(4).reshape(6,2,2)
>>> bsr_matrix((data,indices,indptr), shape=(6,6)).todense()
matrix([[1, 1, 0, 0, 2, 2],
        [1, 1, 0, 0, 2, 2],
        [0, 0, 0, 0, 3, 3],
        [0, 0, 0, 0, 3, 3],
        [4, 4, 5, 5, 6, 6],
        [4, 4, 5, 5, 6, 6]])

優(yōu)點(diǎn)：很類似于csr，更適合于矩陣的某些子矩陣很多值，在某些情況下比csr和csc計(jì)算更高效。

5. 稀疏矩陣的存取

5.1 用save_npz保存單個(gè)稀疏矩陣

>>> scipy.sparse.save_npz('sparse_matrix.npz', sparse_matrix)
>>> sparse_matrix = scipy.sparse.load_npz('sparse_matrix.npz')

稀疏矩陣存儲(chǔ)大小比較：

a = np.arange(100000).reshape(1000,100)
a[10: 300] = 0
b = sparse.csr_matrix(a)
# 稀疏矩陣壓縮存儲(chǔ)到npz文件
sparse.save_npz('b_compressed.npz', b, True)  # 文件大小：100KB
# 稀疏矩陣不壓縮存儲(chǔ)到npz文件
sparse.save_npz('b_uncompressed.npz', b, False)  # 文件大小：560KB
# 存儲(chǔ)到普通的npy文件
np.save('a.npy', a)  # 文件大小：391KB
# 存儲(chǔ)到壓縮的npz文件
np.savez_compressed('a_compressed.npz', a=a)  # 文件大小：97KB

對(duì)于存儲(chǔ)到npz文件中的CSR格式的稀疏矩陣，內(nèi)容為：

data.npy
format.npy
indices.npy
indptr.npy
shape.npy

6. 總結(jié)

加載數(shù)據(jù)文件時(shí)使用coo_matrix快速構(gòu)建稀疏矩陣，然后調(diào)用to_csr()、to_csc()、to_dense()把它轉(zhuǎn)換成CSR或稠密矩陣(numpy.matrix)。

coo_matrix格式常用于從文件中進(jìn)行稀疏矩陣的讀寫，而csr_matrix格式常用于讀入數(shù)據(jù)后進(jìn)行稀疏矩陣計(jì)算。

7. 參考

【1】https://blog.csdn.net/winycg/article/details/80967112
【2】https://blog.csdn.net/vor234/article/details/124935384