算法及编程语言 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

风花雪月 发表于 2006-11-12 06:58

拟牛顿法（变尺度法）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 *direction); //0.618法线搜索
double line_search(double (*pf)(double *x), int n,double *start,double *direction); //解析法线搜索
double DFP(double (*pf)(double *x), 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;
for(i=1;i<=n;i++)
{
temp=point;
}
for(i=1;i<=n;i++)
{
temp+=0.5*h;
grad=4*pf(temp)/(3*h);
temp-=h;
grad-=4*pf(temp)/(3*h);
temp+=(3*h/2);
grad-=(pf(temp)/(6*h));
temp-=(2*h);
grad+=(pf(temp)/(6*h));
temp=point;
}
delete[] temp;
}

//一维搜索模块
//参数：指向目标函数的指针，变量个数，出发点，搜索方向
//返回：最优步长
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;

grad=new double;
temp_point=new double;
comput_grad(pf,n,start,grad);
diver_a=0;
for(i=1;i<=n;i++)
diver_a=diver_a+grad*direction;
do
{
b=a+step;
for(i=1;i<=n;i++)
temp_point=start+b*direction;
comput_grad(pf,n,temp_point,grad);
diver_b=0;
for(i=1;i<=n;i++)
diver_b=diver_b+grad*direction;
if( fabs(diver_b)<1E-10 )
{
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=start+a*direction;
value_a=(*pf)(temp_point);
for(i=1;i<=n;i++)
temp_point=start+b*direction;
value_b=(*pf)(temp_point);

do
{
s=3*(value_b-value_a)/(b-a);
z=s-diver_a-diver_b;
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++)
temp_point=start+t*direction;
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*direction;
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) );

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;

double *g0=new double; //梯度
double *g1=new double;
double *dg=new double;

double *p=new double; //搜索方向 =-g
double t; //一维搜索步长

double *x0=new double;
double *x1=new double;
double *dx=new double;

double **H=new double*;
for (i=0; i<n; i++) H = new double;

double **tempH=new double*;
for (i=0; i<n; i++) tempH = new double;

double *gH=new double;
double *Hg=new double;
double num1;
double num2;

for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(i==j) H=1.0; // H0=I
else H=0.0;
tempH=0.0;
}

for(i=0;i<n;i++)
x0=min_point;

comput_grad(pf,n,x0,g0);

g_norm=0.0;
for(i=0;i<n;i++) g_norm=g_norm+g0*g0;
g_norm=sqrt(g_norm);
if (g_norm<e)
{
for(i=0;i<n;i++) min_point=x0;

delete[] g0;
delete[] g1;
delete[] dg;
delete[] p;
delete[] x0;
delete[] x1;
delete[] dx;
for (i=0; i<n; i++) delete[] H;
delete []H;
for (i=0; i<n; i++) delete[] tempH;
delete []tempH;
delete[] gH;
delete[] Hg;

return pf(min_point);
}

for(i=0;i<n;i++) p=-g0;

do
{
t=line_search(pf,n,x0,p);
for(i=0;i<n;i++) x1=x0+t*p;
comput_grad(pf,n,x1,g1);
g_norm=0.0;
for(i=0;i<n;i++) g_norm=g_norm+g1*g1;
g_norm=sqrt(g_norm);
//cout<<k<<" "<<x0<<"    "<<x0<<"    "<<g_norm<<"\n";
if (g_norm<e)
{
for(i=0;i<n;i++) min_point=x1;

delete[] g0;
delete[] g1;
delete[] dg;
delete[] p;
delete[] x0;
delete[] x1;
delete[] dx;
for (i=0; i<n; i++) delete[] H;
delete []H;
for (i=0; i<n; i++) delete[] tempH;
delete []tempH;
delete[] gH;
delete[] Hg;

return pf(min_point);
}
for(i=0;i<n;i++)
{
dx=x1-x0;
dg=g1-g0;
}

//////////////////求Hk+1的矩阵运算

//g*H,H*g
for(i=0;i<n;i++)
{
gH=0.0;
Hg=0.0;
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
gH=gH+dg*H;
//Hg=Hg+H*dg;
Hg=gH;
}
}

