| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/**********************************/
/* 準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法 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;
}
/****************************************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* 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;
}
/*********************************************/
/* 与えられた点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;
}
/********************/
/* 微係数を計算する */
/********************/
// 関数１
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]));
}
/**********************************/
/* 準Newton法によるの最小値の計算 */
/* coded by Y.Suganuma */
/**********************************/
import java.io.*;
import java.util.*;
public class Test {
public static void main(String args[]) throws IOException
{
double eps, H[][], step, x[], dx[], y[];
int fun, i1, max, method, n, opt_1, sw = 0;
StringTokenizer str;
String line;
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
// データの入力
// 1 行目
line = in.readLine();
str = new StringTokenizer(line, " ");
line = str.nextToken();
line = str.nextToken();
fun = Integer.parseInt(line);
line = str.nextToken();
line = str.nextToken();
n = Integer.parseInt(line);
line = str.nextToken();
line = str.nextToken();
max = Integer.parseInt(line);
line = str.nextToken();
line = str.nextToken();
opt_1 = Integer.parseInt(line);
// 2 行目
line = in.readLine();
str = new StringTokenizer(line, " ");
line = str.nextToken();
line = str.nextToken();
method = Integer.parseInt(line);
line = str.nextToken();
line = str.nextToken();
eps = Double.parseDouble(line);
line = str.nextToken();
line = str.nextToken();
step = Double.parseDouble(line);
// 3 行目
x = new double [n];
dx = new double [n];
y = new double [1];
H = new double [n][n];
line = in.readLine();
str = new StringTokenizer(line, " ");
line = str.nextToken();
for (i1 = 0; i1 < n; i1++) {
line = str.nextToken();
x[i1] = Double.parseDouble(line);
}
// 実行
Kansu kn = new Kansu(fun);
sw = kn.DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H);
// 結果の出力
if (sw < 0) {
System.out.print(" 収束しませんでした！");
switch (sw) {
case -1:
System.out.println("（収束回数）");
break;
case -2:
System.out.print("（１次元最適化の区間）");
break;
case -3:
System.out.print("（黄金分割法）");
break;
}
}
else {
System.out.print(" 結果＝");
for (i1 = 0; i1 < n; i1++)
System.out.print(x[i1] + " ");
System.out.println(" 最小値＝" + y[0] + " 回数＝" + sw);
}
}
}
/************************/
/* 関数値，微係数の計算 */
/************************/
class Kansu extends P_Newton
{
private int sw;
// コンストラクタ
Kansu (int s) {sw = s;}
// 与えられた点xから，dx方向へk*dxだけ進んだ
// 点における関数値を計算する
double snx(double k, double x[], double dx[])
{
double x1, y1, z1, y = 0.0, w[] = new double [2];
switch (sw) {
case 1:
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;
break;
case 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;
break;
case 3:
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;
break;
}
return y;
}
// 関数の微分
void dsnx(double x[], double dx[])
{
switch (sw) {
case 1:
dx[0] = 2.0 * (x[0] - 1.0);
dx[1] = 2.0 * (x[1] - 2.0);
break;
case 2:
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]);
break;
case 3:
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]));
break;
}
}
}
abstract class P_Newton {
/********************************************************/
/* DFP法 or BFGS法 */
/* method : =0 : DFP法 */
/* =1 : BFGS法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x1 : x(初期値と答え) */
/* dx : 関数の微分値 */
/* H : Hesse行列の逆行列 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
abstract double snx(double k, double x[], double dx[]);
abstract void dsnx(double x[], double dx[]);
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 f1, f2, g[], g1[], g2[], H1[][], H2[][], I[][], k, s[], sm, sm1, sm2,
sp, w[], x2[], y1[] = new double [1], y2[] = new double [1];
int count = 0, i1, i2, i3, sw = 0, sw1[] = new int [1];
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[0] = snx(0.0, x1, dx);
dsnx(x1, g1);
H1 = new double [n][n];
H2 = new double [n][n];
I = new double [n][n];
for (i1 = 0; i1 < n; i1++) {
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[0] = snx(0.0, x2, dx);
}
// １次元最適化を行う
else {
// 区間を決める
sw1[0] = 0;
f1 = y1[0];
sp = step;
f2 = snx(sp, x1, dx);
if (f2 > f1)
sp = -step;
for (i1 = 0; i1 < max && sw1[0] == 0; i1++) {
f2 = snx(sp, x1, dx);
if (f2 > f1)
sw1[0] = 1;
else {
sp *= 2.0;
f1 = f2;
}
}
// 区間が求まらない
if (sw1[0] == 0)
sw = -2;
// 区間が求まった
else {
// 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x1, dx);
// 黄金分割法が失敗
if (sw1[0] < 0)
sw = -3;
// 黄金分割法が成功
else {
// 新しい点
for (i1 = 0; i1 < n; i1++)
x2[i1] = x1[i1] + k * dx[i1];
}
}
}
if (sw == 0) {
// 収束（関数値の変化＜ｅｐｓ）
if (Math.abs(y2[0]-y1[0]) < eps) {
sw = count;
y[0] = y2[0];
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 = Math.sqrt(sm);
// 収束（傾き＜ｅｐｓ）
if (sm < eps) {
sw = count;
y[0] = y2[0];
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[0] = y2[0];
for (i1 = 0; i1 < n; i1++) {
g1[i1] = g2[i1];
x1[i1] = x2[i1];
}
}
}
}
}
if (sw == 0)
sw = -1;
return sw;
}
/****************************************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* return : 結果(w＋y＊dwのy) */
/****************************************************************/
double gold(double a, double b, double eps, double val[], int ind[], int max, double w[], double dw[])
{
double f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
int count = 0;
// 初期設定
ind[0] = -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] < 0) {
count += 1;
if (f2 > f1) {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x1;
val[0] = f1;
}
else {
b = x2;
x2 = x1;
x1 = a + (1.0 - tau) * (b - a);
f2 = f1;
f1 = snx(x1, w, dw);
}
}
else {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x2;
val[0] = 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] == 0) {
ind[0] = count;
fa = snx(a, w, dw);
fb = snx(b, w, dw);
if (fb < fa) {
a = b;
fa = fb;
}
if (fa < f1) {
x = a;
val[0] = fa;
}
}
return x;
}
}
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>最適化（準 Newton 法）</TITLE>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
<SCRIPT TYPE="text/javascript">
str1 = "";
str2 = new Array();
function main()
{
// データの設定
let n = parseInt(document.getElementById("n").value);
let max = parseInt(document.getElementById("max").value);
let step = parseFloat(document.getElementById("step").value);
let x = document.getElementById("x0").value.split(/ {1,}/)
for (let i1 = 0; i1 < n; i1++)
x[i1] = parseFloat(x[i1]);
let opt_1 = 0;
if (document.getElementById("line").options[1].selected)
opt_1 = 1;
let method = 0;
if (document.getElementById("method").options[1].selected)
method = 1;
str1 = document.getElementById("func").value + ";";
str2 = document.getElementById("dfunc").value.split("\n")
for (let i1 = 0; i1 < n; i1++)
str2[i1] = str2[i1] + ";";
let eps = 1.0e-10;
let y = new Array();
let dx = new Array();
let H = new Array();
for (let i1 = 0; i1 < n; i1++) {
dx[i1] = 0.0;
H[i1] = new Array();
}
// 実行
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx, dsnx);
// 結果
let res;
if (sw < 0) {
res = " 収束しませんでした！";
switch (sw) {
case -1:
res += "（収束回数）";
break;
case -2:
res += "（１次元最適化の区間）";
break;
case -3:
res += "（黄金分割法）";
break;
}
document.getElementById("res").value = res;
}
else {
let c = 100000;
res = " x ＝";
for (let i1 = 0; i1 < n; i1++) {
let s1 = Math.round(x[i1] * c) / c;
res = res + " " + s1;
}
res += "\n";
let s2 = Math.round(y[0] * c) / c;
res = res + " 最小値 ＝ " + s2 + " 回数 ＝ " + sw;
document.getElementById("res").value = res;
}
}
/****************/
/* 関数値の計算 */
/****************/
function snx(k, x, dx)
{
let y = eval(str1);
return y;
}
/************************/
/* 関数値の微分値の計算 */
/************************/
function dsnx(n, x, dx)
{
for (let i1 = 0; i1 < n; i1++)
dx[i1] = eval(str2[i1]);
}
/******************************************************/
/* DFP法 or BFGS法 */
/* method : =0 : DFP法 */
/* =1 : BFGS法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x1 : x(初期値と答え) */
/* dx : 関数の微分値 */
/* H : Hesse行列の逆行列 */
/* fn : 関数値を計算する関数 */
/* dfn : 関数の微分値を計算する関数 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/******************************************************/
function DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, H, fn, dfn)
{
let f1, f2, k, sm, sm1, sm2, sp, count = 0, sw = 0;
let g = new Array();
let s = new Array();
let w = new Array();
let x2 = new Array();
let y1 = new Array();
let y2 = new Array();
let g1 = new Array();
let g2 = new Array();
let sw1 = new Array();
let H1 = new Array();
let H2 = new Array();
let I = new Array();
for (let i1 = 0; i1 < n; i1++) {
H1[i1] = new Array();
H2[i1] = new Array();
I[i1] = new Array();
}
y1[0] = fn(0.0, x1, dx);
dfn(n, x1, g1);
for (let i1 = 0; i1 < n; i1++) {
for (let 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 (let i1 = 0; i1 < n; i1++) {
dx[i1] = 0.0;
for (let i2 = 0; i2 < n; i2++)
dx[i1] -= H[i1][i2] * g1[i2];
}
// １次元最適化を行わない
if (opt_1 == 0) {
// 新しい点
for (let i1 = 0; i1 < n; i1++)
x2[i1] = x1[i1] + step * dx[i1];
// 新しい関数値
y2[0] = fn(0.0, x2, dx);
}
// １次元最適化を行う
else {
// 区間を決める
sw1[0] = 0;
f1 = y1[0];
sp = step;
f2 = fn(sp, x1, dx);
if (f2 > f1)
sp = -step;
for (let i1 = 0; i1 < max && sw1[0] == 0; i1++) {
f2 = fn(sp, x1, dx);
if (f2 > f1)
sw1[0] = 1;
else {
sp *= 2.0;
f1 = f2;
}
}
// 区間が求まらない
if (sw1[0] == 0)
sw = -2;
// 区間が求まった
else {
// 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x1, dx, snx);
// 黄金分割法が失敗
if (sw1[0] < 0)
sw = -3;
// 黄金分割法が成功
else {
// 新しい点
for (let i1 = 0; i1 < n; i1++)
x2[i1] = x1[i1] + k * dx[i1];
}
}
}
if (sw == 0) {
// 収束（関数値の変化＜ｅｐｓ）
if (Math.abs(y2[0]-y1[0]) < eps) {
sw = count;
y[0] = y2[0];
for (let i1 = 0; i1 < n; i1++)
x1[i1] = x2[i1];
}
// 関数値の変化が大きい
else {
// 傾きの計算
dfn(n, x2, g2);
sm = 0.0;
for (let i1 = 0; i1 < n; i1++)
sm += g2[i1] * g2[i1];
sm = Math.sqrt(sm);
// 収束（傾き＜ｅｐｓ）
if (sm < eps) {
sw = count;
y[0] = y2[0];
for (let i1 = 0; i1 < n; i1++)
x1[i1] = x2[i1];
}
// 収束していない
else {
for (let i1 = 0; i1 < n; i1++) {
g[i1] = g2[i1] - g1[i1];
s[i1] = x2[i1] - x1[i1];
}
sm1 = 0.0;
for (let i1 = 0; i1 < n; i1++)
sm1 += s[i1] * g[i1];
// DFP法
if (method == 0) {
for (let i1 = 0; i1 < n; i1++) {
w[i1] = 0.0;
for (let i2 = 0; i2 < n; i2++)
w[i1] += g[i2] * H[i2][i1];
}
sm2 = 0.0;
for (let i1 = 0; i1 < n; i1++)
sm2 += w[i1] * g[i1];
for (let i1 = 0; i1 < n; i1++) {
w[i1] = 0.0;
for (let i2 = 0; i2 < n; i2++)
w[i1] += H[i1][i2] * g[i2];
}
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++)
H1[i1][i2] = w[i1] * g[i2];
}
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++) {
H2[i1][i2] = 0.0;
for (let i3 = 0; i3 < n; i3++)
H2[i1][i2] += H1[i1][i3] * H[i3][i2];
}
}
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++)
H[i1][i2] = H[i1][i2] + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2;
}
}
// BFGS法
else {
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++)
H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1;
}
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++) {
H2[i1][i2] = 0.0;
for (let i3 = 0; i3 < n; i3++)
H2[i1][i2] += H1[i1][i3] * H[i3][i2];
}
}
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++) {
H[i1][i2] = 0.0;
for (let i3 = 0; i3 < n; i3++)
H[i1][i2] += H2[i1][i3] * H1[i3][i2];
}
}
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++)
H[i1][i2] += s[i1] * s[i2] / sm1;
}
}
y1[0] = y2[0];
for (let i1 = 0; i1 < n; i1++) {
g1[i1] = g2[i1];
x1[i1] = x2[i1];
}
}
}
}
}
if (sw == 0)
sw = -1;
return sw;
}
/******************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* fn : 関数値を計算する関数 */
/* return : 結果(w＋y＊dwのy) */
/******************************************/
function gold(a, b, eps, val, ind, max, w, dw, fn)
{
let f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
let count = 0;
// 初期設定
ind[0] = -1;
x1 = b - tau * (b - a);
x2 = a + tau * (b - a);
f1 = fn(x1, w, dw);
f2 = fn(x2, w, dw);
// 計算
while (count < max && ind[0] < 0) {
count += 1;
if (f2 > f1) {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x1;
val[0] = f1;
}
else {
b = x2;
x2 = x1;
x1 = a + (1.0 - tau) * (b - a);
f2 = f1;
f1 = fn(x1, w, dw);
}
}
else {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x2;
val[0] = f2;
f1 = f2;
}
else {
a = x1;
x1 = x2;
x2 = b - (1.0 - tau) * (b - a);
f1 = f2;
f2 = fn(x2, w, dw);
}
}
}
// 収束した場合の処理
if (ind[0] == 0) {
ind[0] = count;
fa = fn(a, w, dw);
fb = fn(b, w, dw);
if (fb < fa) {
a = b;
fa = fb;
}
if (fa < f1) {
x = a;
val[0] = fa;
}
}
return x;
}
</SCRIPT>
</HEAD>
<BODY STYLE="font-size: 130%; background-color: #eeffee;">
<H2 STYLE="text-align:center"><B>最適化（準 Newton 法）</B></H2>
<DL>
<DT> テキストフィールドおよびテキストエリアには，例として，以下に示すような関数の最小値を求める場合に対する値（ BFGS 法において，一次元最適化を行わない場合は刻み幅を 0.02，行う場合は 0.002 にしてみて下さい）が設定されています（他の問題を実行する場合は，それらを適切に修正してください）．
<P STYLE="text-align:center"><IMG SRC="s_newton.gif"></P>
<DT>n 変数関数 f(<B>x</B>) に対する準 Newton 法においては，現在の点 <B>x</B> における傾き，
<P STYLE="text-align:center"><IMG SRC="s_newton1.gif"></P>
<DT>から方向 <B>p</B> を求め，その方向にαだけ進んだ点の関数値 f( <B>x</B> + α<B>p</B>) を計算する，という処理を繰り返します．目的関数の箇所には，この例に対する f( <B>x</B> + α<B>p</B>) の計算式，目的関数の微分の箇所には，<B>g</B> の計算式が与えられています．ただし，プログラム内では，変数 x は，この例の x と y からなる配列，変数 dx は <B>g</B> に対応し，関数 f を x 及び y で微分したものからなる配列，また，k はαに対応する変数として表現されています．なお，目的関数，及び，目的関数の微分は，JavaScript の仕様に適合した形式で記述してあることに注意してください．
</DL>
<DIV STYLE="text-align:center">
変数の数(n)：<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="2">
最大繰り返し回数：<INPUT ID="max" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="1000">
刻み幅：<INPUT ID="step" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.2"><BR>
初期値(n個)：<INPUT ID="x0" STYLE="font-size: 100%" TYPE="text" SIZE="20" VALUE="0.0 0.0">
１次元最適化：<SELECT STYLE="font-size: 100%" SIZE="1" ID="line">
<OPTION VALUE="no"> 行わない </OPTION><BR>
<OPTION VALUE="yes"> 行う </OPTION><BR>
</SELECT>
方法：<SELECT STYLE="font-size: 100%" SIZE="1" ID="method">
<OPTION VALUE="no"> DFP法 </OPTION><BR>
<OPTION VALUE="yes"> BFGS法 </OPTION><BR>
</SELECT><BR><BR>
目的関数： f(x+αp) = <INPUT ID="func" STYLE="font-size: 100%" TYPE="text" SIZE="70" VALUE="100.0 * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) + (1.0 - (x[0] + k * dx[0])) * (1.0 - (x[0] + k * dx[0]))"><BR><BR>
目的関数の微分： f'(x) = <TEXTAREA ID="dfunc" STYLE="font-size: 110%" COLS="70" ROWS="10">
400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
-200.0 * (x[1] - x[0] * x[0])
</TEXTAREA><BR><BR>
<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR><BR>
結果：<TEXTAREA ID="res" STYLE="font-size: 110%" COLS="40" ROWS="2"></TEXTAREA>
</DIV>
</BODY>
</HTML>
<?php
/**********************************/
/* 準Newton法によるの最小値の計算 */
/* coded by Y.Suganuma */
/**********************************/
$sw = 0;
// データの入力
fscanf(STDIN, "%*s %d %*s %d %*s %d %*s %d", $fun, $n, $max, $opt_1);
fscanf(STDIN, "%*s %d %*s %lf %*s %lf", $method, $eps, $step);
$x = array($n);
$dx = array(0.0, 0.0);
$H = array($n);
$str = trim(fgets(STDIN));
strtok($str, " ");
for ($i1 = 0; $i1 < $n; $i1++) {
$x[$i1] = floatval(strtok(" "));
$H[$i1] = array($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);
}
/********************************************************/
/* 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 : 黄金分割法が失敗 */
/********************************************************/
function DFP_BFGS($method, $opt_1, $max, $n, $eps, $step, &$y, &$x1, $dx, $H, $snx, $dsnx)
{
$count = 0;
$sw = 0;
$x2 = array($n);
$g = array($n);
$g1 = array($n);
$g2 = array($n);
$s = array($n);
$w = array($n);
$y1 = $snx(0.0, $x1, $dx);
$dsnx($x1, $g1);
$H1 = array($n);
$H2 = array($n);
$I = array($n);
for ($i1 = 0; $i1 < $n; $i1++) {
$H1[$i1] = array($n);
$H2[$i1] = array($n);
$I[$i1] = array($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 (abs($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;
return $sw;
}
/*********************************************/
/* 与えられた点xから，dx方向へk*dxだけ進んだ */
/* 点における関数値を計算する */
/*********************************************/
// 関数１
function snx1($k, $x, $dx)
{
$w = array(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;
}
// 関数２
function snx2($k, $x, $dx)
{
$w = array(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;
}
// 関数３
function snx3($k, $x, $dx)
{
$w = array(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;
}
/********************/
/* 微係数を計算する */
/********************/
// 関数１
function dsnx1($x, &$dx)
{
$dx[0] = -2.0 * ($x[0] - 1.0);
$dx[1] = -2.0 * ($x[1] - 2.0);
}
// 関数２
function dsnx2($x, &$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]);
}
// 関数３
function dsnx3($x, &$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]));
}
/****************************************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* snx : 関数値を計算する関数の名前 */
/* return : 結果(w＋y＊dwのy) */
/****************************************************************/
function gold($a, $b, $eps, &$val, &$ind, $max, $w, $dw, $snx)
{
// 初期設定
$tau = (sqrt(5.0) - 1.0) / 2.0;
$x = 0.0;
$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 (abs($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 (abs($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;
}
?>
#*********************************/
# 準Newton法によるの最小値の計算 */
# coded by Y.Suganuma */
#*********************************/
#********************************************/
# 与えられた点xから，dx方向へk*dxだけ進んだ */
# 点における関数値及び微係数を計算する */
#********************************************/
# 関数１
snx1 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(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
x1 * x1 + y1 * y1
else
dx[0] = -2.0 * (x[0] - 1.0)
dx[1] = -2.0 * (x[1] - 2.0)
end
}
# 関数２
snx2 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(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]
100.0 * x1 * x1 + y1 * y1
else
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])
end
}
# 関数３
snx3 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(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])
x1 * x1 + y1 * y1 + z1 * z1
else
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]))
end
}
#***************************************************************/
# 黄金分割法(与えられた方向での最小値) */
# a,b : 初期区間 a < b */
# eps : 許容誤差 */
# val : 関数値 */
# ind : 計算状況 */
# >= 0 : 正常終了(収束回数) */
# = -1 : 収束せず */
# max : 最大試行回数 */
# w : 位置 */
# dw : 傾きの成分 */
# snx : 関数値を計算する関数の名前 */
# return : 結果(w＋y＊dwのy) */
#***************************************************************/
def gold(a, b, eps, val, ind, max, w, dw, &snx)
# 初期設定
tau = (Math.sqrt(5.0) - 1.0) / 2.0
x = 0.0
count = 0
ind[0] = -1
x1 = b - tau * (b - a)
x2 = a + tau * (b - a)
f1 = snx.call(0, x1, w, dw)
f2 = snx.call(0, x2, w, dw)
# 計算
while count < max && ind[0] < 0
count += 1
if f2 > f1
if (b-a).abs() < eps
ind[0] = 0
x = x1
val[0] = f1
else
b = x2
x2 = x1
x1 = a + (1.0 - tau) * (b - a)
f2 = f1
f1 = snx.call(0, x1, w, dw)
end
else
if (b-a).abs() < eps
ind[0] = 0
x = x2
val[0] = f2
f1 = f2
else
a = x1
x1 = x2
x2 = b - (1.0 - tau) * (b - a)
f1 = f2
f2 = snx.call(0, x2, w, dw)
end
end
end
# 収束した場合の処理
if ind[0] == 0
ind[0] = count
fa = snx.call(0, a, w, dw)
fb = snx.call(0, b, w, dw)
if fb < fa
a = b
fa = fb
end
if fa < f1
x = a
val[0] = fa
end
end
return x
end
#*******************************************************/
# DFP法 or BFGS法 */
# method : =0 : DFP法 */
# =1 : BFGS法 */
# opt_1 : =0 : １次元最適化を行わない */
# =1 : １次元最適化を行う */
# max : 最大繰り返し回数 */
# n : 次元 */
# eps : 収束判定条件 */
# step : きざみ幅 */
# y : 最小値 */
# x1 : x(初期値と答え) */
# dx : 関数の微分値 */
# H : Hesse行列の逆行列 */
# snx : 関数値とその微分を計算する関数名 */
# （微分は，その符号を変える) */
# return : >=0 : 正常終了(収束回数) */
# =-1 : 収束せず */
# =-2 : １次元最適化の区間が求まらない */
# =-3 : 黄金分割法が失敗 */
#*******************************************************/
def DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, h, &snx)
count = 0
sw = 0
y2 = Array.new(1)
sw1 = Array.new(1)
x2 = Array.new(n)
g = Array.new(n)
g1 = Array.new(n)
g2 = Array.new(n)
s = Array.new(n)
w = Array.new(n)
y1 = snx.call(0, 0.0, x1, dx)
snx.call(1, 0.0, x1, g1)
h1 = Array.new(n)
h2 = Array.new(n)
i = Array.new(n)
for i1 in 0 ... n
h1[i1] = Array.new(n)
h2[i1] = Array.new(n)
i[i1] = Array.new(n)
for i2 in 0 ... n
h[i1][i2] = 0.0
i[i1][i2] = 0.0
end
h[i1][i1] = 1.0
i[i1][i1] = 1.0
end
while count < max && sw == 0
count += 1
# 方向の計算
for i1 in 0 ... n
dx[i1] = 0.0
for i2 in 0 ... n
dx[i1] -= h[i1][i2] * g1[i2]
end
end
# １次元最適化を行わない
if opt_1 == 0
# 新しい点
for i1 in 0 ... n
x2[i1] = x1[i1] + step * dx[i1]
end
# 新しい関数値
y2[0] = snx.call(0, 0.0, x2, dx)
# １次元最適化を行う
else
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx.call(0, sp, x1, dx)
if f2 > f1
sp = -step
end
for i1 in 0 ... max
f2 = snx.call(0, sp, x1, dx)
if f2 > f1
sw1[0] = 1
break
else
sp *= 2.0
f1 = f2
end
end
# 区間が求まらない
if sw1[0] == 0
sw = -2
# 区間が求まった
else
# 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x1, dx, &snx)
# 黄金分割法が失敗
if sw1[0] < 0
sw = -3
# 黄金分割法が成功
else
# 新しい点
for i1 in 0 ... n
x2[i1] = x1[i1] + k * dx[i1]
end
end
end
end
if sw == 0
# 収束（関数値の変化＜ｅｐｓ）
if (y2[0]-y1).abs() < eps
sw = count
y[0] = y2[0]
for i1 in 0 ... n
x1[i1] = x2[i1]
end
# 関数値の変化が大きい
else
# 傾きの計算
snx.call(1, 0.0, x2, g2)
sm = 0.0
for i1 in 0 ... n
sm += g2[i1] * g2[i1]
end
sm = Math.sqrt(sm)
# 収束（傾き＜ｅｐｓ）
if sm < eps
sw = count
y[0] = y2[0]
for i1 in 0 ... n
x1[i1] = x2[i1]
end
# 収束していない
else
for i1 in 0 ... n
g[i1] = g2[i1] - g1[i1]
s[i1] = x2[i1] - x1[i1]
end
sm1 = 0.0
for i1 in 0 ... n
sm1 += s[i1] * g[i1]
end
# DFP法
if method == 0
for i1 in 0 ... n
w[i1] = 0.0
for i2 in 0 ... n
w[i1] += g[i2] * h[i2][i1]
end
end
sm2 = 0.0
for i1 in 0 ... n
sm2 += w[i1] * g[i1]
end
for i1 in 0 ... n
w[i1] = 0.0
for i2 in 0 ... n
w[i1] += h[i1][i2] * g[i2]
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h1[i1][i2] = w[i1] * g[i2]
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h2[i1][i2] = 0.0
for i3 in 0 ... n
h2[i1][i2] += h1[i1][i3] * h[i3][i2]
end
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h[i1][i2] = h[i1][i2] + s[i1] * s[i2] / sm1 - h2[i1][i2] / sm2
end
end
# BFGS法
else
for i1 in 0 ... n
for i2 in 0 ... n
h1[i1][i2] = i[i1][i2] - s[i1] * g[i2] / sm1
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h2[i1][i2] = 0.0
for i3 in 0 ... n
h2[i1][i2] += h1[i1][i3] * h[i3][i2]
end
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h[i1][i2] = 0.0
for i3 in 0 ... n
h[i1][i2] += h2[i1][i3] * h1[i3][i2]
end
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h[i1][i2] += s[i1] * s[i2] / sm1
end
end
end
y1 = y2[0]
for i1 in 0 ... n
g1[i1] = g2[i1]
x1[i1] = x2[i1]
end
end
end
end
end
if sw == 0
sw = -1
end
return sw
end
# データの入力
str = gets()
a = str.split(" ")
fun = Integer(a[1])
n = Integer(a[3])
max = Integer(a[5])
opt_1 = Integer(a[7])
str = gets();
a = str.split(" ");
method = Integer(a[1])
eps = Float(a[3])
step = Float(a[5])
x = Array.new(n)
dx = Array.new(n)
h = Array.new(n)
y = Array.new(1)
sw = 0
str = gets();
a = str.split(" ");
for i1 in 0 ... n
x[i1] = Float(a[i1+1])
dx[i1] = 0.0
h[i1] = Array.new(n)
end
# 実行
case fun
when 1
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx1)
when 2
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx2)
when 3
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx3)
end
# 結果の出力
if sw < 0
printf(" 収束しませんでした！")
case sw
when -1
printf("（収束回数）\n")
when -2
printf("（１次元最適化の区間）\n")
when -3
printf("（黄金分割法）\n")
end
else
printf(" 結果＝")
for i1 in 0 ... n
printf("%f ", x[i1])
end
printf(" 最小値＝%f 回数＝%d\n", y[0], sw)
end
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
############################################
# 黄金分割法(与えられた方向での最小値)
# a,b : 初期区間 a < b
# eps : 許容誤差
# val : 関数値
# ind : 計算状況
# >= 0 : 正常終了(収束回数)
# = -1 : 収束せず
# max : 最大試行回数
# w : 位置
# dw : 傾きの成分
# snx : 関数値を計算する関数の名前
# return : 結果(w＋y＊dwのy)
# coded by Y.Suganuma
############################################
def gold(a, b, eps, val, ind, max, w, dw, snx) :
# 初期設定
tau = (sqrt(5.0) - 1.0) / 2.0
x = 0.0
count = 0
ind[0] = -1
x1 = b - tau * (b - a)
x2 = a + tau * (b - a)
f1 = snx(x1, w, dw)
f2 = snx(x2, w, dw)
# 計算
while count < max and ind[0] < 0 :
count += 1
if f2 > f1 :
if abs(b-a) < eps :
ind[0] = 0
x = x1
val[0] = f1
else :
b = x2
x2 = x1
x1 = a + (1.0 - tau) * (b - a)
f2 = f1
f1 = snx(x1, w, dw)
else :
if abs(b-a) < eps :
ind[0] = 0
x = x2
val[0] = 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] == 0 :
ind[0] = count
fa = snx(a, w, dw)
fb = snx(b, w, dw)
if fb < fa :
a = b
fa = fb
if fa < f1 :
x = a
val[0] = fa
return x
############################################
# 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 : 黄金分割法が失敗
# coded by Y.Suganuma
############################################
def DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, H, snx, dsnx) :
count = 0
sw = 0
sw1 = np.empty(1, np.int)
y2 = np.empty(1, np.float)
x2 = np.empty(n, np.float)
g = np.empty(n, np.float)
g1 = np.empty(n, np.float)
g2 = np.empty(n, np.float)
s = np.empty(n, np.float)
w = np.empty(n, np.float)
I = np.identity(n, np.float)
H1 = np.empty((n, n), np.float)
H2 = np.empty((n, n), np.float)
y1 = snx(0.0, x1, dx)
dsnx(x1, g1)
while count < max and sw == 0 :
count += 1
# 方向の計算
for i1 in range(0, n) :
dx[i1] = 0.0
for i2 in range(0, n) :
dx[i1] -= H[i1][i2] * g1[i2]
# １次元最適化を行わない
if opt_1 == 0 :
# 新しい点
for i1 in range(0, n) :
x2[i1] = x1[i1] + step * dx[i1]
# 新しい関数値
y2[0] = snx(0.0, x2, dx)
# １次元最適化を行う
else :
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx(sp, x1, dx)
if f2 > f1 :
sp = -step
for i1 in range(0, max) :
f2 = snx(sp, x1, dx)
if f2 > f1 :
sw1[0] = 1
break
else :
sp *= 2.0
f1 = f2
# 区間が求まらない
if sw1[0] == 0 :
sw = -2
# 区間が求まった
else :
# 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x1, dx, snx)
# 黄金分割法が失敗
if sw1[0] < 0 :
sw = -3
# 黄金分割法が成功
else :
# 新しい点
for i1 in range(0, n) :
x2[i1] = x1[i1] + k * dx[i1]
if sw == 0 :
# 収束（関数値の変化＜ｅｐｓ）
if abs(y2[0]-y1) < eps :
sw = count
y[0] = y2[0]
for i1 in range(0, n) :
x1[i1] = x2[i1]
# 関数値の変化が大きい
else :
# 傾きの計算
dsnx(x2, g2)
sm = 0.0
for i1 in range(0, n) :
sm += g2[i1] * g2[i1]
sm = sqrt(sm)
# 収束（傾き＜ｅｐｓ）
if sm < eps :
sw = count
y[0] = y2[0]
for i1 in range(0, n) :
x1[i1] = x2[i1]
# 収束していない
else :
for i1 in range(0, n) :
g[i1] = g2[i1] - g1[i1]
s[i1] = x2[i1] - x1[i1]
sm1 = 0.0
for i1 in range(0, n) :
sm1 += s[i1] * g[i1]
# DFP法
if method == 0 :
for i1 in range(0, n) :
w[i1] = 0.0
for i2 in range(0, n) :
w[i1] += g[i2] * H[i2][i1]
sm2 = 0.0
for i1 in range(0, n) :
sm2 += w[i1] * g[i1]
for i1 in range(0, n) :
w[i1] = 0.0
for i2 in range(0, n) :
w[i1] += H[i1][i2] * g[i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H1[i1][i2] = w[i1] * g[i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H2[i1][i2] = 0.0
for i3 in range(0, n) :
H2[i1][i2] += H1[i1][i3] * H[i3][i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H[i1][i2] = H[i1][i2] + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2
# BFGS法
else :
for i1 in range(0, n) :
for i2 in range(0, n) :
H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1
for i1 in range(0, n) :
for i2 in range(0, n) :
H2[i1][i2] = 0.0
for i3 in range(0, n) :
H2[i1][i2] += H1[i1][i3] * H[i3][i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H[i1][i2] = 0.0
for i3 in range(0, n) :
H[i1][i2] += H2[i1][i3] * H1[i3][i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H[i1][i2] += s[i1] * s[i2] / sm1
y1 = y2[0]
for i1 in range(0, n) :
g1[i1] = g2[i1]
x1[i1] = x2[i1]
if sw == 0 :
sw = -1
return sw
############################################
# 準Newton法によるの最小値の計算
# coded by Y.Suganuma
############################################
############################################
# 与えられた点xから，dx方向へk*dxだけ進んだ
# 点における関数値を計算する
############################################
# 関数１
def snx1(k, x, dx) :
w = np.empty(2, np.float)
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
# 関数２
def snx2(k, x, dx) :
w = np.empty(2, np.float)
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
# 関数３
def snx3(k, x, dx) :
w = np.empty(2, np.float)
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
############################################
# 微係数を計算する
############################################
# 関数１
def dsnx1(x, dx) :
dx[0] = -2.0 * (x[0] - 1.0)
dx[1] = -2.0 * (x[1] - 2.0)
# 関数２
def dsnx2(x, 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])
# 関数３
def dsnx3(x, 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]))
# データの入力
sw = 0
line = sys.stdin.readline()
s = line.split()
fun = int(s[1])
n = int(s[3])
max = int(s[5])
opt_1 = int(s[7])
line = sys.stdin.readline()
s = line.split()
method = int(s[1])
eps = float(s[3])
step = float(s[5])
x = np.empty(n, np.float)
dx = np.empty(n, np.float)
H = np.identity(n, np.float)
y = np.empty(1, np.float)
line = sys.stdin.readline()
s = line.split()
for i1 in range(0, n) :
x[i1] = float(s[i1+1])
# 実行
if fun == 1 :
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx1, dsnx1)
elif fun == 2 :
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx2, dsnx2)
else :
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx3, dsnx3)
# 結果の出力
if sw < 0 :
print(" 収束しませんでした！", end="")
if sw == -1 :
print("（収束回数）")
elif sw == -2 :
print("（１次元最適化の区間）")
else :
print("（黄金分割法）")
else :
print(" 結果＝", end="")
for i1 in range(0, n) :
print(str(x[i1]) + " ", end="")
print(" 最小値＝" + str(y[0]) + " 回数＝" + str(sw))
/**********************************/
/* 準Newton法によるの最小値の計算 */
/* coded by Y.Suganuma */
/**********************************/
using System;
class Program
{
static void Main()
{
Test1 ts = new Test1();
}
}
class Test1
{
public Test1()
{
// データの入力
// 1 行目
string[] str = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
int fun = int.Parse(str[1]);
int n = int.Parse(str[3]);
int max = int.Parse(str[5]);
int opt_1 = int.Parse(str[7]);
// 2 行目
str = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
int method = int.Parse(str[1]);
double eps = double.Parse(str[3]);
double step = double.Parse(str[5]);
// 3 行目
double[][] H = new double [n][];
double[] x = new double [n];
double[] dx = new double [n];
double y = 0.0;
str = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
for (int i1 = 0; i1 < n; i1++) {
x[i1] = double.Parse(str[i1+1]);
H[i1] = new double [n];
}
// 実行
int sw;
if (fun == 1)
sw = DFP_BFGS(method, opt_1, max, n, eps, step, ref y, x, dx, H, snx1, dsnx1);
else if (fun == 2)
sw = DFP_BFGS(method, opt_1, max, n, eps, step, ref y, x, dx, H, snx2, dsnx2);
else
sw = DFP_BFGS(method, opt_1, max, n, eps, step, ref y, x, dx, H, snx3, dsnx3);
// 結果の出力
if (sw < 0) {
Console.Write(" 収束しませんでした！");
switch (sw) {
case -1:
Console.WriteLine("（収束回数）");
break;
case -2:
Console.WriteLine("（１次元最適化の区間）");
break;
case -3:
Console.WriteLine("（黄金分割法）");
break;
}
}
else {
Console.Write(" 結果＝");
for (int i1 = 0; i1 < n; i1++)
Console.Write(x[i1] + " ");
Console.WriteLine(" 最小値＝" + y + " 回数＝" + sw);
}
}
// 関数１の値の計算
double snx1(double k, double[] x, double[] dx)
{
double[] w = new double [2];
w[0] = x[0] + k * dx[0];
w[1] = x[1] + k * dx[1];
double x1 = w[0] - 1.0;
double y1 = w[1] - 2.0;
double y = x1 * x1 + y1 * y1;
return y;
}
// 関数２の値の計算
double snx2(double k, double[] x, double[] dx)
{
double[] w = new double [2];
w[0] = x[0] + k * dx[0];
w[1] = x[1] + k * dx[1];
double x1 = w[1] - w[0] * w[0];
double y1 = 1.0 - w[0];
double y = 100.0 * x1 * x1 + y1 * y1;
return y;
}
// 関数３の値の計算
double snx3(double k, double[] x, double[] dx)
{
double[] w = new double [2];
w[0] = x[0] + k * dx[0];
w[1] = x[1] + k * dx[1];
double x1 = 1.5 - w[0] * (1.0 - w[1]);
double y1 = 2.25 - w[0] * (1.0 - w[1] * w[1]);
double z1 = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
double y = x1 * x1 + y1 * y1 + z1 * z1;
return y;
}
// 関数１の微分
int dsnx1(double[] x, double[] dx)
{
dx[0] = -2.0 * (x[0] - 1.0);
dx[1] = -2.0 * (x[1] - 2.0);
return 0;
}
// 関数２の微分
int 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]);
return 0;
}
// 関数３の微分
int 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]));
return 0;
}
/********************************************************/
/* DFP法 or BFGS法 */
/* method : =0 : DFP法 */
/* =1 : BFGS法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x1 : x(初期値と答え) */
/* dx : 関数の微分値 */
/* H : Hesse行列の逆行列 */
/* fn : 関数値を計算する関数 */
/* dfn : 関数の微分値を計算する関数 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
int DFP_BFGS(int method, int opt_1, int max, int n, double eps, double step,
ref double y, double[] x1, double[] dx, double[][] H,
Func<double, double[], double[], double> fn,
Func<double[], double[], int> dfn)
{
double[] x2 = new double [n];
double[] g = new double [n];
double[] g1 = new double [n];
double[] g2 = new double [n];
double[] s = new double [n];
double[] w = new double [n];
double y1 = fn(0.0, x1, dx);
double y2 = 0.0;
dfn(x1, g1);
double[][] H1 = new double [n][];
double[][] H2 = new double [n][];
double[][] I = new double [n][];
for (int i1 = 0; i1 < n; i1++) {
H1[i1] = new double [n];
H2[i1] = new double [n];
I[i1] = new double [n];
for (int 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;
}
int count = 0;
int sw = 0;
while (count < max && sw == 0) {
count++;
// 方向の計算
for (int i1 = 0; i1 < n; i1++) {
dx[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++)
dx[i1] -= H[i1][i2] * g1[i2];
}
// １次元最適化を行わない
if (opt_1 == 0) {
// 新しい点
for (int i1 = 0; i1 < n; i1++)
x2[i1] = x1[i1] + step * dx[i1];
// 新しい関数値
y2 = fn(0.0, x2, dx);
}
// １次元最適化を行う
else {
// 区間を決める
int sw1 = 0;
double f1 = y1;
double sp = step;
double f2 = fn(sp, x1, dx);
if (f2 > f1)
sp = -step;
for (int i1 = 0; i1 < max && sw1 == 0; i1++) {
f2 = fn(sp, x1, dx);
if (f2 > f1)
sw1 = 1;
else {
sp *= 2.0;
f1 = f2;
}
}
// 区間が求まらない
if (sw1 == 0)
sw = -2;
// 区間が求まった
else {
// 黄金分割法
double k = gold(0.0, sp, eps, ref y2, ref sw1, max, x1, dx, fn);
// 黄金分割法が失敗
if (sw1 < 0)
sw = -3;
// 黄金分割法が成功
else {
// 新しい点
for (int i1 = 0; i1 < n; i1++)
x2[i1] = x1[i1] + k * dx[i1];
}
}
}
if (sw == 0) {
// 収束（関数値の変化＜ｅｐｓ）
if (Math.Abs(y2-y1) < eps) {
sw = count;
y = y2;
for (int i1 = 0; i1 < n; i1++)
x1[i1] = x2[i1];
}
// 関数値の変化が大きい
else {
// 傾きの計算
dfn(x2, g2);
double sm = 0.0;
for (int i1 = 0; i1 < n; i1++)
sm += g2[i1] * g2[i1];
sm = Math.Sqrt(sm);
// 収束（傾き＜ｅｐｓ）
if (sm < eps) {
sw = count;
y = y2;
for (int i1 = 0; i1 < n; i1++)
x1[i1] = x2[i1];
}
// 収束していない
else {
for (int i1 = 0; i1 < n; i1++) {
g[i1] = g2[i1] - g1[i1];
s[i1] = x2[i1] - x1[i1];
}
double sm1 = 0.0;
for (int i1 = 0; i1 < n; i1++)
sm1 += s[i1] * g[i1];
// DFP法
if (method == 0) {
for (int i1 = 0; i1 < n; i1++) {
w[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++)
w[i1] += g[i2] * H[i2][i1];
}
double sm2 = 0.0;
for (int i1 = 0; i1 < n; i1++)
sm2 += w[i1] * g[i1];
for (int i1 = 0; i1 < n; i1++) {
w[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++)
w[i1] += H[i1][i2] * g[i2];
}
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
H1[i1][i2] = w[i1] * g[i2];
}
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
H2[i1][i2] = 0.0;
for (int i3 = 0; i3 < n; i3++)
H2[i1][i2] += H1[i1][i3] * H[i3][i2];
}
}
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
H[i1][i2] = H[i1][i2] + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2;
}
}
// BFGS法
else {
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1;
}
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
H2[i1][i2] = 0.0;
for (int i3 = 0; i3 < n; i3++)
H2[i1][i2] += H1[i1][i3] * H[i3][i2];
}
}
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
H[i1][i2] = 0.0;
for (int i3 = 0; i3 < n; i3++)
H[i1][i2] += H2[i1][i3] * H1[i3][i2];
}
}
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
H[i1][i2] += s[i1] * s[i2] / sm1;
}
}
y1 = y2;
for (int i1 = 0; i1 < n; i1++) {
g1[i1] = g2[i1];
x1[i1] = x2[i1];
}
}
}
}
}
if (sw == 0)
sw = -1;
return sw;
}
/******************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* fn : 関数値を計算する関数 */
/* return : 結果(w＋y＊dwのy) */
/******************************************/
double gold(double a, double b, double eps, ref double val, ref int ind,
int max, double[] w, double[] dw,
Func<double, double[], double[], double> fn)
{
// 初期設定
double tau = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0;
double x1 = b - tau * (b - a);
double x2 = a + tau * (b - a);
double f1 = fn(x1, w, dw);
double f2 = fn(x2, w, dw);
int count = 0;
ind = -1;
// 計算
while (count < max && ind < 0) {
count += 1;
if (f2 > f1) {
if (Math.Abs(b-a) < eps) {
ind = 0;
x = x1;
val = f1;
}
else {
b = x2;
x2 = x1;
x1 = a + (1.0 - tau) * (b - a);
f2 = f1;
f1 = fn(x1, w, dw);
}
}
else {
if (Math.Abs(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 = fn(x2, w, dw);
}
}
}
// 収束した場合の処理
if (ind == 0) {
ind = count;
double fa = fn(a, w, dw);
double fb = fn(b, w, dw);
if (fb < fa) {
a = b;
fa = fb;
}
if (fa < f1) {
x = a;
val = fa;
}
}
return x;
}
}
''''''''''''''''''''''''''''''''''
' 準Newton法によるの最小値の計算 '
' coded by Y.Suganuma '
''''''''''''''''''''''''''''''''''
Imports System.Text.RegularExpressions
Module Test
Sub Main()
' データの入力
' 1 行目
Dim MS As Regex = New Regex("\s+")
Dim str() As String = MS.Split(Console.ReadLine().Trim())
Dim fun As Integer = Integer.Parse(str(1))
Dim n As Integer = Integer.Parse(str(3))
Dim max As Integer = Integer.Parse(str(5))
Dim opt_1 As Integer = Integer.Parse(str(7))
' 2 行目
str = MS.Split(Console.ReadLine().Trim())
Dim method As Integer = Integer.Parse(str(1))
Dim eps As Double = Double.Parse(str(3))
Dim step1 As Double = Double.Parse(str(5))
' 3 行目
Dim H(n,n) As Double
Dim x(n) As Double
Dim dx(n) As Double
Dim y As Double = 0.0
str = MS.Split(Console.ReadLine().Trim())
For i1 As Integer = 0 To n-1
x(i1) = Double.Parse(str(i1+1))
Next
' 関数１の値の計算（ラムダ式）
Dim snx1 = Function(k0, x0, dx0) As Double
Dim w(2) As Double
w(0) = x0(0) + k0 * dx0(0)
w(1) = x0(1) + k0 * dx0(1)
Dim x1 As Double = w(0) - 1.0
Dim y1 As Double = w(1) - 2.0
Dim y0 As Double = x1 * x1 + y1 * y1
Return y0
End Function
' 関数２の値の計算（ラムダ式）
Dim snx2 = Function(k0, x0, dx0) As Double
Dim w(2) As Double
w(0) = x0(0) + k0 * dx0(0)
w(1) = x0(1) + k0 * dx0(1)
Dim x1 As Double = w(1) - w(0) * w(0)
Dim y1 As Double = 1.0 - w(0)
Dim y0 As Double = 100.0 * x1 * x1 + y1 * y1
Return y0
End Function
' 関数３の値の計算（ラムダ式）
Dim snx3 = Function(k0, x0, dx0) As Double
Dim w(2) As Double
w(0) = x0(0) + k0 * dx0(0)
w(1) = x0(1) + k0 * dx0(1)
Dim x1 As Double = 1.5 - w(0) * (1.0 - w(1))
Dim y1 As Double = 2.25 - w(0) * (1.0 - w(1) * w(1))
Dim z1 As Double = 2.625 - w(0) * (1.0 - w(1) * w(1) * w(1))
Dim y0 As Double = x1 * x1 + y1 * y1 + z1 * z1
Return y0
End Function
' 関数１の微分（ラムダ式）
Dim dsnx1 = Function(x0, dx0) As Integer
dx0(0) = -2.0 * (x0(0) - 1.0)
dx0(1) = -2.0 * (x0(1) - 2.0)
Return 0
End Function
' 関数２の微分（ラムダ式）
Dim dsnx2 = Function(x0, dx0) As Integer
dx0(0) = 400.0 * x0(0) * (x0(1) - x0(0) * x0(0)) + 2.0 * (1.0 - x0(0))
dx0(1) = -200.0 * (x0(1) - x0(0) * x0(0))
Return 0
End Function
' 関数３の微分（ラムダ式）
Dim dsnx3 = Function(x0, dx0) As Integer
dx0(0) = 2.0 * (1.0 - x0(1)) * (1.5 - x0(0) * (1.0 - x0(1))) +
2.0 * (1.0 - x0(1) * x0(1)) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) +
2.0 * (1.0 - x0(1) * x0(1) * x0(1)) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
dx0(1) = -2.0 * x0(0) * (1.5 - x0(0) * (1.0 - x0(1))) -
4.0 * x0(0) * x0(1) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) -
6.0 * x0(0) * x0(1) * x0(1) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
Return 0
End Function
' 実行
Dim sw As Integer
If fun = 1
sw = DFP_BFGS(method, opt_1, max, n, eps, step1, y, x, dx, H, snx1, dsnx1)
ElseIf fun = 2
sw = DFP_BFGS(method, opt_1, max, n, eps, step1, y, x, dx, H, snx2, dsnx2)
Else
sw = DFP_BFGS(method, opt_1, max, n, eps, step1, y, x, dx, H, snx3, dsnx3)
End If
' 結果の出力
If sw < 0
Console.Write(" 収束しませんでした！")
If sw = -1
Console.WriteLine("（収束回数）")
ElseIf sw = -2
Console.WriteLine("（１次元最適化の区間）")
ElseIf sw = -3
Console.WriteLine("（黄金分割法）")
End If
Else
Console.Write(" 結果＝")
For i1 As Integer = 0 To n-1
Console.Write(x(i1) & " ")
Next
Console.WriteLine(" 最小値＝" & y & " 回数＝" & sw)
End If
End Sub
''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' DFP法 or BFGS法 '
' method : =0 : DFP法 '
' =1 : BFGS法 '
' opt_1 : =0 : １次元最適化を行わない '
' =1 : １次元最適化を行う '
' max : 最大繰り返し回数 '
' n : 次元 '
' eps : 収束判定条件 '
' step1 : きざみ幅 '
' y : 最小値 '
' x1 : x(初期値と答え) '
' dx : 関数の微分値 '
' H : Hesse行列の逆行列 '
' fn : 関数値を計算する関数 '
' dfn : 関数の微分値を計算する関数 '
' return : >=0 : 正常終了(収束回数) '
' =-1 : 収束せず '
' =-2 : １次元最適化の区間が求まらない '
' =-3 : 黄金分割法が失敗 '
''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function DFP_BFGS(method As Integer, opt_1 As Integer, max As Integer, n As Integer,
eps As Double, step1 As Double, ByRef y As Double, x1() As Double,
dx() As Double, H(,) As Double,
fn As Func(Of Double, Double(), Double(), Double),
dfn As Func(Of Double(), Double(), Integer))
Dim x2(n) As Double
Dim g(n) As Double
Dim g1(n) As Double
Dim g2(n) As Double
Dim s(n) As Double
Dim w(n) As Double
Dim y1 As Double = fn(0.0, x1, dx)
Dim y2 As Double = 0.0
dfn(x1, g1)
Dim H1(n,n) As Double
Dim H2(n,n) As Double
Dim I(n,n) As Double
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H(i1,i2) = 0.0
I(i1,i2) = 0.0
Next
H(i1,i1) = 1.0
I(i1,i1) = 1.0
Next
Dim count As Integer = 0
Dim sw As Integer = 0
Do While count < max and sw = 0
count += 1
' 方向の計算
For i1 As Integer = 0 To n-1
dx(i1) = 0.0
For i2 As Integer = 0 To n-1
dx(i1) -= H(i1,i2) * g1(i2)
Next
Next
' １次元最適化を行わない
If opt_1 = 0
' 新しい点
For i1 As Integer = 0 To n-1
x2(i1) = x1(i1) + step1 * dx(i1)
Next
' 新しい関数値
y2 = fn(0.0, x2, dx)
' １次元最適化を行う
Else
' 区間を決める
Dim sw1 As Integer = 0
Dim f1 As Double = y1
Dim sp As Double = step1
Dim f2 As Double = fn(sp, x1, dx)
If f2 > f1
sp = -step1
End If
Dim ii As Integer = 0
Do While ii < max and sw1 = 0
f2 = fn(sp, x1, dx)
If f2 > f1
sw1 = 1
Else
sp *= 2.0
f1 = f2
End If
ii += 1
Loop
' 区間が求まらない
If sw1 = 0
sw = -2
' 区間が求まった
Else
' 黄金分割法
Dim k As Double = gold(0.0, sp, eps, y2, sw1, max, x1, dx, fn)
' 黄金分割法が失敗
If sw1 < 0
sw = -3
' 黄金分割法が成功
Else
' 新しい点
For i1 As Integer = 0 To n-1
x2(i1) = x1(i1) + k * dx(i1)
Next
End If
End If
End If
If sw = 0
' 収束（関数値の変化＜ｅｐｓ）
If Math.Abs(y2-y1) < eps
sw = count
y = y2
For i1 As Integer = 0 To n-1
x1(i1) = x2(i1)
Next
' 関数値の変化が大きい
Else
' 傾きの計算
dfn(x2, g2)
Dim sm As Double = 0.0
For i1 As Integer = 0 To n-1
sm += g2(i1) * g2(i1)
Next
sm = Math.Sqrt(sm)
' 収束（傾き＜ｅｐｓ）
If sm < eps
sw = count
y = y2
For i1 As Integer = 0 To n-1
x1(i1) = x2(i1)
Next
' 収束していない
Else
For i1 As Integer = 0 To n-1
g(i1) = g2(i1) - g1(i1)
s(i1) = x2(i1) - x1(i1)
Next
Dim sm1 As Double = 0.0
For i1 As Integer = 0 To n-1
sm1 += s(i1) * g(i1)
Next
' DFP法
If method = 0
For i1 As Integer = 0 To n-1
w(i1) = 0.0
For i2 As Integer = 0 To n-1
w(i1) += g(i2) * H(i2,i1)
Next
Next
Dim sm2 As Double = 0.0
For i1 As Integer = 0 To n-1
sm2 += w(i1) * g(i1)
Next
For i1 As Integer = 0 To n-1
w(i1) = 0.0
For i2 As Integer = 0 To n-1
w(i1) += H(i1,i2) * g(i2)
Next
Next
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H1(i1,i2) = w(i1) * g(i2)
Next
Next
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H2(i1,i2) = 0.0
For i3 As Integer = 0 To n-1
H2(i1,i2) += H1(i1,i3) * H(i3,i2)
Next
Next
Next
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H(i1,i2) = H(i1,i2) + s(i1) * s(i2) / sm1 - H2(i1,i2) / sm2
Next
Next
' BFGS法
Else
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H1(i1,i2) = I(i1,i2) - s(i1) * g(i2) / sm1
Next
Next
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H2(i1,i2) = 0.0
For i3 As Integer = 0 To n-1
H2(i1,i2) += H1(i1,i3) * H(i3,i2)
Next
Next
Next
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H(i1,i2) = 0.0
For i3 As Integer = 0 To n-1
H(i1,i2) += H2(i1,i3) * H1(i3,i2)
Next
Next
Next
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H(i1,i2) += s(i1) * s(i2) / sm1
Next
Next
End If
y1 = y2
For i1 As Integer = 0 To n-1
g1(i1) = g2(i1)
x1(i1) = x2(i1)
Next
End If
End If
End If
Loop
If sw = 0
sw = -1
End If
Return sw
End Function
''''''''''''''''''''''''''''''''''''''''''
' 黄金分割法(与えられた方向での最小値) '
' a,b : 初期区間 a < b '
' eps : 許容誤差 '
' val : 関数値 '
' ind : 計算状況 '
' >= 0 : 正常終了(収束回数) '
' = -1 : 収束せず '
' max : 最大試行回数 '
' w : 位置 '
' dw : 傾きの成分 '
' fn : 関数値を計算する関数 '
' return : 結果(w＋y＊dwのy) '
''''''''''''''''''''''''''''''''''''''''''
Function gold(a As Double, b As Double, eps As Double, ByRef val As Double,
ByRef ind As Integer, max As Integer, w() As Double, dw() As Double,
fn As Func(Of Double, Double(), Double(), Double))
' 初期設定
Dim tau As Double = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0
Dim x1 As Double = b - tau * (b - a)
Dim x2 As Double = a + tau * (b - a)
Dim f1 As Double = fn(x1, w, dw)
Dim f2 As Double = fn(x2, w, dw)
Dim count As Integer = 0
ind = -1
' 計算
Do While count < max and ind < 0
count += 1
If f2 > f1
If Math.Abs(b-a) < eps
ind = 0
x = x1
val = f1
Else
b = x2
x2 = x1
x1 = a + (1.0 - tau) * (b - a)
f2 = f1
f1 = fn(x1, w, dw)
End If
Else
If Math.Abs(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 = fn(x2, w, dw)
End If
End If
Loop
' 収束した場合の処理
If ind = 0
ind = count
Dim fa As Double = fn(a, w, dw)
Dim fb As Double = fn(b, w, dw)
If fb < fa
a = b
fa = fb
End If
If fa < f1
x = a
val = fa
End If
End If
Return x
End Function
End Module
| 情報学部 | 菅沼ホーム | 目次 | 索引 |