基于GPU的LLE算法加速及性能優(yōu)化

李 繁，嚴(yán) 星，張曉宇

(1.新疆財(cái)經(jīng)大學(xué) 網(wǎng)絡(luò)與實(shí)驗(yàn)教學(xué)中心，新疆 烏魯木齊 830012；2.新疆財(cái)經(jīng)大學(xué) 信息管理學(xué)院，新疆 烏魯木齊 830012)

0 引 言

LLE數(shù)據(jù)降維算法是利用局部的信息計(jì)算相關(guān)權(quán)重，從而產(chǎn)生一個(gè)大型稀疏矩陣，通過(guò)解大型稀疏矩陣的特征值來(lái)解決非線性數(shù)據(jù)降維的問(wèn)題。在LLE算法提出后又有一些改進(jìn)的LLE算法提出：如HLLE算法結(jié)合了LLE和拉普拉斯特征映射來(lái)修改LLE的第二個(gè)步驟，并采用Hessain來(lái)評(píng)估測(cè)量權(quán)重，從而避免LLE算法存在的一些out-of-sample問(wèn)題[1]。MLLE算法也是修改LLE的重建權(quán)重部分，通過(guò)引進(jìn)對(duì)每一個(gè)鄰點(diǎn)計(jì)算多個(gè)獨(dú)立的線性權(quán)重來(lái)修改LLE的第二個(gè)步驟[2]，這個(gè)方法一樣是在改進(jìn)LLE計(jì)算權(quán)重時(shí)的正確性或可能會(huì)發(fā)生的一些問(wèn)題。所以上面兩個(gè)算法都是修改LLE的第二個(gè)步驟，來(lái)提高LLE算法的準(zhǔn)確性。

在探討LLE算法的準(zhǔn)確性問(wèn)題時(shí)，更改的通常都只有在重建權(quán)重的那個(gè)部分[3]，由于在任何LLE變種算法中都要計(jì)算KNN和解決稀疏矩陣特征值這兩個(gè)問(wèn)題，所以在整體的計(jì)算效率上，LLE已經(jīng)通過(guò)KNN的這個(gè)步驟來(lái)選擇與數(shù)據(jù)點(diǎn)接近的K個(gè)點(diǎn)來(lái)保持局部的信息，并將不在鄰近的數(shù)據(jù)點(diǎn)權(quán)重設(shè)為0，這樣一來(lái)通過(guò)重建權(quán)重的時(shí)候就可以產(chǎn)生一個(gè)大型的稀疏矩陣，產(chǎn)生大型的稀疏矩陣可降低后面特征值的計(jì)算量，所以通過(guò)KNN的這個(gè)步驟可以大量減少后面特征值的計(jì)算量，這也就是為何LLE的數(shù)據(jù)降維技術(shù)比其它的方法計(jì)算速度還快，不過(guò)就算KNN這步降低了特征值的計(jì)算量，其大型稀疏矩陣求特征值也存在效率上的問(wèn)題，所以要加快整體的LLE計(jì)算效率就要先從KNN及稀疏矩陣這兩部分來(lái)加速。

1 KNN算法優(yōu)化

KNN算法這部分在LLE數(shù)據(jù)降維技術(shù)中占很重要的地位，就是因?yàn)橛肒NN可以使LLE后面的計(jì)算量大幅降低，只需根據(jù)KNN選出來(lái)的鄰點(diǎn)，并計(jì)算權(quán)重就會(huì)產(chǎn)生一個(gè)大型稀疏矩陣來(lái)降低運(yùn)算元素[4-6]。LLE算法保存鄰點(diǎn)的權(quán)重，最主要的目的是只需要保持局部的信息而不用計(jì)算整個(gè)全局?jǐn)?shù)據(jù)點(diǎn)，利用保持局部狀態(tài)的數(shù)據(jù)點(diǎn)投射到低維度空間。因?yàn)槊恳粋€(gè)數(shù)據(jù)點(diǎn)有權(quán)重的在高維空間都是鄰點(diǎn)，所以此方法不會(huì)破壞非線性特性。接下來(lái)我們將探討如何在通用計(jì)算圖形處理器(GPGPU)上加速KNN算法，以及加速后的GPGPU算法可以比正常的KNN算法加快多少。

