三種GPU并行的自適應(yīng)鄰域模擬退火算法

林 敏，鐘一文

福建農(nóng)林大學(xué) 計(jì)算機(jī)與信息學(xué)院，福州 350002

1 引言

模擬退火算法是一種啟發(fā)式的全局最優(yōu)化算法，但其有兩個(gè)缺陷：（1）容易陷入局部最優(yōu)解；（2）需要花費(fèi)較多的時(shí)間以求得一個(gè)較優(yōu)的解[1]。為克服這些缺陷，學(xué)者們提出了各種方法以加速退火算法的收斂速率。這些算法主要有：依賴溫度的快速退火方案，如Guassion分布退火算法，Cauthy分布的退火算法[2]，ASA[3]；與其他智能算法相結(jié)合的退火方案，如遺傳算法[4]，粒子群算法[5]，神經(jīng)網(wǎng)絡(luò)[6]；并行的退火算法，如 PSA[7-9]。

這些算法在不同方面提高了退火算法的效率，但在實(shí)際應(yīng)用中仍需較大的計(jì)算量。近年來，隨著GPU技術(shù)的發(fā)展，GPU并行計(jì)算的研究掀起了一個(gè)熱點(diǎn)。文中提出了新的基于GPU并行的自適應(yīng)鄰域的模擬退火算法，并借助了nonu-SA[10]算法的狀態(tài)生成函數(shù)，結(jié)合GPU并行的特點(diǎn)，提出了三種GPU并行的算法。算法通過產(chǎn)生大量的并行線程加快搜索進(jìn)程，并使鄰域的選擇可以從多個(gè)位置出發(fā)，即每個(gè)位置被限制在局部范圍搜索下一個(gè)解，而多個(gè)位置的組合使其在尋找下一個(gè)解時(shí)可以有選擇地“長跳”，不易陷入局部最優(yōu)解。

2 自適應(yīng)鄰域的模擬退火算法

鄰域的分布一直是模擬退火算法的一個(gè)研究重點(diǎn)，Gaussian分布較適于局部搜索，而產(chǎn)生大擾動(dòng)的概率幾乎為零，難以脫離平坦區(qū)和多峰區(qū)[11]，所以Gaussian分布極易陷入局部最優(yōu)解。Cauchy分布在原點(diǎn)處的峰值比Gaussian分布小，而兩端長扁平形狀趨于零的速度比Gaussian分布慢，因此更容易產(chǎn)生大步長擾動(dòng)，使得SA更容易跳出局部極小，但其產(chǎn)生小幅度擾動(dòng)的能力相對減弱，從而局部趨化性搜索的能力下降。Cauchy分布的特點(diǎn)使其維持在半局部的范圍內(nèi)進(jìn)行搜索，并時(shí)不時(shí)來一個(gè)“長跳”，“長跳”使其不易陷入局部極小值，但當(dāng)最小值與局部極小值的距離很近的時(shí)候，“長跳”會(huì)使其容易遠(yuǎn)離極小值。

文獻(xiàn)[10]提出了nonu-SA算法（自適應(yīng)鄰域的退火算法），其性能優(yōu)于Guassion的鄰域半徑隨著迭代次數(shù)的增加，而呈現(xiàn)遞減狀態(tài)，當(dāng)溫度很高的時(shí)候鄰域半徑可覆蓋整個(gè)解空間，隨著迭代次數(shù)的增加搜索范圍將隨之收窄到很小的范圍。假設(shè)當(dāng)前解X={x1,x2,…,xk,…,xm}，元素xk隨機(jī)擾動(dòng)得到xk'和新解X′={x1,x2,…,xk',…,xm}，nonu-SA采用如下擾動(dòng)公式：

由遞減函數(shù)的性質(zhì)可知鄰域半徑隨迭代次數(shù)n的增加，而逐漸縮窄，不像Guassion分布和Cauthy分布鄰域依賴于溫度。

nonu-SA算法收斂速度與迭代次數(shù)n有關(guān)，在數(shù)據(jù)量大的時(shí)候需要較多的迭代次數(shù)，這樣降低了其收斂速度，且隨著迭代次數(shù)的增加，鄰域半徑收斂到很小的范圍，這時(shí)如果全局最小值與局部最小值相隔較遠(yuǎn)的距離，那么算法極易陷入局部最小值。

3 基于GPU并行的自適應(yīng)鄰域退火算法

