こんにちは。
今回は「差分進化」を扱います。
時間が取れそうだったので何かつくろうかと思っていたところ、Nature of Code の遺伝的アルゴリズム関連はまだやっていなかったなと。
だいぶ前ですがエッチな画像をつくるみたいな企画もあって、「遺伝的アルゴリズム」という名前だけでも知っている方がいるかもしれません。
ただ Nature of Code は割と分量が多く、まずはあっさりできそうなものからと思って関連する Wikipedia のページを眺めていました。
そこで「差分進化」のページを見つけまして、最適化問題云々と書いてあったので応用 (?) もやりやすそうだと思い取り組んでみることにしました。
目次
差分進化法とは
Wikipedia には、
与えられた評価尺度に関する候補解を反復的に改良していき、問題を最適化する手法である。
とあります。
最急降下法のような最適化手法と似た概念ではありますが、最急降下法ではある関数から漸化式をつくります。
しかし差分進化法ではそのようなものは無く、適当な初期値からある操作によって動かしてそれがより目標に近ければ採用、というようなやり方です。
また、最急降下法の漸化式では考えている関数の微分形が必要ですが、差分進化法では必要ないので、数式で表現できなかったり微分できなかったりしても最適化できるということになります。
ざっくり、最急降下法はガチガチの数式から答えを出し、差分進化法ではシミュレーション的に答えを求めていく、と表現できると思います。
アルゴリズム
詳しい説明は調べればいくらでも出てきますが、それほど難しいものでもなさそうなのでここにも書いておきます。
ここでは、高低差のある地形のようなものに対してその中で一番低いところを探すことを考えてみます。
(下の図で一番内側の楕円の中心が最低地点ということにします。)
初期化
まずは初期化を行います。
これから考える候補解をランダムに生成します。
突然変異・交叉
続いて、全ての候補解について子を決めます。
候補解 $\boldsymbol{x}_{i}$ を考えるとき、他の候補解からランダムに3つ $\boldsymbol{x}_{r1}$、$\boldsymbol{x}_{r2}$、$\boldsymbol{x}_{r3}$ を選びます。
(太字になっている文字はベクトルであることを表しています。)
選んだ3つから変異ベクトル $\boldsymbol{v}$ をつくります。
$$ \boldsymbol{v} = \boldsymbol{x}_{r1} + F ( \boldsymbol{x}_{r2} - \boldsymbol{x}_{r3} ) $$
ここで、$F$ は適当な定数です。
例えば、$\boldsymbol{x}_{3}$ に注目して $\boldsymbol{x}_{4}$、$\boldsymbol{x}_{2}$、$\boldsymbol{x}_{1}$ を選んだとすると、
この変異ベクトル $\boldsymbol{v}$ と今考えている候補解 $\boldsymbol{x}_{i}$ から、どのように変化させて子とするかを決めます。
決め方としては、候補解 $\boldsymbol{x}$ の各成分 $x_{j}$ について、ある確率 (交叉率) $P$ で変異ベクトル $\boldsymbol{v}$ の成分 $v_{j}$ を引き継ぐ、という感じです。
それとは別に、候補解の次元 $D$ を最大値とする乱数 $j_{\mathrm{rand}}$ を選んでおき、成分 $v_{j_{\mathrm{rand}}}$ は必ず引き継ぐようにします。
子を $\boldsymbol{u}$ とすると、その成分 $u_{j}$ はこんな感じに表せそうです。
$$ u_{j} = \begin{cases} v_{j} & \mathrm{if} ~~ \mathrm{rand} \le P ~~ \mathrm{or} ~~ j = j_{\mathrm{rand}} \\ x_{j} & \mathrm{otherwise} \end{cases} $$
今は二次元を考えているので $D$ は 2 で、$j_{\mathrm{rand}}$ は 1 または 2 となります。
例えば先程の図の状態を考えるとすると、$\boldsymbol{x}_{3}$ に注目しているので、このベクトルの成分がどちらかは必ず $\boldsymbol{v}$ のものを引き継ぎ、もう一方は一定の確率で引き継ぎます。
つまり、$\boldsymbol{x}_{3}$ の子は、$(v, x)$、$(x, v)$、$(v, v)$ の3つのいずれかになる可能性があるということになります。
(上の図の赤丸と橙丸がそれです。)
淘汰
最後に、次の世代を確定させます。
より目標に近い種に進化させていきたいので、何か評価方法を用意しておき $\boldsymbol{x}$ と $\boldsymbol{u}$ で評価がより高い方を選びます。
今考えている例では、楕円の中心により近い方を残すことになります。
もし $\boldsymbol{u}$ が下の橙丸になったら、$\boldsymbol{x}_{3}$ と $\boldsymbol{u}$ で評価を比較し、より楕円に近い $\boldsymbol{x}_{3}$ が次の世代に選ばれます。
同様にして、全ての候補解に対して終了条件を満たすまで操作を繰り返します。
終了条件は何世代まで計算するかとか目標にどのくらい近づいたらとかそんな感じです。
最終的には中心に近づいていくようになり、この中で一番低い位置にいるものが最適解ということになります。
注意点 他
これで差分進化法による最適化のアルゴリズムが説明できましたが、これが必ずしも最も低い、真の解になっているとは限らないことに注意が必要です。
若干ズレているかもしれないというのもそうですが、他にくぼみがあってそちらに集まってしまった場合もありえるからです。
この場合はより浅いくぼみに集まったことになるので、真の解とはほど遠い結果となってしまいます。
これを避けるためいろいろな初期値から最適化する、という考え方もあります。
(最も深いくぼみは global minimum、より浅いくぼみは local minimum とよばれたりします。)
またこれは実装時の注意点かもしれませんが、$\boldsymbol{u}$ の計算時にもし枠からはみ出すようだったら、
$$ u_{j} = \begin{cases} x_{r1, j} + \mathrm{rand} ( \underline{x}_{j} - x_{r1, j} ) & \mathrm{if} ~ u_{j} \lt \underline{x}_{j} \\ x_{r1, j} + \mathrm{rand} ( \overline{x}_{j} - x_{r1, j} ) & \mathrm{if} ~ u_{j} \gt \overline{x}_{j} \end{cases} $$
のようにして枠内に収まるように修正すればいいようです。
ここで、$\underline{x}_{j}$ と $\overline{x}_{j}$ はそれぞれ $j$ 番目の成分の最小値と最大値を表しています。
実装
理論的な部分は説明できましたので、実際に processing で動くようにしてみます。
sphere関数とよばれている、
$$ f(\boldsymbol{x}) = \sum _{i = 1} ^{n} x_i ^2 $$
二次関数の三次元バージョンみたいなものを考えることにします。
ここでは、この図の下にある平面のような、高さを白黒で表した図の上で候補解を動かしていき、最も低い地点を探してみることにします。
白くなるにしたがって低くなるように描画しているので、最も明るいところを目指すという感じになります。
まずはもろもろの変数を設定します。
// 候補解の数 n = 20; // 探索範囲 min = -3; max = 3; // 1マスの大きさ len = 5; // 1辺の点数 points = int(width/len); // 初期値 for(int i = 0; i < n; i ++){ agents.add(new PVector(random(0, width), random(0, width))); }
ここで、「探索範囲」はグラフの話をしています。つまり $x$ は $[-3, 3]$、$y$ は $[-3,3]$ を考える、という意味です。
また、候補解を入れた配列 agents はベクトルの ArrayList として宣言しておき、後で出し入れしやすいようにしておきます。
盤面と候補解の描画は、
// 盤面の描画 for(int i = 0; i < points; i ++){ for(int j = 0; j < points; j ++){ float c = map(evaluate(new PVector(i*len, j*len)), 0, 5, 255, 0); fill(c); rect(i*len, j*len, len, len); } } // 候補解の描画 for(int i = 0; i < agents.size(); i ++){ fill(0); ellipse(agents.get(i).x, agents.get(i).y, 5, 5); }
塗る色 c の設定において、考えている関数によって第2,3変数は変わると思います。
evaluate() は評価関数で、単純に agents の座標を代入しているだけです。
// 評価用の関数 float evaluate(PVector v){ float x = Xtox(v.x); float y = Xtox(v.y); float f = sqrt(x*x + y*y); return f; }
ここで、$(X, Y)$ を描画用の座標系、$(x, y)$ を計算用の座標系としています。
agents の座標は描画用の座標 (0 ~ width の間のランダム) で決めているので、これを計算用の座標系に変換 (Xtox() を適用) してから評価しています。
// 描画用の (X, Y) から計算用の (x, y) への変換 float Xtox(float v){ return map(v, 0, width, min, max); }
そして、肝である最適化の処理が続きます。
// 最適化 for(int j = 0; j < agents.size(); j ++){ // ランダムで被らないように配列をつくる ArrayList<PVector> list = new ArrayList<PVector>(); for(int i = 0; i < agents.size(); i ++){ list.add(agents.get(i)); } // 次の世代を決めたいもの PVector agent = list.get(j); list.remove(j); // 計算に使う3つをランダムに選ぶ PVector[] r = new PVector[3]; for(int i = 0; i < 3; i ++){ int a = floor(random(0, list.size())); r[i] = list.get(a); list.remove(a); } // 変異ベクトルの計算 float F = 0.5; // 係数 PVector v = new PVector(r[0].x + F * (r[1].x - r[2].x), r[0].y + F * (r[1].y - r[2].y)); //// 以下、子を決める処理 //// PVector u = new PVector(agent.x, agent.y); // 必ず引き継ぐもの int j_rand = int(random(1, 3)); // 二次元のため 1 or 2 if(j_rand == 1){ u.x = v.x; } // 1 ならx成分を引き継ぐ if(j_rand == 2){ u.y = v.y; } // 2 ならy成分を引き継ぐ // 交叉率に基づいて引き継ぐ float Cr = 0.1; // 交叉率 float rand = random(0, 1); if(random(0, 1) <= Cr){ u.x = v.x; } if(random(0, 1) <= Cr){ u.y = v.y; } // 枠外に出たら修正 if(u.x < 0) { u.x = r[0].x + random(0, 1) * (0 - r[0].x); } if(u.x > width){ u.x = r[0].x + random(0, 1) * (width - r[0].x); } if(u.y < 0) { u.y = r[0].y + random(0, 1) * (0 - r[0].y); } if(u.y > width){ u.y = r[0].y + random(0, 1) * (width - r[0].y); } //// 子を決める処理はここまで //// // 評価が高い方を次の世代に if(evaluate(agent) >= evaluate(u)){ agents.set(j, u); } }
流れを説明すると…
まず agents のコピー list をつくっておきます。
(そこから1つ選んでそれを消去、を繰り返すことで被りなく選べるようにしています。)
次の世代を決めたい候補解を選んだ後、変位ベクトルの計算用に3つ選びます。
そして、変位ベクトルを計算します。
// ランダムで被らないように配列をつくる ArrayList<PVector> list = new ArrayList<PVector>(); for(int i = 0; i < agents.size(); i ++){ list.add(agents.get(i)); } // 次の世代を決めたいもの PVector agent = list.get(j); list.remove(j); // 計算に使う3つをランダムに選ぶ PVector[] r = new PVector[3]; for(int i = 0; i < 3; i ++){ int a = floor(random(0, list.size())); r[i] = list.get(a); list.remove(a); } // 変異ベクトルの計算 float F = 0.5; // 係数 PVector v = new PVector(r[0].x + F * (r[1].x - r[2].x), r[0].y + F * (r[1].y - r[2].y));
続いて、子を決めます。
必ず引き継ぐ要素と確率 (交叉率) で引き継ぐものがあるので、それぞれ決めて、枠外に出たら修正しておしまいです。
PVector u = new PVector(agent.x, agent.y); // 必ず引き継ぐもの int j_rand = int(random(1, 3)); // 二次元のため 1 or 2 if(j_rand == 1){ u.x = v.x; } // 1 ならx成分を引き継ぐ if(j_rand == 2){ u.y = v.y; } // 2 ならy成分を引き継ぐ // 交叉率に基づいて引き継ぐ float Cr = 0.1; // 交叉率 if(random(0, 1) <= Cr){ u.x = v.x; } if(random(0, 1) <= Cr){ u.y = v.y; } // 枠外に出たら修正 if(u.x < 0) { u.x = r[0].x + random(0, 1) * (0 - r[0].x); } if(u.x > width){ u.x = r[0].x + random(0, 1) * (width - r[0].x); } if(u.y < 0) { u.y = r[0].y + random(0, 1) * (0 - r[0].y); } if(u.y > width){ u.y = r[0].y + random(0, 1) * (width - r[0].y); }
最後に、次の世代に親を残すか子を残すか決めます。
evaluate() を使い、結果がより低い方に置き換えます。
// 評価が高い方を次の世代に if(evaluate(agent) >= evaluate(u)){ agents.set(j, u); }
ちなみに、候補解の数は $10D$、係数 $F$ は $0.5$、交叉率 $CR$ は $0.1$ とするのがいいようです。
($D$ は次元数を表しており、今回は 2 になります。)
今のところ終了条件は無しで、結果を出力しておいて収束していそうだったら止める、という感じにしています。
世代数を表示するようにしても良いかもしれません。
以上で差分進化法の実装ができました。
ぶつ切りだったので、全体をここに載せておきます。
int n; ArrayList<PVector> agents = new ArrayList<PVector>(); int generation = 0; int points; float len; float max, min; void setup(){ size(500, 500); noStroke(); // 候補解の数 n = 20; // 探索範囲 min = -3; max = 3; // 1マスの大きさ len = 5; // 1辺の点数 points = int(width/len); // 初期値 for(int i = 0; i < n; i ++){ agents.add(new PVector(random(0, width), random(0, width))); } } void draw(){ // 盤面の描画 for(int i = 0; i < points; i ++){ for(int j = 0; j < points; j ++){ float c = map(evaluate(new PVector(i*len, j*len)), 0, 5, 255, 0); fill(c); rect(i*len, j*len, len, len); } } // 候補解の描画 for(int i = 0; i < agents.size(); i ++){ fill(0); ellipse(agents.get(i).x, agents.get(i).y, 5, 5); } // 最適化 for(int j = 0; j < agents.size(); j ++){ // ランダムで被らないように配列をつくる ArrayList<PVector> list = new ArrayList<PVector>(); for(int i = 0; i < agents.size(); i ++){ list.add(agents.get(i)); } // 次の世代を決めたいもの PVector agent = list.get(j); list.remove(j); // 計算に使う3つをランダムに選ぶ PVector[] r = new PVector[3]; for(int i = 0; i < 3; i ++){ int a = floor(random(0, list.size())); r[i] = list.get(a); list.remove(a); } // 変異ベクトルの計算 float F = 0.5; // 係数 PVector v = new PVector(r[0].x + F * (r[1].x - r[2].x), r[0].y + F * (r[1].y - r[2].y)); //// 以下、子を決める処理 //// PVector u = new PVector(agent.x, agent.y); // 必ず引き継ぐもの int j_rand = int(random(1, 3)); // 二次元のため 1 or 2 if(j_rand == 1){ u.x = v.x; } // 1 ならx成分を引き継ぐ if(j_rand == 2){ u.y = v.y; } // 2 ならy成分を引き継ぐ // 交叉率に基づいて引き継ぐ float Cr = 0.1; // 交叉率 if(random(0, 1) <= Cr){ u.x = v.x; } if(random(0, 1) <= Cr){ u.y = v.y; } // 枠外に出たら修正 if(u.x < 0) { u.x = r[0].x + random(0, 1) * (0 - r[0].x); } if(u.x > width){ u.x = r[0].x + random(0, 1) * (width - r[0].x); } if(u.y < 0) { u.y = r[0].y + random(0, 1) * (0 - r[0].y); } if(u.y > width){ u.y = r[0].y + random(0, 1) * (width - r[0].y); } //// 子を決める処理はここまで //// // 評価が高い方を次の世代に if(evaluate(agent) >= evaluate(u)){ agents.set(j, u); } // 候補解の出力 println("generation =", generation); for(int i = 0; i < agents.size(); i ++){ println("[", Xtox(agents.get(i).x), Xtox(agents.get(i).y), "]"); } } generation += 1; } // 評価用の関数 float evaluate(PVector v){ float x = Xtox(v.x); float y = Xtox(v.y); float f = sqrt(x*x + y*y); return f; } // 描画用の (X, Y) から計算用の (x, y) への変換 float Xtox(float v){ return map(v, 0, width, min, max); }
実行してみると、
徐々に中央に集まっていく様子が見えます。
50世代くらいで全ての候補解がほとんど $(0, 0)$ になり、最低値をとる座標が $(0, 0)$ だとわかりました。
いろいろなベンチマーク関数
グラフに特徴的な構造を持った、最大値や最小値 (最適解) がわかっている関数をベンチマーク関数とよび、最適化アルゴリズムの評価に使われるようです。
検索するとたくさん出てきます。
(どれも割と複雑な関数なので実装が大変ですが。)
例えば…
three-hump関数
Easom関数
Michalewicz関数
なお、今回の実装では$y$軸が逆転したグラフになっているということに注意が必要です。
(気になる方は Ytoy() みたいな関数をつくって調節すると良いかと思います。)
まとめ
以上、processing に差分進化を実装してみました。
最適化のアルゴリズムと基本的な評価方法は触れましたが、実際の問題を解くのには使っていません。
ただの数式の最大や最小を求める分には良いのですが、例えば $x+y=10$ を満たしながら、みたいな条件が入ってくるとうまく書けなかったのでまた今度にしたいと思います。
最後まで読んでいただいてありがとうございました。また次回。
