天天看點
傳回

sopt:一個簡單的python最優化庫引言

引言

    最近有些朋友總來問我有關遺傳算法的東西，我是在大學搞數學模組化的時候接觸過一些最優化和進化算法方面的東西，以前也寫過幾篇部落格記錄過,比如

遺傳算法的C語言實作(一):以非線性函數求極值為例

和

C語言實作粒子群算法(PSO)一

等，如果對原理有興趣的話可以去我的部落格具體檢視:

Lyrichu's Blog

。是以突發奇想，幹脆把以前寫的一些進化算法比如遺傳算法(GA),粒子群算法(PSO),模拟退火算法(SA)以及最近看的基于梯度的一些優化算法比如Gradient Descent,SGD,Momentum等整理一下，寫成一個python庫，友善那些有需要的朋友使用。斷斷續續花了幾天的時間，初步完成了一些基本功能，這裡簡單介紹一下。

sopt 簡介

    sopt是simple optimization的簡稱，目前我已經将代碼托管到pypi,位址是

sopt

,可以直接通過pip指令下載下傳安裝使用,由于項目隻依賴numpy,是以在windows和linux環境下安裝都很友善，直接

pip install sopt

即可。項目的github位址是

。目前sopt包含的優化方法如下:

  • 遺傳算法(Genetic Algorithm,GA)
  • 粒子群算法(Particle Swarm Optimization,PSO)
  • 模拟退火算法(Simulated Anealing,SA)
  • 随機遊走算法(Random Walk):
  • 梯度下降法(Gradient Descent,GD)
  • 動量優化算法(Momentum)
  • 自适應梯度算法(AdaGrad)
  • RMSProp

  • Adam

    由于隻是一個初步的版本，後續如果有時間的話，會加上更多的優化算法進去。目前所有的優化算法暫時隻支援連續實函數的優化;除了基于梯度的幾個優化算法,GA、PSO、SA以及Random Walk都同時支援無限制優化，線性限制優化以及非線性限制優化。具體我會在下面詳細說明。

    sopt 使用詳解以及執行個體示範

    1.SGA 使用

        SGA是Simple Genetic Algorithm的簡稱，是最基本的遺傳算法。其編碼方式采用二進制編碼,選擇方法采用輪盤賭法,交叉方法采用單點交叉,變異方式采用均勻變異,預設是求函數的最小值。下面是sopt中SGA的一個簡單使用執行個體: 
    from sopt.SGA import SGA
from math import sin
def func1(x):
return (x[0]-1)**2 + (sin(x[1])-0.5)**4 + 2
if __name__ == '__main__':
sga = SGA.SGA(func = func1,func_type = 'min',variables_num = 2,
    lower_bound = 0,upper_bound = 2,generations = 20,
    binary_code_length = 10)
# run SGA
sga.run()
# show the SGA optimization result in figure
sga.save_plot()
# print the result
sga.show_result()
    運作結果如下: 
    -------------------- SGA config is: --------------------
lower_bound:[0, 0]
generations:20
cross_rate:0.7
variables_num:2
mutation_rate:0.1
func_type:min
upper_bound:[2, 2]
population_size:100
func:<function func1 at 0x7f3d2311b158>
binary_code_length:10
-------------------- SGA caculation result is: --------------------
global best generation index/total generations:3/20
global best point:[1.00488759 0.45356794]
global best target:2.00003849823336

    用圖像展示為圖1所示:

    圖1 SGA 運作結果

    上面定義的目标函數為\(f(x_1,x_2)=(x_1-1)^2+(sin(x_2)-0.5)^4+2\),其中\(0<x_1,x_2<2\),函數的最小值點為(1,0.5236),最小值為2,與尋找到的最小值點(1.00488759 0.45356794)以及最小值2.00003849823336是非常接近的。

    SGA類的全部參數定義如下:

  • lower_bound:一個int或者float的數字，或者是長度為

    variables_num

     的ndarray,表示每個變量的下界,如果是一個數字的話，我們認為所有的下界都是一樣的(必填);
  • upper_bound:一個int或者float的數字，或者是長度為

    variables_num

     的ndarray,表示每個變量的上界,如果是一個數字的話，我們認為所有的上界都是一樣的(必填);
  • variables_num:表示目标函數變量的個數(必填);
  • func(必填):表示目标函數;
  • cross_rate:GA的交叉率,預設為0.7;
  • mutation_rate:變異率,預設為0.1;
  • population_size:種群數目,預設為100;
  • generations:進化代數,預設是200;
  • binary_code_length:GA二進制編碼的長度,預設是10;
  • func_type:函數優化類型,求最小值取值為'min',求最大值取值為'max',預設是'min';

    2. GA 使用

        GA是相比SGA更加一般的遺傳算法實作，其對于編碼方式，選擇方式，交叉方式以變異方式等會有更多的支援。首先還是看一個例子: 
    from sopt.GA.GA import GA