并行算法主要有幾種實(shí)現(xiàn)方式，一是基于多核CPU并行的，二是基于集群架構(gòu)的，三是基于GPU并行[12]的。基于CPU并行的方式由于CPU核心數(shù)量的限制，并行的規(guī)模非常有限；而基于集群的方式其成本較高，且數(shù)據(jù)的傳輸較慢；GPU具有大量的核心，可以同時(shí)運(yùn)行幾千幾萬個(gè)線程，GPU擁有最高的浮點(diǎn)計(jì)算能力，甚至超越了最先進(jìn)的高性能計(jì)算（HPC）中心的設(shè)備。許多前沿研究才能申請?jiān)囉眠@種每秒萬億次浮點(diǎn)計(jì)算的超級(jí)計(jì)算機(jī)，如今由帶有若干GPGPU與快速冗余獨(dú)立磁盤陣列的工作站即可實(shí)現(xiàn)。本文提出的GPU并行的退火算法是根據(jù)GPU及CUDA框架的特點(diǎn)進(jìn)行設(shè)計(jì)的，并利用了nonu-SA生成鄰域的方法，使算法的精度和收斂性大大提高。

3.1 GPU與CUDA

CUDA是一種基于新的并行編程模型和指令集架構(gòu)的通用計(jì)算架構(gòu)，CUDA能夠利用Nvidia GPU的并行計(jì)算引擎比CPU更高效地解決許多復(fù)雜計(jì)算任務(wù)，它包含一個(gè)讓開發(fā)者能夠使用C作為高級(jí)編程語言的軟件環(huán)境。

運(yùn)行在GPU上的CUDA并行計(jì)算函數(shù)稱為kernel（內(nèi)核函數(shù)），內(nèi)核函數(shù)并不是一個(gè)完整的程序，而是整個(gè)CUDA程序中的一個(gè)可以被并行執(zhí)行的步驟。每一個(gè)內(nèi)核會(huì)被N個(gè)CUDA線程執(zhí)行N次。每個(gè)kernel函數(shù)以線程網(wǎng)格（Grid）的形式組成，每個(gè)Grid有許多同樣大小的線程塊Block組成，每個(gè)線程塊又由若干個(gè)線程（Thread）組成，如圖1所示。例如BLOCK為16，thread為256，則總的線程數(shù)為16×256，在進(jìn)行CUDA程序設(shè)計(jì)時(shí)，需要考慮合理的BLOCK數(shù)和塊內(nèi)線程數(shù)Thread。

圖1 CUDA的kernel結(jié)構(gòu)

3.2 解空間的存儲(chǔ)結(jié)構(gòu)

GPU及CUDA框架的特點(diǎn)，決定了GPU并行程序計(jì)算中需要在GPU端產(chǎn)生大量的解，選擇合適的存儲(chǔ)結(jié)構(gòu)將大大減少GPU端訪問存儲(chǔ)器的開銷。假設(shè)要求解的函數(shù)的其中一個(gè)解表示為Xi=(x1,x2,…,xm)，而并行計(jì)算中有n個(gè)線程則有X1,X2,X3,…,Xi,…,Xn個(gè)解，如圖2所示，每一列表示一個(gè)解，之所以選擇這種結(jié)構(gòu)將在4.1節(jié)中解釋。圖2中所示的結(jié)構(gòu)可以用二維數(shù)組或一維數(shù)組存儲(chǔ)，本文中使用一維數(shù)組存儲(chǔ)。同時(shí)定義變量X_Valuei表示每個(gè)解Xi的函數(shù)值，即X_Valuei=f(Xi)，用一維數(shù)組X_Value[n]存儲(chǔ)所有f(Xi)的值。為了區(qū)分GPU端存儲(chǔ)空間，下文中將用GPU_Xi表示存儲(chǔ)在GPU端的解向量，GPU_X_Valuei表示對應(yīng)的函數(shù)值。

圖2 解空間的存儲(chǔ)結(jié)構(gòu)

運(yùn)算中需要產(chǎn)生大量的隨機(jī)數(shù)，傳統(tǒng)的做法是在主機(jī)端先生成大量的隨機(jī)數(shù)，并將其復(fù)制到GPU的全局內(nèi)存中，這種做法雖然簡單，但需要占用大量的GPU存儲(chǔ)空間。本文的作法是在GPU端動(dòng)態(tài)的并行生成隨機(jī)數(shù)，既可以節(jié)省大量的存儲(chǔ)空間，又可以提高執(zhí)行效率。

3.3 基于多條馬爾可夫鏈的GPU并行算法

