焼きなまし法で遊ぶ - Negative/Positive Thinking

はじめに

焼きなまし法について、どんな挙動で、どこを気をつけないといけないのかよくわかってないので、遊んでみる。

焼きなまし法メモ : http://d.hatena.ne.jp/jetbead/20111014/1318598381

評価関数

簡単のために、4次関数を使って、初期位置は小さい山の外側(x0=2.0)に配置してみる。
$f(x)=-3x^4+3x^3+10x^2-10x+15$

【関数の形】
(google)

(wolfram alpha)
http://www.wolframalpha.com/input/?i=-3x%5E4%2B3x%5E3%2B10x%5E2-10x%2B15

best_x = argmax f(x)を求めたい。

コード

#include <iostream>
#include <cmath>

//xorshift
// 注意: longではなくint(32bit)にすべき
unsigned long xor128(){
  static unsigned long x=123456789, y=362436069, z=521288629, w=88675123;
  unsigned long t;
  t=(x^(x<<11));
  x=y; y=z; z=w;
  return w=(w^(w>>19))^(t^(t>>8));
}
//[0,1)の一様乱数
// 注意: int_maxぐらいで割るべき
double frand(){
  return xor128()%ULONG_MAX/static_cast<double>(ULONG_MAX); 
}

//焼きなまし
class SA {
  double x; //探索中の解
  double z; //暫定解
  double best; //暫定解の評価値
  double T; //温度
  int R; //反復回数

  //評価関数
  //f(x) = -3x^4+3x^3+10x^2-10x+15
  //http://www.wolframalpha.com/input/?i=-3x%5E4%2B3x%5E3%2B10x%5E2-10x%2B15
  double f(double x){
    return -3*x*x*x*x + 3*x*x*x + 10*x*x -10*x + 15;
  }

  //近傍からランダムに選ぶ
  double get_random_neighbor(double xx){
    //左右どちらかにちょっとだけ動かしたものを近傍N(x)としてみる
    if(frand()<0.5) return xx+0.01;
    return xx-0.01;
  }

  //温度の更新
  double next_T(double t){
    return t*0.995;
  }

  //反復回数の更新
  double next_R(double r){
    return r;
  }

public:
  //xの初期値、温度の初期値、反復回数の初期値
  SA(double init_x, double t, int r){
    x = init_x;
    z = init_x;
    best = f(z);
    T = t;
    R = r;
  }

  //探索
  void calc(){

    while(T>1.0){//十分冷えるまで
      for(int r=0; r<R; r++){
	//xの近傍からランダムに選ぶ
	double y = get_random_neighbor(x);
	//変化量
	double delta = f(x) - f(y);
	
	if(delta<0.0){ //よりよい解が見つかった場合
	  x = y;
	}else{ //悪くなる解の場合
	  if(exp(-delta/T) > frand()){ //ある程度の確率で探索を許す
	    x = y;
	  }
	}
	
	//最適解の更新
	if(f(x) > best){
	  best = f(x);
	  z = x;
	}

        //探索を出力	
	std::cout << x << std::endl;

      }
      
      std::cout << "#" << z << " " << best << std::endl;

      //温度の更新
      T = next_T(T);
      //反復回数の更新
      R = next_R(R);   
    }

  }

  //最適解
  double get_best(){ return z; }

};

int main(){
  SA sa(2.0, 100, 1000);

  sa.calc();

  std::cout << sa.get_best() << std::endl;

  return 0;
}