Go to Qiita Advent Calendar 2024 Top
LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@rustcat

勾配法を可視化してみた

Last updated at Posted at 2024-08-29

前回は「最急降下法（勾配降下法）」を可視化してみたので、今回は「これなら分かる最適化数学」という本で紹介されていた勾配法を可視化してみました。
アルゴリズムは下記サイトで詳しく解説されていたので、参考にさせていただきました。
前回の勾配降下法よりは収束が早い感じです。
なおグラフは都合により前回と比べて上下逆になってしまっていますが、基本的に同じものです。

4本のプログラムは関数部分以外は全く同じなので、2本目からは関数部分のみを掲載しました。

間違いなどありましたら、指摘していただけるとありがたいです。

参考にさせていただいたサイト

参考にさせていただいた本

  • これなら分かる最適化数学　基礎原理から計算手法まで

参考にさせていただいたグラフ表示に関するサイト（実際に参考にしたのは「Visual C++の易しい使い方（23） ―三次元グラフの陰線処理―」の方ですが、今は公開されていないようです）

図1～図4に各関数の探索の様子をグラフに表示しました。

図1 $-5(x-1)^2-2(y-2)^2$
GradientMethod0003.gif

図2 $-x^2-y^2-20\sin^2 x$
GradientMethod0004.gif

図3 $-x^3-y^3+9xy-27$
GradientMethod0005.gif

図4 $-0.5x^4-0.1y^4+0.01x^3y^3+x+0.3y$
GradientMethod0006.gif

GradientMethod1.cpp
#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>
#include <fstream>
using namespace std;

// 0以上1以下の実数乱数
#define RAND_01 ((double)rand() / RAND_MAX)

const int ITERATION = 100; // 試行回数

const int SWARM_SIZE = 20; // 群の大きさ

double X_MAX = 10; // graphのX軸の表示範囲の上限
double X_MIN = -10; // graphのX軸の表示範囲の下限

double Y_MAX = 10; // graphのY軸の表示範囲の上限
double Y_MIN = -10; // graphのY軸の表示範囲の下限

double Z_MAX = 200; //  graphのZ軸の表示範囲の上限
double Z_MIN = 0; // graphのZ軸の表示範囲の下限

string str_formula = "x**2 + y**2"; // グラフに表示する数式

double x[ITERATION][SWARM_SIZE];
double y[ITERATION][SWARM_SIZE];
double z[ITERATION][SWARM_SIZE];

double function(double x, double y)
{
	double z;

	z = -5 * (x - 1) * (x - 1) - 2 * (y - 2) * (y - 2);

	return z;
}

void set_function_conditions()
{
	X_MAX = 4;
	X_MIN = -6;
	Y_MAX = 4;
	Y_MIN = -6;
	Z_MAX = 0;
	Z_MIN = -400;
	str_formula = "-5 * (x - 1)**2 - 2 * (y - 2)**2";
}

void write_data()
{
	ofstream fout;
	fout.open("graph_data.txt");

	fout << str_formula << endl;
	cout << str_formula << endl ;

	fout << "ITERATION " << ITERATION << endl;
	cout << "ITERATION " << ITERATION << endl;

	fout << "SWARM_SIZE " << SWARM_SIZE << endl;
	cout << "SWARM_SIZE " << SWARM_SIZE << endl;

	fout << "X_MAX " << X_MAX << endl;
	cout << "X_MAX " << X_MAX << endl;

	fout << "X_MIN " << X_MIN << endl;
	cout << "X_MIN " << X_MIN << endl;

	fout << "Y_MAX " << Y_MAX << endl;
	cout << "Y_MAX " << Y_MAX << endl;

	fout << "Y_MIN " << Y_MIN << endl;
	cout << "Y_MIN " << Y_MIN << endl;

	fout << "Z_MAX " << Z_MAX << endl;
	cout << "Z_MAX " << Z_MAX << endl;

	fout << "Z_MIN " << Z_MIN << endl;
	cout << "Z_MIN " << Z_MIN << endl;

	fout << "DATA" << endl;
	cout << "DATA" << endl;

	for (int i = 0; i < ITERATION; i++) {
		cout << i << endl;
		for (int j = 0; j < SWARM_SIZE; j++) {
			cout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
			fout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
		}
		cout << endl;
		fout << endl;
	}

	fout.close();
}

void initialize()
{
	for (int i = 0; i < SWARM_SIZE; i++) {
		x[0][i] = RAND_01 * (X_MAX - X_MIN) + X_MIN;
		y[0][i] = RAND_01 * (Y_MAX - Y_MIN) + Y_MIN;
		z[0][i] = function(x[0][i], y[0][i]);
	}
}
	
double __delta_h = 0.000001;

double diff_x_func(double x1, double y1) // x軸方向への偏微分
{
	double f1 = function(x1 + __delta_h, y1);
	double f2 = function(x1 - __delta_h, y1);
	double f3 = (f1 - f2) / (2 * __delta_h);
	return f3;
}

double diff_y_func(double x1, double y1) // y軸方向への偏微分
{
	double f1 = function(x1, y1 + __delta_h);
	double f2 = function(x1, y1 - __delta_h);
	double f3 = (f1 - f2) / (2 * __delta_h);
	return f3;
}