GPU的每個(gè)線程獨(dú)立運(yùn)行一條馬爾可夫鏈[13]，直到滿足停止條件。這時(shí)在GPU端并行計(jì)算出每個(gè)BLOCK內(nèi)的最優(yōu)解，然后再從這些最優(yōu)解并行計(jì)算出最優(yōu)解，最后將最優(yōu)解傳送給CPU端。以下是本算法的偽代碼表示：

GPU端產(chǎn)生了BLOCK_NUM*THREAD_NUM個(gè)線程，每個(gè)線程都執(zhí)行相同的Markov_SA（）過程，每個(gè)Markov_SA（）過程相當(dāng)于負(fù)責(zé)一條獨(dú)立的馬爾可夫鏈，Markov_SA（）過程的偽代碼如下：

m個(gè)線程將產(chǎn)生m個(gè)解決方案，最終要從這m個(gè)解決方案中找出最優(yōu)解，一般的做法是由CPU端完成這個(gè)過程，但當(dāng)數(shù)據(jù)量較大的時(shí)候，將花費(fèi)比較多的時(shí)間。所以本算法采用歸約法在GPU端完成最優(yōu)值的比較查找。SeleBest過程以塊作為單位，塊內(nèi)每個(gè)線程負(fù)責(zé)兩個(gè)數(shù)的比較，每個(gè)塊的最優(yōu)解保存到塊的第一個(gè)元素，然后將結(jié)果傳輸給CPU端，由CPU端在剩余的最優(yōu)解中找出最優(yōu)解。

3.4 基于GPU并行的遺傳-模擬退火算法

遺傳算法[14]是一種基于自然選擇和群體遺傳機(jī)理的尋優(yōu)算法，它具有三個(gè)基本算子：選擇、交叉和變異。遺傳算法具有天生的并行性，但它極易陷入“早熟”，本算法利用遺傳算法的交叉算子與退火算法相結(jié)合，既可避免其過早進(jìn)入“早熟”，又可加快退火算法的搜索進(jìn)程。

算法的主過程描述如下：

算法Cross采用遺傳算法的交叉算子，按一定的交叉概率從所有解中選出部分解，對這些解隨機(jī)地進(jìn)行兩兩交叉。經(jīng)過交叉會(huì)產(chǎn)生兩個(gè)新解，加上原來的兩個(gè)解總共有四個(gè)解，在這四個(gè)解中選擇最好的一個(gè)解，并將其復(fù)制給另外一個(gè)解做為新解以進(jìn)行下一步的計(jì)算。其算法描述如下。交叉運(yùn)算的過程表示如下：

3.5 基于BLOCK分塊的并行算法

GPU中的線程以BLOCK進(jìn)行劃分，塊內(nèi)線程可通過共享存儲(chǔ)器和同步執(zhí)行協(xié)作。每個(gè)BLOCK內(nèi)的所有線程計(jì)算出下一個(gè)解后，通過歸約法選出該BLOCK中的最優(yōu)解，將其作為此BLOCK內(nèi)所有線程的下一個(gè)解。其主過程與4.2節(jié)中類似，每個(gè)線程的執(zhí)行過程表示如下：

算法中__shared__表示存儲(chǔ)類型是共享內(nèi)存，共享內(nèi)存是接近核心的內(nèi)存，所以訪問速度很快。算法中在計(jì)算塊內(nèi)最優(yōu)解時(shí)，需要比較最優(yōu)解，并將最優(yōu)解所處的線程號(hào)記錄下來，以便計(jì)算最優(yōu)解在全局內(nèi)存中的對應(yīng)位置。所以定義了mysharenode這個(gè)結(jié)構(gòu)體類型，mysharenode包含min和tx兩個(gè)數(shù)據(jù)成員，min表示最優(yōu)解，tx表示塊內(nèi)線程號(hào)。

4CUDA程序的優(yōu)化

GPU并行計(jì)算，要在多種存儲(chǔ)器之間進(jìn)行數(shù)據(jù)傳輸，不同存儲(chǔ)器的使用方法和傳輸速率有很大差異，在理想的情況下，GPU端運(yùn)行足夠多的線程，在所有存儲(chǔ)器傳輸進(jìn)行的時(shí)候，GPU的各個(gè)核心始終在計(jì)算著，這樣能很好地隱藏各種訪存延遲。

GPU中的存儲(chǔ)結(jié)構(gòu)主要有全局內(nèi)存、共享內(nèi)存、常量存儲(chǔ)器、寄存器、紋理存儲(chǔ)器[15]，不同的存儲(chǔ)器有各自的優(yōu)勢，需要采用不同的訪問方式和存儲(chǔ)結(jié)構(gòu)。本文中主要用到了全局存儲(chǔ)器，共享內(nèi)存，對于使用全局存儲(chǔ)器要盡量滿足合并內(nèi)存訪問的要求，對于共享存儲(chǔ)器要盡量避免bank conflict。

