基于GPU的高效并行l(wèi)1最小化算法

高家全,李澤界

(浙江工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，浙江 杭州 310023)

?

基于GPU的高效并行l(wèi)1最小化算法

高家全,李澤界

(浙江工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，浙江 杭州 310023)

多數(shù)l1最小化算法主要由稠密矩陣矢量乘(如Ax和ATx)和矢量運(yùn)算組成.為使其適應(yīng)大數(shù)據(jù)環(huán)境下的性能需求，基于GPU，利用其新的特征,提出了兩個(gè)新穎的基于GPU的并行矩陣矢量乘.這兩個(gè)算法實(shí)現(xiàn)了全局內(nèi)存的合并訪問(wèn)，對(duì)任意給定矩陣，通過(guò)所使用的自適應(yīng)分配線程數(shù)或warp數(shù)的策略，增加了魯棒性.基于這兩個(gè)算法，并以兩個(gè)流行的l1最小化算法為例：快速迭代收縮閾值算法(FISTA)和增廣拉格朗日乘子法(ALM)，提出了兩個(gè)高效基于GPU的并行l(wèi)1最小化算法.實(shí)驗(yàn)結(jié)果驗(yàn)證了提出的算法是高效的，并有良好的性能.

l1最小化；GPU；矩陣矢量乘；FISTA；ALM

l1最小化問(wèn)題為min‖x‖1，須滿足Ax=b，其中A∈(Rm×n(m?n)是一個(gè)滿秩的稠密矩陣，b(Rm是預(yù)先設(shè)定的矢量，x∈(Rn是未知解.l1最小化問(wèn)題的解，也稱為稀疏表達(dá)，已經(jīng)廣泛應(yīng)用到多個(gè)領(lǐng)域，例如信號(hào)處理[1-2]，機(jī)器學(xué)習(xí)[3]和統(tǒng)計(jì)推理[4]等.為了解決l1最小化問(wèn)題，研究者已開發(fā)出許多有效的算法.例如，梯度投影法[5]，截?cái)嗯ｎD內(nèi)點(diǎn)法[6]，同倫法[7]，迭代收縮閾值法[8]，增廣拉格朗日乘子法(ALM)[9]等.隨著問(wèn)題規(guī)模的增長(zhǎng)，算法的執(zhí)行效率很大程度上被消減.為提高其效率，一個(gè)有效的途徑就是將這些算法移植到分布式或多核的架構(gòu)上，例如目前流行的圖形處理單元(GPU).自從NVIDIA公司在2007年[10]介紹了CUDA編程模型以來(lái)，基于GPU加速數(shù)據(jù)處理的研究已成為人們研究的熱點(diǎn)[11-12].

大多數(shù)l1最小化算法，其運(yùn)算主要由稠密的矩陣矢量乘和矢量操作組成.由于CUBLAS庫(kù)[13]中已包含這些運(yùn)算的高效實(shí)現(xiàn)，所以現(xiàn)有基于GPU加速的l1最小化算法[14]主要基于CUBLAS庫(kù).然而，在GTX980顯卡上，對(duì)CUBLAS實(shí)現(xiàn)的矩陣矢量乘進(jìn)行了測(cè)試，發(fā)現(xiàn)其會(huì)產(chǎn)生性能波動(dòng)，隨著行數(shù)或者列數(shù)的增長(zhǎng)，且最大最小性能差距顯著.因此，本研究基于GPU，針對(duì)Ax(GEMV)和ATx(GEMV-T)，分別提出一個(gè)自適應(yīng)warp分配策略的GEMV核和thread分配策略的GEMV-T核.實(shí)驗(yàn)結(jié)果顯示：該方法比CUBLAS更加魯棒，而且總能保持高性能.基于GEMV核和GEMV-T核，以快速迭代收縮閾值算法FISTA和增廣拉格朗日乘子法ALM為例，利用內(nèi)核融合和新的GPU特性，例如洗牌指令(The shuffle instruction)和只讀數(shù)據(jù)內(nèi)存，提出了兩個(gè)高效的基于GPU的并行l(wèi)1最小化算法.該研究工作的主要貢獻(xiàn)：1) 提出了兩種新穎的基于GPU的并行稠密矩陣矢量乘算法；2)利用內(nèi)核融合技術(shù)和新的GPU特性，并基于所提出的并行稠密矩陣矢量算法思想，實(shí)現(xiàn)了兩個(gè)高效的并行l(wèi)1最小化算法；3)通過(guò)實(shí)驗(yàn)，驗(yàn)證了方法的有效性.