1.1 基于GPU的Brute force KNN

Brute force KNN大致可分成3個(gè)部分：

(1)計(jì)算每一點(diǎn)的平面距離；

(2)利用一個(gè)排序算法對(duì)平面距離進(jìn)行排序；

(3)挑出排序完成后的前K個(gè)數(shù)據(jù)點(diǎn)并回傳。

整個(gè)計(jì)算流程的時(shí)間復(fù)雜度是O(NlogN)， 其中N是我們所輸入數(shù)據(jù)的筆數(shù)。可以看出，在完整的計(jì)算過(guò)程中最耗費(fèi)時(shí)間的就是第一和第二個(gè)步驟，而這兩個(gè)步驟也是并行度最高的兩個(gè)步驟，我們將針對(duì)這兩個(gè)步驟進(jìn)行高度的并行計(jì)算。

首先，先定義平面距離公式，假設(shè)我們有兩組Data set，一個(gè)是參照數(shù)據(jù)點(diǎn)R，另一組數(shù)據(jù)集合為搜索數(shù)據(jù)點(diǎn)Q。而平面距離公式是計(jì)算Q到R矩陣?yán)锩恳恍兄g的距離

(1)

因此根據(jù)我們的平面距離公式來(lái)計(jì)算兩個(gè)數(shù)據(jù)集合內(nèi)數(shù)據(jù)點(diǎn)的平面距離。由于LLE是要查找我們輸入的Data set內(nèi)每一點(diǎn)的KNN，所以R和Q都是一樣的二維矩陣，并且Q矩陣內(nèi)的每一行元素都要對(duì)R矩陣的每一行元素做平面距離計(jì)算，因此我們?cè)谶@的時(shí)間復(fù)雜度是O(N2M)，M是數(shù)據(jù)維度。

ComputeKNNdistanceonCPU(Pseudocode)

Input：RandQmatrix(DATA， SIZE， Dimension)

Output： DistanceDmatrix

For i=1 to SIZE

For j=1 to SIZE

For k=1 to Dimension

Dist=∑((ri-qj)×(ri-qj));

End

D(i,j)=sqrt(Distance);

End

End

由此可以看出數(shù)據(jù)點(diǎn)的距離計(jì)算只是多個(gè)循環(huán)迭代重復(fù)計(jì)算，因此在計(jì)算距離這個(gè)計(jì)算函數(shù)具有高度的可并行度，所以我們將這個(gè)計(jì)算函數(shù)放置到GPU的計(jì)算單元上，利用GPU高度并行計(jì)算的特性來(lái)大幅降低此計(jì)算函數(shù)所需要的執(zhí)行時(shí)間，以加快KNN算法。

ComputeKNNdistanceonGPU(Pseudocode)

Input：RandQmatrix(DATA， SIZE， Dimension)

Output： DistanceDmatrix

__shared__ float shared_R[BLOCK_DIM][ BLOCK_DIM];

__shared__ float shared_Q[BLOCK_DIM][ BLOCK_DIM];

int threadx=threadIdx.x;

int thready=threadIdx.y;

begin_R=BLOCK_DIM*blockIdx.y;

begin_Q=BLOCK_DIM*blockIdx.x;

step_R=BLOCK_DIM*R;

step_Q=BLOCK_DIM*Q;

end_R=begin_R+(dim-1)*R;

