機(jī)器學(xué)習(xí)優(yōu)化算法之爬山算法小結(jié)

 機(jī)器學(xué)習(xí)的項(xiàng)目,不可避免的需要補(bǔ)充一些優(yōu)化算法,對(duì)于優(yōu)化算法,爬山算法還是比較重要的.鑒于此,花了些時(shí)間仔細(xì)閱讀了些爬山算法的paper.基于這些,做一些總結(jié).
 目錄
　　1. 爬山算法簡(jiǎn)單描述
 　  2. 爬山算法的主要算法
        2.1 首選爬山算法
        2.2 最陡爬山算法
        2.3 隨機(jī)重新開(kāi)始爬山算法
        2.4 模擬退火算法(也是爬山算法)
      3. 實(shí)例求解
 正文
    爬山算法,是一種局部貪心的最優(yōu)算法. 該算法的主要思想是:每次拿相鄰點(diǎn)與當(dāng)前點(diǎn)進(jìn)行比對(duì),取兩者中較優(yōu)者,作為爬坡的下一步.
舉一個(gè)例子,求解下面表達(dá)式的最大值. 且假設(shè) x,y均按為0.1間隔遞增.
為了更好的描述,我們先使用pyhton畫(huà)出該函數(shù)的圖像:

圖像的python代碼:

 1 # encoding:utf8
 2 from matplotlib import pyplot as plt
 3 import numpy as np
 4 from mpl_toolkits.mplot3d import Axes3D
 5
 6
 7 def func(X, Y, x_move=0, y_move=0):
 8     def mul(X, Y, alis=1):
 9         return alis * np.exp(-(X * X + Y * Y))
10
11     return mul(X, Y) + mul(X - x_move, Y - y_move, 2)
12
13
14 def show(X, Y):
15     fig = plt.figure()
16     ax = Axes3D(fig)
17     X, Y = np.meshgrid(X, Y)
18     Z = func(X, Y, 1.7, 1.7)
19     plt.title("demo_hill_climbing")
20     ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', )
21     ax.set_xlabel('x label', color='r')
22     ax.set_ylabel('y label', color='g')
23     ax.set_zlabel('z label', color='b')
24     # 具體函數(shù)方法可用 help(function) 查看，如：help(ax.plot_surface)
25     # ax.scatter(X,Y,Z,c='r') #繪點(diǎn)
26     plt.show()
27
28 if __name__ == '__main__':
29     X = np.arange(-2, 4, 0.1)
30     Y = np.arange(-2, 4, 0.1)
31
32     show(X,Y)
View Code
     對(duì)于上面這個(gè)問(wèn)題,我們使用爬山算法該如何求解呢? 下面我們從爬山算法中的幾種方式分別求解一下這個(gè)小題.
  1. 首選爬山算法
　　依次尋找該點(diǎn)X的鄰近點(diǎn)中首次出現(xiàn)的比點(diǎn)X價(jià)值高的點(diǎn),并將該點(diǎn)作為爬山的點(diǎn)(此處說(shuō)的價(jià)值高,在該題中是指Z或f(x,y)值較大). 依次循環(huán),直至該點(diǎn)的鄰近點(diǎn)中不再有比其大的點(diǎn). 我們成為該點(diǎn)就是山的頂點(diǎn),又稱(chēng)為最優(yōu)點(diǎn).
     那么解題思路就有:
　　　　　1.  隨機(jī)選擇一個(gè)登山的起點(diǎn)S(x0,y0,z0),并以此為起點(diǎn)開(kāi)始登山.直至"登頂".
   下面是我們實(shí)現(xiàn)的代碼:

 1 # encoding:utf8
 2 from random import random, randint
 3
 4 from matplotlib import pyplot as plt
 5 import numpy as np
 6 from mpl_toolkits.mplot3d import Axes3D
 7
 8
 9 def func(X, Y, x_move=1.7, y_move=1.7):
