Differential Evolutionで大域的最適化

はじめに

Differential Evolutionという手法を使うと、少ないパラメーターで関数の最小化を行うことができます。いわゆる遺伝的アルゴリズムと呼ばれるものです。

アルゴリズム

例としてAckley Functionと呼ばれる関数の最小点を探索してみます。
この関数は(x1, x2)=(0,0)で最小値0を持つので、最小値探索の結果が(0,0)となれば正しく動いていることになります。
Ackley Function

double Ackley(double x1, double x2)
{
  const double n = 2.0;
  const double c = cos(2.0*M_PI*x1) + cos(2.0*M_PI*x2);
  return 20.0 + exp(1.0) - 20.0*exp(-0.2*sqrt(1.0/n*(x1*x1 + x2*x2))) 
                    - exp(1.0/n*c);
}

元ネタ

[1] A comparative study of differential evolution variants for global optimization, 2006.
[2]http://www.cybernet.co.jp/mds/solution/sol3.html

初期状態

パラメーターが存在しそうな範囲[min,max]に対してランダムに点を生成します。

// // [min, max]でランダムに点を選ぶ
// MatrixXd x(N_, dim_);
// Initialization(x);

void DifferentialEvolution::Initialization(MatrixXd& x)
{
  std::uniform_real_distribution<double> dist(min_, max_);

  x.noalias() = x.unaryExpr([&](double _)  { return dist(rnd_); });

  std::cout << x << std::endl;
  // rowは生成したデータ点、colはパラメーター
  // -0.82978 -4.07661
  //  4.97185 -1.03419
  // ...
}

ランダムに３点選んで１点つくる

さっき生成した点から3点(r1,r2,r3)を選びます。遺伝的アルゴリズムなので、これが親の役割となり、この３点から新たな点を生成します。

Fは定数で0.5にしています。
これ以外にも、もっと多くの点を使ったり前回よかった点を使うなどの種類があります。

点を混ぜる

３点から作った点vと現在注目している点xを混ぜます。

CRは定数0.5で、この場合だと50%の確率で新しい点を採用することになります。jrandはパラメーターのインデックスjの乱数で、新しい点uがxとまるっきり同じになることを防ぐ値です。

採用の決定

新しく作った点uがより悪い点、つまり関数値を増やす場合には採用しません。

以上の処理を繰り返すと収束します。

ソースコード

// [0,n)の数列を混ぜる
// n=4で1,0,2,3とか
std::vector<int> DifferentialEvolution::RandPerm(int n)
{
  std::vector<int> r(n);
  std::iota(begin(r), end(r), 0);
  std::shuffle(begin(r), end(r), rnd_);

  return r;
}

void DifferentialEvolution::Mutation(MatrixXd& x)
{
  std::uniform_int_distribution<int> distInt(0, dim_);
  std::uniform_real_distribution<double> distReal(0.0, 1.0);

  for (int i = 0; i < N_; i++) // データ点
  {
    // 3点選ぶ
    std::vector<int> r(RandPerm(N_));
    VectorXd r1 = x.row(r[0]);
    VectorXd r2 = x.row(r[1]);
    VectorXd r3 = x.row(r[2]);

    int jrand = distInt(rnd_);
    VectorXd ui = VectorXd::Zero(dim_);

    for (int j = 0; j < dim_; j++) // パラメータ次元
    {
	double uj = 0.0;
	if (distReal(rnd_) <= CR_ || j == jrand)
	{
	  VectorXd vi = r1 + F_ * (r2 - r3);
	  uj = vi(j);
	}
	else
	{
	  uj = x(i, j);
	}

	ui(j) = uj;
        // 追記：こっちの方がよかった
        //ui(j) = (distReal(rnd_) <= CR_ || j == jrand) 
	//          ? r1(j) + F_ * (r2(j) - r3(j)) 
	//          : x(i, j);
    }

    if (Ackley(ui(0), ui(1)) <=  Ackley(x(i, 0), x(i, 1)))
    {
      x.row(i) = ui;
    }
  }
}

曲線フィッティング

具体例として点が与えられたときに、曲線をフィッティングする問題を考えてみます。だいぶ前にレーベンバーグ・マーカート(Levenberg-Marquardt)法の記事で使った誤差関数を使ってみます。wildpie.hatenablog.com
観測された14点から下式のパラメーターb1, b2を求めることが目的となります。

// 観測データ
double xa[] = {77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0};
double ya[] = {10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0};

// 誤差関数
double Misra1a(double b1, double b2) 
{
  double sum = 0.0;
  for (int i = 0; i < 14; i++)
  {
    double temp = b1*(1.0 - exp(b2*xa[i])) - ya[i];
    sum += temp*temp;
  }
  return sum;
}

これを実行すると数百回目の試行で238.942, -0.000550157のパラメーターを得ることができました。
微分も使わず力わざで最小値を求めるのも悪くないですね。