1.3.3　稀疏矩陣向量乘法

稀疏矩陣向量乘法（SPMV）可在很多情況下代替稠密矩陣運算，可以大量節省內存占用，減少計算開銷。矩陣向量乘法不同于矩陣和矩陣的乘法，這是完全訪存密集型的計算，我們主要的優化方向是提升訪存效率或減少訪存開銷。

稀疏矩陣一般只存儲非零元的信息，非零元的存儲格式決定了訪存的模式，這需要根據非零元的分布模式和要做的計算類型來設計。我們假設分布模式并非對角線分布，整體分布較均勻，局部可能會有聚集，計算類型是稀疏矩陣乘以稠密向量，結果為稠密向量。

標準的稀疏矩陣存儲格式主要有：COO（Coordinate Format）和CSR（Compressed Sparse Row）等。COO很簡單，就是使用3個數組，分別存儲全部非零元的行下標（row index）、列下標（column index）和值（value）；CSR稍復雜，對行下標進行了壓縮，假設矩陣行數是m，則壓縮后的數組長度為m+1，記作（row ptr），其中第i個元素（0-base）表示矩陣前i行的非零元個數。

圖1-1和圖1-2展示了COO和CSR格式存儲稀疏矩陣的一個例子。

圖1-1　一個稀疏矩陣的例子

我們來考慮矩陣向量乘法計算y=Ax，其中A是稀疏矩陣，維度是m和n，非零元個數是k；x和y是稠密向量，維度分別是n和m，m×n>>k>>max（m，n）。做這個稀疏矩陣向量乘法就要遍歷A的每一行，和x對應位置相乘，把結果累加到y的對應位置。這個過程對A的k個非零元全部訪問了一遍，對x也訪問了k個元素（重疊），對y訪問了一遍，所以優化重點在于減少訪問A的冗余，并提升訪問x的效率。下面這幾個優化標準稀疏矩陣存儲格式的方法，可以提升訪存效率，減少冗余。

圖1-2　COO和CSR格式的稀疏矩陣存儲方法

（1）對矩陣A做行列分塊處理

對x的訪問每次總是從左到右進行稀疏的遍歷，如果n很大（比如上百萬甚至更多），則訪問x的空間局部性較差。所以我們首先改進矩陣A的訪問順序，將矩陣A分解成多個方形的子矩陣。子矩陣的維度適應較高層CPU硬件cache的大小，這樣在遍歷每一個子矩陣時，對x的訪問相對集中于一個較小的區間，這個區間內的x會被cache緩存，這樣能夠大大提高訪問效率。分塊方式如圖1-3所示。

圖1-3　稀疏矩陣的分塊存儲

（2）自適應分塊存儲結構

由于稀疏矩陣的非零元分布不一定均勻，有的分塊會非常稀疏，有的則會相對稠密。對于極稀疏的分塊（非零元數量遠小于行數），如果用和CSR相似的壓縮行存儲策略，則會浪費空間，所以用COO的方式反而更能節省存儲空間，提高訪問效率。

對于哪些分塊使用CSR，哪些使用COO方式，可以通過實驗的方式確定一個非零元的數量和分塊大小的比值。高于該值的用CSR方式存儲，否則用COO方式存儲。

如圖1-4所示，一共使用5個數組存儲自適應分塊信息的稀疏矩陣，灰色的部分是CSR的相關信息，白色的部分是COO的相關信息。col_idx和vals的意義不變；types存儲分塊類型，標識當前分塊是CSR還是COO；如果當前分塊是CSR，則row_info存儲類似row_ptr的信息（第k個元素表示分塊內第k行的非零元個數），否則存儲COO的row_idx的信息；row_id存儲每個分塊在row_info上的起始地址。

圖1-4　自適應分塊的稀疏矩陣存儲格式

（3）減少下標存儲的冗余

矩陣分塊后，分塊內間址的下標并不需要4字節int型整數存儲，比如分塊維度在64K以內，可以用2字節的unsigned short來存儲。這樣，無論是CSR或COO的row_idx、row_ptr，還是col_idx，都可以減少50%的存儲空間，并同時提升訪存效率。

（4）多線程和NUMA特性

單處理器多核多線程并行計算稀疏矩陣向量乘的過程比較簡單，只需把矩陣劃分成線程數量的子矩陣。這里采用橫切的方法，計算結果不用合并。

但是對于多處理器非一致內存訪問（NUMA），就需要對數據在內存中的分布做特殊處理，才能最大程度地利用全部的內存帶寬。

一個典型的Intel X86雙路服務器的拓撲架構如圖1-5所示。

Memory#0是CPU#0的本地內存，Memory#1是CPU#1的本地內存，它們有各自獨立的內存帶寬。CPU#0訪問Memory#1需要經過內部總線（在Intel的架構中叫QPI總線），這個總線的帶寬一般小于內存帶寬。另外如果要訪問的數據只集中在一顆CPU的本地內存中，那么只能利用一個NUMA node的內存帶寬，這就限制了系統的總體吞吐。

所以需要把稀疏矩陣的存儲均勻地分配到兩顆處理器各自的本地內存中。對于一個雙CPU，每顆CPU一共4核的系統，需要開8個線程，并把這8個線程分別綁定到8顆CPU核上，使線程的上下文不會在核間遷移。對于每個線程要處理的稀疏矩陣數據，也通過系統調用（在Linux中是mbind），綁定到所在CPU核的本地內存中。這樣每個核處理的數據一定是從本地內存中獲得的，不會經過QPI總線。這就最大程度地利用了系統內存的帶寬。經過實測，這個優化方法可以提升70%左右的內存帶寬。

圖1-5　雙路NUMA服務器系統互連拓撲結構

對于我們測試的一個維度大約1M、稀疏度0.0001的稀疏矩陣來說，所有優化加起來，相對Intel MKL庫中CSR矩陣的SpMV API加速了2.5x左右。學術界還有很多針對稀疏矩陣存儲格式的討論和研究，其中有些還利用了SIMD向量指令，這里介紹的稀疏矩陣乘法方法，更多是為了討論內存和cache優化的一些基本原理。稀疏矩陣根據稀疏度和非零元分布的不同，需要使用不同的存儲策略，所以遇到實際的稀疏矩陣問題，需要根據實際情況開發不同的存儲格式，不必局限于本節描述的方法。