10     def mul(X, Y, alis=1):
11         return alis * np.exp(-(X * X + Y * Y))
12
13     return mul(X, Y) + mul(X - x_move, Y - y_move, 2)
14
15
16 def show(X, Y, Z):
17     fig = plt.figure()
18     ax = Axes3D(fig)
19     plt.title("demo_hill_climbing")
20     ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', )
21     ax.set_xlabel('x label', color='r')
22     ax.set_ylabel('y label', color='g')
23     ax.set_zlabel('z label', color='b')
24     # ax.scatter(X,Y,Z,c='r') #繪點(diǎn)
25     plt.show()
26
27
28 def drawPaht(X, Y, Z,px,py,pz):
29     fig = plt.figure()
30     ax = Axes3D(fig)
31     plt.title("demo_hill_climbing")
32     ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', )
33     ax.set_xlabel('x label', color='r')
34     ax.set_ylabel('y label', color='g')
35     ax.set_zlabel('z label', color='b')
36     ax.plot(px,py,pz,'r.') #繪點(diǎn)
37     plt.show()
38
39
40 def hill_climb(X, Y):
41     global_X = []
42     global_Y = []
43
44     len_x = len(X)
45     len_y = len(Y)
46     # 隨機(jī)登山點(diǎn)
47     st_x = randint(0, len_x-1)
48     st_y = randint(0, len_y-1)
49
50     def argmax(stx, sty, alisx=0, alisy=0):
51         cur = func(X[0][st_x], Y[st_y][0])
52         next = func(X[0][st_x + alisx], Y[st_y + alisy][0])
53
54         return cur < next and True or False
55
56     while (len_x > st_x >= 0) or (len_y > st_y >= 0):
57         if st_x + 1 < len_x and argmax(st_x, st_y, 1):
58             st_x += 1
59         elif st_y + 1 < len_x and argmax(st_x, st_y, 0, 1):
60             st_y += 1
61         elif st_x >= 1 and argmax(st_x, st_y, -1):
62             st_x -= 1
63         elif st_y >= 1 and argmax(st_x, st_y, 0, -1):
64             st_y -= 1
65         else:
66             break
67         global_X.append(X[0][st_x])
68         global_Y.append(Y[st_y][0])
69     return global_X, global_Y, func(X[0][st_x], Y[st_y][0])
70
71
72 if __name__ == '__main__':
73     X = np.arange(-2, 4, 0.1)
74     Y = np.arange(-2, 4, 0.1)
75     X, Y = np.meshgrid(X, Y)
76     Z = func(X, Y, 1.7, 1.7)
77     px, py, maxhill = hill_climb(X, Y)
78     print px,py,maxhill
79     drawPaht(X, Y, Z,px,py,func(np.array(px), np.array(py), 1.7, 1.7))
View Code
對(duì)比幾次運(yùn)行的結(jié)果:

從上圖中,我們可以比較清楚的觀察到,首選爬山算法的缺陷.
2.那么最陡爬山算法呢?
   簡(jiǎn)單描述:
              最陡爬山算法是在首選爬山算法上的一種改良,它規(guī)定每次選取鄰近點(diǎn)價(jià)值最大的那個(gè)點(diǎn)作為爬上的點(diǎn).
   下面我們來(lái)實(shí)現(xiàn)一下它:

 1 # encoding:utf8
 2 from random import random, randint
 3
 4 from matplotlib import pyplot as plt
 5 import numpy as np
 6 from mpl_toolkits.mplot3d import Axes3D
 7
 8
 9 def func(X, Y, x_move=1.7, y_move=1.7):
