基于GPU的LLE算法加速及性能優化

李 繁，嚴 星，張曉宇

(1.新疆財經大學 網絡與實驗教學中心，新疆 烏魯木齊 830012；2.新疆財經大學 信息管理學院，新疆 烏魯木齊 830012)

0 引 言

LLE數據降維算法是利用局部的信息計算相關權重，從而產生一個大型稀疏矩陣，通過解大型稀疏矩陣的特征值來解決非線性數據降維的問題。在LLE算法提出后又有一些改進的LLE算法提出：如HLLE算法結合了LLE和拉普拉斯特征映射來修改LLE的第二個步驟，并采用Hessain來評估測量權重，從而避免LLE算法存在的一些out-of-sample問題[1]。MLLE算法也是修改LLE的重建權重部分，通過引進對每一個鄰點計算多個獨立的線性權重來修改LLE的第二個步驟[2]，這個方法一樣是在改進LLE計算權重時的正確性或可能會發生的一些問題。所以上面兩個算法都是修改LLE的第二個步驟，來提高LLE算法的準確性。

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

1 KNN算法優化

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

1.1 基于GPU的Brute force KNN

Brute force KNN大致可分成3個部分：

(1)計算每一點的平面距離；

(2)利用一個排序算法對平面距離進行排序；

(3)挑出排序完成后的前K個數據點并回傳。

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

首先，先定義平面距離公式，假設我們有兩組Data set，一個是參照數據點R，另一組數據集合為搜索數據點Q。而平面距離公式是計算Q到R矩陣里每一行之間的距離

(1)

因此根據我們的平面距離公式來計算兩個數據集合內數據點的平面距離。由于LLE是要查找我們輸入的Data set內每一點的KNN，所以R和Q都是一樣的二維矩陣，并且Q矩陣內的每一行元素都要對R矩陣的每一行元素做平面距離計算，因此我們在這的時間復雜度是O(N2M)，M是數據維度。

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

由此可以看出數據點的距離計算只是多個循環迭代重復計算，因此在計算距離這個計算函數具有高度的可并行度，所以我們將這個計算函數放置到GPU的計算單元上，利用GPU高度并行計算的特性來大幅降低此計算函數所需要的執行時間，以加快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);

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

圖1 GPU計算架構

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

接下來KNN算法除了之前的距離計算公式上有很大的計算負擔之外，還有另一個需要大量計算效能的函數，就是利用一個排序算法來對平面距離進行排序，而排序有很多種方法，像插入排序(Insertion sort)、泡沫排序(Bubble sort)、基數排序(Radix sort)、快速排序(Quick sort)等等，因此我們必須要選擇一個適合用來并行的排序算法來進行KNN的排序。根據上述的排序算法，Quick sort這個方法在CPU的計算是快速排序的平均時間復雜度比較低O(NlogN)， 但因為它是用遞歸的方式來運作，所以沒辦法進行有效率的并行運算，因此我們將利用其它的排序方法來加快我們的排序。相比較來說，雖然插入算法在CPU的運算上具有比較低的運算效率O(N2)， 但因為Insertion sort是用iteration排序的方法并且計算之間是獨立的(Independent)，所以在并行度上來講插入算法是比較適合用GPU來進行運算加速的，所以我們在暴力KNN算法中用GPU加快了數據點的距離計算和排序函數，接下來我們要探討另外一種方法，也就是KD tree on GPU。

1.2 基于GPU的KD tree

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

圖2 GPU KD tree構造

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是用遞歸的方式來建構，在GPU上實現并行處理相對比較難，但在分割的時候我們還是可以用GPU來實現Radix sort的排序算法來降低檢索中位數的時間復雜度，然后在第(2)步會用一個循環進行KD tree分割的部分，所以在此步會用比較元素和中位數來建立左右List，這個部分我們也可以用GPU來分類，將整個Data set切成左右兩個Sub data set，所以在KD tree的建構上我們還是可以用GPU來找尋中位數和分割元素兩個步驟的加速運算，而接下來我們將探討如何利用GPU進行KD tree的 KNN搜索，如圖3所示。

圖3 GPU KNN搜索

KDTreeSearchGPUAlgorithm

