bfgs.cpp

共轭梯度的几种方法对标准函数的测试结果比较 附比较结果

#include<iostream.h>
#include<time.h>
#include<math.h>
#define N 2
/*---subroutine calf(n,x,f)---*/
void CALF(int n,double x[],double &f){
	int i;
	double c=100.0;
	f=0.0;
	for(i=0;i<n-1;i++)
        f=f+c*(x[i+1]-x[i]*x[i])*(x[i+1]-x[i]*x[i])+(1.0-x[i])*(1.0-x[i]);
}//calf
/*------Gradient subroutine calg(n,x,g)--------*/
void CALG(int n,double x[],double g[]){
	int i;
	double c=100.0;
	g[0]=-4.0*c*x[0]*(x[1]-x[0]*x[0])-2.0*(1.0-x[0]);
	for(i=1;i<n-1;i++)
        g[i]=2.0*c*(x[i]-x[i-1]*x[i-1])-4.0*c*x[i]*(x[i+1]-x[i]*x[i])-2.0*(1.0-x[i]);
	g[n-1]=2.0*c*(x[n-1]-x[n-2]*x[n-2]);
}//calg
void DOIT1(double H[N][N],double g[N],double p[N]){
	//DOIT1
	int i,j;
	for(i=0;i<N;i++)
		p[i]=0;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			p[i]+=H[i][j]*g[j];

}//DOIT1
double DOIT(double m1[N],double m2[N]){
	//DOIT
	double sum=0;
	for(int i=0;i<N;i++)
		sum+=m1[i]*m2[i];
	return sum;
}//DOIT
void DOIT2(double y1[N],double y2[N],double temp[N][N]){
	//DOIT2
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp[i][j]=y1[i]*y2[j];
}//DOIT2
void DOIT3(double temp1[N][N],double temp2[N][N],double temp3[N][N]){
	//DOIT3
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp3[i][j]=0;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			for(int k=0;k<N;k++)
			temp3[i][j]+=temp1[i][k]*temp2[k][j];
}//DOIT3
void DOIT4(double temp1[N],double temp2[N][N],double temp3[N]){
	//DOIT4
	int i,j;
	for(i=0;i<N;i++)
		temp3[i]=0;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp3[i]+=temp1[j]*temp2[i][j];
}
int if_g(double x[N]){
	double g[N];
	CALG(N,x,g);
	int i=0;
	double e=0.001;
	for(i=0;i<N;i++)
		if(fabs(g[i])>e)return 0;
		return 1;
}//if_g
double DOIT1(int n,double m1[N],double m2[N]){
	double sum=0;
	for(int i=0;i<n;i++)
		sum+=m1[i]*m2[i];
	return sum;
}//DOIT
//
//      SUBROUTINE FOR INEXACT LINE SEARCH
//
//           PROGRAMMED BY Y. YUAN (1985), modified by Q. NI.
//
double DDOT1(int n,double m1[N],double m2[N]){
	double sum=0;
	for(int i=0;i<n;i++)
		sum+=m1[i]*m2[i];
	return sum;
}//DOIT
//
//      SUBROUTINE FOR INEXACT LINE SEARCH
//
//           
//
void INEXLS(int n,double X[N],double D[N],double &ALPHA,int &count1,int &count2){
	//INEXLS
//      int  A21;
//	  double a=0.1,b=1000;
//	  double XX;
	  double F,F1,g1,g;
	  double alpha1=0,alpha2=1000;
	  double alpha;
//	  double TRY;
//	  double DG,DGO,DG1,DG2;
//	  double T1,T2,T3,T4;
//	  double A1,A2;
//	  double RATIO1,RATIO2;
	  double G[N],G1[N];
	  double p=0.1,t=0.6;
//	  int INFO;
//	  int NCALF=0,NCALG=0,NCALFN;
	  int i;
////////////////////step1//////////////////////////
      ALPHA=0.5;
      CALF(N,X,F1);
	  count1++;
	  CALG(N,X,G1);
	  count2++;
      g1=DOIT(G1,D);
	  //step2//////
	  for(i=0;i<N;i++)
		  X[i]=X[i]+ALPHA*D[i];
	  CALF(N,X,F);
	  count1++;
	  if(F-F1<=p*ALPHA*g1)goto step3;
	  else {alpha=alpha1+(ALPHA-alpha1)/(2*(1+(F1-F)/((ALPHA-alpha1)*g1)));
      alpha2=ALPHA;
	  ALPHA=alpha;}
step3: //for(i=0;i<N;i++)
		  //X[i]=X[i]+ALPHA*D[i]; 
	   CALG(N,X,G);
	   count2++;
	   g=DOIT(G,D);
	   if(g>=t*g1)goto end1;
	   else {alpha=ALPHA+((ALPHA-alpha1)*g)/(g1-g);
		   alpha=ALPHA;
	   F1=F;
	   g1=g;
	   ALPHA=alpha;
	   }
end1:;

}