10     def mul(X, Y, alis=1):
11         return alis * np.exp(-(X * X + Y * Y))
12
13     return mul(X, Y) + mul(X - x_move, Y - y_move, 2)
14
15
16 def show(X, Y, Z):
17     fig = plt.figure()
18     ax = Axes3D(fig)
19     plt.title("demo_hill_climbing")
20     ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', )
21     ax.set_xlabel('x label', color='r')
22     ax.set_ylabel('y label', color='g')
23     ax.set_zlabel('z label', color='b')
24     # ax.scatter(X,Y,Z,c='r') #繪點(diǎn)
25     plt.show()
26
27
28 def drawPaht(X, Y, Z, px, py, pz):
29     fig = plt.figure()
30     ax = Axes3D(fig)
31     plt.title("demo_hill_climbing")
32     ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', )
33     ax.set_xlabel('x label', color='r')
34     ax.set_ylabel('y label', color='g')
35     ax.set_zlabel('z label', color='b')
36     ax.plot(px, py, pz, 'r.')  # 繪點(diǎn)
37     plt.show()
38
39
40 def hill_climb(X, Y):
41     global_X = []
42     global_Y = []
43
44     len_x = len(X)
45     len_y = len(Y)
46     # 隨機(jī)登山點(diǎn)
47     st_x = randint(0, len_x - 1)
48     st_y = randint(0, len_y - 1)
49
50     def argmax(stx, sty, alisx, alisy):
51         cur = func(X[0][stx], Y[sty][0])
52         next = func(X[0][alisx], Y[alisy][0])
53         if cur < next:
54             return alisx, alisy
55         return stx, sty
56         #return cur < next and alisx, alisy or stx, sty
57
58     tmp_x = st_x
59     tmp_y = st_y
60     while (len_x > st_x >= 0) or (len_y > st_y >= 0):
61         if st_x + 1 < len_x:
62             tmp_x, tmp_y = argmax(tmp_x, tmp_y, (st_x + 1), st_y)
63
64         if st_x >= 1:
65             tmp_x, tmp_y = argmax(tmp_x, tmp_y, st_x - 1, st_y)
66
67         if st_y + 1 < len_x:
68             tmp_x, tmp_y = argmax(tmp_x, tmp_y, st_x, st_y + 1)
69
70         if st_y >= 1:
71             tmp_x, tmp_y = argmax(tmp_x, tmp_y, st_x, st_y - 1)
72
73         if tmp_x != st_x or tmp_y != st_y:
74             st_x = tmp_x
75             st_y = tmp_y
76         else:
77             break
78         global_X.append(X[0][st_x])
79         global_Y.append(Y[st_y][0])
80     return global_X, global_Y, func(X[0][st_x], Y[st_y][0])
81
82
83 if __name__ == '__main__':
84     X = np.arange(-2, 4, 0.1)
85     Y = np.arange(-2, 4, 0.1)
86     X, Y = np.meshgrid(X, Y)
87     Z = func(X, Y, 1.7, 1.7)
88     px, py, maxhill = hill_climb(X, Y)
89     print px, py, maxhill
90     drawPaht(X, Y, Z, px, py, func(np.array(px), np.array(py), 1.7, 1.7))
View Code

從這個(gè)結(jié)果來(lái)看,因?yàn)榉秶鷶U(kuò)大了一點(diǎn),所以效果會(huì)好一點(diǎn)點(diǎn),當(dāng)依舊是一個(gè)局部最優(yōu)算法.
3.隨機(jī)重新開(kāi)始爬山算法呢?
   簡(jiǎn)單的描述:
   　　　　隨機(jī)重新開(kāi)始爬山算法是基于最陡爬山算法,其實(shí)就是加一個(gè)達(dá)到全局最優(yōu)解的條件,如果滿足該條件,就結(jié)束運(yùn)算,反之則無(wú)限次重復(fù)運(yùn)算最陡爬山算法.
  由于此題,并沒(méi)有結(jié)束的特征條件,我們這里就不給予實(shí)現(xiàn).
4.模擬退火算法
   簡(jiǎn)單描述:
(1)隨機(jī)挑選一個(gè)單元k，并給它一個(gè)隨機(jī)的位移，求出系統(tǒng)因此而產(chǎn)生的能量變化ΔEk。
(2)若ΔEk?0，該位移可采納，而變化后的系統(tǒng)狀態(tài)可作為下次變化的起點(diǎn)；
若ΔEk>0，位移后的狀態(tài)可采納的概率為
式中T為溫度，然后從(0,1)區(qū)間均勻分布的隨機(jī)數(shù)中挑選一個(gè)數(shù)R，若R<Pk，則將變化后的狀態(tài)作為下次的起點(diǎn)；否則，將變化前的狀態(tài)作為下次的起點(diǎn)。 數(shù)據(jù)分析培訓(xùn)
(3)轉(zhuǎn)第(1)步繼續(xù)執(zhí)行，知道達(dá)到平衡狀態(tài)為止。
代碼實(shí)現(xiàn)為:

 1 # encoding:utf8
 2 from random import random, randint
 3
 4 from matplotlib import pyplot as plt
 5 import numpy as np
 6 from mpl_toolkits.mplot3d import Axes3D
 7
 8
 9 def func(X, Y, x_move=1.7, y_move=1.7):