1 兩個(gè)l1最小化算法

1.1 快速迭代收縮閾值算法

l1最小化問(wèn)題也被稱為基追蹤問(wèn)題(The basis pursuit problem).在實(shí)際中，一個(gè)測(cè)量數(shù)據(jù)b經(jīng)常含有噪聲(例如測(cè)量誤差ε)，這個(gè)被稱為BPDN問(wèn)題.它的一個(gè)變異被稱為攜帶標(biāo)量λ的無(wú)約束BPDN問(wèn)題：

(1)

快速迭代收縮閾值算法[8](FISTA)是通過(guò)結(jié)合Nesterovs最優(yōu)梯度算法實(shí)現(xiàn)的一個(gè)加速版本，擁有非漸進(jìn)的收斂率O(k2).對(duì)于FISTA，它添加了一個(gè)新的序列{yk,k=1,2,…,n}，即

(2)

式中：soft(u,a)=sign(u)·max{|u|-a,0}為軟閾值算子；y1=x0；t1=1；Lf為▽f(·)關(guān)聯(lián)的李普希茨常數(shù)，可由計(jì)算ATA的特征譜得到(‖ATA‖2).

1.2 增廣拉格朗日乘子法

增廣拉格朗日乘子法[9](ALM)結(jié)合罰函數(shù)法和拉格朗日乘子法，對(duì)應(yīng)增廣拉格朗日函數(shù)為

Lρ(x*,λ)=‖x‖1+λT(Ax-b)+

(3)

式(3)可以通過(guò)迭代求解，即

(4)

式中{ρk}為單調(diào)遞增的正序列.通過(guò)式(4)，可以同時(shí)求解最優(yōu)解x*和λ*.式(4)的第一步稱為x最小化，它是一個(gè)凸優(yōu)化，可以通過(guò)FISTA求解.

2 CUDA架構(gòu)

CUDA是一個(gè)通用的并行計(jì)算平臺(tái)和編程模型.開發(fā)者可以使用CUDA C/C++定義函數(shù)(稱為內(nèi)核)，并通過(guò)CUDA線程并行執(zhí)行這些內(nèi)核.一個(gè)GPU擁有眾多的計(jì)算核心，稱為CUDA核心；這些核心被組織成一個(gè)可擴(kuò)展的流多處理器(SM)的陣列.一個(gè)SM可以并行地執(zhí)行成千上百個(gè)線程.當(dāng)一個(gè)線程塊被分配給一個(gè)SM，它又被劃分為多個(gè)warp，每個(gè)warp擁有32個(gè)線程.最理想的情況，當(dāng)這32個(gè)線程具有相同的執(zhí)行路徑時(shí)，指令執(zhí)行效率最高.

GPU內(nèi)存包括片上內(nèi)存，如共享內(nèi)存和L1緩存，和板載內(nèi)存，如全局內(nèi)存.全局內(nèi)存容量最大，并被所有SM共享，然而由于全局內(nèi)存位于片外的板載DRAM，其低帶寬和高延遲經(jīng)常導(dǎo)致存取訪問(wèn)慢.共享內(nèi)存被同一個(gè)線程塊內(nèi)的所有線程共享，提供高帶寬和低延遲，但容量小.L1緩存和只讀數(shù)據(jù)緩存被同一個(gè)SM內(nèi)的所有線程共享，且存取速度能與共享內(nèi)存一樣快.相比不可控的L1緩存，只讀數(shù)據(jù)緩存是可控的，容易被利用.

3 GPU實(shí)現(xiàn)

3.1 魯棒的矩陣矢量乘