void BFGS_design(double H[N][N],double s[N],double y[N]){
	//DFP_design
	double t1,t2;
	double temp1[N],temp2[N];
	double temp3[N][N],temp4[N][N],temp5[N][N],temp_w[N][N];
	double temp_w1[N],temp_w2[N],temp_w3[N];
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp3[i][j]=H[i][j];
		DOIT1(H,y,temp1);
		DOIT2(temp1,y,temp4);
		DOIT3(temp4,H,temp5);
		DOIT4(y,temp3,temp2);
		t1=DOIT(temp2,y);
		DOIT2(s,s,temp3);
		t2=DOIT(y,s);
		for(i=0;i<N;i++)
			temp_w1[i]=s[i]/t2;
		for(i=0;i<N;i++)
			temp_w2[i]=temp1[i]/t1;
		for(i=0;i<N;i++)
			temp_w3[i]=sqrt(t1)*(temp_w1[i]-temp_w2[i]);
		DOIT2(temp_w3,temp_w3,temp_w);
	for(i=0;i<N;i++)
		for(j=0;j<N;j++){
			temp5[i][j]/=t1;
			temp3[i][j]/=t2;
		}
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			H[i][j]=H[i][j]-temp5[i][j]+temp3[i][j]+temp_w[i][j];
}
void main()
{
	int i,j;
	int k=0,k1=0;
	int count1=0,count2=0,count3=0;
	double f;
	double alpha;
	double p[N];//search 
	double g[N];
	double s[N],y[N],tempg[N],tempx[N];
//step1:
	double x[N];
	double H[N][N];
	clock_t start,finish;
	for(i=0;i<N;i++)
		if(i%2==0)x[i]=-1.2;
			else x[i]=1;
	CALG(N,x,g);
	CALF(N,x,f);
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			if(i==j)H[i][j]=1;
			else H[i][j]=0;
   start=clock();   //时间起点
step2:
	DOIT1(H,g,p);
	for(i=0;i<N;i++)
			p[i]=-p[i];
//step3:
	for(i=0;i<N;i++)
		tempx[i]=x[i];
	INEXLS(N,tempx,p,alpha,count1,count2);
	k1+=count1;
	count3+=count2;
//	for(i=0;i<N;i++)
//	cout<<x[i]<<'\t';
	cout<<"alpha:"<<alpha<<'\n';
//step4:
     for(i=0;i<N;i++)
		x[i]=x[i]+alpha*p[i];
//step5:
	if(if_g(x))goto end;
	CALG(N,x,tempg);
	count3++;
	for(i=0;i<N;i++){
		s[i]=alpha*p[i];
		y[i]=tempg[i]-g[i];
	}
	for(i=0;i<N;i++)
		g[i]=tempg[i];
//step6:
	BFGS_design(H,s,y);
	k++;
	goto step2;
end:
  finish=clock();  //时间终点
  CALF(N,x,f);
  cout<<"the opt solution:\n";
  for(i=0;i<N;i++)
	cout<<x[i]<<'\t';
  cout<<endl;
  cout<<"The opt value:"<<f<<endl;
  cout<<"number of iter:"<<k<<endl;
  k1=k1+k;
  cout<<"number of elavation:"<<k1<<endl;
  cout<<"NUMBER OF GRADIENT EVALUATION:"<<count3<<endl;
  cout<<"The duration is :"<<double(finish-start)/CLOCKS_PER_SEC<<'m';
}