目次
はじめに
コスト関数の微分不要なソルバーには、
Nelder Meadが有名ですが、
このソルバーは基本的に局所最適化用なので、
複雑な問題には局所解に落ち込んでしまうという問題があります。
そこで微分不要な大域最適化ソルバーとして、
差分進化法(Differential Evolution)があります。
この記事では差分進化法による大域的最適化の概要と
Pythonサンプルコードを紹介したいと思います。
差分進化法とは?
差分進化法(Differential Evolution)は、
大域的最適化アルゴリズムの一つです。
進化的アルゴリズム(Evolutionary Algorithm)の一種であり、
確率的な探索によって、解候補個体を用いた多点探索を行う手法です。
DEは非線形モデルや、微分不可能な問題などの
現実的な実行時間では厳密な解を求めることが困難な問題に対して、
近似解を求めることが可能であり、
様々な最適化問題に利用されています。
SciPyのoptimizeモジュールでも大域最適化ソルバーとして利用されています。
DEの特徴としてランダムに選択した個体間の差分ベクトルを用いて
新たな個体を生成するため、差分進化法(Differential Evolution)と呼ばれます。
差分ベクトルを使うことで、探索の最初は広い範囲に探索点を持ちますが、
計算が進むにつれて、解の近くに探索点が集まってくる形になります。
各候補点の差分ベクトルを使うというのはNelder-Mead法に似ているため、
差分進化法はNelder-Mead法に進化的アルゴリズム(EA)を応用したものと説明する人もいます。
差分進化法のアルゴリズム概要
差分進化法では、初期探索点を一定範囲にばらまいたあと、
下記のプロセスをそれぞれの探索点に対して繰り返して、
探索点を真の解の近くに近づけていきます。
1.Mutation(突然変異)
ある候補点に着目して、それ以外の候補点を3つ選び、
その内、一つの候補点を、残り2つの候補点の差分ベクトルで
移動した点(変異点)を生成します。
これはNelder-meadのReflectionに似た処理です。
2.Crossover(交叉)
着目した候補点と、1で生成された変異点の両方の一部を
合成した点(交叉点)を生成します。
お互いにどれだけのデータを合成するかは、
必ず一部は変異点の分が選ばれるように、乱数などで選択されます。
3.Selection(生存選択)
最後に生成された交叉点が、注目点よりもコストが良ければ、
注目点を交叉点に置き換えます。
差分進化法のPythonサンプルコード
下記のコードを実行すると、
冒頭のgifアニメーションのような最適化を実行できます。
import numpy as np def differential_evolution(func, bounds, *, pop_size=15, max_iter=1000, F=0.5, cr=0.7, ftol=10**-8, callback=None, rng=None): """ differential evolution solver :param func: The objective function to be minimized. :param bounds: Bounds for variables. (min, max) pairs for each element in x, defining the finite lower and upper bounds for the optimizing argument of func. :param pop_size: :param max_iter: The maximum number of iterations over which the entire population is evolved :param F: The mutation constant. it should be in the range [0, 2] :param cr: The crossover probability constant, should be in the range [ 0, 1]. Increasing this value allows a larger number of mutants to progress into the next generation, but at the risk of population stability. :param ftol: Relative tolerance for convergence. :param callback: A function to follow the progress of the minimization. Arguments are i: iteration count best_x: current best x best_obj: current best objective populations: current populations :param rng: Random generator :return: best_x: found best x, best_obj: found best objective """ if rng is None: rng = np.random.default_rng() xdim = len(bounds) populations = rng.uniform(low=bounds[:, 0], high=bounds[:, 1], size=(pop_size, xdim)) obj_list = [func(pop) for pop in populations] best_x = populations[np.argmin(obj_list)] best_obj = min(obj_list) prev_obj = best_obj for i in range(max_iter): for j in range(pop_size): # (1) mutation mutated = mutation(F, bounds, j, pop_size, populations, rng) # (2) crossover trial = crossover(mutated, populations[j], xdim, cr, rng) # (3) selection obj_list, populations = selection(func, j, obj_list, populations, trial) best_obj = min(obj_list) best_x = populations[np.argmin(obj_list)] if best_obj < prev_obj: if abs(prev_obj - best_obj) <= ftol: break # converge prev_obj = best_obj if callback is not None: callback(i, best_x, best_obj, populations) return best_x, best_obj def mutation(F, bounds, j, pop_size, populations, rng): indexes = [i for i in range(pop_size) if i != j] # pick up 3 points. a, b, c = populations[rng.choice(indexes, 3, replace=False)] mutated = a + F * (b - c) mutated = np.clip(mutated, bounds[:, 0], bounds[:, 1]) return mutated def crossover(mutated, target, dims, cr, rng): p = rng.random(dims) # At least one element of mutated should be selected. p[rng.choice([i for i in np.arange(len(p))], 1)] = 0.0 trial = [mutated[i] if p[i] < cr else target[i] for i in range(dims)] return trial def selection(func, j, obj_list, populations, trial): obj_trial = func(trial) if obj_trial < obj_list[j]: # better point is found populations[j] = trial obj_list[j] = obj_trial return obj_list, populations def himmelblau(x, y): a = x ** 2 + y - 11 b = x + y ** 2 - 7 return a ** 2 + b ** 2 def main(): import matplotlib.pyplot as plt def plot_callback(i, best_vector, best_obj, pop): if i % 5 == 0: print(f"Iter: {i} f([{best_vector}]) = {best_obj}") plt.cla() x = np.linspace(-5.0, 5.0, 100) y = np.linspace(-5.0, 5.0, 100) XX, YY = np.meshgrid(x, y) plt.contour(x, y, himmelblau(XX, YY), levels=30) plt.plot(pop[:, 0], pop[:, 1], ".r") plt.plot(best_vector[0], best_vector[1], "xr") plt.pause(0.1) # define lower and upper bounds for every dimension bounds = np.asarray([(-5.0, 5.0), (-5.0, 5.0)]) xopt, fopt = differential_evolution(lambda x: himmelblau(x[0], x[1]), bounds, callback=plot_callback, rng=np.random.default_rng(12345) ) print(f"Solution: f({xopt}) = {fopt}") plt.show() if __name__ == '__main__': main()
参考資料
MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。