三種GPU并行的自適應鄰域模擬退火算法

林 敏，鐘一文

福建農林大學 計算機與信息學院，福州 350002

1 引言

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

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

2 自適應鄰域的模擬退火算法

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

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

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

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

3 基于GPU并行的自適應鄰域退火算法

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

3.1 GPU與CUDA

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

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

圖1 CUDA的kernel結構

3.2 解空間的存儲結構

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

圖2 解空間的存儲結構

運算中需要產生大量的隨機數，傳統的做法是在主機端先生成大量的隨機數，并將其復制到GPU的全局內存中，這種做法雖然簡單，但需要占用大量的GPU存儲空間。本文的作法是在GPU端動態的并行生成隨機數，既可以節省大量的存儲空間，又可以提高執行效率。

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

GPU的每個線程獨立運行一條馬爾可夫鏈[13]，直到滿足停止條件。這時在GPU端并行計算出每個BLOCK內的最優解，然后再從這些最優解并行計算出最優解，最后將最優解傳送給CPU端。以下是本算法的偽代碼表示：

GPU端產生了BLOCK_NUM*THREAD_NUM個線程，每個線程都執行相同的Markov_SA（）過程，每個Markov_SA（）過程相當于負責一條獨立的馬爾可夫鏈，Markov_SA（）過程的偽代碼如下：

m個線程將產生m個解決方案，最終要從這m個解決方案中找出最優解，一般的做法是由CPU端完成這個過程，但當數據量較大的時候，將花費比較多的時間。所以本算法采用歸約法在GPU端完成最優值的比較查找。SeleBest過程以塊作為單位，塊內每個線程負責兩個數的比較，每個塊的最優解保存到塊的第一個元素，然后將結果傳輸給CPU端，由CPU端在剩余的最優解中找出最優解。

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

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

算法的主過程描述如下：

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

3.5 基于BLOCK分塊的并行算法

GPU中的線程以BLOCK進行劃分，塊內線程可通過共享存儲器和同步執行協作。每個BLOCK內的所有線程計算出下一個解后，通過歸約法選出該BLOCK中的最優解，將其作為此BLOCK內所有線程的下一個解。其主過程與4.2節中類似，每個線程的執行過程表示如下：

算法中__shared__表示存儲類型是共享內存，共享內存是接近核心的內存，所以訪問速度很快。算法中在計算塊內最優解時，需要比較最優解，并將最優解所處的線程號記錄下來，以便計算最優解在全局內存中的對應位置。所以定義了mysharenode這個結構體類型，mysharenode包含min和tx兩個數據成員，min表示最優解，tx表示塊內線程號。

4CUDA程序的優化

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

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

4.1 全局內存的合并訪問

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

圖3 全局內存的合并訪問

本文設計的存儲結構一列為一個解向量，當計算函數值時，GPU可以通過二次合并訪問即可讀出所需的數據，并可通過大量的線程隱藏延時，如果數據是float型，則僅需一次合并訪問，但double型數據量更大，所以其吞吐量是相同的。

4.2 處理bank conflict

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

在計算最優解時，分別先計算每個block內的最優解，然后將數據從全局內存加載到共享內存，計算完最優解再將數據傳輸給全局內存。這里共享內存中的數據可能會產生bank conflict，為了避免這個問題，算法中專門設計了相應的存儲結構。

5 實驗與分析

算法中選用了11個基準函數，如表5所示，在GTX650TI BOOST，Intel G630平臺上進行了實驗，b為5，馬爾可夫鏈長k為3，BLOCK數為16，Thread數為256，其中迭代次數N為2 000，每一次迭代表示運行的時間為單位1，則2 000次迭代運行的時間為2 000。對每個函數分別運行30次，其結果如表1至4所示。

表1 nonu-SA算法

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

表3 基于BLOCK的GPU并行算法

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

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

表5 基準函數

為了進一步檢驗算法的收斂速度，分別在迭代次數為100，500，1 000，2 000，4 000的情況下進行了實驗，結果如表6至9所示。從實驗數據看并行算法表現出了較好的穩定性，并且隨著迭代次數的增加，多數函數的收斂速度大大優于nonu-SA算法。

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

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

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

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

表9 nonu-SA

6 結束語

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

表10 算法的運行時間比較

[1]張霖斌，姚振興，紀晨，等.快速模擬退火算法及應用[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分布狀態發生器的模擬退火算法[J].清華大學學報：自然科學版，2000，40（9）：109-112.

[12]張慶科，楊波.基于GPU的現代并行優化算法[J].計算機科學，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應用設計與開發：方法與最佳實踐[M].于玉龍，譯.北京：機械工業出版社，2013：94-114.