4.1 全局內(nèi)存的合并訪問

訪問一次全局存儲(chǔ)器需要幾百個(gè)時(shí)鐘周期，對全局內(nèi)存的訪問是否滿足合并訪問條件時(shí)對CUDA程序性能影響最明顯的因素之一。如果warp內(nèi)的32線程所有訪問地址都位于一個(gè)128 Byte的段內(nèi)，則只會(huì)進(jìn)行一次合并訪問，否則要分成多次訪問，如圖3所示。

圖3 全局內(nèi)存的合并訪問

本文設(shè)計(jì)的存儲(chǔ)結(jié)構(gòu)一列為一個(gè)解向量，當(dāng)計(jì)算函數(shù)值時(shí)，GPU可以通過二次合并訪問即可讀出所需的數(shù)據(jù)，并可通過大量的線程隱藏延時(shí)，如果數(shù)據(jù)是float型，則僅需一次合并訪問，但double型數(shù)據(jù)量更大，所以其吞吐量是相同的。

4.2 處理bank conflict

共享存儲(chǔ)器位于GPU內(nèi)，速度比全局內(nèi)存快很多，在不發(fā)生bank conflict的情況下，共享內(nèi)存的延遲幾乎只有局部存儲(chǔ)器和共享存儲(chǔ)器的1/100。

在計(jì)算最優(yōu)解時(shí)，分別先計(jì)算每個(gè)block內(nèi)的最優(yōu)解，然后將數(shù)據(jù)從全局內(nèi)存加載到共享內(nèi)存，計(jì)算完最優(yōu)解再將數(shù)據(jù)傳輸給全局內(nèi)存。這里共享內(nèi)存中的數(shù)據(jù)可能會(huì)產(chǎn)生bank conflict，為了避免這個(gè)問題，算法中專門設(shè)計(jì)了相應(yīng)的存儲(chǔ)結(jié)構(gòu)。

5 實(shí)驗(yàn)與分析

算法中選用了11個(gè)基準(zhǔn)函數(shù)，如表5所示，在GTX650TI BOOST，Intel G630平臺(tái)上進(jìn)行了實(shí)驗(yàn)，b為5，馬爾可夫鏈長k為3，BLOCK數(shù)為16，Thread數(shù)為256，其中迭代次數(shù)N為2 000，每一次迭代表示運(yùn)行的時(shí)間為單位1，則2 000次迭代運(yùn)行的時(shí)間為2 000。對每個(gè)函數(shù)分別運(yùn)行30次，其結(jié)果如表1至4所示。

表1 nonu-SA算法

表2 基于GPU并行的遺傳-模擬退火算法

表3 基于BLOCK的GPU并行算法

表4 基于多條馬爾可夫鏈的GPU并行算法

從實(shí)驗(yàn)可以看到基于BLOCK的GPU并行算法收斂速度最快，基于GPU并行的遺傳-模擬退火算法在某些方面優(yōu)于基于BLOCK的GPU并行算法，如函數(shù)4和函數(shù)7。這是由于遺傳算法的特性，使其更適合于這兩個(gè)函數(shù)，并且當(dāng)調(diào)高交叉率PC將進(jìn)一步提高該算法的性能。基于GPU的多路并行算法，其算法性能整體上不如其他兩種并行算法，但其設(shè)計(jì)非常簡單占用空間小，大部分情況表現(xiàn)比nonu-SA出色，且隨著線程數(shù)的增加更加明顯。對于函數(shù)9結(jié)果比較特殊，幾種算法的結(jié)果都比較接近，甚至nonu-SA結(jié)果更好。

表5 基準(zhǔn)函數(shù)

為了進(jìn)一步檢驗(yàn)算法的收斂速度，分別在迭代次數(shù)為100，500，1 000，2 000，4 000的情況下進(jìn)行了實(shí)驗(yàn)，結(jié)果如表6至9所示。從實(shí)驗(yàn)數(shù)據(jù)看并行算法表現(xiàn)出了較好的穩(wěn)定性，并且隨著迭代次數(shù)的增加，多數(shù)函數(shù)的收斂速度大大優(yōu)于nonu-SA算法。

表6 基于BLOCK的GPU并行模擬退火算法

