模擬退火演算法

在實際日常中,人們會經常遇到如下問題:在某個給定的定義域X” role=”presentation”>XXX內,求函式f(x)” role=”presentation”>f(x)f(x)f(x)對應的最優值。此處以最小值問題舉例(最大值問題可以等價轉化成最小值問題),形式化為:

minx∈Xf(x).” role=”presentation”>minx∈Xf(x).minx∈Xf(x).

\min_{x \in X} f(x).
如果X” role=”presentation”>XXX是離散有限取值,那麼可以通過窮取法獲得問題的最優解;如果X” role=”presentation”>XXX連續,但f(x)” role=”presentation”>f(x)f(x)f(x)是凸的,那可以通過梯度下降等方法獲得最優解;如果X” role=”presentation”>XXX連續且f(x)” role=”presentation”>f(x)f(x)f(x)非凸,雖說根據已有的近似求解法能夠找到問題解,可解是否是最優的還有待考量,很多時候若初始值選擇的不好,非常容易陷入區域性最優值。

這裡寫圖片描述

隨著日常業務場景的複雜化,第三種問題經常遇見。如何有效地避免區域性最優的困擾?模擬退火演算法應運而生。其實模擬退火也算是啟發式演算法的一種,具體學習的是冶金學中金屬加熱-冷卻的過程。由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所發明的,V.Čern在1985年也獨立發明此演演算法。

不過模擬退火演算法到底是如何模擬金屬退火的原理?主要是將熱力學的理論套用到統計學上,將搜尋空間內每一點想像成空氣內的分子;分子的能量,就是它本身的動能;而搜尋空間內的每一點,也像空氣分子一樣帶有“能量”,以表示該點對命題的合適程度。演演算法先以搜尋空間內一個任意點作起始:每一步先選擇一個“鄰居”,然後再計算從現有位置到達“鄰居”的概率。若概率大於給定的閾值,則跳轉到“鄰居”;若概率較小,則停留在原位置不動。

一、模擬退火演算法基本思想

模擬退火其實也是一種貪心演算法,但是它的搜尋過程引入了隨機因素。在迭代更新可行解時,以一定的概率來接受一個比當前解要差的解,因此有可能會跳出這個區域性的最優解,達到全域性的最優解。以下圖為例,假定初始解為左邊藍色點A,模擬退火演算法會快速搜尋到區域性最優解B,但在搜尋到區域性最優解後,不是就此結束,而是會以一定的概率接受到左邊的移動。也許經過幾次這樣的不是區域性最優的移動後會到達全域性最優點D,於是就跳出了區域性最小值。

這裡寫圖片描述
根據熱力學的原理,在溫度為T” role=”presentation”>TTT時,出現能量差為dE” role=”presentation”>dEdEdE的降溫的概率為p(dE)” role=”presentation” style=”position: relative;”>p(dE)p(dE)p(dE),表示為:

p(dE)=exp⁡(dEkT).” role=”presentation”>p(dE)=exp(dEkT).p(dE)=exp⁡(dEkT).

p(dE) = \exp(\frac{dE}{kT}).
其中k” role=”presentation” style=”position: relative;”>kkk是波爾茲曼常數,值為k=1.3806488(13)×10−23″ role=”presentation” style=”position: relative;”>k=1.3806488(13)×10−23k=1.3806488(13)×10−23k= 1.3806488(13) × 10^{−23},exp” role=”presentation” style=”position: relative;”>expexp\exp表示自然指數,且dE&lt;0″ role=”presentation” style=”position: relative;”>dE<0dE<0dE < 0。因此dE/kT&lt;0″ role=”presentation” style=”position: relative;”>dE/kT<0dE/kT<0dE/kT< 0,所以p(dE)” role=”presentation” style=”position: relative;”>p(dE)p(dE)p(dE)函式的取值範圍是(0,1)。滿足概率密度函式的定義。其實這條公式更直觀意思就是:溫度越高,出現一次能量差為p(dE)” role=”presentation” style=”position: relative;”>p(dE)p(dE)p(dE)的降溫的概率就越大;溫度越低,則出現降溫的概率就越小。

