はじめに

コスト関数の微分不要なソルバーには、

Nelder Meadが有名ですが、

myenigma.hatenablog.com

このソルバーは基本的に局所最適化用なので、

複雑な問題には局所解に落ち込んでしまうという問題があります。

そこで微分不要な大域最適化ソルバーとして、

差分進化法(Differential Evolution)があります。

この記事では差分進化法による大域的最適化の概要と

Pythonサンプルコードを紹介したいと思います。

差分進化法とは？

差分進化法(Differential Evolution)は、

大域的最適化アルゴリズムの一つです。

進化的アルゴリズム(Evolutionary Algorithm)の一種であり、

確率的な探索によって、解候補個体を用いた多点探索を行う手法です。

DEは非線形モデルや、微分不可能な問題などの

現実的な実行時間では厳密な解を求めることが困難な問題に対して、

近似解を求めることが可能であり、

様々な最適化問題に利用されています。

SciPyのoptimizeモジュールでも大域最適化ソルバーとして利用されています。

docs.scipy.org

DEの特徴としてランダムに選択した個体間の差分ベクトルを用いて

新たな個体を生成するため、差分進化法(Differential Evolution)と呼ばれます。

差分ベクトルを使うことで、探索の最初は広い範囲に探索点を持ちますが、

計算が進むにつれて、解の近くに探索点が集まってくる形になります。

各候補点の差分ベクトルを使うというのはNelder-Mead法に似ているため、

myenigma.hatenablog.com

差分進化法はNelder-Mead法に進化的アルゴリズム(EA)を応用したものと説明する人もいます。

差分進化法のアルゴリズム概要

差分進化法では、初期探索点を一定範囲にばらまいたあと、

下記のプロセスをそれぞれの探索点に対して繰り返して、

探索点を真の解の近くに近づけていきます。

1.Mutation(突然変異)

ある候補点に着目して、それ以外の候補点を3つ選び、

その内、一つの候補点を、残り２つの候補点の差分ベクトルで

移動した点(変異点)を生成します。

これは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()