Newton法を可視化してみた

「これなら分かる最適化数学」という本で紹介されていたニュートン法を可視化してみました。
ニュートン法は高速ですが、このプログラムの場合は極大値と極小値の両方に収束していますね。（図2の場合は、極小値というよりは鞍点というべきですね。図3も参考にさせていただいたサイトでは鞍点と言っていますね。鞍点に関連して、参考にさせていただいたサイトを一つ追加しました）
もう少し解説などできればいいんですが、私の理解が不足しすぎているので、詳しく知りたい方は参考にさせていただいたサイトや本を読んでいただくのがいいと思います。

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

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

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

• 【最適化問題の基礎】ニュートン法で停留点を探す（鞍点と固有値の関係）
• ニュートン法はどのように収束していくのか – 1変数の場合で視覚化
• 多変数のニュートン法の収束の様子を可視化
• C++ で連立一次方程式を解く(LU分解法)

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

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

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

• 隠線処理を用いた三次元グラフの作成

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

図1 $\quad-5(x-1)^2-2(y-2)^2$

図2 $\quad-x^2-y^2-20\sin^2 x$

図3 $\quad-x^3-y^3+9xy-27$

図4 $\quad-0.5x^4-0.1y^4+0.01x^3y^3+x+0.3y$

図5 $\quad-x^2+xy-y^2+x+y$

Newton1.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 = 40; // 試行回数

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();

#define DIM 2

double delta_h = 0.00001;

double d_func_dx(double x, double y)
{
	double f1 = function(x + delta_h, y);
	double f2 = function(x - delta_h, y);
	double f3 = (f1 - f2) / (2 * delta_h);
	return f3;
}

double d_func_dy(double x, double y)
{
	double f1 = function(x, y + delta_h);
	double f2 = function(x, y - delta_h);
	double f3 = (f1 -f2) / (2 * delta_h);
	return f3;
}

void calc_nabla_f(double nabla_f[], double x1, double y1)
{
	nabla_f[0] = d_func_dx(x1, y1);
	nabla_f[1] = d_func_dy(x1, y1);
}

double dd_func_dxdx(double x, double y)
{
	double f1 = d_func_dx(x + delta_h, y);
	double f2 = d_func_dx(x - delta_h, y);
	double f3 = (f1 - f2) / (2 * delta_h);
	return f3;
}

double dd_func_dxdy(double x, double y)
{
	double f1 = d_func_dx(x, y + delta_h);
	double f2 = d_func_dx(x, y - delta_h);
	double f3 = (f1 - f2) / (2 * delta_h);
	return f3;
}

double dd_func_dydx(double x, double y)
{
	double f1 = d_func_dy(x + delta_h, y);
	double f2 = d_func_dy(x - delta_h, y);
	double f3 = (f1 - f2) / (2* delta_h);
	return f3;
}

double dd_func_dydy(double x, double y)
{
	double f1 = d_func_dy(x, y + delta_h);
	double f2 = d_func_dy(x, y - delta_h);
	double f3 = (f1 - f2) / (2 * delta_h);
	return f3;
}

void calc_hasse_H(double hasse_H[][DIM], double x1, double y1)
{
	hasse_H[0][0] = dd_func_dxdx(x1, y1);
	hasse_H[0][1] = dd_func_dxdy(x1, y1);
	hasse_H[1][0] = dd_func_dydx(x1, y1);
	hasse_H[1][1] = dd_func_dydy(x1, y1);
}

void decomp(double a[][DIM], double b[], double x3[]);

void search()
{
	double sigma = 0.000001; //収束条件

	for (int i = 1; i < ITERATION; i++) {
		for (int j = 0; j < SWARM_SIZE; j++) {
			double x_bar = x[i - 1][j];
			double y_bar = y[i - 1][j];

			double nabla_f[DIM];
			calc_nabla_f(nabla_f, x_bar, y_bar);

			double F[DIM];
			F[0] = -nabla_f[0];
			F[1] = -nabla_f[1];

			double hasse_H[DIM][DIM];
			calc_hasse_H(hasse_H, x_bar, y_bar);
			
			// LU分解でΔxを求める
			double delta_x[DIM];
			decomp(hasse_H, F, delta_x);

			// x、yの更新
			x[i][j] = x_bar + delta_x[0];
			if (x[i][j] > X_MAX)
				x[i][j] = X_MAX;
			if (x[i][j] < X_MIN)
				x[i][j] = X_MIN;

			y[i][j] = y_bar + delta_x[1];
			if (y[i][j] > Y_MAX)
				y[i][j] = Y_MAX;
			if (y[i][j] < Y_MIN)
				y[i][j] = Y_MIN;

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

//			double norm = sqrt(delta_x[0] * delta_x[0] + delta_x[1] * delta_x[1]);
//			if (norm < sigma)
//				break;
		}
	}
}

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]);
	}
}

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

	set_function_conditions();

	// 初期化
	initialize();

	// ニュートン法による最適値探索
	search();

	write_data();
}