在實際問題中,這裡的“一定的概率”的計算參考了金屬冶煉的退火過程。假定當前可行解為x” role=”presentation” style=”position: relative;”>xxx,迭代更新後的解為x_new” role=”presentation” style=”position: relative;”>x_newx_newx\_new,那麼對應的“能量差”定義為:

Δf=f(x_new)−f(x).” role=”presentation”>Δf=f(x_new)−f(x).Δf=f(x_new)−f(x).

Δf = f(x\_new) − f(x).
其對應的“一定概率”為:

p(Δf)={exp⁡(−ΔfkT)最小值优化问题exp⁡(ΔfkT)最大值优化问题” role=”presentation”>p(Δf)=⎧⎩⎨exp(−ΔfkT)exp(ΔfkT)最小值優化問題最大值優化問題p(Δf)={exp⁡(−ΔfkT)最小值優化問題exp⁡(ΔfkT)最大值優化問題

p(Δf) = \left\{ \begin{array}{ll}
\exp(-\frac{Δf}{kT}) & 最小值優化問題 \\
\exp(\frac{Δf}{kT}) & 最大值優化問題
\end{array} \right.
注:在實際問題中,可以設定k=1″ role=”presentation” style=”position: relative;”>k=1k=1k=1。因為kT” role=”presentation” style=”position: relative;”>kTkTkT可以等價於一個引數 T” role=”presentation” style=”position: relative;”>TTT。如設定k=2″ role=”presentation” style=”position: relative;”>k=2k=2k=2、T=1000″ role=”presentation” style=”position: relative;”>T=1000T=1000T=1000,等於直接設定T=2000″ role=”presentation” style=”position: relative;”>T=2000T=2000T=2000的效果。

二、模擬退火演算法描述

  1. 初始化:初始溫度T” role=”presentation” style=”position: relative;”>TTT(充分大),溫度下限Tmin” role=”presentation” style=”position: relative;”>TminTminT_{min}(充分小),初始解狀態x” role=”presentation” style=”position: relative;”>xxx(是演算法迭代的起點),每個T” role=”presentation” style=”position: relative;”>TTT值的迭代次數L” role=”presentation” style=”position: relative;”>LLL;
  2. 對l=1,2,…,L” role=”presentation” style=”position: relative;”>l=1,2,…,Ll=1,2,…,Ll= 1,2,…,L做第3至第6步;
  3. 產生新解x_new” role=”presentation” style=”position: relative;”>x_newx_newx\_new: (x_new=x Δx” role=”presentation” style=”position: relative;”>x_new=x Δxx_new=x Δxx\_new = x Δx);
  4. 利計算增量Δf=f(x_new)−f(x)” role=”presentation” style=”position: relative;”>Δf=f(x_new)−f(x)Δf=f(x_new)−f(x)Δf = f(x\_new) − f(x),其中f(x)” role=”presentation” style=”position: relative;”>f(x)f(x)f(x)為優化目標;
  5. 若Δf&lt;0″ role=”presentation” style=”position: relative;”>Δf<0Δf<0Δf<0(若尋找最大值,Δf&gt;0″ role=”presentation” style=”position: relative;”>Δf>0Δf>0Δf>0)則接受x_new” role=”presentation” style=”position: relative;”>x_newx_newx\_new作為新的當前解,否則以概率exp⁡(−Δf/(kT))” role=”presentation” style=”position: relative;”>exp(−Δf/(kT))exp⁡(−Δf/(kT))\exp(-Δf/(kT ))接受x_new” role=”presentation” style=”position: relative;”>x_newx_newx\_new作為新的當前解;
  6. 如果滿足終止條件則輸出當前解作為最優解,結束程式。(終止條件通常取為連續若干個新解都沒有被接受時終止演算法。);
  7. T” role=”presentation” style=”position: relative;”>TTT逐漸減少,且T&gt;Tmin” role=”presentation” style=”position: relative;”>T>TminT>TminT> T_{min},然後轉第2步。

這裡寫圖片描述

三、模擬退火演算法的優缺點

模擬退火演算法的應用很廣泛,可以高效地求解NP完全問題,如貨郎擔問題(Travelling Salesman Problem,簡記為TSP)、最大截問題(Max Cut Problem)、0-1揹包問題(Zero One Knapsack Problem)、圖著色問題(Graph Colouring Problem)等等,但其引數難以控制,不能保證一次就收斂到最優值,一般需要多次嘗試才能獲得(大部分情況下還是會陷入區域性最優值)。觀察模擬退火演算法的過程,發現其主要存在如下三個引數問題:

  (1) 溫度T的初始值設定問題
  溫度T” role=”presentation” style=”position: relative;”>TTT的初始值設定是影響模擬退火演算法全域性搜尋效能的重要因素之一、初始溫度高,則搜尋到全域性最優解的可能性大,但因此要花費大量的計算時間;反之,則可節約計算時間,但全域性搜尋效能可能受到影響。

  (2) 退火速度問題,即每個T” role=”presentation” style=”position: relative;”>TTT值的迭代次數
  模擬退火演算法的全域性搜尋效能也與退火速度密切相關。一般來說,同一溫度下的“充分”搜尋是相當必要的,但這也需要計算時間。迴圈次數增加必定帶來計算開銷的增大。

  (3) 溫度管理問題
  溫度管理問題也是模擬退火演算法難以處理的問題之一。實際應用中,由於必須考慮計算複雜度的切實可行性等問題,常採用如下所示的降溫方式:

T=α×T.α∈(0,1).” role=”presentation”>T=α×T.α∈(0,1).T=α×T.α∈(0,1).

T=\alpha \times T.\alpha \in (0,1).
注:為了保證較大的搜尋空間,α” role=”presentation” style=”position: relative;”>αα\alpha一般取接近於1的值,如0.95、0.9。

四、模擬退火演算法Python實戰

經過上面理論知識的薰陶,相信大家已經對模擬退火演算法有了較深入的理解,接下來通過實戰再強化一下大家的認識,此處利用模擬退火演算法求解如下優化問題:

minf(x)=(x−2)∗(x 3)∗(x 8)∗(x−9)” role=”presentation”>minf(x)=(x−2)∗(x 3)∗(x 8)∗(x−9)minf(x)=(x−2)∗(x 3)∗(x 8)∗(x−9)

\min f(x)=(x-2)*(x 3)*(x 8)*(x-9)

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
def inputfun(x):
return (x-2)*(x 3)*(x 8)*(x-9)
initT = 1000 #初始溫度
minT = 1 #溫度下限
iterL = 1000 #每個T值的迭代次數
delta = 0.95 #溫度衰減係數
k = 1
initx = 10*(2*np.random.rand()-1)
nowt = initT
print "初始解:",initx
xx = np.linspace(-10,10,300)
yy = inputfun(xx)
plt.figure()
plt.plot(xx,yy)
plt.plot(initx,inputfun(initx),'o')
#模擬退火演算法尋找最小值過程
while nowt>minT:
for i in np.arange(1,iterL,1):
funVal = inputfun(initx)
xnew = initx (2*np.random.rand()-1)
if xnew>=-10 and xnew<=10:
funnew = inputfun(xnew)
res = funnew-funVal
if res<0:
initx = xnew
else:
p = np.exp(-(res)/(k*nowt))
if np.random.rand()<p:
initx = xnew
#            print initx-xnew
#    print initx
#    print nowt
nowt = nowt*delta
print "最優解:",initx
print "最優值:",inputfun(initx)
plt.plot(initx,inputfun(initx),'*r')
plt.show()

這裡寫圖片描述
這裡寫圖片描述
可以看到,即使初始解選取的有風險,模擬退火演算法經過迭代,也可以成功跳出區域性最優值,最終收斂到問題的全域性最優。

參考資料

  1. http://wiki.mbalib.com/wiki/模擬退火演算法 模擬退火演算法