本節(jié)提出兩個(gè)矩陣矢量乘內(nèi)核，包括Ax(GEMV)和ATx(GEMV-T)，其中A∈Rm×n.這里使用以0開始索引的按行存儲(chǔ)方式來(lái)存儲(chǔ)矩陣A，并利用填充策略優(yōu)化全局內(nèi)存訪問(wèn)性能，從而提高了訪問(wèn)效率.

3.1.1 GEMV核

Ax(GEMV)運(yùn)算由m個(gè)內(nèi)積組成(每個(gè)內(nèi)積為A的一行和x相乘)，且每個(gè)內(nèi)積是獨(dú)立的.因此，在設(shè)計(jì)GEMV核時(shí)，可分配一個(gè)warp或者多個(gè)warp去計(jì)算一個(gè)內(nèi)積.為優(yōu)化GEMV核的性能，本研究提出如下的自適應(yīng)warp分配策略，選取k個(gè)warp計(jì)算一個(gè)內(nèi)積為

minw=sm×2 048/(k×32)m≤w

(5)

式中：sm為流多處理器的數(shù)目；m為矩陣的行數(shù).GEMV內(nèi)核偽代碼描述為

輸入:A,x,k,h,m,n;

輸出:b;

1. __shared__xP[],bP[];

2. bVAL←0.0f;

3. lane←threadIdx.x&31; wid←threadIdx.x>>5;

4. wgSize←32×k; wgid←threadIdx.x/wgSize;

5. hub←threadIdx.x&(wgSize-1);

6. row←blockIdx.x×h+wgid;//行索引分配

7. for row tom-1 with row+=h×gridDim.x

//第一個(gè)階段

8. bVAL←0; Aid←row×n+hub;

9. for xid=threadIdx.x ton-1 with xid+=blockDim.x

//x-load步

10.xP[threadIdx.x]←x[xid];

11. __syncthreads();

//partial-reduction步

12. fori=hub to blockDim.x-1 withi+=wgSize

13. bVAL←bVAL+A[Aid]×xP[i];

14.Aid←Aid+wgSize;

15.end

16.__syncthreads();

17.end

//warp-reduction步

18. 使用洗牌指令完成warp內(nèi)的歸約操作

19.iflane==0 bP[wid]←bVAL;

20.__syncthreads();

//第二個(gè)階段

21.ifwid==0

22.bVAL←bP[lane];

23. 使用洗牌指令完成warp內(nèi)的歸約操作

24.iftid%k==0 b[row+threadIdx.x/k]←bVAL;

25.end

26.end

GEMV內(nèi)核的主要步驟包括兩個(gè)階段，第一個(gè)階段有三個(gè)步驟：x-load步,partial-reduction步和warp-reduction步.在x-load步中，一個(gè)線程塊的所有線程并行地讀取矢量x的每個(gè)元素到共享內(nèi)存xP.因?yàn)閤的尺寸很大，所以x被分塊讀進(jìn)共享內(nèi)存，大小為線程塊尺寸.通過(guò)上述方式，保證了x的訪問(wèn)是合并的，并且通過(guò)共享x的片段，減少了訪問(wèn)次數(shù).在partial-reduction步中，每次完成x片段的加載后，同一個(gè)warp組(k個(gè)warp組成一個(gè)warp組)內(nèi)所有線程并行地執(zhí)行x的部分歸約操作(詳見程序描述的12～15).顯然，每個(gè)線程最多執(zhí)行|n/wgSize|(n為矩陣A的列數(shù)，wgSize為warp組內(nèi)線程數(shù))次歸約操作，并且對(duì)于訪問(wèn)全局內(nèi)存上的矩陣A也是合并的.然后，使用洗牌指令實(shí)現(xiàn)warp-reduction步(warp內(nèi)歸約操作)，并將結(jié)果存在共享內(nèi)存上.

在第二階段，將warp-reduction步的結(jié)果歸約為最終的輸出.當(dāng)k=1時(shí)，一個(gè)warp組只有一個(gè)warp.只需要第一個(gè)階段，便可將結(jié)果輸出(因?yàn)槠耸÷?.同時(shí)，對(duì)于k=32的情況，由于不需要加載矢量x到共享內(nèi)存，可直接將x加載到寄存器，所以需一個(gè)新核(因?yàn)槠摵耸÷?.

3.1.2 GEMV-T內(nèi)核