from sopt.util.functions import *
from sopt.util.ga_config import *
from sopt.util.constraints import *
class TestGA:
def __init__(self):
    self.func = quadratic11
    self.func_type = quadratic11_func_type
    self.variables_num = quadratic11_variables_num
    self.lower_bound = quadratic11_lower_bound
    self.upper_bound = quadratic11_upper_bound
    self.cross_rate = 0.8
    self.mutation_rate = 0.05
    self.generations = 200
    self.population_size = 100
    self.binary_code_length = 20
    self.cross_rate_exp = 1
    self.mutation_rate_exp = 1
    self.code_type = code_type.binary
    self.cross_code = False
    self.select_method = select_method.proportion
    self.rank_select_probs = None
    self.tournament_num = 2
    self.cross_method = cross_method.uniform
    self.arithmetic_cross_alpha = 0.1
    self.arithmetic_cross_exp = 1
    self.mutation_method = mutation_method.uniform
    self.none_uniform_mutation_rate = 1
    #self.complex_constraints = [constraints1,constraints2,constraints3]
    self.complex_constraints = None
    self.complex_constraints_method = complex_constraints_method.penalty
    self.complex_constraints_C = 1e6
    self.M = 1e8
    self.GA = GA(**self.__dict__)
def test(self):
    start_time = time()
    self.GA.run()
    print("GA costs %.4f seconds!" % (time()-start_time))
    self.GA.save_plot()
    self.GA.show_result()
if __name__ == '__main__':
TestGA().test()
    上面代碼的運作結果為: 
    GA costs 6.8320 seconds!
