拟牛顿法(变尺度法)DFP算法的cc 源码

*direction);
//0.618 法线搜索
double line_search(double (*pf)(double *x), int n, double *start, double
*direction);
//解析法线Leabharlann Baidu索
double DFP(double (*pf)(double *x),
g_norm=0.0;
for(i=0;i<n;i++)
g_norm=g_norm+g1[i]*g1[i];
g_norm=sqrt(g_norm);
//cout<<k<<"
"<<x0[0]<<"
"<<g_norm<<"\n";
if (g_norm<e)
{
for(i=0;i<n;i++)
min_point[i]=x1[i];
int n,
double *min_point);
//无约束变尺度法
//梯度计算模块 //参数：指向目标函数的指针，变量个数，求梯度的点，结果 void comput_grad(double (*pf)(double *x),
int n, double *point, double *grad) { double h=1E-3; int i; double *temp; temp = new double[n]; for(i=1;i<=n;i++) { temp[i-1]=point[i-1]; } for(i=1;i<=n;i++) { temp[i-1]+=0.5*h; grad[i-1]=4*pf(temp)/(3*h); temp[i-1]-=h; grad[i-1]-=4*pf(temp)/(3*h); temp[i-1]+=(3*h/2); grad[i-1]-=(pf(temp)/(6*h)); temp[i-1]-=(2*h); grad[i-1]+=(pf(temp)/(6*h)); temp[i-1]=point[i-1]; } delete[] temp; }
delete[] H[i]; delete[] tempH[i];
return pf(min_point); } for(i=0;i<n;i++) {
dx[i]=x1[i]-x0[i]; dg[i]=g1[i]-g0[i]; }
//////////////////求 Hk+1 的矩阵运算
//g*H,H*g for(i=0;i<n;i++) {
delete[] grad; delete[] temp_point; return(b); } if( diver_b<-1E-15 ) { a=b; diver_a=diver_b; step=2*step;
} }while( diver_b<=1E-15 );
for(i=1;i<=n;i++) temp_point[i-1]=start[i-1]+a*direction[i-1];
"<<x0[1]<<"
delete[] g0; delete[] g1; delete[] dg; delete[] p; delete[] x0; delete[] x1; delete[] dx; for (i=0; i<n; i++) delete []H; for (i=0; i<n; i++) delete []tempH; delete[] gH; delete[] Hg;
temp_point[i-1]=start[i-1]+t*direction[i-1];
value_t=(*pf)(temp_point);
comput_grad(pf,n,temp_point,grad);
diver_t=0;
for(i=1;i<=n;i++)
diver_t=diver_t+grad[i-1]*direction[i-1];
delete[] grad; delete[] temp_point; return(t);
}
//无约束变尺度法 DFP 函数声明 // //参数：pf 指向目标函数的指针，n 变量个数，min_point 接受初始点、存放结果 //返回：极小点处函数值 // double DFP(
double (*pf)(double *x), int n, double *min_point ) { int i,j; int k=0; double e=1E-5; double g_norm;
num1=num1+dx[i]*dg[i]; num2=num2+gH[i]*dg[i]; }
//tempH[i][j] for(i=0;i<n;i++)
for(j=0;j<n;j++) tempH[i][j]=0.0;
for(i=0;i<n;i++) {
for(j=0;j<n;j++) {
tempH[i][j]=tempH[i][j]+H[i][j]; tempH[i][j]=tempH[i][j]+dx[i]*dx[j]/num1; tempH[i][j]=tempH[i][j]-Hg[i]*gH[j]/num2; } }
//一维搜索模块 //参数：指向目标函数的指针，变量个数，出发点，搜索方向 //返回：最优步长 double line_search(
double (*pf)(double *x), int n, double *start, double *direction) { int i; double step=0.001; double a=0,value_a,diver_a; double b,value_b,diver_b; double t,value_t,diver_t; double s,z,w; double *grad,*temp_point;
delete[] H[i]; delete[] tempH[i];
return pf(min_point); }
for(i=0;i<n;i++)
p[i]=-g0[i];
do
{
t=line_search(pf,n,x0,p);
for(i=0;i<n;i++)
x1[i]=x0[i]+t*p[i];
comput_grad(pf,n,x1,g1);
拟牛顿法（变尺度法）DFP 算法的 c/c++源码 #include "iostream.h" #include "math.h"
void comput_grad(double (*pf)(double *x), int n, double *point, double *grad);
//计算梯度
double line_search1(double (*pf)(double *x), int n, double *start, double
{
p[i]=p[i]-H[i][j]*g1[j];
}
}
for(i=0;i<n;i++) {
g0[i]=g1[i]; x0[i]=x1[i]; } k=k+1; }while(g_norm>e);
for(i=0;i<n;i++)
min_point[i]=x1[i];
delete[] g0; delete[] g1; delete[] dg; delete[] p; delete[] x0; delete[] x1; delete[] dx; for (i=0; i<n; i++) delete []H; for (i=0; i<n; i++) delete []tempH; delete[] gH; delete[] Hg;
b=a+step; for(i=1;i<=n;i++)
temp_point[i-1]=start[i-1]+b*direction[i-1]; comput_grad(pf,n,temp_point,grad); diver_b=0; for(i=1;i<=n;i++)
diver_b=diver_b+grad[i-1]*direction[i-1]; if( fabs(diver_b)<1E-10 ) {
if(diver_t>1E-6)
{
b=t;
value_b=value_t;
diver_b=diver_t;
}
else if(diver_t<-1E-6)
{
a=t;
value_a=value_t;
diver_a=diver_t;
}
else break;
}while( (fabs(diver_t)>=1E-6) && (fabs(b-a)>1E-6) );
grad=new double[n]; temp_point=new double[n]; comput_grad(pf,n,start,grad); diver_a=0; for(i=1;i<=n;i++)
diver_a=diver_a+grad[i-1]*direction[i-1]; do {
for(i=0;i<n;i++) for(j=0;j<n;j++)
{
if(i==j)
H[i][j]=1.0;
else
H[i][j]=0.0;
tempH[i][j]=0.0;
}
// H0=I
for(i=0;i<n;i++) x0[i]=min_point[i];
comput_grad(pf,n,x0,g0);
for(i=0;i<n;i++) {
for(j=0;j<n;j++) {
H[i][j]=tempH[i][j]; } } /////////////////////////////
//P
for(i=0;i<n;i++)
p[i]=0.0;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
gH[i]=0.0; Hg[i]=0.0; } for(i=0;i<n;i++) { for(j=0;j<n;j++) {
gH[i]=gH[i]+dg[j]*H[j][i]; //Hg[i]=Hg[i]+H[i][j]*dg[j];
Hg[i]=gH[i]; } }
//num1,num2 num1=0.0; num2=0.0; for(i=0;i<n;i++) {
double *g0=new double[n]; double *g1=new double[n]; double *dg=new double[n];
//梯度
double *p=new double[n]; double t;
//搜索方向 =-g //一维搜索步长
double *x0=new double[n]; double *x1=new double[n]; double *dx=new double[n];
double **H=new double*[n]; for (i=0; i<n; i++)
H[i] = new double[n];
double **tempH=new double*[n]; for (i=0; i<n; i++)
tempH[i] = new double[n];
double *gH=new double[n]; double *Hg=new double[n]; double num1; double num2;
w=sqrt(
fabs(z*z-diver_a*diver_b)
);
//////////////////!!!!!!!!!!!!!!!!!!!!!!
t=a+(w-z-diver_a)*(b-a)/(diver_b-diver_a+2*w);
value_b=(*pf)(temp_point);
for(i=1;i<=n;i++)
value_a=(*pf)(temp_point); for(i=1;i<=n;i++)
temp_point[i-1]=start[i-1]+b*direction[i-1]; value_b=(*pf)(temp_point);
do
{
s=3*(value_b-value_a)/(b-a);
z=s-diver_a-diver_b;
g_norm=0.0;
for(i=0;i<n;i++)
g_norm=g_norm+g0[i]*g0[i];
g_norm=sqrt(g_norm);
if (g_norm<e)
{
for(i=0;i<n;i++)
min_point[i]=x0[i];
delete[] g0; delete[] g1; delete[] dg; delete[] p; delete[] x0; delete[] x1; delete[] dx; for (i=0; i<n; i++) delete []H; for (i=0; i<n; i++) delete []tempH; delete[] gH; delete[] Hg;