newton_dfp.c

拟牛顿法求函数极小值

double newton_DFP(double (*pf)(double *x),int n,double *min_point)
{
        int i,j;
        int k=0;
        double e=1E-5;
        double g_norm;
        double *g0;            
        double *g1;
        double *dg;
        double *p;              
        double t;                                               
        double *x0;
        double *x1;
        double *dx;
        double **H;
        double **tempH;
        double *gH;
        double *Hg;
        double num1;
        double num2;
        g0=dvector(0,n-1);
        g1=dvector(0,n-1);
        dg=dvector(0,n-1);
        p=dvector(0,n-1);
        x0=dvector(0,n-1);
        x1=dvector(0,n-1);
       dx=dvector(0,n-1);
        H=dmatrix(0,n-1,0,n-1);
        tempH=dmatrix(0,n-1,0,n-1);
       gH=dvector(0,n-1);
       Hg=dvector(0,n-1);
        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;
                }
        for(i=0;i<n;i++)
                x0[i]=min_point[i];  
        for(i=0;i<n;i++)         printf("initial x0=%f\n",x0[i]);
        grad(pf,n,x0,g0);
     printf("pf=%f\n",*pf);
        g_norm=0.0;
        for(i=0;i<n;i++)        g_norm=g_norm+g0[i]*g0[i];        
        g_norm=sqrt(g_norm);
        printf("g_norm=%f\n",g_norm);

        if (g_norm<e) 
        {
                for(i=0;i<n;i++)        min_point[i]=x0[i];
                free_dvector(g0,0,n-1);   
                free_dvector(g1,0,n-1);                
                free_dvector(dg,0,n-1);   
                free_dvector(p,0,n-1);               
                free_dvector(x0,0,n-1);                
                free_dvector(x1,0,n-1);     
                free_dvector(dx,0,n-1);          
                free_dmatrix(H,0,n-1,0,n-1);
                free_dmatrix(tempH,0,n-1,0,n-1);
                free_dvector(gH,0,n-1);                 
                free_dvector(Hg,0,n-1);     
                return pf(min_point);
        }


        for(i=0;i<n;i++)        p[i]=-g0[i];        
  printf("DFP00OK\n") ;
     do
   {

                 t=line_search0618(pf,n,x0,p);   
printf("t=%f\n",t) ;                           
                for(i=0;i<n;i++)        x1[i]=x0[i]+t*p[i];
                grad(pf,n,x1,g1);
                g_norm=0.0;
                for(i=0;i<n;i++)        g_norm=g_norm+g1[i]*g1[i];
                g_norm=sqrt(g_norm);
printf("g_norm=%f\n",g_norm);
                if (g_norm<e) 
                {
                        for(i=0;i<n;i++)        min_point[i]=x1[i];
                free_dvector(g0,0,n-1);   
                free_dvector(g1,0,n-1);                
                free_dvector(dg,0,n-1);   
                free_dvector(p,0,n-1);               
                free_dvector(x0,0,n-1);                
                free_dvector(x1,0,n-1);     
                free_dvector(dx,0,n-1);          
                free_dmatrix(H,0,n-1,0,n-1);
                free_dmatrix(tempH,0,n-1,0,n-1);
                free_dvector(gH,0,n-1);                 
                free_dvector(Hg,0,n-1);  
printf("hehe in do g_norm<e\n");
                return pf(min_point);
                }
                for(i=0;i<n;i++)
                {
                        dx[i]=x1[i]-x0[i];
                        dg[i]=g1[i]-g0[i];
                }

                for(i=0;i<n;i++)
                {
                        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++)
                {
                        num1=num1+dx[i]*dg[i];
                        num2=num2+gH[i]*dg[i];
                }
                 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;
                        }
                }
                for(i=0;i<n;i++)
                {
                        for(j=0;j<n;j++)
                        {
                                H[i][j]=tempH[i][j];
                        }
                }
                for(i=0;i<n;i++)        p[i]=0.0;
                for(i=0;i<n;i++)
                {
                        for(j=0;j<n;j++)
                        {
                                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;
       if (k>6)
        {printf("k>6\n");
        break;}
     }while(g_norm>e);
printf("k=%d\n",k);
        for(i=0;i<n;i++)        min_point[i]=x1[i];
                free_dvector(g0,0,n-1);   
                free_dvector(g1,0,n-1);                
                free_dvector(dg,0,n-1);   
                free_dvector(p,0,n-1);               
                free_dvector(x0,0,n-1);                
                free_dvector(x1,0,n-1);     
                free_dvector(dx,0,n-1);          
                free_dmatrix(H,0,n-1,0,n-1);
                free_dmatrix(tempH,0,n-1,0,n-1);
                free_dvector(gH,0,n-1);                 
                free_dvector(Hg,0,n-1);  
        return pf(min_point);
}