最適化(共役勾配法)

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

-------------------------makefile-------------------------
#
#     リンク
#
CFLAGS = -c -Wall -O2
OBJECT = test.o conjugate.o dsnx.o gold.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
conjugate.o: conjugate.cpp
	g++ $(CFLAGS) conjugate.cpp
gold.o: gold.cpp
	g++ $(CFLAGS) gold.cpp
snx.o: snx.cpp
	g++ $(CFLAGS) snx.cpp

-------------------------i_data-------------------------
					// 関数 a,一次元最適化を使用しない
関数 1 変数の数 2 最大試行回数 100 一次元最適化 0
許容誤差 1.0e-10 刻み幅 0.5
初期値 0.0 0.0
					// 関数 a,一次元最適化を使用する
関数 1 変数の数 2 最大試行回数 100 一次元最適化 1
許容誤差 1.0e-10 刻み幅 0.5
初期値 0.0 0.0
					// 関数 b,一次元最適化を使用しない
関数 2 変数の数 2 最大試行回数 5000 一次元最適化 0
許容誤差 1.0e-10 刻み幅 0.003
初期値 0.0 0.0
					// 関数 b,一次元最適化を使用する
関数 2 変数の数 2 最大試行回数 100 一次元最適化 1
許容誤差 1.0e-10 刻み幅 0.001
初期値 0.0 0.0
					// 関数 c,一次元最適化を使用しない
関数 3 変数の数 2 最大試行回数 500 一次元最適化 0
許容誤差 1.0e-10 刻み幅 0.04
初期値 0.0 0.0
					// 関数 c,一次元最適化を使用する
関数 3 変数の数 2 最大試行回数 100 一次元最適化 1
許容誤差 1.0e-10 刻み幅 0.01
初期値 0.0 0.0

-------------------------main-------------------------
/********************************/
/* 共役勾配法による最小値の計算 */
/*      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 conjugate(int, int, int, double, double, double *, double *, double *,
          double (*)(double, double *, double *),
          void (*)(double *, double *));

int main()
{
	double eps, step, *x, *dx, y;
	int fun, i1, max, n, opt_1, sw = 0;
					// データの入力
	scanf("%*s %d %*s %d %*s %d %*s %d", &fun, &n, &max, &opt_1);
	scanf("%*s %lf %*s %lf", &eps, &step);
	x  = new double [n];
	dx = new double [n];
	scanf("%*s");
	for (i1 = 0; i1 < n; i1++)
		scanf("%lf", &x[i1]);
					// 実行
	switch (fun) {
		case 1:
			sw = conjugate(opt_1, max, n, eps, step, &y, x, dx, snx1, dsnx1);
			break;
		case 2:
			sw = conjugate(opt_1, max, n, eps, step, &y, x, dx, snx2, dsnx2);
			break;
		case 3:
			sw = conjugate(opt_1, max, n, eps, step, &y, x, dx, 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;
}

-------------------------conjugate.cpp-------------------------
/********************************************************/
/* 共役勾配法                                           */
/*      opt_1 : =0 : 1次元最適化を行わない             */
/*              =1 : 1次元最適化を行う                 */
/*      max : 最大繰り返し回数                          */
/*      n : 次元                                        */
/*      eps : 収束判定条件                              */
/*      step : きざみ幅                                 */
/*      y : 最小値                                      */
/*      x : x(初期値と答え)                             */
/*      dx : 関数の微分値                               */
/*      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 conjugate(int opt_1, int max, int n, double eps, double step, double *y,
              double *x, double *dx, double (*snx)(double, double *, double *), 
              void (*dsnx)(double *, double *))
{
	double b = 0.0, f1, f2, *g, k, s1 = 1.0, s2, sp, y1, y2;
	int count = 0, i1, sw = 0, sw1;

	g  = new double [n];

	y1 = snx(0.0, x, dx);

	while (count < max && sw == 0) {
					// 傾きの計算
		dsnx(x, g);
					// 傾きのチェック
		s2 = 0.0;
		for (i1 = 0; i1 < n; i1++)
			s2 += g[i1] * g[i1];
					// 収束していない
		if (sqrt(s2) > eps) {
						// 方向の計算
			count++;
			if (count%n == 1)
				b = 0.0;
			else
				b = s2 / s1;
			for (i1 = 0; i1 < n; i1++)
				dx[i1] = g[i1] + b * dx[i1];
						// 1次元最適化を行わない
			if (opt_1 == 0) {
							// 新しい点
				for (i1 = 0; i1 < n; i1++)
					x[i1] += step * dx[i1];
							// 新しい関数値
				y2 = snx(0.0, x, dx);
								// 関数値の変化が大きい
				if (fabs(y2-y1) > eps) {
					y1 = y2;
					s1 = s2;
				}
								// 収束(関数値の変化<eps)
				else {
					sw = count;
					*y = y2;
				}
			}
						// 1次元最適化を行う
			else {
							// 区間を決める
				sw1 = 0;
				f1  = y1;
				sp  = step;
				f2  = snx(sp, x, dx);
				if (f2 > f1)
					sp = -step;
				for (i1 = 0; i1 < max && sw1 == 0; i1++) {
					f2 = snx(sp, x, 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, x, dx, snx);
								// 黄金分割法が失敗
					if (sw1 < 0)
						sw = -3;
								// 黄金分割法が成功
					else {
									// 新しい点
						for (i1 = 0; i1 < n; i1++)
							x[i1] += k * dx[i1];
										// 関数値の変化が大きい
						if (fabs(y1-y2) > eps) {
							y1 = y2;
							s1 = s2;
						}
										// 収束(関数値の変化<eps)
						else {
							sw = count;
							*y = y2;
						}
					}
				}
			}
		}
					// 収束(傾き<eps)
		else {
			sw = count;
			*y = y1;
		}
	}

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

	delete [] g;

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