//num1,num2
num1=0.0;
num2=0.0;
for(i=0;i<n;i++)
{
num1=num1+dx*dg;
num2=num2+gH*dg;
}

//tempH
for(i=0;i<n;i++)
for(j=0;j<n;j++)
tempH=0.0;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
tempH=tempH+H;
tempH=tempH+dx*dx/num1;
tempH=tempH-Hg*gH/num2;
}
}

for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
H=tempH;
}
}
/////////////////////////////

//P
for(i=0;i<n;i++) p=0.0;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
p=p-H*g1;
}
}

for(i=0;i<n;i++)
{
g0=g1;
x0=x1;
}
k=k+1;
}while(g_norm>e);

for(i=0;i<n;i++) min_point=x1;

delete[] g0;
delete[] g1;
delete[] dg;
delete[] p;
delete[] x0;
delete[] x1;
delete[] dx;
for (i=0; i<n; i++) delete[] H;
delete []H;
for (i=0; i<n; i++) delete[] tempH;
delete []tempH;
delete[] gH;
delete[] Hg;

return pf(min_point);
}

/////////////
double fun(double *x)
{
return 100*( x-x*x )*( x-x*x ) + (1-x)*(1-x);
}

void main()
{
int n=2;
double min_point={-5,10};
double min_value=DFP(fun,n,min_point);
cout<<min_point<<" and "<<min_point<<" "<<min_value;

}

//0.618法线搜索
//
//参数：指向目标函数的指针，变量个数，出发点，搜索方向
//返回：最优步长
//
double line_search1(
   double (*pf)(double *x),
   int n,
   double *start,
   double *direction)
{
int i;
int k;

double l,a,b,c,u,lamda,t,e;

double *xa=new double;
double *xb=new double;
double *xc=new double;
double *xl=new double;
double *xu=new double;

//确定初始搜索区间
l=0.001;
a=0;

k=0;
do
{
k++;
c=a+l;
for(i=0;i<n;i++)
{
xc=start+c*direction;
xa=start+a*direction;
}
l=l/3;
}while( pf(xc) >= pf(xa) ); // ???

k=0;
do
{
k++;
l=2*l;
b=c+l;
for(i=0;i<n;i++)
{
xc=start+c*direction;
xb=start+b*direction;
}
a=c;
c=b;
}while( pf(xb) <= pf(xc) );

a=0;
b=0.1;

//寻优
t=0.618;
e=0.000001;

lamda=b-t*(b-a);
u=a+t*(b-a);

for(i=0;i<n;i++)
{
xl=start+lamda*direction;
xu=start+u*direction;
}

k=0;
do
{
k++;
if( pf(xl)<pf(xu) )
{
b=u;
u=lamda;
lamda=b-t*(b-a);
}
else
{
a=lamda;
lamda=u;
u=t*(b-a);
}
}while( b-a>=e );

lamda=(a+b)/2;

delete[] xa;
delete[] xb;
delete[] xc;
delete[] xl;
delete[] xu;

return lamda ;

}

haixu2000 发表于 2007-5-26 21:41

回复 #1 风花雪月的帖子

好贴，帮了我大忙，程序很好懂阿:lol :lol

hitcjh 发表于 2007-11-26 23:51

非常感谢，可以参考的资料

yeting0504@hotm 发表于 2008-7-8 10:04

非常感谢，看这个就是舒服啊！谢谢

20031374 发表于 2008-12-14 17:50

非常感谢，看这个就是舒服啊:@P

siping4840 发表于 2009-8-26 11:03

不错，受用了！

不错，受用了！不错，受用了！

modest_xiu 发表于 2012-5-7 17:36

参考~~谢谢

程瑞岩 发表于 2012-6-17 10:07

谢谢学习了。 {:{27}:}

nari 发表于 2012-7-18 22:57

拟牛顿法主要是什么思想呀，可否有个简介呀

www.zdscj.com| 发表于 2012-7-21 15:55

新人来学习一下先

页: [1]

声振论坛's Archiver

拟牛顿法（变尺度法）DFP算法的c/c++源码

回复 #1 风花雪月 的帖子

不错 ，受用了！

回复 #1 风花雪月的帖子

不错，受用了！