最適化(準 Newton 法)

mskefile,データ例,プログラム

-------------------------makefile-------------------------
#
#     リンク
#
CFLAGS = -c -Wall -O2
OBJECT = test.o dsnx.o gold.o DFP_BFGS.o snx.o
pgm: $(OBJECT)
	g++ $(OBJECT) -o test -lm
#
#     コンパイル
#
test.o: test.cpp
	g++ $(CFLAGS) test.cpp
dsnx.o: dsnx.cpp
	g++ $(CFLAGS) dsnx.cpp
gold.o: gold.cpp
	g++ $(CFLAGS) gold.cpp
DFP_BFGS.o: DFP_BFGS.cpp
	g++ $(CFLAGS) DFP_BFGS.cpp
snx.o: snx.cpp
	g++ $(CFLAGS) snx.cpp

-------------------------i_data-------------------------
					// 関数 a,DFP法,一次元最適化を使用しない
関数 1 変数の数 2 最大試行回数 100 一次元最適化 0
方法 0 許容誤差 1.0e-10 刻み幅 0.5
初期値 0.0 0.0
					// 関数 a,DFP法,一次元最適化を使用する
関数 1 変数の数 2 最大試行回数 100 一次元最適化 1
方法 0 許容誤差 1.0e-10 刻み幅 0.5
初期値 0.0 0.0
					// 関数 a,BFGS法,一次元最適化を使用しない
関数 1 変数の数 2 最大試行回数 100 一次元最適化 0
方法 1 許容誤差 1.0e-10 刻み幅 0.5
初期値 0.0 0.0
					// 関数 a,BFGS法,一次元最適化を使用する
関数 1 変数の数 2 最大試行回数 100 一次元最適化 1
方法 1 許容誤差 1.0e-10 刻み幅 0.5
初期値 0.0 0.0
					// 関数 b,DFP法,一次元最適化を使用しない
関数 2 変数の数 2 最大試行回数 100 一次元最適化 0
方法 0 許容誤差 1.0e-10 刻み幅 0.2
初期値 0.0 0.0
					// 関数 b,DFP法,一次元最適化を使用する
関数 2 変数の数 2 最大試行回数 100 一次元最適化 1
方法 0 許容誤差 1.0e-10 刻み幅 0.1
初期値 0.0 0.0
					// 関数 b,BFGS法,一次元最適化を使用しない
関数 2 変数の数 2 最大試行回数 1000 一次元最適化 0
方法 1 許容誤差 1.0e-10 刻み幅 0.02
初期値 0.0 0.0
					// 関数 b,BFGS法,一次元最適化を使用する
関数 2 変数の数 2 最大試行回数 100 一次元最適化 1
方法 1 許容誤差 1.0e-10 刻み幅 0.002
初期値 0.0 0.0
					// 関数 c,DFP法,一次元最適化を使用しない
関数 3 変数の数 2 最大試行回数 200 一次元最適化 0
方法 0 許容誤差 1.0e-10 刻み幅 0.1
初期値 0.0 0.0
					// 関数 c,DFP法,一次元最適化を使用する
関数 3 変数の数 2 最大試行回数 100 一次元最適化 1
方法 0 許容誤差 1.0e-10 刻み幅 0.1
初期値 0.0 0.0
					// 関数 c,BFGS法,一次元最適化を使用しない
関数 3 変数の数 2 最大試行回数 200 一次元最適化 0
方法 1 許容誤差 1.0e-10 刻み幅 0.09
初期値 0.0 0.0
					// 関数 c,BFGS法,一次元最適化を使用する
関数 3 変数の数 2 最大試行回数 200 一次元最適化 1
方法 1 許容誤差 1.0e-10 刻み幅 0.09
初期値 0.0 0.0

-------------------------main-------------------------
/**********************************/
/* 準Newton法によるの最小値の計算 */
/*      coded by Y.Suganuma       */
/**********************************/
#include <stdio.h>

void dsnx1(double *, double *);
double snx1(double, double *, double *);
void dsnx2(double *, double *);
double snx2(double, double *, double *);
void dsnx3(double *, double *);
double snx3(double, double *, double *);
int DFP_BFGS(int, int, int, int, double, double, double *, double *, double *, double **,
          double (*)(double, double *, double *),
          void (*)(double *, double *));

