拟牛顿法.docx
《拟牛顿法.docx》由会员分享，可在线阅读，更多相关《拟牛顿法.docx（17页珍藏版）》请在冰点文库上搜索。
拟牛顿法
拟牛顿法（变尺度法）DFP算法的c/c++源码
#include"iostream.h"
#include"math.h"
voidcomput_grad（double（*pf）（double*x）,intn,double*point,double*grad）;//计算梯度
doubleline_search1（double（*pf）（double*x）,intn,double*start,double*direction）;//0.618法线搜索
doubleline_search（double（*pf）（double*x）,intn,double*start,double*direction）;//解析法线搜索
doubleDFP（double（*pf）（double*x）,intn,double*min_point）;//无约束变尺度法
//梯度计算模块
//参数：
指向目标函数的指针，变量个数，求梯度的点，结果
voidcomput_grad（double（*pf）（double*x）,
intn,
double*point,
double*grad）
{
doubleh=1E-3;
inti;
double*temp;
temp=newdouble[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;
}
//一维搜索模块
//参数：
指向目标函数的指针，变量个数，出发点，搜索方向
//返回：
最优步长
doubleline_search（
double（*pf）（double*x）,
intn,
double*start,
double*direction）
{
inti;
doublestep=0.001;
doublea=0,value_a,diver_a;
doubleb,value_b,diver_b;
doublet,value_t,diver_t;
doubles,z,w;
double*grad,*temp_point;
grad=newdouble[n];
temp_point=newdouble[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
{
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）
{
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];
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;
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[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];
if（diver_t>1E-6）
{
b=t;
value_b=value_t;
diver_b=diver_t;
}
elseif（diver_t<-1E-6）
{
a=t;
value_a=value_t;
diver_a=diver_t;
}
elsebreak;
}while（（fabs（diver_t）>=1E-6）&&（fabs（b-a）>1E-6））;
delete[]grad;
delete[]temp_point;
return（t）;
}
//无约束变尺度法DFP函数声明
//
//参数：
pf指向目标函数的指针，n变量个数，min_point接受初始点、存放结果
//返回：
极小点处函数值
//
doubleDFP（
double（*pf）（double*x）,
intn,
double*min_point
）
{
inti,j;
intk=0;
doublee=1E-5;
doubleg_norm;
double*g0=newdouble[n];//梯度
double*g1=newdouble[n];
double*dg=newdouble[n];
double*p=newdouble[n];//搜索方向=-g
doublet;//一维搜索步长
double*x0=newdouble[n];
double*x1=newdouble[n];
double*dx=newdouble[n];
double**H=newdouble*[n];
for（i=0;idouble**tempH=newdouble*[n];
for（i=0;idouble*gH=newdouble[n];
double*Hg=newdouble[n];
doublenum1;
doublenum2;
for（i=0;ifor（j=0;j{
if（i==j）H[i][j]=1.0;//H0=I
elseH[i][j]=0.0;
tempH[i][j]=0.0;
}
for（i=0;ix0[i]=min_point[i];
comput_grad（pf,n,x0,g0）;
g_norm=0.0;
for（i=0;ig_norm=sqrt（g_norm）;
if（g_norm{
for（i=0;idelete[]g0;
delete[]g1;
delete[]dg;
delete[]p;
delete[]x0;
delete[]x1;
delete[]dx;
for（i=0;idelete[]H;
for（i=0;idelete[]tempH;
delete[]gH;
delete[]Hg;
returnpf（min_point）;
}
for（i=0;ido
{
t=line_search（pf,n,x0,p）;
for（i=0;icomput_grad（pf,n,x1,g1）;
g_norm=0.0;
for（i=0;ig_norm=sqrt（g_norm）;
//cout<if（g_norm{
for（i=0;idelete[]g0;
delete[]g1;
delete[]dg;
delete[]p;
delete[]x0;
delete[]x1;
delete[]dx;
for（i=0;idelete[]H;
for（i=0;idelete[]tempH;
delete[]gH;
delete[]Hg;
returnpf（min_point）;
}
for（i=0;i{
dx[i]=x1[i]-x0[i];
dg[i]=g1[i]-g0[i];
}
//////////////////求Hk+1的矩阵运算
//g*H,H*g
for（i=0;i{
gH[i]=0.0;
Hg[i]=0.0;
}
for（i=0;i{
for（j=0;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{
num1=num1+dx[i]*dg[i];
num2=num2+gH[i]*dg[i];
}
//tempH[i][j]
for（i=0;ifor（j=0;jtempH[i][j]=0.0;
for（i=0;i{
for（j=0;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{
for（j=0;j{
H[i][j]=tempH[i][j];
}
}
/////////////////////////////
//P
for（i=0;ifor（i=0;i{
for（j=0;j{
p[i]=p[i]-H[i][j]*g1[j];
}
}
for（i=0;i{
g0[i]=g1[i];
x0[i]=x1[i];
}
k=k+1;
}while（g_norm>e）;
for（i=0;idelete[]g0;
delete[]g1;
delete[]dg;
delete[]p;
delete[]x0;
delete[]x1;
delete[]dx;
for（i=0;idelete[]H;
for（i=0;idelete[]tempH;
delete[]gH;
delete[]Hg;
returnpf（min_point）;
}
/////////////
doublefun（double*x）
{
return100*（x[1]-x[0]*x[0]）*（x[1]-x[0]*x[0]）+（1-x[0]）*（1-x[0]）;
}
voidmain（）
{
intn=2;
doublemin_point[2]={-5,10};
doublemin_value=DFP（fun,n,min_point）;
cout<}
//0.618法线搜索
//
//参数：
指向目标函数的指针，变量个数，出发点，搜索方向
//返回：
最优步长
//
doubleline_search1（
double（*pf）（double*x）,
intn,
double*start,
double*direction）
{
inti;
intk;
doublel,a,b,c,u,lamda,t,e;
double*xa=newdouble[n];
double*xb=newdouble[n];
double*xc=newdouble[n];
double*xl=newdouble[n];
double*xu=newdouble[n];
//确定初始搜索区间
l=0.001;
a=0;
k=0;
do
{
k++;
c=a+l;
for（i=0;i{
xc[i]=start[i]+c*direction[i];
xa[i]=start[i]+a*direction[i];
}
l=l/3;
}while（pf（xc）>=pf（xa））;//?
?
?
k=0;
do
{
k++;
l=2*l;
b=c+l;
for（i=0;i{
xc[i]=start[i]+c*direction[i];
xb[i]=start[i]+b*direction[i];
}
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{
xl[i]=start[i]+lamda*direction[i];
xu[i]=start[i]+u*direction[i];
}
k=0;
do
{
k++;
if（pf（xl）{
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;
returnlamda;
}
升级会员 