一種利用GPU優化大規模小方陣奇異值分解的新方法

李曉敏，鄢社鋒，侯朝煥

(中國科學院 聲學研究所，北京 100190)

傳統意義的 GPU主要針對游戲加速和圖形圖像處理。NVIDIA公司于2007年發布了基于CUDA的GPU，這種意義上的GPU具有很強的并行計算能力和數據吞吐能力，在諸多領域里能使程序性能提高好幾個數量級，對未來信息社會處理海量數據的需求具有很強的適應性。

基于CUDA的GPU自問世4年以來，由于其在性能、成本和開發周期上的顯著優勢，一經推出即在學術界和產業界引起了熱烈反響，已經廣泛應用于金融、石油、天文學、流體力學、信號處理、電磁仿真、模式識別、圖像處理、視頻壓縮等眾多領域[1]。

求解一般矩陣的 SVD最常用算法可以分為兩大類：基于QR分解的算法和基于Jacobi旋轉的分解[2]。前者首先通過一系列Householder變換實現雙對角，然后通過QR變換實現對角。其中，QR分解是一種嚴格串行的過程，只能獲得有限的并行度，因此不適合應用于GPU。后者使用平面旋轉進行正交化。雙邊Jacobi每次正交化矩陣的兩行和兩列，調度時同時涉及到行和列，不適合 GPU的并行架構。而單邊Jacobi每次旋轉只需要正交化矩陣的兩列，很容易實現并行訪問。

目前關于在GPU上進行SVD的研究并不多。最初，文獻[3]采用在GPU上實現雙對角變換和在CPU上實現對角變換的混合編程法完成了SVD。2009年，文獻[5]首次使用 GPU實現了雙對角-對角算法。同年，文獻[6]在 3D重建過程中將VL數據庫里的SVD程序轉換成CUDA語言，從而實現實時重建。2010年，文獻[7]采用 Webb算法在GPU上實現了SVD。然而，這些方法相對于MATLAB和IntelMKL的優勢僅僅體現在針對單個大矩陣SVD，并且當矩陣維數增大到成千上萬時這種優勢更加明顯。到目前為止，僅有文獻[8]實現了對單個512*512小型實矩陣 SVD的優化，而針對多個小型復矩陣SVD的優化幾乎沒有。

本文充分利用了 GPU內眾多并行內核以及共享存儲器和常數存儲器隱藏訪問時間的優勢，借助Brent和Luk提出的Jacobi行列調度法[9]，采用單邊 Jacobi旋轉，迭代實現了大規模小方陣的 SVD優化。

1 單邊Jacobi方法

對M*M維方陣A做SVD分解：

單邊Jacobi方法通過一系列平面旋轉(一般文獻中在介紹Jacobi方法時為了方便也將平面旋轉稱為Jacobi旋轉，即右乘正交陣V)使矩陣A正交化。首先，對矩陣A重復進行Jacobi旋轉，直到旋轉后的矩陣W各列兩兩正交，得到右奇異向量V；然后，對正交矩陣W歸一化，得到奇異值向量和左奇異向量U。

1.1 旋轉正交化過程

在矩陣 A 中任取一組列對(i，j)，計算第 i列內積存入a，計算第j列內積存入b，計算第i列和第j列的差積存入c；根據a，b，c計算旋轉角cs和旋轉因子sn；根據cs和sn對A中(i，j)兩列的每個元素進行代數運算，從而實現這兩列的正交，與此同時，對單位矩陣V中(i，j)兩列也進行同樣的代數運算。循環重復以上步驟，直到|c|/小于一個接近于0的正小數 tol，即表明A中任意兩列均正交。

1.2 歸一化和計算SVD

2 基于GPU的SVD實現

2.1 優化方法

CUDA是一種將GPU作為數據并行計算設備的軟硬件體系，采用了比較容易掌握的類C語言進行開發。它是一個 SIMD(Single Instruction Multiple Data)系統，即一個程序編譯一次以后，CUDA將計算任務映射為大量的可以并行執行的線程，并由擁有大量內核的硬件動態調度和執行這些線程，從而顯著提高運算速度。如圖1所示，將一個可以并行化執行的任務首先分配給若干個線程網格(Grid)，其次將每個Grid內的任務分配給若干個線程塊(Block)，最后再將每個Block內的任務細分給若干個線程(Thread)。Grid中的所有Blocks并行執行，Block中的所有 Threads并行執行，這種兩層并行模型是 CUDA最重要的創新之一。