ATx(GEMV-T)運(yùn)算由n個(gè)內(nèi)積組成(每個(gè)內(nèi)積為A的一列和x相乘)，且每個(gè)內(nèi)積是獨(dú)立的.由于GEMV-T中的矢量x尺寸較小，所以在設(shè)計(jì)GEMV-T核時(shí)，可分配一個(gè)線程或者多個(gè)線程去計(jì)算一個(gè)內(nèi)積.為優(yōu)化GEMV-T核性能，本研究提出如下的自適應(yīng)thread分配策略，選取k個(gè)線程計(jì)算一個(gè)內(nèi)積為

mint=sm×2 048/kn≤t

(6)

式中：sm為流多處理器的數(shù)目；m為矩陣的行數(shù).GEMV-T內(nèi)核偽代碼描述為

輸入:A,x,k,h,m,n;

輸出:b;

1.__shared__xP[], bP[];

2.bVAL;

3.wid←threadIdx.x>>5;lane←threadIdx.x&31;

4.wgSize←32×h;wgid←threadIdx.x/(wgSize);

5.col←blockIdx.x×32+32×wgid+lane; //列索引分配

6.forcolton-1withcol+=k×32×gridDim.x

//第一個(gè)階段

7.bVAL←0;Aid←col+wid%h×n;

8.forxid=threadIdx.xtom-1withxid+=blockDim.x

// x-load步

9. xP[threadIdx.x]←x[xid];

10. __syncthreads();

11. //partial-reduction步

12. fori= wid%hto blockDim.x-1 withi+=h

13. bVAL←bVAL+A[Aid]×xP[i];

14. Aid←Aid+h×n;

15. end

16. __syncthreads();

17. end

18.bP[tid]←bVAL;

19. __syncthreads();

//第二個(gè)階段

20. if wid%h==0

21. for bPid=tid to tid+32×h-1 with bPid+=32

22. bVAL←bVAL+bP[bPid];

23. end

24.b[col]←bVAL;

25. end

26. end

GEMV-T內(nèi)核的主要過(guò)程也包括兩個(gè)階段，第一個(gè)階段有兩個(gè)步驟：x-load步和partial-reduction步.其中，x-load步跟GEMV內(nèi)核執(zhí)行一樣的功能.在partial-reduction步中，由于使用以0開始索引的按行存儲(chǔ)方法存儲(chǔ)矩陣，因此如果線程組(k個(gè)線程組成一個(gè)線程組)的組織不合理，對(duì)于全局內(nèi)存上的矩陣A的訪問(wèn)將是不合并的.

因此，為確保訪問(wèn)是合并的，線程組的創(chuàng)建根據(jù)如下：

定義1 假設(shè)線程塊大小是s，h個(gè)線程一起被分配給ATx中的一個(gè)內(nèi)積，而且z=s/h.那么，線程組組織如下：{0,z,…,(h-1)×z},{1,z+1,…,(h-1)×z+1},…,{z-1,2×z-1,…,2×(h-1)×z-1}.

對(duì)于已經(jīng)完成加載的x片段，同一個(gè)線程組內(nèi)所有線程并行地執(zhí)行歸約操作，與GEMV核相似.由于同一個(gè)線程組內(nèi)線程經(jīng)常不屬于同一個(gè)warp，所以不能使用洗牌指令實(shí)現(xiàn)warp內(nèi)歸約操作.因此，在第二階段，將partial-reduction步結(jié)果保存在共享內(nèi)存上，然后再進(jìn)行歸約，并將歸約結(jié)果輸出.

3.2 并行的FISTA和ALM

在FISTA中，完成每次迭代需要6個(gè)操作，如圖1所示.為了減少對(duì)應(yīng)的內(nèi)核數(shù)目，節(jié)約內(nèi)核調(diào)用開銷，避免雙倍數(shù)據(jù)加載，這里對(duì)這些內(nèi)核進(jìn)行了融合.如圖1所示，內(nèi)核1融合了GEMV和矢量操作1；內(nèi)核3融合了余下的矢量操作.進(jìn)而，通過(guò)融合技術(shù)，核總數(shù)從6個(gè)降到3個(gè).對(duì)于圖1中的內(nèi)核，基于GEMV核和GEMV-T核的思想，很容易實(shí)現(xiàn)內(nèi)核1和內(nèi)核2.CUBLAS雖然有著矢量運(yùn)算的高效實(shí)現(xiàn)，但無(wú)法實(shí)現(xiàn)多核合并.因此，這里采用文獻(xiàn)[15]的方法，通過(guò)多核融合，實(shí)現(xiàn)了內(nèi)核3.