int main()
{
	double eps, **H, step, *x, *dx, y;
	int fun, i1, max, method, n, opt_1, sw = 0;
					// データの入力
	scanf("%*s %d %*s %d %*s %d %*s %d", &fun, &n, &max, &opt_1);
	scanf("%*s %d %*s %lf %*s %lf", &method, &eps, &step);
	x  = new double [n];
	dx = new double [n];
	H  = new double * [n];
	scanf("%*s");
	for (i1 = 0; i1 < n; i1++) {
		scanf("%lf", &x[i1]);
		H[i1] = new double [n];
	}
					// 実行
	switch (fun) {
		case 1:
			sw = DFP_BFGS(method, opt_1, max, n, eps, step, &y, x, dx, H, snx1, dsnx1);
			break;
		case 2:
			sw = DFP_BFGS(method, opt_1, max, n, eps, step, &y, x, dx, H, snx2, dsnx2);
			break;
		case 3:
			sw = DFP_BFGS(method, opt_1, max, n, eps, step, &y, x, dx, H, snx3, dsnx3);
			break;
	}
					// 結果の出力
	if (sw < 0) {
		printf("   収束しませんでした!");
		switch (sw) {
			case -1:
				printf("(収束回数)\n");
				break;
			case -2:
				printf("(1次元最適化の区間)\n");
				break;
			case -3:
				printf("(黄金分割法)\n");
				break;
		}
	}
	else {
		printf("   結果=");
		for (i1 = 0; i1 < n; i1++)
			printf("%f ", x[i1]);
		printf(" 最小値=%f  回数=%d\n", y, sw);
	}

	return 0;
}

-------------------------DFP_BFGS.cpp-------------------------
/********************************************************/
/* DFP法 or BFGS法                                      */
/*      method : =0 : DFP法                             */
/*               =1 : BFGS法                            */
/*      opt_1 : =0 : 1次元最適化を行わない             */
/*              =1 : 1次元最適化を行う                 */
/*      max : 最大繰り返し回数                          */
/*      n : 次元                                        */
/*      eps : 収束判定条件                              */
/*      step : きざみ幅                                 */
/*      y : 最小値                                      */
/*      x1 : x(初期値と答え)                            */
/*      dx : 関数の微分値                               */
/*      H : Hesse行列の逆行列                           */
/*      snx : 関数値を計算する関数名                    */
/*      dsnx : 関数の微分を計算する関数名(符号を変える) */
/*      return : >=0 : 正常終了(収束回数)               */
/*               =-1 : 収束せず                         */
/*               =-2 : 1次元最適化の区間が求まらない   */
/*               =-3 : 黄金分割法が失敗                 */
/********************************************************/
#include <math.h>

double gold(double, double, double, double *, int *, int, double *, double *,
            double (*)(double, double *, double *));