int cond1=(begin_Q+tx

int cond2=(begin_R+ty

For (i=begin_R, j=begin_Q; i to end_R; i+=step_R, Q+=step_Q)

if (cond2 && cond1)

For (k=0; k

tmp=shared_R[k][ty]-shared_Q[k][tx];

ssd+=tmp*tmp;

End

End if

__syncthreads();

End

If (cond2 && cond1)

D[ (begin_R+thready)*R+begin_Q+threadx ]=sqrt(ssd);

以上就是我們將此計(jì)算函數(shù)放到GPGPU計(jì)算單元上的偽碼，我們將R、Q、D矩陣設(shè)為共享的變量，并且利用Block Size和Block ID來(lái)計(jì)算我們?cè)赗矩陣和Q矩陣的起始位和邊界避免內(nèi)存溢出，并利用cond1和cond2來(lái)控制GPU的內(nèi)存存取，如果不符合cond1和cond2的話就無(wú)法對(duì)該內(nèi)存區(qū)塊進(jìn)行存取，整個(gè)GPGPU的計(jì)算處理方法和內(nèi)存的配置方式如圖1所示。

圖1 GPU計(jì)算架構(gòu)

由圖1我們可以得知在GPU的運(yùn)算架構(gòu)上有多個(gè)處理器(Processor)，每一個(gè)處理器都有自己的寄存器(Register)，第一層的高速緩存(Cache)Constant cache是一個(gè)只讀(Read-only)的高速緩存，Constant cache是一個(gè)8 KB的SM芯片，在CUDA的程序架構(gòu)中可以用來(lái)做優(yōu)化的內(nèi)存存取，接著第二層是Texture cache，這是用來(lái)支持二維(two dimension)內(nèi)存的Locality，讓整個(gè)GPU的內(nèi)存存取更有效率，之后最外層才是Device memory，也就是一般GPU的數(shù)據(jù)存放的內(nèi)存區(qū)塊，接下來(lái)GPU在Thread的運(yùn)用中Kernel可以看成一個(gè)grid，每個(gè)grid內(nèi)都有多個(gè)block thread，因此整個(gè)內(nèi)存的配置上要搭配block thread的利用才可以使GPU運(yùn)算更加有效率。

接下來(lái)KNN算法除了之前的距離計(jì)算公式上有很大的計(jì)算負(fù)擔(dān)之外，還有另一個(gè)需要大量計(jì)算效能的函數(shù)，就是利用一個(gè)排序算法來(lái)對(duì)平面距離進(jìn)行排序，而排序有很多種方法，像插入排序(Insertion sort)、泡沫排序(Bubble sort)、基數(shù)排序(Radix sort)、快速排序(Quick sort)等等，因此我們必須要選擇一個(gè)適合用來(lái)并行的排序算法來(lái)進(jìn)行KNN的排序。根據(jù)上述的排序算法，Quick sort這個(gè)方法在CPU的計(jì)算是快速排序的平均時(shí)間復(fù)雜度比較低O(NlogN)， 但因?yàn)樗怯眠f歸的方式來(lái)運(yùn)作，所以沒(méi)辦法進(jìn)行有效率的并行運(yùn)算，因此我們將利用其它的排序方法來(lái)加快我們的排序。相比較來(lái)說(shuō)，雖然插入算法在CPU的運(yùn)算上具有比較低的運(yùn)算效率O(N2)， 但因?yàn)镮nsertion sort是用iteration排序的方法并且計(jì)算之間是獨(dú)立的(Independent)，所以在并行度上來(lái)講插入算法是比較適合用GPU來(lái)進(jìn)行運(yùn)算加速的，所以我們?cè)诒┝NN算法中用GPU加快了數(shù)據(jù)點(diǎn)的距離計(jì)算和排序函數(shù)，接下來(lái)我們要探討另外一種方法，也就是KD tree on GPU。

1.2 基于GPU的KD tree

在前一節(jié)已經(jīng)探討了暴力KNN算法(Brute force KNN)如何用GPGPU來(lái)進(jìn)行加速，還有另一種方法是更有效率的搜索KNN，那就是KD tree的數(shù)據(jù)結(jié)構(gòu)(Data structure)[7]，我們將探討KD tree在GPU的計(jì)算架構(gòu)上是否一樣可以發(fā)揮出卓越的效能，是否會(huì)有不一樣的效果，因此我們研究并探討如何并行KD tree，雖然KD tree是用遞歸(Recursive)的方式來(lái)建構(gòu)，遞歸算法在GPU上并不容易并行，但根據(jù)之前的研究與探討，KD tree還是可以并行的，如圖2所示。接著我們將探討如何在GPU上進(jìn)行KD tree的并行運(yùn)算及是否可以加速。

圖2 GPU KD tree構(gòu)造

KDTreeConstructGPUAlgorithm

Input： KD Tree construct(data_list,number_data_point)