// LU分解法
void decomp(double a[][DIM], double b[], double x3[])
{
	// ピボット選択
	int p[DIM] = {0, 1};
	for(int pivot = 0; pivot < DIM; pivot++) {
		// 各列で 一番値が大きい行を 探す
		int col = pivot;
		int max_row = pivot;
		double max_val = 0;
		for (int i = p[pivot]; i < DIM; i++) {
			int row = p[i];
			if (fabs(a[row][col]) > max_val) {
				// 一番値が大きい行
				max_val = fabs(a[row][col]);
				max_row = row;
			}
		}

		// 一番値が大きい行と入れ替え
		if (max_row != p[pivot]) {
			int tmp	= p[pivot];
			p[pivot] = max_row;
			p[max_row] = tmp;
		}
	}

/*
	// 並べ替え後の状態を確認
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) 
			cout << setw(14) << fixed << setprecision(10) << a[p[i]][j] << "\t";
		cout << setw(14) << fixed << setprecision(10) << b[p[i]] << "\t";
		cout << endl;
	}
	cout << endl;
*/

	// 前進消去
	for (int pivot = 0; pivot < DIM - 1; pivot++) {
		int p_row = p[pivot];
		int p_col = pivot;
		
		for (int i = pivot + 1; i < DIM; i++) {
			int row = p[i];
			double s = a[row][p_col] / a[p_row][p_col];
			
			for (int col = pivot; col < DIM; col++)
				a[row][col] -= a[p_row][col] * s; // これが 上三角行列
			a[row][p_col] = s;					// これが 下三角行列
			// b[row]	-= b[p_row] * s;		 // この値は変更しない
		}
	}

/*
	// 前進消去後の状態を確認
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) 
			cout << setw(14) << fixed << setprecision(10) << a[p[i]][j] << "\t";
		cout << setw(14) << fixed << setprecision(10) << b[p[i]] << "\t";
		cout << endl;
	}
	cout << endl;
*/

	// Ly = b から y を求める
	double y3[DIM] = {0, 0};
	for (int i = 0; i < DIM; i++) {
		int row = p[i];
		for (int col = 0; col < i; col++)
			b[row] -= a[row][col] * y3[p[col]];
		y3[row] = b[row];
	}

/*
	// y の 値
	for (int i = 0; i < 2; i++)
		cout << setw(12) << fixed << setprecision(10) << y3[p[i]] << "\t";
	cout << endl << endl;
*/

	// Ux = y から x を求める
	for (int i = DIM - 1; i >= 0; i--) {
		int row = p[i];
		for (int col = DIM - 1; col >= i + 1; col--)
			y3[row] -= a[row][col] * x3[col];
		x3[i] = y3[row] / a[row][i];
	}
}

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();
}

Newton2.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 = 200;
	Z_MIN = 0;
	str_formula = "x**2 + y**2 + 20.0 * sin(x)**2";
}

Newton3.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 = 1500;
	Z_MIN = -3000;
	str_formula = "x**3 + y**3 - 9 * x * y + 27";
}

Newton4.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 * x * x * x - 0.1 * y * y * y * y + 0.01 * x * x * x * y * y * y + x + 0.3 * y";
}

Newton5.cpp

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

	z = -x * x + x * y - y * y + x + y;

	return z;
}

void set_function_conditions()
{
	X_MAX = 10;
	X_MIN = -10;
	Y_MAX = 10;
	Y_MIN = -10;
	Z_MAX = 0;
	Z_MIN = -300;
	str_formula = "-x**2 + x * y - y**2 + x + y";
}