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("（１次元最適化の区間）\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 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x1 : x(初期値と答え) */
/* dx : 関数の微分値 */
/* H : Hesse行列の逆行列 */
/* snx : 関数値を計算する関数名 */
/* dsnx : 関数の微分を計算する関数名(符号を変える) */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-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];
}
// １次元最適化を行わない
if (opt_1 == 0) {
// 新しい点
for (i1 = 0; i1 < n; i1++)
x2[i1] = x1[i1] + step * dx[i1];
// 新しい関数値
y2 = snx(0.0, x2, dx);
}
// １次元最適化を行う
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) {
// 収束（関数値の変化＜ｅｐｓ）
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);
// 収束（傾き＜ｅｐｓ）
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だけ進んだ */
/* 点における関数値を計算する */
/*********************************************/
// 関数１
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;
}
// 関数２
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;
}
// 関数３
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-------------------------
/********************/
/* 微係数を計算する */
/********************/
// 関数１
void dsnx1(double *x, double *dx)
{
dx[0] = 2.0 * (x[0] - 1.0);
dx[1] = 2.0 * (x[1] - 2.0);
}
// 関数２
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]);
}
// 関数３
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]));
}