圖1 GPU內線程結構Fig.1 Configuration of threads in GPU

CUDA在GPU內定義了8種寄存器，圖2所示展示了常用的4種。其中，共享存儲器(Shared Memory)是GPU片內的高速存儲器，通常用于存儲中間結果或 Block的公共結果；常數存儲器(ConstantMemory)是只讀的地址空間，雖然位于顯存，但是擁有緩存加速，在CUDA程序中用于存儲需要頻繁訪問的只讀參數。

本文的SVD優化方法將任務細分給多個Blocks內的多個 Threads，不同維數的矩陣，任務細分策略不同；同時，將需要頻繁使用的矩陣數據和中間結果存入共享存儲器，有效隱藏存儲器訪問延遲；另外，采用Brent和Luk提出的Jacobi行列調度法，將正交化過程中需要選取的所有列隊序號制成表格，存入常數存儲器，避免多次重復取數，使程序高效化執行。

圖2 GPU內主要存儲器結構Fig.2 Configuration of main memory in GPU

在每次任兩列旋轉正交之前，均執行算法1中的步驟。首先，初始化Lk=2k 1，Rk=2k。

算法1 Jacobi行列調度算法

1：如果 Lk

2：如果k=1，則outRk←Rk

3：否則如果k

4：如果k>1，則outLk← Rk，直到所有Rk賦值給 outLk；

5：如果 k>M/2，則 Rk←inRk，否則 Rk←Lk；6：如果k>1，則Lk←inLk

對于M維方陣，需要(M-1)輪正交化才能保證任意兩列之間均正交化。以M=6為例：

第一步正交的列對：(1，2)，(3，4)，(5，6)；

第二步正交的列對：(1，4)，(2，6)，(3，5)；

第三步正交的列對：(1，6)，(4，5)，(2，3)；

第四步正交的列對：(1，5)，(6，3)，(4，2)；

第五步正交的列對：(1，3)，(5，2)，(6，4)；

在旋轉變換之前，計算每一步需要正交的列對編號，并將其制成表格。

2.2 優化過程

本文顯卡采用NVIDIAGTX460，擁有48KB共享存儲器和16KB常數存儲器。設矩陣維數為M，需要做SVD的矩陣個數為N。

由于單邊 Jacobi方法每次需要正交化一個列對，所以矩陣維數最好為2的整數倍，本文分別選取16，32和64維的情況進行說明。執行正交化之前，文獻[10]將所有復矩陣轉化為具有等價奇異值和奇異向量的實矩陣，矩陣維數由原來的M增加到2*M。為了讓描述簡潔直觀，之后的介紹中 M實際代表2*M。

由于本方法旨在充分利用共享存儲器和常數存儲器，因此要特別注意所用資源不應超過 GPU提供的上限，不同維數的矩陣，任務細分策略不同。當矩陣維數為16時，經過計算可以將A和V兩個矩陣完全從全局存儲器拷貝到共享存儲器。使用N個Blocks，每個Block處理一個小方陣A，Block內有(M/2)*M個Threads，分別處理每次正交化過程中M/2個列對的M個方陣元素；當矩陣維數為32和64時候，經過計算 shared memory已經不能完全存儲A和V。因此，使用N*(M/2)個Blocks，每個Block只讀入一個列對到共享存儲器中進行處理，內有2*M個Threads，正交化該列隊的2*M個矩陣元素。

處理步驟如下：

1.將待處理的數據讀入shared memory。

當矩陣維數為16時。每個Block將要處理的矩陣A從global memory中讀入shared memory存為dA[M*M]，同時將單位陣V也讀入shared memory存為dV[M*M]。每個Block里有(M/2)個列對，因此定義 cs[M/2]、ss[M/2]和 c[M/2]分別存放每個列對里第1列內積、第2列內積以及兩列差積。

當矩陣維數為32和64時。根據當前Block的ID號查表得到所處理的是哪兩列，每個Block將要處理的兩列數據從矩陣A中讀入shared memory存為DA[2*M]，同時將單位陣V的相應位置兩列也讀入 shared memory存為 DV[2*M]。每個Block里有1個列對，因此定義cs、ss和c分別存放第1列內積、第2列內積和兩列差積。

2.求內積、差積和旋轉因子。