int DFP_BFGS(int method, int opt_1, int max, int n, double eps, double step, double *y,
          double *x1, double *dx, double **H, double (*snx)(double, double *, double *), 
          void (*dsnx)(double *, double *))
{
	double f1, f2, *g, *g1, *g2, **H1, **H2, **I, k, *s, sm, sm1, sm2,
           sp, *w, *x2, y1, y2;
	int count = 0, i1, i2, i3, sw = 0, sw1;

	x2 = new double [n];
	g  = new double [n];
	g1 = new double [n];
	g2 = new double [n];
	s  = new double [n];
	w  = new double [n];

	y1 = snx(0.0, x1, dx);
	dsnx(x1, g1);

	H1 = new double * [n];
	H2 = new double * [n];
	I  = new double * [n];
	for (i1 = 0; i1 < n; i1++) {
		H1[i1] = new double [n];
		H2[i1] = new double [n];
		I[i1]  = new double [n];
		for (i2 = 0; i2 < n; i2++) {
			H[i1][i2] = 0.0;
			I[i1][i2] = 0.0;
		}
		H[i1][i1] = 1.0;
		I[i1][i1] = 1.0;
	}

	while (count < max && sw == 0) {
		count++;
					// 方向の計算
		for (i1 = 0; i1 < n; i1++) {
			dx[i1] = 0.0;
			for (i2 = 0; i2 < n; i2++)
				dx[i1] -= H[i1][i2] * g1[i2];
		}
					// 1次元最適化を行わない
		if (opt_1 == 0) {
						// 新しい点
			for (i1 = 0; i1 < n; i1++)
				x2[i1] = x1[i1] + step * dx[i1];
						// 新しい関数値
			y2 = snx(0.0, x2, dx);
		}
					// 1次元最適化を行う
		else {
						// 区間を決める
			sw1 = 0;
			f1  = y1;
			sp  = step;
			f2  = snx(sp, x1, dx);
			if (f2 > f1)
				sp = -step;
			for (i1 = 0; i1 < max && sw1 == 0; i1++) {
				f2 = snx(sp, x1, dx);
				if (f2 > f1)
					sw1 = 1;
				else {
					sp *= 2.0;
					f1  = f2;
				}
			}
						// 区間が求まらない
			if (sw1 == 0)
				sw = -2;
						// 区間が求まった
			else {
							// 黄金分割法
				k = gold(0.0, sp, eps, &y2, &sw1, max, x1, dx, snx);
							// 黄金分割法が失敗
				if (sw1 < 0)
					sw = -3;
							// 黄金分割法が成功
				else {
								// 新しい点
					for (i1 = 0; i1 < n; i1++)
						x2[i1] = x1[i1] + k * dx[i1];
				}
			}
		}

		if (sw == 0) {
					// 収束(関数値の変化<eps)
			if (fabs(y2-y1) < eps) {
				sw = count;
				*y = y2;
				for (i1 = 0; i1 < n; i1++)
					x1[i1] = x2[i1];
			}
					// 関数値の変化が大きい
			else {
						// 傾きの計算
				dsnx(x2, g2);
				sm = 0.0;
				for (i1 = 0; i1 < n; i1++)
					sm += g2[i1] * g2[i1];
				sm = sqrt(sm);
						// 収束(傾き<eps)
				if (sm < eps) {
					sw = count;
					*y = y2;
					for (i1 = 0; i1 < n; i1++)
						x1[i1] = x2[i1];
				}
						// 収束していない
				else {
					for (i1 = 0; i1 < n; i1++) {
						g[i1] = g2[i1] - g1[i1];
						s[i1] = x2[i1] - x1[i1];
					}
					sm1 = 0.0;
					for (i1 = 0; i1 < n; i1++)
						sm1 += s[i1] * g[i1];
							// DFP法
					if (method == 0) {
						for (i1 = 0; i1 < n; i1++) {
							w[i1] = 0.0;
							for (i2 = 0; i2 < n; i2++)
								w[i1] += g[i2] * H[i2][i1];
						}
						sm2 = 0.0;
						for (i1 = 0; i1 < n; i1++)
							sm2 += w[i1] * g[i1];
						for (i1 = 0; i1 < n; i1++) {
							w[i1] = 0.0;
							for (i2 = 0; i2 < n; i2++)
								w[i1] += H[i1][i2] * g[i2];
						}
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++)
								H1[i1][i2] = w[i1] * g[i2];
						}
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++) {
								H2[i1][i2] = 0.0;
								for (i3 = 0; i3 < n; i3++)
									H2[i1][i2] += H1[i1][i3] * H[i3][i2];
							}
						}
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++)
								H[i1][i2] = H[i1][i2]  + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2;
						}
					}
							// BFGS法
					else {
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++)
								H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1;
						}
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++) {
								H2[i1][i2] = 0.0;
								for (i3 = 0; i3 < n; i3++)
									H2[i1][i2] += H1[i1][i3] * H[i3][i2];
							}
						}
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++) {
								H[i1][i2] = 0.0;
								for (i3 = 0; i3 < n; i3++)
									H[i1][i2] += H2[i1][i3] * H1[i3][i2];
							}
						}
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++)
								H[i1][i2] += s[i1] * s[i2] / sm1;
						}
					}
					y1 = y2;
					for (i1 = 0; i1 < n; i1++) {
						g1[i1] = g2[i1];
						x1[i1] = x2[i1];
					}
				}
			}
		}
	}

	if (sw == 0)
		sw = -1;

	delete [] x2;
	delete [] g;
	delete [] g1;
	delete [] g2;
	delete [] s;
	delete [] w;

	for (i1 = 0; i1 < n; i1++) {
		delete [] H1[i1];
		delete [] H2[i1];
		delete [] I[i1];
	}
	delete [] H1;
	delete [] H2;
	delete [] I;

	return sw;
}

