稀疏矩陣存儲格式總結+存儲效率對比:COO,CSR,DIA,ELL,HYB

稀疏矩陣是指矩陣中的元素大部分是0的矩陣，事實上，實際問題中大規模矩陣基本上都是稀疏矩陣，不少稀疏度在90%甚至99%以上。所以咱們須要有高效的稀疏矩陣存儲格式。本文總結幾種典型的格式：COO,CSR,DIA,ELL,HYB。

 

（1）Coordinate（COO）

這是最簡單的一種格式，每個元素須要用一個三元組來表示，分別是（行號，列號，數值），對應上圖右邊的一列。這種方式簡單，可是記錄單信息多（行列），每一個三元組本身能夠定位，所以空間不是最優。

 

（2）Compressed Sparse Row (CSR)

CSR是比較標準的一種，也須要三類數據來表達：數值，列號，以及行偏移。CSR不是三元組，而是總體的編碼方式。數值和列號與COO一致，表示一個元素以及其列號，行偏移表示某一行的第一個元素在values裏面的起始偏移位置。如上圖中，第一行元素1是0偏移，第二行元素2是2偏移，第三行元素5是4偏移，第4行元素6是7偏移。在行偏移的最後補上矩陣總的元素個數，本例中是9。

 

CSC是和CSR相對應的一種方式，即按列壓縮的意思。

以上圖中矩陣爲例：

Values：        [1 5 7 2 6 8 3 9 4]

Row Indices：[0 2 0 1 3 1 2 2 3]

Column Offsets：[0 2 5 7 9]

 

再來看一個CSR的例子[4]：

 

（3）ELLPACK (ELL)

用兩個和原始矩陣相同行數的矩陣來存：第一個矩陣存的是列號，第二個矩陣存的是數值，行號就不存了，用自身所在的行來表示；這兩個矩陣每一行都是從頭開始放，若是沒有元素了就用個標誌好比*結束。上圖中間矩陣有誤，第三行應該是  0 2 3。

 

注：這樣若是某一行不少元素，那麼後面兩個矩陣就會很胖，其餘行結尾*不少，浪費。能夠存成數組，好比上面兩個矩陣就是：

0 1 * 1 2 * 0 2 3 * 1 3 *

1 7 * 2 8 * 5 3 9 * 6 4 *

可是這樣要取一行就比較不方便了

（4）Diagonal (DIA)

對角線存儲法，按對角線方式存，列表明對角線，行表明行。省略全零的對角線。(從左下往右上開始：第一個對角線是零忽略，第二個對角線是5，6，第三個對角線是零忽略，第四個對角線是1，2，3，4，第五個對角線是7，8，9，第六第七個對角線忽略)。[3]

這裏行對應行，因此5和6是分別在第三行第四行的，前面補上無效元素*。若是對角線中間有0，存的時候也須要補0，因此若是原始矩陣就是一個對角性很好的矩陣那壓縮率會很是高，好比下圖，可是若是是隨機的那效率會很是糟糕。

 

（5）Hybrid (HYB) ELL + COO

爲了解決（3）ELL中提到的，若是某一行特別多，形成其餘行的浪費，那麼把這些多出來的元素（好比第三行的9，其餘每一行最大都是2個元素）用COO單獨存儲。

 

選擇稀疏矩陣存儲格式的一些經驗[2]：

1. DIA和ELL格式在進行稀疏矩陣-矢量乘積(sparse matrix-vector products)時效率最高，因此它們是應用迭代法(如共軛梯度法)解稀疏線性系統最快的格式；
2. COO和CSR格式比起DIA和ELL來，更加靈活，易於操做；
3. ELL的優勢是快速，而COO優勢是靈活，兩者結合後的HYB格式是一種不錯的稀疏矩陣表示格式；
4. 根據Nathan Bell的工做，CSR格式在存儲稀疏矩陣時非零元素平均使用的字節數(Bytes per Nonzero Entry)最爲穩定（float類型約爲8.5，double類型約爲12.5），而DIA格式存儲數據的非零元素平均使用的字節數與矩陣類型有較大關係，適合於StructuredMesh結構的稀疏矩陣（float類型約爲4.05，double類型約爲8.10），對於Unstructured Mesh以及Random Matrix,DIA格式使用的字節數是CSR格式的十幾倍；
5. 從我使用過的一些線性代數計算庫來講，COO格式經常使用於從文件中進行稀疏矩陣的讀寫，如matrix market即採用COO格式，而CSR格式經常使用於讀入數據後進行稀疏矩陣計算。

 

一些特殊類型矩陣的存儲效率（數值越小說明壓縮率越高，即存儲效率越高）:

Structured Mesh

Unstructured Mesh

Random matrix

Power-Law Graph

 

格式適用性總結：

 

 

下面摘自[2]

6. Skyline Storage Format

The skyline storage format is important for the direct sparse solvers, and it is well suited for Cholesky or LU decomposition when no pivoting is required.

The skyline storage format accepted in Intel MKL can store only triangular matrix or triangular part of a matrix. This format is specified by two arrays:values andpointers. The following table describes these arrays:

values

A scalar array. For a lower triangular matrix it contains the set of elements from each row of the matrix starting from the first non-zero element to and including the diagonal element. For an upper triangular matrix it contains the set of elements from each column of the matrix starting with the first non-zero element down to and including the diagonal element. Encountered zero elements are included in the sets.

pointers

An integer array with dimension(m+1), where m is the number of rows for lower triangle (columns for the upper triangle).pointers(i) -pointers(1)+1gives the index of element invalues that is first non-zero element in row (column)i. The value ofpointers(m+1)is set tonnz+pointers(1), wherennz is the number of elements in the arrayvalues.

7. Block Compressed Sparse Row Format (BSR)

The Intel MKL block compressed sparse row (BSR) format for sparse matrices is specified by four arrays:values,columns,pointerB, andpointerE. The following table describes these arrays.

values

A real array that contains the elements of the non-zero blocks of a sparse matrix. The elements are stored block-by-block in row-major order. A non-zero block is the block that contains at least one non-zero element. All elements of non-zero blocks are stored, even if some of them is equal to zero. Within each non-zero block elements are stored in column-major order in the case of one-based indexing, and in row-major order in the case of the zero-based indexing.

columns

Element i of the integer array columns is the number of the column in the block matrix that contains thei-th non-zero block.

pointerB

Element j of this integer array gives the index of the element in thecolumns array that is first non-zero block in a rowj of the block matrix.

pointerE

Element j of this integer array gives the index of the element in thecolumns array that contains the last non-zero block in a rowj of the block matrix plus 1.

 

 

 

 

[1] Sparse Matrix Representations & Iterative Solvers, Lesson 1 by Nathan Bell. http://www.bu.edu/pasi/files/2011/01/NathanBell1-10-1000.pdf

[2] http://blog.csdn.net/anshan1984/article/details/8580952

[3] http://zhangjunhd.github.io/2014/09/29/sparse-matrix.html

[4] http://www.360doc.com/content/09/0204/17/96202_2458312.shtml

[5] Implementing Sparse Matrix-Vector Multiplication on Throughput-Oriented Processors, Nathan Bell and Michael Garland, Proceedings of Supercomputing '09

[6] Efficient Sparse Matrix-Vector Multiplication on CUDA, Nathan Bell and Michael Garland, NVIDIA Technical Report NVR-2008-004, December 2008