search
LoginSignup
8

More than 1 year has passed since last update.

posted at

最適化アルゴリズムを実装していくぞ(差分進化)

はじめに

最適化アルゴリズムの実装シリーズです。
まずは概要を見てください。

コードはgithubにあります。

差分進化

差分進化(Differential evolution:DE)は、生物の進化をもとに考案された遺伝的アルゴリズムと似ている手法です。

大きく、Mutation(突然変異)、Crossover(交叉)、Selection(生存選択)の3つのフェーズがあります。
遺伝的アルゴリズムみたいにフェーズ毎にアルゴリズムがいろいろあるようですが、本記事では一番簡単なアルゴリズムで実装しています。

参考
差分進化法でハイパーパラメータチューニング
差分進化(Wikipedia)

  • アルゴリズムのフロー

draw2-DE.png

  • 用語の対応
問題 差分進化
入力値の配列 エージェントの位置
入力値 エージェント
評価値 エージェントの評価値
  • ハイパーパラメータに関して
変数名 意味 所感
crossover_rate 交叉する確率 変異ベクトルになる確率
scaling 変異ベクトルのスケール率(長さ) 大きいほど移動範囲が大きい

Mutation(突然変異)

ランダムに選んだ3人のエージェントから変異ベクトルを作成します。

import random

i = 自分を表すindex

# i番目を除いた3エージェントをランダムに選択
r1, r2, r3 = random.sample([ j for j in range(len(agents)) if j != i ], 3)
pos1 = agents[r1].getArray()
pos2 = agents[r2].getArray()
pos3 = agents[r3].getArray()

変異ベクトルは以下の式で出します。

$$ \vec{v} = \vec{x_1} + F (\vec{x_2} - \vec{x_3}) $$

$F$ は変異差分の適用率(スケール因子:scaling factor)を表す0~2の実数です。
pythonコードだと以下になります。

# ベクトル計算しやすいようにnumpy化
import numpy as np
pos1 = np.asarray(pos1)
pos2 = np.asarray(pos2)
pos3 = np.asarray(pos3)

m_pos = pos1 + scaling * (pos2 - pos3)

Crossover(交叉)

一様交叉(binary crossover)で交叉します。
ある確率(crossover_rate)で成分を入れ替える交叉となります。

また交叉は、変異ベクトル側の要素を1成分だけ必ず取り入れます。

import random

# 変異ベクトルで交叉させる(一様交叉)
pos = agent.getArray()
ri = random.randint(0, len(pos))  # 1成分は必ず変異ベクトル
for j in range(len(pos)):
    if  ri == j or random.random() < self.crossover_rate:
        pos[j] = m_pos[j]
    else:
        pass  # 更新しない

# 新しい位置のエージェントを作成
new_agent = problem.create(pos)

生存選択

交叉でできた新しいエージェントと今のエージェントを比べ、更新されていれば置き換えます。

# 優れている個体なら置き換える
if agents[i].getScore() < new_agent.getScore():
    agents[i] = new_agent

コード全体

コード全体です。
クラス化しているので上記コードと少し違います。

import math
import random

import numpy as np

class DE():
    def __init__(self,
        agent_max,           # エージェント数
        crossover_rate=0.5,  # 交叉率
        scaling=0.5,         # 差分の適用率
    ):
        self.agent_max = agent_max
        self.crossover_rate = crossover_rate
        self.scaling = scaling

    def init(self, problem):
        self.problem = problem

        # 初期位置の生成
        self.agents = []
        for _ in range(self.agent_max):
            self.agents.append(problem.create())


    def step(self):

        for i, agent in enumerate(self.agents):

            # iを含まない3個体をランダムに選択
            r1, r2, r3 = random.sample([ j for j in range(len(self.agents)) if j != i ], 3)
            pos1 = self.agents[r1].getArray()
            pos2 = self.agents[r2].getArray()
            pos3 = self.agents[r3].getArray()

            # 3個体から変異ベクトルをだす
            pos1 = np.asarray(pos1)
            pos2 = np.asarray(pos2)
            pos3 = np.asarray(pos3)
            m_pos = pos1 + self.scaling * (pos2 - pos3)

            # 変異ベクトルで交叉させる(一様交叉)
            pos = agent.getArray()
            ri = random.randint(0, len(pos))  # 1成分は必ず変異ベクトル
            for j in range(len(pos)):
                if  ri == j or random.random() < self.crossover_rate:
                    pos[j] = m_pos[j]
                else:
                    pass  # 更新しない

            # 優れている個体なら置き換える
            new_agent = self.problem.create(pos)
            self.count += 1
            if agent.getScore() < new_agent.getScore():
                self.agents[i] = new_agent


ハイパーパラメータ例

各問題に対して optuna でハイパーパラメータを最適化した結果です。
最適化の1回の試行は、探索時間を2秒間として結果を出しています。
これを100回実行し、最適なハイパーパラメータを optuna に探してもらいました。

問題 agent_max crossover_rate scaling
EightQueen 5 0.008405098138137779 1.7482804860765253
function_Ackley 36 0.4076390525351224 0.2908895854800526
function_Griewank 14 0.27752386128521395 0.4629100940098222
function_Michalewicz 12 0.1532879607238835 0.0742830755371933
function_Rastrigin 28 0.33513859646880306 0.0754225020709786
function_Schwefel 13 0.00032331965923372563 0.13153649005308807
function_StyblinskiTang 39 0.21247741932099348 0.08732185323441227
function_XinSheYang 33 0.0955103914325307 0.008270969294347359
LifeGame 39 0.6612227467897149 1.136453380180552
OneMax 4 0.1190487045395953 1.1581036102901494
TSP 23 0.41212989299137665 0.014644735558753091

実際の動きの可視化

1次元は6個体、2次元は20個体で50step実行した結果です。
赤い丸がそのstepでの最高スコアを持っている個体となります。

パラメータは以下で実行しました。

DE(N, crossover_rate=0, scaling=0.4)

function_Ackley

  • 1次元

function_Ackley_DE_2.gif

  • 2次元

function_Ackley_DE_3.gif

function_Rastrigin

  • 1次元

ffunction_Rastrigin_DE_2.gif

  • 2次元

function_Rastrigin_DE_3.gif

function_Schwefel

  • 1次元

function_Schwefel_DE_2.gif

  • 2次元

function_Schwefel_DE_3.gif

function_StyblinskiTang

  • 1次元

function_StyblinskiTang_DE_2.gif

  • 2次元

function_StyblinskiTang_DE_3.gif

function_XinSheYang

  • 1次元

function_XinSheYang_DE_2.gif

  • 2次元

function_XinSheYang_DE_3.gif

あとがき

いい感じに収束していきますね。
ただランダム移動がないので局所解に陥りやすい気はします。

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
8