10     def mul(X, Y, alis=1):
11         return alis * np.exp(-(X * X + Y * Y))
12
13     return mul(X, Y) + mul(X - x_move, Y - y_move, 2)
14
15
16 def show(X, Y, Z):
17     fig = plt.figure()
18     ax = Axes3D(fig)
19     plt.title("demo_hill_climbing")
20     ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', )
21     ax.set_xlabel('x label', color='r')
22     ax.set_ylabel('y label', color='g')
23     ax.set_zlabel('z label', color='b')
24     # ax.scatter(X,Y,Z,c='r') #繪點(diǎn)
25     plt.show()
26
27
28 def drawPaht(X, Y, Z, px, py, pz):
29     fig = plt.figure()
30     ax = Axes3D(fig)
31     plt.title("demo_hill_climbing")
32     ax.plot_surface(X, Y, Z, rstride=1, cstride=1, color='b' )
33     ax.set_xlabel('x label', color='r')
34     ax.set_ylabel('y label', color='g')
35     ax.set_zlabel('z label', color='b')
36     ax.plot(px, py, pz, 'r.')  # 繪點(diǎn)
37     plt.show()
38
39
40 def hill_climb(X, Y):
41     global_X = []
42     global_Y = []
43     # 初始溫度
44     temperature = 105.5
45     # 溫度下降的比率
46     delta = 0.98
47     # 溫度精確度
48     tmin = 1e-10
49
50     len_x = len(X)
51     len_y = len(Y)
52
53     # 隨機(jī)登山點(diǎn)
54     st_x = X[0][randint(0, len_x - 1)]
55     st_y = Y[randint(0, len_y - 1)][0]
56     st_z = func(st_x, st_y)
57
58     def argmax(stx, sty, alisx, alisy):
59         cur = func(st_x, st_y)
60         next = func(alisx, alisy)
61
62         return cur < next and True or False
63
64     while (temperature > tmin):
65         # 隨機(jī)產(chǎn)生一個(gè)新的鄰近點(diǎn)
66         # 說(shuō)明: 溫度越高幅度鄰近點(diǎn)跳躍的幅度越大
67         tmp_x = st_x + (random() * 2 - 1) * temperature
68         tmp_y = st_y + + (random() * 2 - 1) * temperature
69         if 4 > tmp_x >= -2 and 4 > tmp_y >= -2:
70             if argmax(st_x, st_y, tmp_x, tmp_y):
71                 st_x = tmp_x
72                 st_y = tmp_y
73             else:  # 有機(jī)會(huì)跳出局域最優(yōu)解
74                 pp = 1.0 / (1.0 + np.exp(-(func(tmp_x, tmp_y) - func(st_x, st_y)) / temperature))
75                 if random() < pp:
76                     st_x = tmp_x
77                     st_y = tmp_y
78         temperature *= delta  # 以一定的速率下降
79         global_X.append(st_x)
80         global_Y.append(st_y)
81     return global_X, global_Y, func(st_x, st_y)
82
83
84 if __name__ == '__main__':
85     X = np.arange(-2, 4, 0.1)
86     Y = np.arange(-2, 4, 0.1)
87     X, Y = np.meshgrid(X, Y)
88     Z = func(X, Y, 1.7, 1.7)
89     px, py, maxhill = hill_climb(X, Y)
90     print px, py, maxhill
91     drawPaht(X, Y, Z, px, py, func(np.array(px), np.array(py), 1.7, 1.7))
View Code
效果: