「差分進化」(Differential Evolution, DE)はヒューリスティックな（大域的）最適化問題をとく標準的なアルゴリズムの一つでscipyにも実装されている。
ここではscipy.optimize.differential_evolutionモジュールの使い方をまとめておく。
アルゴリズム自体の説明は行わないが、Wikipediaの記事がよくまとまっているのでそちらを参照していただきたい。
ここでは基本的な使い方に加えて、特に並列に評価関数を評価する方法についてまとめる。
前提条件
- python 3
- scipy 1.2.0
サンプルコード
もっとも単純なコード
もっとも単純なサンプルコードは以下の通りである。
import pprint
import numpy as np
from scipy.optimize import differential_evolution
bounds = [(0,2), (0, 2), (0, 2)] # 探索するxの定義域範囲
def func(x):
return np.sum( (x-1) * (x-1) ) # 最適化する関数
result = differential_evolution(func, bounds)
pprint.pprint(result) # 出力を見やすくするためにpprintを使って出力している
この例では $f(x_1,x_2,x_3) = \sum_i (x_i-1)^2$ という関数の最小値を求めている。実行すると
{'fun': 0.0,
'message': 'Optimization terminated successfully.',
'nfev': 6214,
'nit': 137,
'success': True,
'x': array([1., 1., 1.])}
となり、最適解$(1,1,1)$がもとまっていることがわかる。xが最適解の座標、funが最適解での評価関数の値、nfevが評価関数を読んだ回数、nitが反復の回数（世代数）である。
disp=True
disp=True を追加すると実行時に世代ごとにログ出力されるようになる。
result = differential_evolution(func, bounds, disp=True)
differential_evolution step 1: f(x)= 0.0363989
differential_evolution step 2: f(x)= 0.0363989
differential_evolution step 3: f(x)= 0.0186435
differential_evolution step 4: f(x)= 0.00906504
differential_evolution step 5: f(x)= 0.00906504
...
maxiter={n}
最大の世代数を指定する。評価関数の評価回数が最大で(maxiter + 1) * popsize * len(x)になる。
result = differential_evolution(func, bounds, maxiter=10, disp=True)
differential_evolution step 1: f(x)= 0.0868267
differential_evolution step 2: f(x)= 0.00339551
differential_evolution step 3: f(x)= 0.00339551
differential_evolution step 4: f(x)= 0.00279434
differential_evolution step 5: f(x)= 0.00175642
differential_evolution step 6: f(x)= 0.000346342
differential_evolution step 7: f(x)= 0.000346342
differential_evolution step 8: f(x)= 0.00020965
differential_evolution step 9: f(x)= 0.000123287
differential_evolution step 10: f(x)= 6.17462e-05
{'fun': 7.500000685994094e-17,
'jac': array([ 0.00000000e+00, -1.11022303e-16, -1.33226762e-15]),
'message': 'Maximum number of iterations has been exceeded.',
'nfev': 511,
'nit': 10,
'success': False,
'x': array([1. , 0.99999999, 0.99999999])}
DEのパラメータの変更
世代ごとの個体の更新時のパラメータとしてstrategy, popsize, mutation, recombination を指定できる。
試してみたところデフォルトの値でも大概うまくいくようだが、より速く大域最適解に収束させるには評価関数の関数形にあった適切な値を選ぶ必要がある。
updating='immediate'
個体の更新を行うタイミングを指定できる。immediate, deferredのどちらかから選ぶ。
immediateの場合には各個体の評価が終わるごとに次々に個体が入れ替えられていく。
一方でdeferredの場合は各世代の個体すべての評価が終わるまでまってから一斉に個体の入れ替えを行う。
一般的にimmediateの方がより少ない世代数で収束するが、deferredの場合は並列に各個体の評価を実施することができる。
ためしに上記の関数で最適化を実施した場合、immediateの場合は123世代、deferredの場合は177世代で最適化が完了した。
workers={n}
評価関数の評価を並列に行うためのオプションも実装されている。multiprocessing.Poolモジュールを用いて複数プロセスで並列に動作させることが可能。
オプションに並列プロセス数を指定する。-1を指定した場合は動作しているマシンのコア数に応じて決まる。
これにより評価関数の評価に時間がかかる場合に大幅に実行時間を短くすることができる。
workersが１より大きい場合には、updatingのオプションは指定しなくてもdeferredになる。
import pprint
import numpy as np
from time import sleep
from scipy.optimize import differential_evolution
bounds = [(0,2), (0, 2), (0, 2)]
def func(x):
sleep(0.05) # 評価関数の実行に0.05secを要する
return np.sum( (x-1) * (x-1) )
result = differential_evolution(func, bounds, updating='deferred', workers=-1, disp=True)
pprint.pprint(result)
例えば12coreのマシンで実行するとトータル30secで完了するが、シリアルに実行すると桁違いに遅い。
workersに整数値を渡した場合にはmultiprocessing.Pool.mapで並列に実行されるが、この並列処理の部分を自分のロジックで置き換えることもできる。
例えば、以下のようなコードを実行すると出力は下のようになる。
import pprint
import numpy as np
from time import sleep
from scipy.optimize import differential_evolution
bounds = [(0,2), (0, 2), (0, 2)]
def func(x):
return np.sum( (x-1) * (x-1) )
def my_map(f, xs):
print('my_map is called')
sleep(0.05)
return map(f, xs)
result = differential_evolution(func, bounds, maxiter=10, updating='deferred', workers=my_map, disp=True)
pprint.pprint(result)
my_map is called
my_map is called
differential_evolution step 1: f(x)= 0.12497
my_map is called
differential_evolution step 2: f(x)= 0.0522403
my_map is called
differential_evolution step 3: f(x)= 0.0114089
my_map is called
differential_evolution step 4: f(x)= 0.0114089
my_map is called
differential_evolution step 5: f(x)= 0.0114089
my_map is called
differential_evolution step 6: f(x)= 0.00742643
my_map is called
differential_evolution step 7: f(x)= 0.00742643
my_map is called
differential_evolution step 8: f(x)= 0.00742643
my_map is called
differential_evolution step 9: f(x)= 0.00335309
my_map is called
differential_evolution step 10: f(x)= 0.00110949
{'fun': 3.052203109076959e-17,
'jac': array([3.92205022e-09, 1.90604633e-08, 1.17477565e-08]),
'message': 'Maximum number of iterations has been exceeded.',
'nfev': 507,
'nit': 10,
'success': False,
'x': array([1., 1., 1.])}