當矩陣維數為16時。根據當前Thread的ID號查表得到所處理的是哪兩列，求得內積和差積放入 shared memory中定義的臨時區域 temp[M*(M/2)]，并對temp進行規約求和將結果放入cs、ss和c中。

當矩陣維數為32和64時。得內積和差積放入shared memory中定義的臨時區域temp[2*M]，并對temp進行規約求和將結果放入cs、ss和c。

3.根據cs、ss和c求旋轉角后對矩陣進行旋轉，并將旋轉后的 dA和 dV(M=16)及 DA和 DV(M=32，64)放回原矩陣A和V中。

4.重復2和3，直到 c小于一定的值，說明當前一系列列對的正交化已經基本完成。

5.在主機端重復調用1-4共(N-1)次，遍歷完Jacobi行列調度表，保證經過旋轉后的A任意兩列之間均正交。

所有Blocks里的全部Threads同時執行以上1-4步，實現了最大并行化。所有4步均在sharedmemory中進行，其中 Jacobi行列調度表放在 constant memory中，大大減少了存儲器訪問時間，提高了數據吞吐量。

3 仿真結果

本節采用基于Intel Dual Core 2.10GHz CPU和NVIDIA GTX460顯卡，分別測試了基于MATLAB 2008a、Intel MKL 10.0.0 LAPACK 和 NVIDIA GTX460上SVD的性能。

隨機生成一批單精度32維方陣，圖4顯示采用MATLAB、MKL和GPU執行SVD的時間隨矩陣個數的變化曲線。

圖4 執行時間隨矩陣個數的變化(M=32)

從圖4可以看出，針對多個小規模方陣，GPU的執行時間明顯小于MATLAB和MKL。同時，基于MATLAB和MKL的SVD執行時間隨矩陣個數的增加成直線增長，而基于GPU的SVD執行時間隨矩陣個數的增加明顯趨于緩慢。

圖5 GPU相對于Matlab和MKL的加速比隨矩陣維數和矩陣個數的變化Fig.5 Speed up of GPU compared with Matlab and MKL mutative with matrix size and matrix number

圖 5展示了針對不同數量不同維數小方陣SVD，GPU相對于MKL和MATLAB的加速比曲線圖。

對于128個16維小方陣的SVD，GPU執行速度相比于 MATLAB提高了17.6倍，相比于 MKL提高了9.2倍。當矩陣維數大于16時，隨著矩陣個數增加，GPU執行SVD相比于這兩種方法的加速比均趨于穩定，分別約為5.1倍和3.4倍。

4 結束語

介紹了一種基于GPU針對大規模小方陣SVD的優化方法，適用于聲納和雷達信號處理算法。該方法充分利用了 GPU眾多并行內核以及共享存儲器和常數存儲器隱藏訪問時間的優勢。在同等條件下，相比于MATLAB和IntelMKL都有一定加速。

[1]張舒，褚艷麗.GPU高性能計算之CUDA[M].北京：中國水利水電出版社，2009：1-3.

[2]張賢達.矩陣分析與應用[M].北京：清華大學出版社，2008：358-365.

[3]Bondhugula V，Govindaraju N，Manocha D.Fast Singular Value Decomposition on Graphics Processors[R].University of North Carolina at Chapel Hill，2006.

[4]Singular Value Decomposition Routines[DB/OL].http://www.culatools.com/features/lapack/.

[5]Sheetal L，Narayanan P J.Singular Value Decomposition on GPU using CUDA[J].IEEE Symposium on Parallel&Distributed Processing，2009：1-10.

[6]Kimikazu K，Tikara H.Singular Value Decomposition for Collaborative Filtering on a GPU[J].IOP Conference Series:Materials Science and Engineering，2010，10(1)：1-5.

[7]Real-time 3D registration of stereo-vision based range images using GPU[C].IEEE Workshop on Applications of Computer Vision(WACV)，2009：1-6.

[8]張舒，竇衡.基于CUDA的矩陣奇異值分解[J].計算機應用研究，2007，24(6)：104-107.

[9]Brent R P，Luk F.T.The solution of Singular Value and Symmetric Eigen problems on Multiprocessor Arrays[J].Sct Stat Comput，1985，6：69-84.

[10]李小波，薛王偉，孫志勇.一種求解復Hermite矩陣特征值的方法[J].數據采集與處理，2005，20(4)：403-406.