林 敏,鐘一文
福建農(nóng)林大學(xué) 計(jì)算機(jī)與信息學(xué)院,福州 350002
模擬退火算法是一種啟發(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ì)算量。近年來(lái),隨著GPU技術(shù)的發(fā)展,GPU并行計(jì)算的研究掀起了一個(gè)熱點(diǎn)。文中提出了新的基于GPU并行的自適應(yīng)鄰域的模擬退火算法,并借助了nonu-SA[10]算法的狀態(tài)生成函數(shù),結(jié)合GPU并行的特點(diǎn),提出了三種GPU并行的算法。算法通過(guò)產(chǎn)生大量的并行線程加快搜索進(jìn)程,并使鄰域的選擇可以從多個(gè)位置出發(fā),即每個(gè)位置被限制在局部范圍搜索下一個(gè)解,而多個(gè)位置的組合使其在尋找下一個(gè)解時(shí)可以有選擇地“長(zhǎng)跳”,不易陷入局部最優(yōu)解。
鄰域的分布一直是模擬退火算法的一個(gè)研究重點(diǎn),Gaussian分布較適于局部搜索,而產(chǎn)生大擾動(dòng)的概率幾乎為零,難以脫離平坦區(qū)和多峰區(qū)[11],所以Gaussian分布極易陷入局部最優(yōu)解。Cauchy分布在原點(diǎn)處的峰值比Gaussian分布小,而兩端長(zhǎng)扁平形狀趨于零的速度比Gaussian分布慢,因此更容易產(chǎn)生大步長(zhǎng)擾動(dòng),使得SA更容易跳出局部極小,但其產(chǎn)生小幅度擾動(dòng)的能力相對(duì)減弱,從而局部趨化性搜索的能力下降。Cauchy分布的特點(diǎn)使其維持在半局部的范圍內(nèi)進(jìn)行搜索,并時(shí)不時(shí)來(lái)一個(gè)“長(zhǎng)跳”,“長(zhǎng)跳”使其不易陷入局部極小值,但當(dāng)最小值與局部極小值的距離很近的時(shí)候,“長(zhǎng)跳”會(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)的距離,那么算法極易陷入局部最小值。
并行算法主要有幾種實(shí)現(xiàn)方式,一是基于多核CPU并行的,二是基于集群架構(gòu)的,三是基于GPU并行[12]的。基于CPU并行的方式由于CPU核心數(shù)量的限制,并行的規(guī)模非常有限;而基于集群的方式其成本較高,且數(shù)據(jù)的傳輸較慢;GPU具有大量的核心,可以同時(shí)運(yùn)行幾千幾萬(wàn)個(gè)線程,GPU擁有最高的浮點(diǎn)計(jì)算能力,甚至超越了最先進(jìn)的高性能計(jì)算(HPC)中心的設(shè)備。許多前沿研究才能申請(qǐng)?jiān)囉眠@種每秒萬(wàn)億次浮點(diǎn)計(jì)算的超級(jí)計(jì)算機(jī),如今由帶有若干GPGPU與快速冗余獨(dú)立磁盤陣列的工作站即可實(shí)現(xiàn)。本文提出的GPU并行的退火算法是根據(jù)GPU及CUDA框架的特點(diǎn)進(jìn)行設(shè)計(jì)的,并利用了nonu-SA生成鄰域的方法,使算法的精度和收斂性大大提高。
CUDA是一種基于新的并行編程模型和指令集架構(gòu)的通用計(jì)算架構(gòu),CUDA能夠利用Nvidia GPU的并行計(jì)算引擎比CPU更高效地解決許多復(fù)雜計(jì)算任務(wù),它包含一個(gè)讓開(kāi)發(fā)者能夠使用C作為高級(jí)編程語(yǔ)言的軟件環(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)
GPU及CUDA框架的特點(diǎn),決定了GPU并行程序計(jì)算中需要在GPU端產(chǎn)生大量的解,選擇合適的存儲(chǔ)結(jié)構(gòu)將大大減少GPU端訪問(wèn)存儲(chǔ)器的開(kāi)銷。假設(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表示對(duì)應(yīng)的函數(shù)值。
圖2 解空間的存儲(chǔ)結(jié)構(gòu)
運(yùn)算中需要產(chǎn)生大量的隨機(jī)數(shù),傳統(tǒng)的做法是在主機(jī)端先生成大量的隨機(jī)數(shù),并將其復(fù)制到GPU的全局內(nèi)存中,這種做法雖然簡(jiǎn)單,但需要占用大量的GPU存儲(chǔ)空間。本文的作法是在GPU端動(dòng)態(tài)的并行生成隨機(jī)數(shù),既可以節(jié)省大量的存儲(chǔ)空間,又可以提高執(zhí)行效率。
GPU的每個(gè)線程獨(dú)立運(yùn)行一條馬爾可夫鏈[13],直到滿足停止條件。這時(shí)在GPU端并行計(jì)算出每個(gè)BLOCK內(nèi)的最優(yōu)解,然后再?gòu)倪@些最優(yōu)解并行計(jì)算出最優(yōu)解,最后將最優(yōu)解傳送給CPU端。以下是本算法的偽代碼表示:
GPU端產(chǎn)生了BLOCK_NUM*THREAD_NUM個(gè)線程,每個(gè)線程都執(zhí)行相同的Markov_SA()過(guò)程,每個(gè)Markov_SA()過(guò)程相當(dāng)于負(fù)責(zé)一條獨(dú)立的馬爾可夫鏈,Markov_SA()過(guò)程的偽代碼如下:
m個(gè)線程將產(chǎn)生m個(gè)解決方案,最終要從這m個(gè)解決方案中找出最優(yōu)解,一般的做法是由CPU端完成這個(gè)過(guò)程,但當(dāng)數(shù)據(jù)量較大的時(shí)候,將花費(fèi)比較多的時(shí)間。所以本算法采用歸約法在GPU端完成最優(yōu)值的比較查找。SeleBest過(guò)程以塊作為單位,塊內(nèi)每個(gè)線程負(fù)責(zé)兩個(gè)數(shù)的比較,每個(gè)塊的最優(yōu)解保存到塊的第一個(gè)元素,然后將結(jié)果傳輸給CPU端,由CPU端在剩余的最優(yōu)解中找出最優(yōu)解。
遺傳算法[14]是一種基于自然選擇和群體遺傳機(jī)理的尋優(yōu)算法,它具有三個(gè)基本算子:選擇、交叉和變異。遺傳算法具有天生的并行性,但它極易陷入“早熟”,本算法利用遺傳算法的交叉算子與退火算法相結(jié)合,既可避免其過(guò)早進(jìn)入“早熟”,又可加快退火算法的搜索進(jìn)程。
算法的主過(guò)程描述如下:
算法Cross采用遺傳算法的交叉算子,按一定的交叉概率從所有解中選出部分解,對(duì)這些解隨機(jī)地進(jìn)行兩兩交叉。經(jīng)過(guò)交叉會(huì)產(chǎn)生兩個(gè)新解,加上原來(lái)的兩個(gè)解總共有四個(gè)解,在這四個(gè)解中選擇最好的一個(gè)解,并將其復(fù)制給另外一個(gè)解做為新解以進(jìn)行下一步的計(jì)算。其算法描述如下。交叉運(yùn)算的過(guò)程表示如下:
GPU中的線程以BLOCK進(jìn)行劃分,塊內(nèi)線程可通過(guò)共享存儲(chǔ)器和同步執(zhí)行協(xié)作。每個(gè)BLOCK內(nèi)的所有線程計(jì)算出下一個(gè)解后,通過(guò)歸約法選出該BLOCK中的最優(yōu)解,將其作為此BLOCK內(nèi)所有線程的下一個(gè)解。其主過(guò)程與4.2節(jié)中類似,每個(gè)線程的執(zhí)行過(guò)程表示如下:
算法中__shared__表示存儲(chǔ)類型是共享內(nèi)存,共享內(nèi)存是接近核心的內(nèi)存,所以訪問(wèn)速度很快。算法中在計(jì)算塊內(nèi)最優(yōu)解時(shí),需要比較最優(yōu)解,并將最優(yōu)解所處的線程號(hào)記錄下來(lái),以便計(jì)算最優(yōu)解在全局內(nèi)存中的對(duì)應(yīng)位置。所以定義了mysharenode這個(gè)結(jié)構(gòu)體類型,mysharenode包含min和tx兩個(gè)數(shù)據(jù)成員,min表示最優(yōu)解,tx表示塊內(nèi)線程號(hào)。
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)勢(shì),需要采用不同的訪問(wèn)方式和存儲(chǔ)結(jié)構(gòu)。本文中主要用到了全局存儲(chǔ)器,共享內(nèi)存,對(duì)于使用全局存儲(chǔ)器要盡量滿足合并內(nèi)存訪問(wèn)的要求,對(duì)于共享存儲(chǔ)器要盡量避免bank conflict。
訪問(wèn)一次全局存儲(chǔ)器需要幾百個(gè)時(shí)鐘周期,對(duì)全局內(nèi)存的訪問(wèn)是否滿足合并訪問(wèn)條件時(shí)對(duì)CUDA程序性能影響最明顯的因素之一。如果warp內(nèi)的32線程所有訪問(wèn)地址都位于一個(gè)128 Byte的段內(nèi),則只會(huì)進(jìn)行一次合并訪問(wèn),否則要分成多次訪問(wèn),如圖3所示。
圖3 全局內(nèi)存的合并訪問(wèn)
本文設(shè)計(jì)的存儲(chǔ)結(jié)構(gòu)一列為一個(gè)解向量,當(dāng)計(jì)算函數(shù)值時(shí),GPU可以通過(guò)二次合并訪問(wèn)即可讀出所需的數(shù)據(jù),并可通過(guò)大量的線程隱藏延時(shí),如果數(shù)據(jù)是float型,則僅需一次合并訪問(wèn),但double型數(shù)據(jù)量更大,所以其吞吐量是相同的。
共享存儲(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è)問(wèn)題,算法中專門設(shè)計(jì)了相應(yīng)的存儲(chǔ)結(jié)構(gòu)。
算法中選用了11個(gè)基準(zhǔn)函數(shù),如表5所示,在GTX650TI BOOST,Intel G630平臺(tái)上進(jìn)行了實(shí)驗(yàn),b為5,馬爾可夫鏈長(zhǎng)k為3,BLOCK數(shù)為16,Thread數(shù)為256,其中迭代次數(shù)N為2 000,每一次迭代表示運(yùn)行的時(shí)間為單位1,則2 000次迭代運(yùn)行的時(shí)間為2 000。對(duì)每個(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)一步提高該算法的性能?;贕PU的多路并行算法,其算法性能整體上不如其他兩種并行算法,但其設(shè)計(jì)非常簡(jiǎn)單占用空間小,大部分情況表現(xiàn)比nonu-SA出色,且隨著線程數(shù)的增加更加明顯。對(duì)于函數(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算法相對(duì)比較接近,但還是優(yōu)于nonu-SA算法幾倍,其他兩種并行算法大大優(yōu)于nonu-SA算法。
表7 GPU并行的遺傳模擬退火算法
表8 多條馬爾可夫鏈GPU并行的退火算法
表9 nonu-SA
根據(jù)GPU并行的特點(diǎn)和優(yōu)勢(shì),提出了三種GPU并行的自適應(yīng)鄰域的模擬退火算法,并通過(guò)11個(gè)基準(zhǔn)函數(shù)對(duì)這三種并行算法加以比較。這三種算法各有優(yōu)勢(shì),基于GPU多路并行算法實(shí)現(xiàn)簡(jiǎn)單,所需存儲(chǔ)空間最小,并且對(duì)個(gè)別函數(shù)效果比較好;基于GPU并行的遺傳-模擬退火算法結(jié)合了遺傳算法的特點(diǎn),而基于BLOCK分塊的GPU并行模擬退火算法,利用了CUDA中的BLOCK概念使線程以BLOCK為單位進(jìn)行共享,對(duì)大部分函數(shù)來(lái)說(shuō)此算法收斂速度最快。這三種算法具有較好的通用性,不僅僅可以應(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ì)與開(kāi)發(fā):方法與最佳實(shí)踐[M].于玉龍,譯.北京:機(jī)械工業(yè)出版社,2013:94-114.