在ALM中，其x最小化可以通過(guò)FISTA求解，因此只需設(shè)計(jì)新的核函數(shù)實(shí)現(xiàn)dual update操作λk+1=λk+ρk(Axk+1-b).顯而易見，基于GEMV核的思想，很容易并行實(shí)現(xiàn)這個(gè)dualupdate操作.這里分別稱并行的FISTA和ALM為GFISTA和GALM.

圖1 內(nèi)核融合Fig.1 Kernel fusion

3.3 優(yōu) 化

隨著l1最小化算法的不斷迭代，軟閾值操作使得矢量x變得越來(lái)越稀疏.因此對(duì)于Ax，當(dāng)輸出向量x的第i個(gè)元素為0，那么A的第i行就不需要被訪問(wèn)了.利用這個(gè)稀疏性可減少對(duì)于全局內(nèi)存上矩陣A的訪問(wèn)次數(shù).

4 實(shí) 驗(yàn)

這一節(jié)測(cè)試所提出來(lái)的GEMV核、GEMV-T核、GFISTA和GALM的有效性.實(shí)驗(yàn)環(huán)境如下：一臺(tái)IntelXeon雙核CPU配一個(gè)NVIDIAGTX980顯卡，使用CUDA7.5的編譯運(yùn)行環(huán)境.所有測(cè)試使用同一個(gè)測(cè)試矩陣集，如表1所示，這些矩陣按照正態(tài)分布隨機(jī)生成.

表1 測(cè)試矩陣

首先，使用上述測(cè)試矩陣與CUBLAS庫(kù)進(jìn)行比較，來(lái)驗(yàn)證所提出的GEMV核和GEMV-T核的有效性和魯棒性.圖2，3分別顯示了兩個(gè)并行稠密矩陣矢量乘內(nèi)核與CUBLAS的性能比較.從圖2可見：GEMV核相比CUBLAS，對(duì)于不同規(guī)模的矩陣，始終能夠維持較高的性能，具有更佳的性能及穩(wěn)定性.類似于GEMV核，GEMV-T核相比CUBLAS有著相似的結(jié)論.

圖2 GEMV核與CUBLAS的性能比較Fig.2 Performance comparison of GEMV and CUBLAS

圖3 GEMV-T核與CUBLAS的性能比較Fig.3 Performance comparison of GEMV-T and CUBLAS

其次，測(cè)試GFISTA和GALM的性能.為每一個(gè)測(cè)試用例，初始x0總是含有1 024個(gè)非零元素且b=Ax0.對(duì)于GFISTA，50次迭代后終止；對(duì)于GALM，10次外迭代后終止，并且每一次外迭代，內(nèi)迭代都執(zhí)行50次.表2列出了所有算法的執(zhí)行時(shí)間，時(shí)間單位是秒.其中CFISTA和CALM分別是對(duì)應(yīng)GFISTA和GALM在CPU上實(shí)現(xiàn)的串行FISTA和ALM算法.從表2可以看出：對(duì)所有用例，相比CFISTA，GFISTA能夠獲得范圍從37.68到53.66倍的加速比，平均加速比是48.22.相比CALM，GALM獲得的最大、最小和平均加速比分別為51.21，35.98，44.0.這些結(jié)果顯示2個(gè)并行算法的性能提升是顯著的.

表2 GFISTA和GALM的性能

5 結(jié) 論

提出了兩個(gè)新穎的GPU并行的矩陣矢量乘，進(jìn)而基于這兩個(gè)并行的矩陣矢量乘算法，利用內(nèi)核融合技術(shù)和新的GPU特性，實(shí)現(xiàn)了兩個(gè)高度并行的l1最小化算法(GFISTA和GALM)，并利用解的稀疏性進(jìn)一步優(yōu)化了性能.實(shí)驗(yàn)結(jié)果表明：提出的方法是有效的和魯棒的，具有較高的并行性.