(1) According data_list find median number

(2) According Median_number to spilt data_list

Spilt_value=Median_number;

Right_list=all element in data_list smaller than Spilt_value;

right_ptr=Right_list;

Left_list=all element in data_list larger than Spilt_value;

left_ptr=Left_list;

(3) Recursive build subtree

KD Tree construct(Right_list, number_data_point_Right_list);

KD Tree construct(Left_list, number_data_point_Left_list);

由此可以看出，雖然KD tree是用遞歸的方式來(lái)建構(gòu)，在GPU上實(shí)現(xiàn)并行處理相對(duì)比較難，但在分割的時(shí)候我們還是可以用GPU來(lái)實(shí)現(xiàn)Radix sort的排序算法來(lái)降低檢索中位數(shù)的時(shí)間復(fù)雜度，然后在第(2)步會(huì)用一個(gè)循環(huán)進(jìn)行KD tree分割的部分，所以在此步會(huì)用比較元素和中位數(shù)來(lái)建立左右List，這個(gè)部分我們也可以用GPU來(lái)分類(lèi)，將整個(gè)Data set切成左右兩個(gè)Sub data set，所以在KD tree的建構(gòu)上我們還是可以用GPU來(lái)找尋中位數(shù)和分割元素兩個(gè)步驟的加速運(yùn)算，而接下來(lái)我們將探討如何利用GPU進(jìn)行KD tree的 KNN搜索，如圖3所示。

圖3 GPU KNN搜索

KDTreeSearchGPUAlgorithm

Input： KD tree結(jié)構(gòu)， 搜索點(diǎn)(Query point)

Output： Nearest Neighbor

If tree_node.left=tree_node.right=Null {leaf case}

Use GPU parallel compute temp_distance

If temp_distance

NN=tree.point;shortest=temp_distance

else

if query_point.axis tree_nod.value then

search_first=left

else

search_first=right

if search_first==left

if query_point.axis-shortest tree_nod.value then

Recursive(query_point, tree_nod.left, w)

if query_point.axis-shortest>tree_nod.value then

Recursive(query_point, tree_nod.right, w)

Else

if query_point.axis-shortest>tree_nod.value then

Recursive(query_point, tree_nod.right, w)

if query_point.axis-shortest tree_nod.value then

Recursive(query_point, tree_nod.left, w)

在GPU KNN的搜索當(dāng)中由于需要通過(guò)Recursive的方式進(jìn)行搜索，所以后面的Search需要用到前面的數(shù)據(jù)進(jìn)行搜索，所以我們?cè)谶@個(gè)部分無(wú)法進(jìn)行并行運(yùn)算，因此我們只能將GPU的運(yùn)算處理加在數(shù)據(jù)點(diǎn)之間的計(jì)算距離上，而在將要進(jìn)行的KNN搜索中，我們要將此數(shù)據(jù)點(diǎn)當(dāng)作最近的數(shù)據(jù)點(diǎn)時(shí)也要先查看在已搜索到的KNN當(dāng)中是否已存在，如果不存在就將此點(diǎn)視為最短的數(shù)據(jù)點(diǎn)，否則將此數(shù)據(jù)點(diǎn)忽略掉。

1.3 KNN算法優(yōu)化結(jié)果

我們將以圖表的方式來(lái)呈現(xiàn)兩個(gè)算法在GPU上執(zhí)行效率相對(duì)于CPU的執(zhí)行效率到底加快了多少。針對(duì)輸入數(shù)據(jù)隨維度增加及效果數(shù)據(jù)點(diǎn)數(shù)的增加判斷算法的加速效果，如圖4、表1所示。

圖4 按維度增加執(zhí)行時(shí)間倍數(shù)關(guān)系

表1 基于維度的執(zhí)行時(shí)間比較(數(shù)據(jù)點(diǎn)=2048，單位/s)

由圖4及表1可以看出，在維度不斷增加的過(guò)程中，在GPU與CPU不同的計(jì)算平臺(tái)上執(zhí)行時(shí)間的距離差距就越來(lái)越遠(yuǎn)，時(shí)間倍數(shù)甚至差到好幾十倍，因此在GPU的運(yùn)算上KNN的算法大幅提升了效能。接下來(lái)我們將看如果增加數(shù)據(jù)點(diǎn)數(shù)的情況下效能會(huì)如何提升。

