Zenn
🦕

最適化アルゴリズムのJulia実装 | 差分進化 (Differential Evolution)

2021/02/13に公開
Julia
最適化
tech

関数 f の最適化技法として，差分進化 (Differential Evolution) という手法が使われている論文を読む機会がありましたので，差分進化法の計算についてまとめ，Juliaで実装を行いました．ここでは関数 f: \mathbb{R}^m \to \mathbb{R} の最小化問題を考えることにします．以下の例で利用するAckley関数は2変数関数なので，m=2です．

差分進化法

概要

参考URLに記したWikipediaやQiitaの記事の中では，どの記事も分かりやすいです．概ね，次のような計算を繰り返すことで，関数値をできるだけ小さくするような入力を求めたいです．

  1. n個のエージェント \mathbf{x}_i を用意する．
  2. 以下を繰り返す: 各エージェント \mathbf{x} について以下を繰り返す．
    1. 3つのエージェント \mathbf{a,b,c} を用意する (全て互いに異なるとする)．
    2. 新しいエージェント \mathbf{z} = \mathbf{a} + \gamma (\mathbf{b} - \mathbf{c}) を作成する．
    3. 最適化する関数の入力次元のうち，1つをランダムに選択し j\in \{1,\dots,m\} とする．
    4. エージェント \mathbf{x} の新世代 \mathbf{x}' を確率 p の変異で作成する
      • \mathbf{x}'_i = \begin{cases}z_i &\text{if }i = j\text{ or prob. }p\\ x_i&\text{otherwise}\end{cases}
    5. 仮に f(\mathbf{x'}) < f(\mathbf{x}) であるとき，\mathbf{x'} で \mathbf{x} を置き換える

例題: Ackley関数

次の関数です．

f(\mathbf{x}) = -a \exp\left(-b\sqrt{\frac{1}{d}\sum_{i=1}^n x_i^2}\right) - \exp\left(\frac{1}{d}\sum_{i=1}^d \cos(cx_i)\right) + a + \exp(1)

ただしここでは a=20, b=0.2, c=2\pi とします．contourを描画したものを次に例で示します．3次元のポヨポヨ図を見たい場合はWikipediaを開いてください．最適解は原点です．

ちなみにこちらのcontour plotをPlots.jlで描画するスクリプトは以下のものです．

function ackley(x; a=20, b=0.2, c=2π)
    d = length(x)
    -a * exp(-b * sqrt(sum(x.^2) / d)) - exp(sum(cos.(c*xi) for xi in x) / d) + a + exp(1)
end
ackley2(x, y) = ackley([x, y])

K = 5
δK = 0.01
xseq = -K:δK:K;
yseq = -K:δK:K;
X = repeat(reshape(xseq, 1, :), length(yseq), 1)
Y = repeat(yseq, 1, length(xseq))
Z = map(ackley2, X, Y);
cs = ColorSchemes.southwest
plot_c = cgrad(cs)
c = contour(xseq, yseq, Z, aspect_ratio=:equal, color=plot_c,
            xlabel="x1", ylabel="x2", size=(400, 400),
            fill=true, xlim=(-4, 4), ylim=(-4, 4))
savefig(c, "ackley.png")

差分進化法の適用結果

作成する個体数を n=20 とし，ランダムで初期化したエージェントを8世代分作成した結果を並べたものを示します．原点の方向に向かっていっている様子が見て取れますね．

実装

以下に実装を説明します（とはいえ擬似コードは上に書いた日本語の説明をそのまま実装するだけですが…）．

エージェント

各エージェントの位置をどのように保存するか，にいろいろと個性が出そうですが，今回はAckley関数だけを想定して書いたので，エージェント n 個分をまとめて行列として持つことにしました．初期化は世の中の実装例からパクってきたものを使いました．

struct Population
    array::Array{Float64, 2}
end
function Population(;population_size=20)
    """
    ランダムな初期化
    """
    array = (rand(population_size, 2) .- 0.5) .* 9
    Population(array)
end

呼び出すと次のように初期化されます．各行がエージェントに相当するため，この行列のeachrowでiterationして使います．

p = Population()
p.array

> 20×2 Array{Float64,2}:
 -3.35705   -0.350142
 -3.8098    -2.66255
 -4.28793   -2.12509
 -0.727233  -4.15515
  1.53554   -0.43424
  1.66609   -4.28848
 -4.14815   -2.34891
  3.09688    3.7733
 -1.2779    -3.74454
  2.48023   -2.20534
 -3.55172   -0.328718
 -3.29623   -2.36474
  3.87191    3.2092
 -2.93712   -0.401041
 -0.81925   -0.0456878
 -1.58764    3.98948
  3.78726    0.836414
  4.39724    1.61375
  2.29069   -0.25602
  4.0606    -1.38312

差分進化法の実装

本体の実装です．差分進化法において各エージェントとは異なる別の3エージェント \mathbf{a, b, c} をサンプリングしてくる必要がありますが，ここではStatsBase.jlからsampleを借りてきて実装しています．do-whileみたいな仕組みで順番に決めてもいいです．StatsBase.jl#sampleの重みweightに，自分自身のインデクスを除いてreplace=falseの形でサンプルを3つ取ってくれば良いです．新しいエージェントの集団 (population，この単語は危ないですね？) を作成する部分については，擬似コードというか日本語の説明の通りです．

using StatsBase
function differential_evolution(f, population; P=0.5, w=0.2)
    n, m = size(population.array)
    array = population.array
    new_array = copy(array)
    
    # 各個体について
    for (id, x) in enumerate(eachrow(array))
        # Obtain three distinct individuals a, b, c (all of x, a, b, c are distinct)
        weight = Weights([j != id for j in 1:n])
        ia, ib, ic = sample(1:n, weight, 3, replace=false)
        a, b, c = array[ia, :], array[ib, :], array[ic, :]
        
        z = a + w * (b - c) # compute a new position z
        j = rand(1:m) # choose a random dimension j in [1, m]
        xnew = zeros(m) # new agent
        for i in 1:m
            if i == j || rand() < P
                xnew[i] = z[i]
            else
                xnew[i] = x[i]
            end
        end
        
        # replace
        if f(xnew) < f(x)
            new_array[id, :] = xnew
        end
    end
    
    Population(new_array)
end

関数の適用と描画

最後に作成した関数を使って最初の適用例を描画する部分です．

# 初期状態 + 7個で8世代分
history = Population[p]
for i in 1:7
    p = differential_evolution(ackley, p)
    push!(history, p)
end

# 8個分保存
list_s = []
for i in 1:length(history)
    pi = history[i]
    c = contour(xseq, yseq, Z, aspect_ratio=:equal, color=plot_c, size=(400, 400),
                label=nothing, fill=true, xlim=(-5, 5), ylim=(-5, 5), colorbar=false)
    s = scatter!(c, pi.array[:, 1], pi.array[:, 2], size=(200, 200), label=nothing,
                 xlim=(-5, 5), ylim=(-5, 5))
    push!(list_s, c)
end

# 2行4列のレイアウトで出力する
p = plot(list_s..., layout=(2, 4), size=(800, 400))
savefig(p, "de.png")

こうして出力されたde.pngが上に示した適用例になります．

参考URL

Discussion