-------------------------gold.cpp-------------------------
/****************************************************************/
/* 黄金分割法(与えられた方向での最小値)                         */
/*      a,b : 初期区間 a < b                                    */
/*      eps : 許容誤差                                          */
/*      val : 関数値                                            */
/*      ind : 計算状況                                          */
/*              >= 0 : 正常終了(収束回数)                       */
/*              = -1 : 収束せず                                 */
/*      max : 最大試行回数                                      */
/*      w : 位置                                                */
/*      dw : 傾きの成分                                         */
/*      snx : 関数値を計算する関数の名前                        */
/*      return : 結果(w+y*dwのy)                              */
/****************************************************************/
#include <math.h>

double gold(double a, double b, double eps, double *val, int *ind, int max,
            double *w, double *dw, double (*snx)(double, double *, double *))
{
	double f1, f2, fa, fb, tau = (sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
	int count = 0;
					// 初期設定
	*ind  = -1;
	x1    = b - tau * (b - a);
	x2    = a + tau * (b - a);
	f1    = snx(x1, w, dw);
	f2    = snx(x2, w, dw);
					// 計算
	while (count < max && *ind < 0) {
		count += 1;
		if (f2 > f1) {
			if (fabs(b-a) < eps) {
				*ind = 0;
				x    = x1;
				*val = f1;
			}
			else {
				b  = x2;
				x2 = x1;
				x1 = a + (1.0 - tau) * (b - a);
				f2 = f1;
				f1 = snx(x1, w, dw);
			}
		}
		else {
			if (fabs(b-a) < eps) {
				*ind = 0;
				x    = x2;
				*val = f2;
				f1   = f2;
			}
			else {
				a  = x1;
				x1 = x2;
				x2 = b - (1.0 - tau) * (b - a);
				f1 = f2;
				f2 = snx(x2, w, dw);
			}
		}
	}
					// 収束した場合の処理
	if (*ind == 0) {
		*ind = count;
		fa   = snx(a, w, dw);
		fb   = snx(b, w, dw);
		if (fb < fa) {
			a  = b;
			fa = fb;
		}
		if (fa < f1) {
			x    = a;
			*val = fa;
		}
	}

	return x;
}

-------------------------snx.cpp-------------------------
/*********************************************/
/* 与えられた点xから,dx方向へk*dxだけ進んだ */
/* 点における関数値を計算する                */
/*********************************************/
					// 関数1
double snx1(double k, double *x, double *dx)
{
	double x1, y1, y, w[2];

	w[0] = x[0] + k * dx[0];
	w[1] = x[1] + k * dx[1];
	x1   = w[0] - 1.0;
	y1   = w[1] - 2.0;
	y    = x1 * x1 + y1 * y1;

	return y;
}
					// 関数2
double snx2(double k, double *x, double *dx)
{
	double x1, y1, y, w[2];

	w[0] = x[0] + k * dx[0];
	w[1] = x[1] + k * dx[1];
	x1   = w[1] - w[0] * w[0];
	y1   = 1.0 - w[0];
	y    = 100.0 * x1 * x1 + y1 * y1;

	return y;
}
					// 関数3
double snx3(double k, double *x, double *dx)
{
	double x1, y1, z1, y, w[2];

	w[0] = x[0] + k * dx[0];
	w[1] = x[1] + k * dx[1];
	x1   = 1.5 - w[0] * (1.0 - w[1]);
	y1   = 2.25 - w[0] * (1.0 - w[1] * w[1]);
	z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
	y    = x1 * x1 + y1 * y1 + z1 * z1;

	return y;
}

-------------------------dsnx.cpp-------------------------
/********************/
/* 微係数を計算する */
/********************/
					// 関数1
void dsnx1(double *x, double *dx)
{
	dx[0] = 2.0 * (x[0] - 1.0);
	dx[1] = 2.0 * (x[1] - 2.0);
}
					// 関数2
void dsnx2(double *x, double *dx)
{
	dx[0] = -400.0 * x[0] * (x[1] - x[0] * x[0]) - 2.0 * (1.0 - x[0]);
	dx[1] = 200.0 * (x[1] - x[0] * x[0]);
}
					// 関数3
void dsnx3(double *x, double *dx)
{
	dx[0] = -2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) -
            2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
            2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
	dx[1] = 2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) +
            4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
            6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
}