如圖5及表2所示，不論是以數(shù)據(jù)點(diǎn)數(shù)增加，還是以數(shù)據(jù)維度增加，在GPU上的算法都大幅提升了KNN的執(zhí)行效率，并且在Data set或維度加大時(shí)，基于GPU的KD tree執(zhí)行都非常高效，因此我們將用基于GPU的KD tree方法來(lái)提高LLE算法在求KNN的部分，相信會(huì)大幅提高LLE的效率。

表2 基于數(shù)據(jù)點(diǎn)數(shù)的執(zhí)行時(shí)間比較(維度=16，單位/s)

圖5 按數(shù)據(jù)點(diǎn)數(shù)增加執(zhí)行時(shí)間倍數(shù)關(guān)系

2 稀疏矩陣特征值求解優(yōu)化

在LLE算法中另一個(gè)需要耗費(fèi)大量運(yùn)算時(shí)間的就是解

大型稀疏矩陣的特征值，由于前一步計(jì)算權(quán)重時(shí)是根據(jù)KNN選出K個(gè)鄰點(diǎn)，同時(shí)將非鄰點(diǎn)之間的權(quán)重設(shè)為0，所以經(jīng)過(guò)權(quán)重重建就會(huì)形成一個(gè)大型的稀疏矩陣，這樣就有效降低了運(yùn)算量。但以大型矩陣要解特征值還是比較耗時(shí)間，因?yàn)樵谇蠼獾倪^(guò)程中需要用到矩陣相乘以及矩陣與向量相乘[8]，這兩個(gè)部分都需要大量的運(yùn)算效能，所以我們將GPU的運(yùn)算平臺(tái)放在這兩個(gè)部分。接下來(lái)，我們將先探討如何配置內(nèi)存可以使GPU并行計(jì)算發(fā)揮最大效果。

2.1 GPU內(nèi)存分配

大型稀疏矩陣的儲(chǔ)存方式，有坐標(biāo)方式(COO)、對(duì)角線壓縮方式(DIA)、行壓縮方式(CSR)等等[9,10]，CSR的儲(chǔ)存方式是最節(jié)省空間且最有效率的，同時(shí)也是最容易用來(lái)直接存取而不需要用任何偏移或計(jì)算公式來(lái)對(duì)應(yīng)回原始位置。但若要將整個(gè)CSR格式放在GPU的內(nèi)存存取方式來(lái)提升GPU的效率還需要做稍微的修改，而在GPU的存取方式有下列幾種方法。

2.1.1 每行分配一個(gè)block

這種方法是一個(gè)很直觀的方法，將CSR的儲(chǔ)存方式每行當(dāng)作一個(gè)block來(lái)計(jì)算，然后block內(nèi)的每一個(gè)thread讀該行數(shù)據(jù)內(nèi)的一個(gè)元素，并且計(jì)算乘積，再利用并行規(guī)約的方法將每個(gè)計(jì)算的結(jié)果加起來(lái)，一直到block計(jì)算完整行的元素，這種方法并行程度較高，并且相鄰的thread就會(huì)存取相鄰的全局內(nèi)存。

雖然有較高的并行度，且存取內(nèi)存也是比較有效率的，但這種方法有兩個(gè)缺點(diǎn)，一是每個(gè)block要有較多的thread才可以得到比較好的性能，但CSR每行的元素個(gè)數(shù)是不確定的，因此在block里有一thread是空運(yùn)算的，所以會(huì)造成計(jì)算能力的浪費(fèi)。二是并行規(guī)約時(shí)需要同步，所以會(huì)延遲時(shí)間。因此有一個(gè)方法可以解決這種問(wèn)題，就是每行分配一個(gè)warp，這可以增加Thread的運(yùn)算能力。

2.1.2 每行分配一個(gè)warp