Input： KD tree結構， 搜索點(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的搜索當中由于需要通過Recursive的方式進行搜索，所以后面的Search需要用到前面的數據進行搜索，所以我們在這個部分無法進行并行運算，因此我們只能將GPU的運算處理加在數據點之間的計算距離上，而在將要進行的KNN搜索中，我們要將此數據點當作最近的數據點時也要先查看在已搜索到的KNN當中是否已存在，如果不存在就將此點視為最短的數據點，否則將此數據點忽略掉。

1.3 KNN算法優化結果

我們將以圖表的方式來呈現兩個算法在GPU上執行效率相對于CPU的執行效率到底加快了多少。針對輸入數據隨維度增加及效果數據點數的增加判斷算法的加速效果，如圖4、表1所示。

圖4 按維度增加執行時間倍數關系

表1 基于維度的執行時間比較(數據點=2048，單位/s)

由圖4及表1可以看出，在維度不斷增加的過程中，在GPU與CPU不同的計算平臺上執行時間的距離差距就越來越遠，時間倍數甚至差到好幾十倍，因此在GPU的運算上KNN的算法大幅提升了效能。接下來我們將看如果增加數據點數的情況下效能會如何提升。

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

表2 基于數據點數的執行時間比較(維度=16，單位/s)

圖5 按數據點數增加執行時間倍數關系

2 稀疏矩陣特征值求解優化

在LLE算法中另一個需要耗費大量運算時間的就是解

大型稀疏矩陣的特征值，由于前一步計算權重時是根據KNN選出K個鄰點，同時將非鄰點之間的權重設為0，所以經過權重重建就會形成一個大型的稀疏矩陣，這樣就有效降低了運算量。但以大型矩陣要解特征值還是比較耗時間，因為在求解的過程中需要用到矩陣相乘以及矩陣與向量相乘[8]，這兩個部分都需要大量的運算效能，所以我們將GPU的運算平臺放在這兩個部分。接下來，我們將先探討如何配置內存可以使GPU并行計算發揮最大效果。

2.1 GPU內存分配

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

2.1.1 每行分配一個block

這種方法是一個很直觀的方法，將CSR的儲存方式每行當作一個block來計算，然后block內的每一個thread讀該行數據內的一個元素，并且計算乘積，再利用并行規約的方法將每個計算的結果加起來，一直到block計算完整行的元素，這種方法并行程度較高，并且相鄰的thread就會存取相鄰的全局內存。

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

2.1.2 每行分配一個warp

每行分配一個warp 可以順利解決block內thread空轉的問題。雖然這種方法可以解決block 內thread 空轉的問題，但到了行尾還是會有運算元素不足的問題，因為前面說過大型的稀疏矩陣不是每一行的元素都相同，所以會有些運算元素不轉的問題，我們將會利用下一個方法來改進相關問題。

2.1.3 對齊計算元素

前兩個方法不管是把每行看成是一個Block還是每行看成一個warp都會有運算元素不足的問題，因此，我們通過在行尾補0對齊的warp_size方法來解決這個問題。我們將行尾補0，這樣一來就不會發生non-coalesced元素運算。因此，我們將利用這個方法來儲存GPU的Sparse matrix。接著，我們來介紹SpMV的算法。

2.2 并行Krylov子空間

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

2.2.1 Span Krylove subspace

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

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函數中我們可以發現它有兩層循環，而且之間的元素相乘都是獨立運算的，所以在這方面很適合GPU來進行并行運算。因此，我們將考慮如何將SpMV的算法放在GPU的平臺上進行運算。

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的計算函數放在GPU的平臺上進行運算，而Thread和Block數據分割的方式如之前內存分配部分所述，利用GPU對內存存取的方式來將GPU的并行度提升至最高。

除了SpMV的大量計算需要用到GPU來進行并行運算之外，還有正交化的計算函數需要用到GPU來并行計算，接下來我們將介紹如何利用GPU進行正交化計算。

2.2.2 正交化

在這部分我們來討論在Arnoldi算法中除了SpMV那個部分花了很多時間計算之外，另一個花大量時間的地方就是在正交化處理的過程中，而在正交化處理所花費時間最多的地方就是矩陣向量相乘，在這個部分我們將利用CUDA內建的Cublas函數庫來進行處理加速，而以下就是在CPU的偽碼，以及用Cublas函數庫寫的程序。

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

接下來我們利用Cublas函數庫重新寫上方的偽碼，我們只簡單調用Cublas的函數庫，來移植這個計算函數到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的函數庫已經有一個Matrix multiplication vector的API函數提供給使用者來調用，并且已經經過了優化，所以在解特征值的問題中就直接套用了CUBLAS的函數。

第三部分是具有高度并行化的標準化部分，因為在標準化的過程中它只是先將整個向量進行一個累加的動作，然后再把向量里的每個元素除以累加的值，接下來我們將探討并行標準化。

2.2.3 標準化

標準化的過程就是將向量里的每個元素標準化后相加起來等于1，因此我們要先定義如下公式來看標準化是如何進行的

(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

以上是標準化在CPU的偽碼，在Normalize有兩個部分來計算，一個是計算向量的元素平方加總，另一個是針對計算向量里每一個元素進行除以加總值開平方根，因此在這兩個部分都可以用GPU來并行運算，接下來是在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上實現在Krylove subspace算法中比較需要高運算效能的計算函數，而這些計算函數實際在GPU上到底加快了多少，下面我們將探討其效能的改進，以及并行Sparse的小結。

2.3 優化結果

先來回顧一下我們的GPU并行Krylov subspace算法，來看我們在那個計算函數進行GPU的并行運算，接著，我們探討將以上3個函數放在GPU的計算平臺上后我們到底加快了多少，并且根據數據資料筆數的增加來考察它們之間的趨勢走勢。

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個運算量負擔量大的函數放在GPU上執行時間效率顯著上升，并且當我們輸入的數據量越大時加速的結果就更為顯著。因此，可以看出GPU可以有效加速Krylov subspace。

圖6 時間倍數Krylov subspace

3 仿真實驗

我們將Swiss roll和S curve兩組Data set作為Input，并根據選取的數據點數增加檢驗用GPU并行LLE的加速效能。

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

圖7 LLE并行算法時間倍數關系

接下來我們依據不同的Data set和數據點數來看GPU LLE算法的輸出圖，如圖8所示。根據GPU并行LLE算法進行運算后的輸出，在LLE的計算過程中是以C++編碼來實現，但在繪圖的時候是將其輸出的Eigen pair導入MATLAB的plot function進行繪制。

圖8 并行LLE算法輸出Swiss roll

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

圖9 并行LLE算法輸出S curve

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

4 結束語

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