為了在運(yùn)行時(shí)間上做進(jìn)一步的論證，將三種并行算法在所取得結(jié)果與nonu-SA相近或者更優(yōu)的情況下所耗費(fèi)的機(jī)器時(shí)間進(jìn)行了比較，得到表10的結(jié)果。從結(jié)果也可以更加證實(shí)三種并行算法的收斂速度優(yōu)于nonu-SA算法，其中3.3節(jié)中算法與nonu-SA算法相對比較接近，但還是優(yōu)于nonu-SA算法幾倍，其他兩種并行算法大大優(yōu)于nonu-SA算法。

表7 GPU并行的遺傳模擬退火算法

表8 多條馬爾可夫鏈GPU并行的退火算法

表9 nonu-SA

6 結(jié)束語

根據(jù)GPU并行的特點(diǎn)和優(yōu)勢，提出了三種GPU并行的自適應(yīng)鄰域的模擬退火算法，并通過11個(gè)基準(zhǔn)函數(shù)對這三種并行算法加以比較。這三種算法各有優(yōu)勢，基于GPU多路并行算法實(shí)現(xiàn)簡單，所需存儲(chǔ)空間最小，并且對個(gè)別函數(shù)效果比較好；基于GPU并行的遺傳-模擬退火算法結(jié)合了遺傳算法的特點(diǎn)，而基于BLOCK分塊的GPU并行模擬退火算法，利用了CUDA中的BLOCK概念使線程以BLOCK為單位進(jìn)行共享，對大部分函數(shù)來說此算法收斂速度最快。這三種算法具有較好的通用性，不僅僅可以應(yīng)用于本文中的鄰域結(jié)構(gòu)，還可以將其移植到采用其他鄰域生成函數(shù)的模擬退火算法中。

表10 算法的運(yùn)行時(shí)間比較

[1]張霖斌，姚振興，紀(jì)晨，等.快速模擬退火算法及應(yīng)用[J].石油地球物理勘探，1997，32（5）：654-660.

[2]Szu H，Hartley R.Fast simulated annealing[J].Physics letters A，1987，122（3）：157-162.

[3]Ingber L，Petraglia A，Petraglia M R，et al.Adaptive simulated annealing[M]//Stochastic Global Optimization and its Applications with Fuzzy Adaptive Simulated Annealing.Berlin Heidelberg：Springer，2012：33-62.

[4]Cordón O，Moya F，Zarco C.A new evolutionary algorithm combining simulated annealing and genetic programming for relevance feedback in fuzzy information retrieval systems[J].Soft Computing，2002，6（5）：308-319.

[5]Soleymani S.Bidding strategy of generation companies using PSO combined with SA method in the pay as bid markets[J].InternationalJournalofElectricalPower&Energy Systems，2011，33（7）：1272-1278.

[6]Salcedo-Sanz S，Yao X.A hybrid Hopfield network-genetic algorithm approach for the terminal assignment problem[J].IEEE Transactions on Systems，Man，and Cybernetics，Part B：Cybernetics，2004，34（6）：2343-2353.

[7]Boissin N，Lutton J L.A parallel simulated annealing algorithm[J].Parallel Computing，1993，19（8）：859-872.

[8]Ram D J，Sreenivas T H，Subramaniam K G.Parallel simulated annealing algorithms[J].Journal of Parallel and Distributed Computing，1996，37（2）：207-212.

[9]Aarts E H L，Korst J H M.Boltzmann machines as a model for parallel annealing[J].Algorithmica，1991，6（1-6）：437-465.

[10]Xinchao Z.Simulated annealing algorithm with adaptive neighborhood[J].Applied Soft Computing，2011，11（2）：1827-1836.

[11]王凌，鄭大鐘.基于Cauchy和Gaussian分布狀態(tài)發(fā)生器的模擬退火算法[J].清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2000，40（9）：109-112.

[12]張慶科，楊波.基于GPU的現(xiàn)代并行優(yōu)化算法[J].計(jì)算機(jī)科學(xué)，2012，39（4）：305-310.

[13]Hastings W K.Monte Carlo sampling methods using Markov chains and their applications[J].Biometrika，1970，57（1）：97-109.

[14]Mahfoud S W，Goldberg D E.Parallel recombinative simulated annealing：a genetic algorithm[J].Parallel Computing，1995，21（1）：1-28.

[15]Farber R.高性能CUDA應(yīng)用設(shè)計(jì)與開發(fā)：方法與最佳實(shí)踐[M].于玉龍，譯.北京：機(jī)械工業(yè)出版社，2013：94-114.