每行分配一個(gè)warp 可以順利解決block內(nèi)thread空轉(zhuǎn)的問(wèn)題。雖然這種方法可以解決block 內(nèi)thread 空轉(zhuǎn)的問(wèn)題，但到了行尾還是會(huì)有運(yùn)算元素不足的問(wèn)題，因?yàn)榍懊嬲f(shuō)過(guò)大型的稀疏矩陣不是每一行的元素都相同，所以會(huì)有些運(yùn)算元素不轉(zhuǎn)的問(wèn)題，我們將會(huì)利用下一個(gè)方法來(lái)改進(jìn)相關(guān)問(wèn)題。

2.1.3 對(duì)齊計(jì)算元素

前兩個(gè)方法不管是把每行看成是一個(gè)Block還是每行看成一個(gè)warp都會(huì)有運(yùn)算元素不足的問(wèn)題，因此，我們通過(guò)在行尾補(bǔ)0對(duì)齊的warp_size方法來(lái)解決這個(gè)問(wèn)題。我們將行尾補(bǔ)0，這樣一來(lái)就不會(huì)發(fā)生non-coalesced元素運(yùn)算。因此，我們將利用這個(gè)方法來(lái)儲(chǔ)存GPU的Sparse matrix。接著，我們來(lái)介紹SpMV的算法。

2.2 并行Krylov子空間

Krylov子空間是解一個(gè)大型矩陣特征值的算法，因此我們利用這個(gè)算法來(lái)解我們的大型矩陣特征值問(wèn)題，并考慮如何使用Arnoldi這個(gè)算法來(lái)進(jìn)行GPU并行，Arnoldi可以并行的部分有Span Krylove subspace、正交化(Orthogonalize)及標(biāo)準(zhǔn)化(Normalize)。

2.2.1 Span Krylove subspace

在擴(kuò)張Krylov subspace這部分其實(shí)就是一個(gè)簡(jiǎn)單的大型稀疏矩陣與向量相乘(SpMV)的方法，因此，我們將研究如何利用GPU提升這個(gè)計(jì)算函數(shù)的計(jì)算效率。

SpMVCPUpseudocode

Input： SpMV (row_ptr,value,index_ptr,vc)

Output： Result

for(i=0; i

dot=0;

for(j=row_ptr[i]; j

dot+=value[j] * vc[index_ptr[j]];

end

(*Result)[i]=dot;

End

在SpMV CPU函數(shù)中我們可以發(fā)現(xiàn)它有兩層循環(huán)，而且之間的元素相乘都是獨(dú)立運(yùn)算的，所以在這方面很適合GPU來(lái)進(jìn)行并行運(yùn)算。因此，我們將考慮如何將SpMV的算法放在GPU的平臺(tái)上進(jìn)行運(yùn)算。

SpMVGPUpseudocode

Input： SpMV (row_ptr,value,index_ptr,vc)

Output： Result

__shared__ float dots[BLOCK_SIZE];

const int tid=threadIdx.x;

const int row=(blockIdx.x*BLOCK_SIZE+tid)/HALF_WARP;

const int lane=tid&HALF_WARP_SB_1;

int begin, end, i;

dots[tid]=0;

if(row

for (i=begin+lane; i

dots[tid]+=val[i]*tex1Dfetch(tex_float, ind[i]);

}

if (lane<8) dots[tid]+=dots[tid+8];

if (lane<4) dots[tid]+=dots[tid+4];

if (lane<2) dots[tid]+=dots[tid+2];

if (lane<1) {dots[tid]+=dots[tid+1]; Result[row]+=dots[tid];}

}

我們將SpMV的計(jì)算函數(shù)放在GPU的平臺(tái)上進(jìn)行運(yùn)算，而Thread和Block數(shù)據(jù)分割的方式如之前內(nèi)存分配部分所述，利用GPU對(duì)內(nèi)存存取的方式來(lái)將GPU的并行度提升至最高。

除了SpMV的大量計(jì)算需要用到GPU來(lái)進(jìn)行并行運(yùn)算之外，還有正交化的計(jì)算函數(shù)需要用到GPU來(lái)并行計(jì)算，接下來(lái)我們將介紹如何利用GPU進(jìn)行正交化計(jì)算。

2.2.2 正交化

在這部分我們來(lái)討論在Arnoldi算法中除了SpMV那個(gè)部分花了很多時(shí)間計(jì)算之外，另一個(gè)花大量時(shí)間的地方就是在正交化處理的過(guò)程中，而在正交化處理所花費(fèi)時(shí)間最多的地方就是矩陣向量相乘，在這個(gè)部分我們將利用CUDA內(nèi)建的Cublas函數(shù)庫(kù)來(lái)進(jìn)行處理加速，而以下就是在CPU的偽碼，以及用Cublas函數(shù)庫(kù)寫(xiě)的程序。

Matrixmultiplicationvectorpseudocode

Input： Matrix A[][]、 Vector V[]

Output： Vector Result[]

for(j=0; j

for(k=0;k

Result[j]+=A [k][j]*V[k];

End

End

接下來(lái)我們利用Cublas函數(shù)庫(kù)重新寫(xiě)上方的偽碼，我們只簡(jiǎn)單調(diào)用Cublas的函數(shù)庫(kù)，來(lái)移植這個(gè)計(jì)算函數(shù)到GPU上。

MatrixmultiplicationvectorCUBLAS

Input： Matrix A[][]、 Vector V[]

Output： Vector Result[]

cublasInit();

cublasSetVector(COL*ROW, sizeof(float), A, 1, CUDA_A, 1);

cublasSetVector(ROW, sizeof(float), V, 1, CUDA_V, 1);

cublasSgemv_(′n′, COL, ROW, α, CUDA_A, COL, CUDA_V, 1, β, CUDA_Result,1);

cublasGetVector(COL, sizeof(float), CUDA_Result, 1, Result, 1);

cublasShutdown();

CUBLAS的函數(shù)庫(kù)已經(jīng)有一個(gè)Matrix multiplication vector的API函數(shù)提供給使用者來(lái)調(diào)用，并且已經(jīng)經(jīng)過(guò)了優(yōu)化，所以在解特征值的問(wèn)題中就直接套用了CUBLAS的函數(shù)。

第三部分是具有高度并行化的標(biāo)準(zhǔn)化部分，因?yàn)樵跇?biāo)準(zhǔn)化的過(guò)程中它只是先將整個(gè)向量進(jìn)行一個(gè)累加的動(dòng)作，然后再把向量里的每個(gè)元素除以累加的值，接下來(lái)我們將探討并行標(biāo)準(zhǔn)化。

2.2.3 標(biāo)準(zhǔn)化

標(biāo)準(zhǔn)化的過(guò)程就是將向量里的每個(gè)元素標(biāo)準(zhǔn)化后相加起來(lái)等于1，因此我們要先定義如下公式來(lái)看標(biāo)準(zhǔn)化是如何進(jìn)行的

(2)

NormalizeCPU

Input： Vector V[]、 DATA_SIZE

Output： Vector Result[]

For (i=0； i

SUM=V[i]*V[i];

End

SUM=sqrt(SUM);

For (i=0； i

V[i]=V[i]/SUM;

End

以上是標(biāo)準(zhǔn)化在CPU的偽碼，在Normalize有兩個(gè)部分來(lái)計(jì)算，一個(gè)是計(jì)算向量的元素平方加總，另一個(gè)是針對(duì)計(jì)算向量里每一個(gè)元素進(jìn)行除以加總值開(kāi)平方根，因此在這兩個(gè)部分都可以用GPU來(lái)并行運(yùn)算，接下來(lái)是在GPU上的偽碼。

NormalizeGPU

Input： Vector V[]、 DATA_SIZE

Output： Vector Result[]

const int tid=threadIdx.x;

const int bid=blockIdx.x;

for(i=bid*THREAD_NUM+tid;i

SUM+=V[i]*V[i];

End

SUM_Array[bid * THREAD_NUM+tid]=SUM;

for(int i=0; i

final_sum+=sum[i];

End

final_sum=sqrt(final_sum);

//Division

int i=blockIdx.x*blockDim.x+threadIdx.x;

Result [i]=V[i]/final_sum;

以上就是在GPU上實(shí)現(xiàn)在Krylove subspace算法中比較需要高運(yùn)算效能的計(jì)算函數(shù)，而這些計(jì)算函數(shù)實(shí)際在GPU上到底加快了多少，下面我們將探討其效能的改進(jìn)，以及并行Sparse的小結(jié)。

2.3 優(yōu)化結(jié)果

先來(lái)回顧一下我們的GPU并行Krylov subspace算法，來(lái)看我們?cè)谀莻€(gè)計(jì)算函數(shù)進(jìn)行GPU的并行運(yùn)算，接著，我們探討將以上3個(gè)函數(shù)放在GPU的計(jì)算平臺(tái)上后我們到底加快了多少，并且根據(jù)數(shù)據(jù)資料筆數(shù)的增加來(lái)考察它們之間的趨勢(shì)走勢(shì)。

KrylovAlgorithm

Input： Sparse matrixA, number of stepsm, and initial vectorv1of norm 1

Output： (Vm,Hm)

(1) Krylov subspace span(Using GPU—SpMV)

(2) Orthogonalizewwith respect toVj(Using CUBLAS)

u=VT×w

(3) SetHj+1,j，并將u加入V集合中(Using GPU—Normalization)

如圖6所示，可看出我們將3個(gè)運(yùn)算量負(fù)擔(dān)量大的函數(shù)放在GPU上執(zhí)行時(shí)間效率顯著上升，并且當(dāng)我們輸入的數(shù)據(jù)量越大時(shí)加速的結(jié)果就更為顯著。因此，可以看出GPU可以有效加速Krylov subspace。

圖6 時(shí)間倍數(shù)Krylov subspace

3 仿真實(shí)驗(yàn)

我們將Swiss roll和S curve兩組Data set作為Input，并根據(jù)選取的數(shù)據(jù)點(diǎn)數(shù)增加檢驗(yàn)用GPU并行LLE的加速效能。

如圖7所示，是根據(jù)不同的數(shù)據(jù)點(diǎn)數(shù)GPU LLE算法加速的情況，我們所設(shè)的參數(shù)K=9、維度(Dimension)=32，根據(jù)時(shí)間倍數(shù)來(lái)看，LLE算法在CPU與GPU不同的平臺(tái)上執(zhí)行效率可達(dá)十多倍的差距，如果數(shù)據(jù)點(diǎn)數(shù)再倍數(shù)增長(zhǎng)，之間的差異就越大。

圖7 LLE并行算法時(shí)間倍數(shù)關(guān)系

接下來(lái)我們依據(jù)不同的Data set和數(shù)據(jù)點(diǎn)數(shù)來(lái)看GPU LLE算法的輸出圖，如圖8所示。根據(jù)GPU并行LLE算法進(jìn)行運(yùn)算后的輸出，在LLE的計(jì)算過(guò)程中是以C++編碼來(lái)實(shí)現(xiàn)，但在繪圖的時(shí)候是將其輸出的Eigen pair導(dǎo)入MATLAB的plot function進(jìn)行繪制。

圖8 并行LLE算法輸出Swiss roll

接著，如圖9所示。是另一個(gè)Data set S curve的輸出圖，讓我們看N分別不同的時(shí)候它的Output的變化。

圖9 并行LLE算法輸出S curve

以上就是根據(jù)不同的Data set和不同的N來(lái)取樣所形成的輸出圖，并根據(jù)特征值的近似值(Approximation rate)是1e-7可知，兩個(gè)個(gè)計(jì)算算法相差在GPU的容許范圍之內(nèi)。

4 結(jié)束語(yǔ)

在傳統(tǒng)的LLE算法中選擇計(jì)算負(fù)擔(dān)較重的KNN算法和大型稀疏矩陣解特征值放在GPU上進(jìn)行并行運(yùn)算，除了利用GPU 的高并行度外還利用GPU對(duì)于浮點(diǎn)數(shù)的特殊處理能力進(jìn)行加速，并在實(shí)驗(yàn)中得到了驗(yàn)證，加速的結(jié)果高達(dá)十幾倍左右，并且輸出的特征值和LLE算法所得特征值的近似值(Approximation rate)也低到1e-7左右。因此，對(duì)于整體的GPU并行LLE算法得到了不錯(cuò)的效果。