本文主要圍繞scipy中的稀疏矩陣展開，也會介紹幾種scipy之外的稀疏矩陣的存儲方式。
dok_matrix
繼承自dict，key是(row,col)構成的二元組，value是非0元素。
優點：
- 非常高效地添加、刪除、查找元素
- 轉換成coo_matrix很快
缺點：
- 繼承了dict的缺點，即內存開銷大
- 不能有重復的(row,col)
適用場景：
- 加載數據文件時使用dok_matrix快速構建稀疏矩陣，然后轉換成其他形式的稀疏矩陣
coo_matrix
如上圖，構造coo_matrix需要3個等長的數組，values數組存放矩陣中的非0元素，row indices存放非0元素的行坐標，column indices存放非0元素的列坐標。
優點：
- 容易構造
- 可以快速地轉換成其他形式的稀疏矩陣
- 支持相同的(row,col)坐標上存放多個值
缺點：
- 構建完成后不允許再插入或刪除元素
- 不能直接進行科學計算和切片操作
適用場景：
- 加載數據文件時使用coo_matrix快速構建稀疏矩陣，然后調用to_csr()、to_csc()、to_dense()把它轉換成CSR或稠密矩陣
csr_matrix
csr_matrix同樣由3個數組組成，values存儲非0元素，column indices存儲非0元素的列坐標，row offsets依次存儲每行的首元素在values中的坐標，如果某行全是0則對應的row offsets值為-1（我猜的）。
優點：
- 高效地按行切片
- 快速地計算矩陣與向量的內積
- 高效地進行矩陣的算術運行，CSR + CSR、CSR * CSR等
缺點：
- 按列切片很慢（考慮CSC）
- 一旦構建完成后，再往里面添加或刪除元素成本很高
csc_matrix
跟csr_matrix剛好反過來。
bsr_matrix
跟CSR/CSC很相近，尤其適用於稀疏矩陣中包含稠密子矩陣的情況。在解決矢量值有限元離散（vector-valued finite element discretizations）這類問題中BSR比CSR/CSC更高效。
dia_matrix
對角線存儲法，按對角線方式存，列代表對角線，行代表行。省略全零的對角線。(從左下往右上開始：第一個對角線是零忽略，第二個對角線是5，6，第三個對角線是零忽略，第四個對角線是1，2，3，4，第五個對角線是7，8，9，第六第七個對角線忽略)。[3]
這里行對應行，所以5和6是分別在第三行第四行的，前面補上無效元素*。如果對角線中間有0，存的時候也需要補0。
適用場景：
- 如果原始矩陣就是一個對角性很好的矩陣那壓縮率會非常高，比如下圖，但是如果是隨機的那效率會非常糟糕。
lil_matrix
內部結構是個二維數組:[[(col,value)]]，第一行對應原矩陣的一行（可以快速地定位到行），行內按列編號排序好（通過折半查找可以快速地定位到列），同樣只存儲非0元素。
優點：
- 快速按行切片
- 高效地添加、刪除、查找元素
缺點：
- 按列切片很慢（考慮CSC）
- 算術運算LIL+LIL很慢（考慮CSR或CSC）
- 矩陣和向量內和解很慢（考慮CSR或CSC）
適用場景：
- 加載數據文件時使用lil_matrix快速構建稀疏矩陣，然后調用to_csr()、to_csc()把它轉換成CSR/CSC進行后續的矩陣運算
- 非0元素非常多時，考慮使用coo_matrix（我個人是這樣理解的，lil_matrix用一個二維數組搞定，二維數組占用的是連續的內存空間，如果非0元素非常多就要申請一塊非常大的連續的內存空間，這樣性能很差。而coo_matrix畢竟是使用的3個一維數組，對連續內存空間的要求沒那么高）
ELLPACK (ELL)
用兩個和原始矩陣相同行數的矩陣來存：第一個矩陣存的是列號，第二個矩陣存的是數值，行號就不存了，用自身所在的行來表示；這兩個矩陣每一行都是從頭開始放，如果沒有元素了就用個標志比如*結束。 上圖中間矩陣有誤，第三行應該是 0 2 3。
注：這樣如果某一行很多元素，那么后面兩個矩陣就會很胖，其他行結尾*很多，浪費。可以存成數組，比如上面兩個矩陣就是：
0 1 * 1 2 * 0 2 3 * 1 3 *
1 7 * 2 8 * 5 3 9 * 6 4 *
但是這樣要取一行就比較不方便了。
Hybrid (HYB) ELL + COO
為了解決ELL中提到的，如果某一行特別多，造成其他行的浪費，那么把這些多出來的元素（比如第三行的9，其他每一行最大都是2個元素）用COO單獨存儲。
skyline matrix storage
沒看明白，自行wiki。
適用場景：
- 非常適合於稀疏矩陣的Cholesky分解或LU分解，這兩種分解都是用來解線性方程組的。
行列雙索引
這是自己實現的一種存儲方式，分別按行和按列建立dict（dict中的key是行號或列號），這樣按下標查找元素很快，但犧牲了空間。為了挽回空間上的犧牲，我們采用二進制來存儲dict中的value。按下標查找元素時，根據行號定位到相應的value，value反序列化后轉成dict，該dict的key是列號。
class DictMatrix():
'''用dict來實現稀疏矩陣
'''
def __init__(self, dft=0.0):
self._data = {}
self._dft = dft # “0元素”的值
self._nums = 0 # 稀疏矩陣中非0元素的個數
def __setitem__(self, index, value):
try:
i, j = index
except:
raise IndexError('invalid index')
# 為了節省內存，我們把j, value打包成字二進制字符串
ik = ('i%d' % i)
ib = struct.pack('if', j, value) # 格式化：i代替integer，f代表float。pack方法返回字符串
jk = ('j%d' % j)
jb = struct.pack('if', i, value)
try:
self._data[ik] += ib # 拼接字符串
except:
self._data[ik] = ib
try:
self._data[jk] += jb
except:
self._data[jk] = jb
self._nums += 1
def __getitem__(self, index):
try:
i, j = index
except:
raise IndexError('invalid index')
if (isinstance(i, int) and isinstance(j, int)):
ik = ('i%d' % i)
if self._data.has_key(ik):
ret = dict(np.fromstring(self._data[ik], dtype='i4,f4'))
return ret.get(j, self._dft)
else:
raise IndexError('invalid index')
return self._dft
def getRow(self, index):
'''獲取某一行的數據，返回dict
'''
rect = dict()
if isinstance(index, int):
ik = ('i%d' % index)
if self._data.has_key(ik):
rect = dict(np.fromstring(self._data[ik], dtype='i4,f4'))
return rect
def getCol(self, index):
'''獲取某一列的數據，返回dict
'''
rect = dict()
if isinstance(index, int):
jk = ('j%d' % index)
if self._data.has_key(jk):
rect = dict(np.fromstring(self._data[jk], dtype='i4,f4'))
return rect
def __len__(self):
return self._nums
def __iter__(self):
pass
def read(self, cache):
'''cache是一個list，其中的每個元素都是個三元組(row,col,value)。
從磁盤中加載稀疏矩陣時，可以先把部分數據加載到cache中，再從cache放到DictMatrix中。
'''
tmpDict = {}
for row, col, value in cache:
if value != self._dft: # 確保添加的是“非0”元素
ik = ('i%d' % row)
ib = struct.pack('if', col, value)
jk = ('j%d' % col)
jb = struct.pack('if', row, value)
try:
tmpDict[ik].write(ib)
except:
# 考慮到字符串拼接性能不太好，我們直接用StringIO的write()來做拼接
tmpDict[ik] = StringIO()
tmpDict[ik].write(ib)
try:
tmpDict[jk].write(jb)
except:
tmpDict[jk] = StringIO()
tmpDict[jk].write(jb)
self._nums += 1
for k, v in tmpDict.items():
v.seek(0)
s = v.read()
try:
self._data[k] += s
except:
self._data[k] = s
if __name__ == '__main__':
dtv = -1.0 # 默認值
matrix = DictMatrix(dtv)
matrix[1, 9] = 58.0
matrix[1, 16] = 20.0
matrix[2, 16] = 9
assert matrix[1, 10] == dtv # 元素不存在，取默認值
for k, v in matrix.getRow(1).items():
print k, v
for k, v in matrix.getCol(16).items():
print k, v
assert matrix[1, 9] == 58.0
assert matrix[1, 16] == 20.0
assert matrix[2, 16] == 9
assert matrix[2, 9] == dtv
print len(matrix)
matrix.read([(3, 3, 15), (9, 3, 100.0)]) # 批量添加數據
assert matrix[1, 9] == 58.0
assert matrix[3, 3] == 15
assert matrix[9, 3] == 100.0
print len(matrix)
上面的代碼中做二進制序列化時用到了struck.pack，來個小例子看下序列化能省多少內存。
# coding=utf-8
__author__ = "orisun"
import struct
from cStringIO import StringIO
import numpy as np
from collections import defaultdict
loop1 = 10000
loop2 = 30
idx1 = 10
idx2 = 10
@profile
def foo1():
position_link = defaultdict(list) # list中每個元素是(int,float)類型的tuple
for i in xrange(loop1):
for j in xrange(loop2):
position_link[i].append((100000, 0.123435465))
print '1', position_link[idx1][idx2]
@profile
def foo2():
position_link = {}
for i in xrange(loop1):
tmp = StringIO()
for j in xrange(loop2):
# 使用StringIO的write()方法做二進制拼接，效率高一些
tmp.write(struct.pack('if', 100000, 0.123435465))
tmp.seek(0)
position_link[i] = tmp.read()
li = list(np.fromstring(position_link[idx1], dtype='i4,f4'))
print '2', li[idx2]
if __name__ == '__main__':
'''foo2比foo1節省16.7%的內存，但是慢了61.7%
使用memory_profiler查看內存使用
使用line_profiler查看耗時
'''
foo1()
foo2()
優點：
- 高效地動態添加元素
- 高效地按下標查找元素
- 高效地按行切片和按列切片
缺點：
- 不支持刪除元素
- 內存占用略大
選擇稀疏矩陣存儲格式的經驗
- DIA和ELL格式在進行稀疏矩陣-矢量乘積(sparse matrix-vector products)時效率最高，所以它們是應用迭代法(如共軛梯度法)解稀疏線性系統最快的格式
- COO格式常用於從文件中進行稀疏矩陣的讀寫，如matrix market即采用COO格式，而CSR格式常用於讀入數據后進行稀疏矩陣計算