// 関数 F(t) の定義
double func1(double x1, double y1, double t)
{
	double f = function(x1 + t * diff_x_func(x1, y1), y1 + t * diff_y_func(x1, y1));
	return f;
}

// 関数の微分を中心差分で求める
double diff_func(double x1, double y1, double t)
{
	double delta_h = 0.000001;
	double f = (func1(x1, y1, t + delta_h) - func1(x1, y1, t - delta_h)) / (2 * delta_h);
	return f;
}

double get_sign(double s)
{
	if (s >= 0.0)
		return 1.0;
	else
		return -1.0;
}

// 1変数の勾配法
double search_max(double x1, double y1, double initial_t, double initial_h, double epsiron)
{
	// step 1 : 初期値の設定
	double t = initial_t;
	double h = initial_h;
	
	while (fabs(diff_func(x1, y1, t)) > epsiron) { // step 5 : 収束のチェック
		// step 2
		h = get_sign(diff_func(x1, y1, t)) * fabs(h);
		double T = t;
		double Td = t + h;

		double fT = func1(x1, y1, T);
		if (fT < func1(x1, y1, Td)) { // step 3
			while (fT <= func1(x1, y1, Td)) { // step 3 (a)
				h = 2 * h;
				T = Td;
				Td = T + h;
				fT = func1(x1, y1, T);
				if (fT > Z_MAX)
					break;
			}
			
			// step 3 (b)
			t = T;
			h = h / 2;
		} else { // step 4
			while (fT >= func1(x1, y1, Td)) { // step 4 (a)
				h = h / 2;
				Td = Td - h;
				fT = func1(x1, y1, T);
				if ((fT > Z_MAX) || (h == 0.0))
					break;
			}
			// step 4 (b)
			t = Td;
			h = 2 * h;
		}
		if ((fT > Z_MAX) || (h == 0.0))
			break;
	}

	return t;
}

double value_of_x(double x1, double y1)
{
	double f = sqrt(x1 * x1 + y1 * y1);
	return f;
}

void hill_climb_max()
{
	double delta = 0.000001; //#収束条件

	for (int i = 1; i < ITERATION; i++) {
		for (int j = 0; j < SWARM_SIZE; j++) {
			double x1 = x[i-1][j];
			double y1 = y[i-1][j];
			// 最大値 t の探索
			double initial_t = 0.1; //## 初期値は常に0.1から始める
			double initial_h = 0.1; //## 初期のステップ幅
			double epsiron = 0.00001; //## 収束判定
			double max_t = search_max(x1, y1, initial_t, initial_h, epsiron);
			// delta_xの計算
			double delta_x_x = max_t * diff_x_func(x1, y1);
			double delta_x_y = max_t * diff_y_func(x1, y1);
			// xの更新
			x1 = x1 + delta_x_x;
			if (x1 > X_MAX)
				x1 = X_MAX;
			else if (x1 < X_MIN)
				x1 = X_MIN;

			y1 = y1 + delta_x_y;
			if (y1 > Y_MAX)
				y1 = Y_MAX;
			else if (y1 < Y_MIN)
				y1 = Y_MIN;

			x[i][j] = x1;
			y[i][j] = y1;
			z[i][j] = function(x1, y1);

//			if (value_of_x(delta_x_x, delta_x_y) < sigma)
//				break;
		}
	}
}

int main()
{
	srand ((unsigned)time(NULL));

	set_function_conditions();

	initialize();

	// 最大値探索
	hill_climb_max();

	write_data();
}
GradientMethod2.cpp
double function(double x, double y)
{
	double z;

	z = -(x * x + y * y + 20.0 * sin(x) * sin(x));

	return z;
}

void set_function_conditions()
{
	X_MAX = 10;
	X_MIN = -10;
	Y_MAX = 10;
	Y_MIN = -10;
	Z_MAX = 0;
	Z_MIN = -200;
	str_formula = "-(x**2 + y**2 + 20.0 * sin(x)**2)";
}
GradientMethod3.cpp
double function(double x, double y)
{
	double z;

	z = -(x * x * x + y * y * y - 9.0 * x * y + 27.0);

	return z;
}

void set_function_conditions()
{
	X_MAX = 10;
	X_MIN = -10;
	Y_MAX = 10;
	Y_MIN = -10;
	Z_MAX = 3000;
	Z_MIN = -1500;
	str_formula = "-(x**3 + y**3 - 9.0 * x * y + 27.0)";
}
GradientMethod4.cpp
double function(double x, double y)
{
	double z;

	z = -0.5 * x * x * x * x - 0.1 * y * y * y * y + 0.01 * x * x * x * y * y * y + x + 0.3 * y;

	return z;
}

void set_function_conditions()
{
	X_MAX = 6;
	X_MIN = -6;
	Y_MAX = 6;
	Y_MIN = -6;
	Z_MAX = 0;
	Z_MIN = -1000;
	str_formula = "-0.5 * x**4 - 0.1 * y**4 + 0.01 * x**3 * y**3 + x + 0.3 * y";
}
0
0
0

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
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?