-------------------- GA config is: --------------------
func:<function quadratic11 at 0x7f998927bd08>
code_type:binary
complex_constraints:None
global_generations_step:200
cross_method:uniform
mutation_method:uniform
cross_rate:0.8
lower_bound:[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
tournament_num:2
variables_num:11
complex_constraints_method:penalty
none_uniform_mutation_rate:1
population_size:100
mutation_rate:0.05
generations:200
arithmetic_cross_alpha:0.1
func_type:min
mutation_rate_exp:1
cross_rate_exp:1
arithmetic_cross_exp:1
M:100000000.0
select_method:proportion
upper_bound:[11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11]
cross_code:False
binary_code_length:20
complex_constraints_C:1000000.0
-------------------- GA caculation result is: --------------------
global best target generation index/total generations:149/200
global best point:[ 1.07431037  2.41401426  2.4807906   4.36291634  4.90653029  6.13753427
  6.58147963  7.7370479   9.42957347 10.46616122 10.87151134]
global best target:2.2685208668743204
        其中

    sopt.util.functions

     預定義了一些測試函數,

    quadratic11

     是一個11元變量的二次函數，函數原型為:\(quadratic11(x_1,...,x_{11})=(x_1-1)^2 +(x_2-2)^2 + ... +(x_{11}-11)^2 + 1\),其中\(1 \le x_1,...,x_{11} \le 11\),函數的最小值點在(1,2,...,11)處取得，最小值為1。另外還定義了其他幾個測試函數為:
  • quadratic50:和quadratic11類似,隻是變量個數變為50,取值範圍變為1-50;
  • quadratic100:和quadratic11類似,隻是變量個數變為100,取值範圍變為1-100;
  • quadratic500:和quadratic11類似,隻是變量個數變為500,取值範圍變為1-500;
  • quadratic1000:和quadratic1000類似,隻是變量個數變為1000,取值範圍變為1-1000;
  • Rosenbrock:函數原型為:\(Rosenbrock(x_1,x_2)=100(x_1^2-x_2)^2+(1-x_1)^2\),其中\(-2.048 \le x_1,x_2 \le 2.048\),這個函數有很多局部極小值點，函數在(-2.048,-2.048)處取得最大值;
  • shubert2函數:函數定義為,\(shubert2(x_1,x_2)=\prod_{i=1}^{2}(\sum_{j=1}^{5}(jcos((j+1)*x_i+j))\),其中\(-10 \le x_1,x_2 \le 10\),函數有760個局部極小值點,全局極小值為-186.731,這裡為了将适應度函數變換為正，我們在原函數的基礎上加上1000,這樣函數的全局極小值點就變為813.269;

  • shubert函數，将shubert2函數中的兩個變量拓展為n個變量(\(\prod_{i=1}^{n}\))就可以得到一般的shubert函數了。

        對于每一個優化方法，我們都預定了一個形如

    xx_config

     的子產品，其中定義了該優化方法的一些常用預設參數，比如

    ga_config

     中就定義了一些GA的一些常用優化參數，

    ga_config.basic_config

     定義了一些基礎參數設定，比如

    basic_config.generations

     是一個預設進化代數,

    basic_config.mutation_rate

     是預設的變異參數；而

    ga_config.cross_method

     則預定義了所有支援的交叉方法,比如

    cross_method.uniform

     表示均勻交叉,

    cross_method.one_point

     表示單點交叉等;

    ga_config.mutation_method

    等也是類似的。有了這些預定義變量，可以免去我們手動輸入很多參數取值以及傳入方法字元串的麻煩(有時候可能會寫錯)。

        觀察上面的運作結果，我們發現對于quadratic11函數，最終找到的全局極小值為2.2685208668743204,和真實的全局極小值點1已經比較接近了；運作耗時6秒多，似乎有些長，這是因為我們種群數目設為100，進化代數設為200,比較大，而且又是采用二進制20位編碼，再加上python腳本語言的運作效率問題，是以稍慢。為了做一個對比,這裡我們将目标函數從quadratic11變為Rosenbrock函數，其他參數設定保持不變，得到的結果如下:

    GA costs 1.7245 seconds!
-------------------- GA config is: --------------------
population_size:100
lower_bound:[-2.048, -2.048]
mutation_rate_exp:1
select_method:proportion
code_type:binary
global_generations_step:200
generations:200
mutation_method:uniform
complex_constraints_method:penalty
binary_code_length:20
cross_method:uniform
arithmetic_cross_alpha:0.1
func:<function Rosenbrock at 0x7fe5fd538a60>
upper_bound:[2.048, 2.048]
cross_code:False
mutation_rate:0.05
cross_rate_exp:1
complex_constraints_C:1000000.0
cross_rate:0.8
variables_num:2
M:100000000.0
complex_constraints:None
none_uniform_mutation_rate:1
tournament_num:2
arithmetic_cross_exp:1
func_type:max
-------------------- GA caculation result is: --------------------
global best target generation index/total generations:75/200
global best point:[-2.04776953 -2.04537109]
global best target:3901.4655271502425

    圖2是每代最優值的計算結果:

    圖2 GA Rosenbrock 函數運作200代尋優結果

    替換成Rosenbrock函數之後，函數的運作時間從6秒多減少到1秒多(函數變量個數明顯減少了),這說明GA的運作時間與目标函數的變量個數是顯著相關的。最後找到的全局極大值點為(-2.04776953,-2.04537109)和真實全局極大值點(-2.048,-2.048)已經很接近了。

        SGA類是GA類的父類，是以GA類和SGA類有一些公共的屬性，比如

    generations

     ,

    population_size

    func_type

     等。和SGA一樣的參數這裡就不再列舉了，如下是GA特有的一些參數:
  • cross_rate_exp:cross_rate指數遞增的參數取值,即變異率按照\(r_t=r_0 \beta^t\),這裡\(r_t\)表示t次疊代的交叉率,\(r_0\)是初始的交叉率，\(\beta\)就是這裡的

    cross_rate_exp

     ,預設取值為1,一般設定為一個比1稍大的數字，比如1.0001等;
  • mutation_rate_exp:mutation_rate指數遞增參數取值,具體含義和

    cross_rate_exp

     類似;
  • code_type:采用什麼樣的編碼方式,有三種取值:'binary'表示二進制編碼(預設),'gray'表示采用格雷編碼,'real'表示采用實數編碼。其中格雷編碼是一種改進的二進制編碼，其優點在于進行交叉、變異等操作時，改變了染色體基因的某幾個位置，二進制編碼對應的實數值可能會發生非常大的變化，而格雷編碼一般不會。二進制編碼與格雷編碼的轉換可以參考附錄;
  • cross_code:這是一個布爾值,表示是否對二進制編碼或格雷碼采用交叉編碼的方式。具體說明參考附錄;
  • select_method:選擇操作方法，有如下的6種取值:
  1. 'proportion':表示比例選擇(輪盤賭法)
  2. 'keep_best':表示保留最優個體
  3. 'determinate_sampling':表示确定性采樣
  4. 'rssr':remainder stochastic sampling with replacement的縮寫，表示無放回餘數随機選擇
  5. 'rank':表示排序選擇
  6. 'stochastic_tournament':表示随機錦标賽法(以上6種選擇方式詳細說明請參考附錄)
  • rank_select_probs:僅當select_method取值為'rank'時，該參數才起作用，表示采用排序選擇時，适應度從低到高每個個體被選擇的機率，預設為None,采用\(p_i=\frac{i}{\sum_{j=1}^{n}j},i=1,2,...,n\)的方式計算，其中\(n\)表示種群個數,\(p_i\)表示第\(i\)個個體被選擇的機率，所有機率之和為1，如果取值不為None,你應該傳入一個size為

    population_size

     的ndarray,所有元素按照遞增排序，數組和為1;
  • tournament_num:如果select_method取值為'stochastic_tournament',該值表示錦标賽競争者數量;
  • cross_method:交叉方法,有如下4種取值:
  1. 'one_point'：表示單點交叉
  2. 'two_point':表示雙點交叉
  3. 'uniform':表示均勻交叉
  4. 'arithmetic':表示算術交叉,隻有當編碼采用實數編碼時，才能使用算術交叉
  • arithmetic_cross_alpha:算術交叉的系數，算術交叉的公式為:\(x_{new1}=\alpha x_1 + (1-\alpha)x_2,x_{new2}=\alpha x_2 + (1-\alpha)x_1\),這裡的\(\alpha\)即表示

    arithmetic_cross_alpha

     ,其預設值為0.1;
  • 'arithmetic_cross_exp':算術交叉系數按照指數遞增，即\(\alpha_t = \alpha r^t\),這裡\(\alpha_t\)表示第t代的算術交叉系數,\(\alpha\)是初始交叉系數,\(r\)即是這裡的

    arithmetic_cross_exp

     ,預設取值為1，一般取一個比1稍大的數，比如1.0001,t是進化代數;
  • 'mutation_method':變異方法，有如下5種取值:
  1. simple:即簡單變異，随機選擇一個染色體上的基因，按照變異機率進行變異
  2. uniform:即均勻變異,每個染色體上的每個基因都按照變異機率進行變異
  3. boundary:邊界變異,每個基因進行變異以後的取值隻能去邊界值,比如說各以0.5的機率取得其上界和下界，這種變異方式僅适用于實數編碼的情況，一般在目标函數的最優點靠近邊界時使用
  4. none_uniform:非均勻變異，我們不是取均勻分布的随機數去替換原來的基因，而是在原來基因的基礎上做一點微小的随機擾動，擾動以後的結果作為新的基因值，具體來說，我們采用的是這樣的變異方式:1)if random(0,1) = 0,\(x_k^{'} = x_k + \bigtriangleup (t,U_{max}^{k}-x_k)\);2)if random(0,1) = 1,\(x_k{'} = x_k - \bigtriangleup (t,x_k -U_{min}{k})\) 其中 \(\bigtriangleup(t,y)=y(1-r^{(1-t/T)b})\),\(r\)是0-1均勻分布的一個随機數，\(T\)是最大進化代數，\(b\)是一個系統參數，表示随機擾動對于進化代數的依賴長度，預設值為1
  5. gaussian:高斯變異,即使用高斯分布去替代均勻分布來進行變異,由高斯分布的特點可以知道，這種變異方式也是在原個體區域附近的某個局部區域進行重點搜尋，這裡高斯分布的均值\(\mu\)定義為\(\mu = \frac{U_{min}^{k}+U_{max}^{k}}{2}\),标準差\(\sigma=\frac{U_{max}^{k}-U_{min}^{k}}{6}\)
  • none_uniform_mutation_rate:即均勻分布

    none_uniform

     中定義的系統參數\(b\);

  • complex_constraints:複雜限制,預設值為None,即沒有複雜限制,隻有簡單的邊界限制,如果有複雜限制，其取值應該是\(\[func_1,func_2,...,func_n\]\),

    其中\(func_i\)表示第\(i\)個複雜限制函數名，比如對于一個複雜限制函數\(func_1:x_1^2+x_2^2 < 3\),其複雜限制函數應該這樣定義:

    def func1(x):
x1 = x[0]
x2 = x[1]
return x1**2 + x2**2 - 3
  • complex_constraints_method:複雜限制求解的方法,預設是

    penalty

     即懲罰函數法，暫時不支援其他的求解方式;
  • complex_constraints_C:采用

    penalty

    求解複雜限制的系數\(C\),比如對于某一個限制\(x_1^2+x_2^2 < 3\),GA在求解過程中，違反了該限制,即解滿足\(x_1^2+x_2^2 \ge 3\),那麼我們對目标函數增加一個懲罰項: \(C(x_1^2+x_2^2-3)\),\(C\)一般取一個很大的正數，預設值為\(10^6\)。

        GA類的參數很多可以調節，是以相比之下使用起來更加麻煩，特别是對于複雜函數的尋優，預設參數未必是最好的，可能需要你根據目标函數的特點自行嘗試，選擇最優的參數組合。

    3. PSO 使用

        PSO算法的原理可以參考我之前的部落格 C語言實作粒子群算法（PSO）二 ,這裡不再贅述。還是先看一個執行個體代碼: 
    from time import time
from sopt.util.functions import *
from sopt.util.pso_config import *
from sopt.PSO.PSO import PSO
from sopt.util.constraints import *
class TestPSO:
def __init__(self):
    self.func = quadratic11
    self.func_type = quadratic11_func_type
    self.variables_num = quadratic11_variables_num
    self.lower_bound = quadratic11_lower_bound
    self.upper_bound = quadratic11_upper_bound
    self.c1 = basic_config.c1
    self.c2 = basic_config.c2
    self.generations = 200
    self.population_size = 100
    self.vmax = 1
    self.vmin = -1
    self.w = 1
    self.w_start = 0.9
    self.w_end = 0.4
    self.w_method = pso_w_method.linear_decrease
    #self.complex_constraints = [constraints1,constraints2,constraints3]
    self.complex_constraints = None
    self.complex_constraints_method = complex_constraints_method.loop
    self.PSO = PSO(**self.__dict__)
def test(self):
    start_time = time()
    self.PSO.run()
    print("PSO costs %.4f seconds!" %(time()-start_time))
    self.PSO.save_plot()
    self.PSO.show_result()
if __name__ == '__main__':
TestPSO().test()
    運作結果為: 
    PSO costs 1.1731 seconds!
-------------------- PSO config is: --------------------
complex_constraints_method:loop
c1:1.49445
lower_bound:[1 1 1 1 1 1 1 1 1 1 1]
w_end:0.4
w_method:linear_decrease
complex_constraints:None
func:<function quadratic11 at 0x7f1ddb81a510>
upper_bound:[11 11 11 11 11 11 11 11 11 11 11]
generations:200
func_type:min
w:1
c2:1.49445
w_start:0.9
population_size:100
vmin:-1
vmax:1
variables_num:11
-------------------- PSO calculation result is: --------------------
global best generation index/total generations: 198/200
global best point: [ 1.          1.99999999  2.99999999  4.          5.          6.
  7.00000001  7.99999999  9.00000001 10.00000001 11.        ]
global best target: 1.0

    圖3 PSO 求解 quadratic11 200代運作結果

    上面的代碼意圖應該是非常明顯的,目标函數是\(quadratic11\),最終求得的最小值點幾乎就是全局最小值點。下面是PSO類中所有參數的具體定義:

  • variables_num:變量個數(必填);

  • variables_num
  • upper_bound:一個int或者float的數字,或者是長度為

    variables_num
  • func:目标函數名(必填);
  • c1:PSO 參數c1,預設值為1.49445;
  • c2:PSO 參數c2,預設值為1.49445;
  • generations:進化代數,預設值為100;
  • population_size:種群數量,預設為50;
  • vmax:粒子最大速度,預設值為1;
  • vmin:粒子最小速度,預設值為-1;
  • w:粒子慣性權重,預設為1;
  • w_start:如果不是使用恒定的慣性權重話,比如使用權重遞減政策,w_start表示初始權重;
  • w_end:相應的,w_end表示末尾權重;
  • w_method:權重遞減的方式,有如下5種取值:
  1. 'constant':權重恒定
  2. 'linear_decrease':線性遞減,即\(w_t = w_{end} + (w_{start} - w_{end})\frac{(T-t)}{T}\),其中\(T\)表示最大進化代數,\(w_t\)表示第\(t\)代權重
  3. 'square1_decrease':第一種平方遞減,即:\(w_t = w_{start} - (w_{start}-w_{end})(\frac{t}{T})^2\)
  4. 'square2_decrease':第二種平方遞減,即：\(w_t = w_{start} - (w_{start}-w_{end})(\frac{2t}{T}-(\frac{t}{T})^2)\)
  5. 'exp_decrease':指數遞減,即:\(w_t = w_{end}(\frac{w_{start}}{w_{end}})^{\frac{1}{1+\frac{10t}{T}}}\)
  • complex_constraints:複雜限制，預設為None,具體含義參考GA類的complex_constraints
  • complex_constraints_method:求解複雜限制的方法，預設為'loop',即如果解不滿足複雜限制，則再次随機産生解，直到滿足限制,暫時不支援其他的求解方式。

    4. SA 使用

        SA,即模拟退火算法，是基于機率的一種尋優方法，具體原理可以參考我的部落格 模拟退火算法（SA）求解TSP 問題（C語言實作） 。sopt中SA的使用執行個體如下: 
    from time import time
from sopt.util.functions import *
from sopt.Optimizers.SA import SA
from sopt.util.sa_config import *
from sopt.util.constraints import *
class TestSA:
def __init__(self):
    self.func = Rosenbrock
    self.func_type = Rosenbrock_func_type
    self.variables_num = Rosenbrock_variables_num
    self.lower_bound = Rosenbrock_lower_bound
    self.upper_bound = Rosenbrock_upper_bound
    self.T_start = 100
    self.T_end = 1e-6
    self.q = 0.9
    self.L = 100
    self.init_pos = None
    #self.complex_constraints = [constraints1,constraints2,constraints3]
    self.complex_constraints_method = complex_constraints_method.loop
    self.SA = SA(**self.__dict__)
def test(self):
    start_time = time()
    self.SA.run()
    print("SA costs %.4f seconds!" %(time()-start_time))
    self.SA.save_plot()
    self.SA.show_result()
if __name__ == '__main__':
TestSA().test()           
    SA costs 0.2039 seconds!
-------------------- SA config is: --------------------
func_type:max
complex_constraints:None
q:0.9
complex_constraints_method:loop
T_start:100
steps:17500
T_end:1e-06
L:100
func:<function Rosenbrock at 0x7f8261799048>
variables_num:2
lower_bound:[-2.048 -2.048]
init_pos:[-2.03887265 -2.02503927]
upper_bound:[2.048 2.048]
-------------------- SA calculation result is: --------------------
global best generation index/total generations:2126/17500
global best point: [-2.03887265 -2.02503927]
global best target: 3830.997799328349

    圖4 SA求解Rosenbrock函數

    SA類的具體參數含義如下:

  • variables_num

  • variables_num
  • T_start:初始溫度,預設值為100;
  • T_end:末尾溫度,預設值為\(10^{-6}\);
  • q:退火系數,預設值為0.9,一般取一個在[0.9,1)之間的數字;
  • L:每個溫度時的疊代次數，即鍊長,預設值為100;
  • init_pos:初始解的取值，預設值為None,此時初始解是随機生成的,或者你也可以指定一個用ndarray表示的初始解;
  • complex_constraints:複雜限制，預設為None,表示沒有限制，具體定義同GA 的 complex_constraints;

  • 5. Random Walk 使用

        Random Walk 即随機遊走算法，是一種全局随機搜尋算法，具體原理可以參考我的部落格 介紹一個全局最優化的方法：随機遊走算法(Random Walk) 。sopt中Random Walk的使用執行個體如下: 
    from time import time
from sopt.util.functions import *
from sopt.util.constraints import *
from sopt.util.random_walk_config import *
from sopt.Optimizers.RandomWalk import RandomWalk
class TestRandomWalk:
def __init__(self):
    self.func = quadratic50
    self.func_type = quadratic50_func_type
    self.variables_num = quadratic50_variables_num
    self.lower_bound = quadratic50_lower_bound
    self.upper_bound = quadratic50_upper_bound
    self.generations = 100
    self.init_step = 10
    self.eps = 1e-4
    self.vectors_num = 10
    self.init_pos = None
    # self.complex_constraints = [constraints1,constraints2,constraints3]
    self.complex_constraints = None
    self.complex_constraints_method = complex_constraints_method.loop
    self.RandomWalk = RandomWalk(**self.__dict__)
def test(self):
    start_time = time()
    self.RandomWalk.random_walk()
    print("random walk costs %.4f seconds!" %(time() - start_time))
    self.RandomWalk.save_plot()
    self.RandomWalk.show_result()
if __name__ == '__main__':
TestRandomWalk().test()           
    Finish 1 random walk!
Finish 2 random walk!
Finish 3 random walk!
Finish 4 random walk!
Finish 5 random walk!
Finish 6 random walk!
Finish 7 random walk!
Finish 8 random walk!
Finish 9 random walk!
Finish 10 random walk!
Finish 11 random walk!
Finish 12 random walk!
Finish 13 random walk!
Finish 14 random walk!
Finish 15 random walk!
Finish 16 random walk!
Finish 17 random walk!
random walk costs 1.0647 seconds!
-------------------- random walk config is: --------------------
init_step:10
eps:0.0001
generations_nums:9042
lower_bound:[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1]
complex_constraints_method:loop
walk_nums:17
complex_constraints:None
vectors_num:10
upper_bound:[50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
 50 50]
variables_num:50
func:<function quadratic50 at 0x7f51555a1840>
generations:100
func_type:min
-------------------- random walk caculation result is: --------------------
global best generation index/total generations:8942/9042
global best point is: [ 1.00004803  1.99999419  3.00000569  3.99998558  5.00002455  5.99999255
  6.99992476  7.99992864  9.00000401  9.99994717 10.99998155 12.00002429
 13.0000035  13.99998567 15.00000421 16.00001454 16.99997252 17.99998041
 19.00002491 20.00003141 21.00004182 21.99998565 22.99997668 23.99999821
 24.99995881 25.99999359 27.00000443 28.00005117 28.99998132 30.00004136
 31.00002021 32.00000616 33.00000678 34.00005423 35.00001799 36.00000051
 37.00002749 38.00000203 39.00007087 39.9999964  41.00004432 42.0000158
 42.99992991 43.99995352 44.99997267 46.00003533 46.9999834  47.99996778
 49.00002904 50.        ]
global best target is: 1.0000000528527013

    圖5 Random Walk 求解quadratic50

    經過實驗發現,Random Walk 具有非常強的全局尋優能力，對于quadratic50這種具有50個變量的複雜目标函數，它也可以很快找到其全局最優點，而且運作速度也很快。RandomWalk類的具體參數含義如下:

  • variables_num

  • variables_num
  • generations:每個step的最大疊代次數,預設是100;
  • init_step:初始步長(step)，預設值為10.0;
  • eps:終止疊代的步長,預設為\(10^{-4}\);
  • vectors_num:使用 improved_random_walk時随機産生向量的個數,預設為10;
  • RandomWalk 除了提供基本的random_walk函數之外，還提供了一個更加強大的improved_random_walk函數，後者的全局尋優能力要更強。

    6. 求解帶複雜限制的目标函數

        上面所述的各種優化方法求解的都是變量僅有簡單邊界限制(形如\(a \le x_i \le b\)),下面介紹如何使用各種優化方法求解帶有複雜限制條件的目标函數。其實，求解方法也非常簡單，以GA為例，下面的例子即對Rosenbrock函數求解了帶有三個複雜限制條件的最優值: 
    from time import time
from sopt.GA.GA import GA
from sopt.util.functions import *
from sopt.util.ga_config import *
from sopt.util.constraints import *
class TestGA:
def __init__(self):
    self.func = Rosenbrock
    self.func_type = Rosenbrock_func_type
    self.variables_num = Rosenbrock_variables_num
    self.lower_bound = Rosenbrock_lower_bound
    self.upper_bound = Rosenbrock_upper_bound
    self.cross_rate = 0.8
    self.mutation_rate = 0.1
    self.generations = 300
    self.population_size = 200
    self.binary_code_length = 20
    self.cross_rate_exp = 1
    self.mutation_rate_exp = 1
    self.code_type = code_type.real 
    self.cross_code = False
    self.select_method = select_method.proportion
    self.rank_select_probs = None
    self.tournament_num = 2
    self.cross_method = cross_method.uniform
    self.arithmetic_cross_alpha = 0.1
    self.arithmetic_cross_exp = 1
    self.mutation_method = mutation_method.uniform
    self.none_uniform_mutation_rate = 1
    self.complex_constraints = [constraints1,constraints2,constraints3]
    self.complex_constraints_method = complex_constraints_method.penalty
    self.complex_constraints_C = 1e8
    self.M = 1e8
    self.GA = GA(**self.__dict__)
def test(self):
    start_time = time()
    self.GA.run()
    print("GA costs %.4f seconds!" % (time()-start_time))
    self.GA.save_plot()
    self.GA.show_result()
if __name__ == '__main__':
TestGA().test()           
    GA costs 1.9957 seconds!
-------------------- GA config is: --------------------
lower_bound:[-2.048, -2.048]
cross_code:False
complex_constraints_method:penalty
mutation_method:uniform
mutation_rate:0.1
mutation_rate_exp:1
cross_rate:0.8
upper_bound:[2.048, 2.048]
arithmetic_cross_exp:1
variables_num:2
generations:300
tournament_num:2
select_method:proportion
func_type:max
complex_constraints_C:100000000.0
cross_method:uniform
complex_constraints:[<function constraints1 at 0x7f5efe2e8d08>, <function constraints2 at 0x7f5efe2e8d90>, <function constraints3 at 0x7f5efe2e8e18>]
func:<function Rosenbrock at 0x7f5efe2e87b8>
none_uniform_mutation_rate:1
cross_rate_exp:1
code_type:real
M:100000000.0
binary_code_length:20
global_generations_step:300
population_size:200
arithmetic_cross_alpha:0.1
-------------------- GA caculation result is: --------------------
global best target generation index/total generations:226/300
global best point:[ 1.7182846  -1.74504313]
global best target:2207.2089435117955

    圖6 GA求解帶有三個複雜限制的Rosenbrock函數

    上面的constraints1,constraints2,constraints3是三個預定義的限制條件函數，其定義分别為:\(constraints1:x_1^2 + x_2^2 - 6 \le 0\);\(constraints2:x_1 + x_2 \le 0\);\(constraints3:-2-x_1 - x_2 \le 0\),函數原型為:

    def constraints1(x):
x1 = x[0]
x2 = x[1]
return x1**2 + x2**2 -3
def constraints2(x):
x1 = x[0]
x2 = x[1]
return x1+x2
def constraints3(x):
x1 = x[0]
x2 = x[1]
return -2 -x1 -x2
    其實觀察可以發現，上面的代碼和原始的GA執行個體代碼唯一的差別，就是其增加了

    self.complex_constraints = [constraints1,constraints2,constraints3]

     這樣一句，對于其他的優化方法，其都定義了

    complex_constraints

    complex_constraints_method

     這兩個屬性，隻要傳入相應的限制條件函數清單以及求解限制條件的方法就可以求解帶複雜限制的目标函數了。比如我們再用Random Walk求解和上面一樣的帶三個限制的Rosenbrock函數，代碼及運作結果如下: 
    from time import time
from sopt.util.functions import *
from sopt.util.constraints import *
from sopt.util.random_walk_config import *
from sopt.Optimizers.RandomWalk import RandomWalk
class TestRandomWalk:
def __init__(self):
    self.func = Rosenbrock
    self.func_type = Rosenbrock_func_type
    self.variables_num = Rosenbrock_variables_num
    self.lower_bound = Rosenbrock_lower_bound
    self.upper_bound = Rosenbrock_upper_bound
    self.generations = 100
    self.init_step = 10
    self.eps = 1e-4
    self.vectors_num = 10
    self.init_pos = None
    self.complex_constraints = [constraints1,constraints2,constraints3]
    self.complex_constraints_method = complex_constraints_method.loop
    self.RandomWalk = RandomWalk(**self.__dict__)
def test(self):
    start_time = time()
    self.RandomWalk.random_walk()
    print("random walk costs %.4f seconds!" %(time() - start_time))
    self.RandomWalk.save_plot()
    self.RandomWalk.show_result()
if __name__ == '__main__':
TestRandomWalk().test()
    運作結果: 
    Finish 1 random walk!
Finish 2 random walk!
Finish 3 random walk!
Finish 4 random walk!
Finish 5 random walk!
Finish 6 random walk!
Finish 7 random walk!
Finish 8 random walk!
Finish 9 random walk!
Finish 10 random walk!
Finish 11 random walk!
Finish 12 random walk!
Finish 13 random walk!
Finish 14 random walk!
Finish 15 random walk!
Finish 16 random walk!
Finish 17 random walk!
random walk costs 0.1543 seconds!
-------------------- random walk config is: --------------------
eps:0.0001
func_type:max
lower_bound:[-2.048 -2.048]
upper_bound:[2.048 2.048]
init_step:10
vectors_num:10
func:<function Rosenbrock at 0x7f547fc952f0>
variables_num:2
walk_nums:17
complex_constraints_method:loop
generations:100
generations_nums:2191
complex_constraints:[<function constraints1 at 0x7f547fc95bf8>, <function constraints2 at 0x7f547fc95c80>, <function constraints3 at 0x7f547fc95d08>]
-------------------- random walk caculation result is: --------------------
global best generation index/total generations:2091/2191
global best point is: [-2.41416736  0.41430367]
global best target is: 2942.6882849234585

    圖7 Random Walk求解帶有三個複雜限制的Rosenbrock函數

    可以發現Random Walk 求解得到的最優解要比GA好，而且運作時間更快，經過實驗發現，在所有的優化方法中，不論是求解帶複雜限制還是不帶複雜限制條件的目标函數，求解效果大體上排序是:Random Walk > PSO > GA > SA 。是以當你在求解具體問題時，不妨多試幾種優化方法，然後擇優選擇。

    7. 基于梯度的系列優化方法

        上面所述的各種優化方法，比如GA,PSO,SA等都是基于随機搜尋的優化算法，其計算是不依賴于目标函數的具體形式，也不需要知道其梯度的，更加傳統的優化算法是基于梯度的算法，比如經典的梯度下降(上升)法(Gradient Descent)以及其一系列變種。下面就簡要介紹sopt中GD,Momentum,AdaGrad,RMSProp以及Adam的實作。關于這些基于梯度的優化算法的具體原理，可以參考我之前的一篇博文 深度學習中常用的優化方法 。另外需要注意，以下所述的基于梯度的各種優化算法，一般都是用在無限制優化問題裡面的，如果是有限制的問題，請選擇上面其他的優化算法。下面是GradientDescent的一個使用執行個體: 
    from time import time
from sopt.util.gradients_config import *
from sopt.util.functions import *
from sopt.Optimizers.Gradients import GradientDescent
class TestGradientDescent:
def __init__(self):
    self.func = quadratic50
    self.func_type = quadratic50_func_type
    self.variables_num = quadratic50_variables_num
    self.init_variables = None
    self.lr = 1e-3
    self.epochs = 5000
    self.GradientDescent = GradientDescent(**self.__dict__)
def test(self):
    start_time = time()
    self.GradientDescent.run()
    print("Gradient Descent costs %.4f seconds!" %(time()-start_time))
    self.GradientDescent.save_plot()
    self.GradientDescent.show_result()
if __name__ == '__main__':
TestGradientDescent().test()           
    Gradient Descent costs 14.3231 seconds!
-------------------- Gradient Descent config is: --------------------
func_type:min
variables_num:50
func:<function quadratic50 at 0x7f74e737b620>
epochs:5000
lr:0.001
-------------------- Gradient Descent caculation result is: --------------------
global best epoch/total epochs:4999/5000
global best point: [ 0.9999524   1.99991045  2.99984898  3.9998496   4.99977767  5.9997246
  6.99967516  7.99964102  8.99958143  9.99951782 10.99947879 11.99944665
 12.99942492 13.99935192 14.99932708 15.99925856 16.99923686 17.99921689
 18.99911527 19.9991255  20.99908968 21.99899699 22.99899622 23.99887832
 24.99883597 25.99885616 26.99881394 27.99869772 28.99869349 29.9986766
 30.99861142 31.99851987 32.998556   33.99849351 34.99845985 35.99836731
 36.99832444 37.99831792 38.99821067 39.99816567 40.99814951 41.99808199
 42.99808161 43.99806655 44.99801207 45.99794449 46.99788003 47.99785468
 48.99780825 49.99771656]
global best target: 1.0000867498727912

    圖7 GradientDescent 求解quadratic50

    下面簡要說明以下GradientDescent,Momentum等類中的主要參數,像

    func

    variables_num

     等含義已經解釋很多次了，不再贅述，這裡主要介紹各類特有的一些參數。
  1. GradientDescent類:
  • lr:學習率,預設是\(10^{-3}\);
  • epochs:疊代次數，預設是100；
  1. Momentum類:
  • beta:Momentum的\(\beta\)參數，預設是0.9;
  1. AdaGrad類:
  • eps:極小的一個正數(防止分母為0),預設為\(10^{-8}\);
  1. RMSProp類:
  • beta:RMSProp的\(\beta\)參數，預設是0.9;
  1. Adam類:
  • beta1:Adam的\(\beta_1\)參數，預設是0.5;
  • beta2:Adam的\(\beta_2\)參數，預設是0.9;
  • 這裡額外說一句,Adam原文作者給出的最優超參數為\(\beta_1 = 0.9,\beta_2 = 0.999\),但是我在實際調試過程中發現,\(\beta_1 = 0.5,\beta_2 = 0.9\)才能達到比較好的效果，具體原因目前不是很清楚。如果想要學習各個優化函數類的具體用法，還可以使用

    from sopt.test import *

     導入預定義的示例測試，來運作觀察結果。比如下面的代碼就運作了PSO的一個示例測試: 
    from sopt.test import test_PSO
test_PSO.TestPSO().test()
    結果: 
    PSO costs 3.4806 seconds!
-------------------- PSO config is: --------------------
lower_bound:[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1]
generations:200
vmin:-1
func:<function quadratic50 at 0x7fd3bf37d8c8>
w_method:linear_decrease
func_type:min
population_size:100
w_start:0.9
complex_constraints_method:loop
vmax:1
c2:1.49445
upper_bound:[50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50]
variables_num:50
c1:1.49445
w:1
complex_constraints:None
w_end:0.4
-------------------- PSO calculation result is: --------------------
global best generation index/total generations: 199/200
global best point: [  1.           1.           2.94484328   4.13875216   5.00293498
   6.13124759   6.99713025   7.92116383   8.87648843  10.02066994
  11.0758768   12.02240279  13.01125368  13.98010373  14.98063168
  15.97776149  17.11878537  18.00246112  18.14780887  20.00637617
  21.00223704  22.00689373  23.14823218  24.0002456   24.98672157
  25.99141686  27.02112321  28.01540506  29.05403155  30.07304888
  31.00414822  32.00982867  32.99444884  33.9114213   34.96631157
  36.22871824  37.0015616   37.98907918  39.01245751  40.1371835
  41.0182043   42.07768102  42.87178292  43.93687997  45.05786395
  46.03778693  47.07913415  50.          48.9964866   50.        ]
global best target: 6.95906097685

    附錄

        附錄會簡要說明上文提到的一些概念，待有空更新。。。

    To do ...

  1. GA的并行化計算
  2. 小生境GA
  3. GA,PSO等求解離散函數最優解
  4. GA,PSO等實際應用舉例
  5. 其他的一些最優化方法,比如牛頓法，拟牛頓法等
  6. 其他

熱愛程式設計，熱愛機器學習！

github:http://www.github.com/Lyrichu

github blog:http://Lyrichu.github.io

個人部落格站點:http://www.movieb2b.com(不再維護)

算法 c語言 Python python自然語言 python算法 python語言開發 python腳本語言 執行python script
上一篇: Jenkins修改管理者密碼.
下一篇: Hibernate緩存和狀态

繼續閱讀