[1] DONOHO D L, ELAD M. On the stability of the basis pursuit in the presence of noise[J]. Signal processing,2006,86(3):511-532.

[2] TROPP J A. Just relax: convex programming methods for identifying sparse signals in noise[J]. IEEE transactions on information theory,2006,52(3):1030-1051.

[3] ELHAMIFAR E,VIDAL R. Sparse subspace clustering: algo-

rithm, theory, and applications[J]. IEEE transactions on software engineering,2013,35(11):2765-2781.

[4] WRIGHT J, MA Y, MAIRAL J, et al. Sparse representation for computer vision and pattern recognition[J]. Proceedings of the IEEE,2010,98(6):1031-1044.

[5] FIGUEIREDO M A T, NOWAK R D, WRIGHT S J. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems[J]. IEEE journal of selected topics in signal processing,2007,1(4):586-597.

[6] KIM S J, KOH K, LUSTIG M, et al. An interior-point method for large-scalel1-regularized least squares[J]. IEEE journal of selected topics in signal processing,2007,1(4):606-617.

[7] DONOHO D L, TSAIG Y. Fast solution of-norm minimization problems when the solution may be sparse[J]. IEEE transactions on information theory,2008,54(11):4789-4812.

[8] BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. Siam journal on imaging sciences,2009,2(1):183-202.

[9] YANG A Y, ZHOU Z H, Ganesh A, et al. Fastl1-minimization algorithms for robust face recognition[J]. IEEE transactions on image processing,2013,22(8):3234-3246.

[10] NVIDIA Corporation. CUDA C Programming guide[EB/OL].[2015-12-05].http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html.

[11] 范菁,吳佳敏.帶顯著性區(qū)域約束的高效視頻全景拼接方法[J].浙江工業(yè)大學(xué)學(xué)報(bào),2015,43(5):479-486.

[12] 丁維龍,陳淑嬌,吳福理.一種基于物理特性的雨滴飛濺模擬算法[J].浙江工業(yè)大學(xué)學(xué)報(bào),2014,42(5):586-590.

[13] NVIDIA Corporation. CUBLAS library[EB/OL].[2015-12-05]. http://docs.nvidia.com/cuda/cublas/index.html.

[14] NAGESH P, GOWDA R, LI B. Fast GPU implementation of large scale dictionary and sparse representation based vision problems[C]//Acoustics Speech and Signal Processing. Dallas: IEEE International Conference,2010:1570-1573.

[15] GAO J Q, LIANG R H, WANG J. Research on the conjugate gradient algorithm with a modified incomplete Cholesky preconditioner on GPU[J]. Journal of parallel and distributed computing,2014,74(2):2088-2098.

(責(zé)任編輯：陳石平)

Highly parallell1-minimization algorithms on GPU

GAO Jiaquan, LI Zejie

(College of Computer Science and Technology, Zhejiang University of Technology, Hangzhou 310023, China)

Thel1-minimization (l1-min) algorithms for solving thel1-min problem have been widely developed. For most ofl1-min algorithms, their components are the dense matrix-vector multiplications such as Ax and ATx, and vector operations. Utilizing the new features of the GPU, two novel matrix-vector multiplications are proposed. The two algorithms access the global memory in a fully coalesced manner, and enhance their robustness by designing a self-adaptive thread (warp) allocation strategy. Besides, based on the two algorithms, two highly parallell1-min algorithms (FISTA and ALM) are presented on the GPU. The experimental results show that the methods proposed are high efficent and have good performance.

l1-minimization; GPU; matrix-vector multiplication; FISTA; ALM

2015-12-17

國(guó)家自然科學(xué)基金資助項(xiàng)目(61379017)

高家全(1972—)，男，河南固始人，副教授，博士，研究方向?yàn)楦咝阅苡?jì)算、大數(shù)據(jù)分析與可視化、智能算法及應(yīng)用，E-mail：gaojq@zjut.edu.cn.

TP338.6

A

1006-4303(2016)05-0495-06