变尺度法BFGS算法的C--源码

1、/viewthread.php?tid=73这个是我那时候交的作业（很粗了。）。嗯。应该标注很清楚了。不过如果你用matlab编写的话，很多大段的代码，基本上一句话就代替了。& A5 n) g, 6 U; E* M3 R3 ! C( W/ c/变尺度法BFGS算法的C+源码4 , F5 # U9 A. n8 I I7 m& G7 v#include iostream.h% Y1 H* O5 t; u6 f4 Y#include math.h& e6 J4 e! a: A9 n% M#define A85 V* k; O- 0 f! s#defin

2、e B10) X3 J6 u) L2 s2 t% k#define C75 E- j+ L4 6 x! R#define D0, ?/ Bb, ku3 O5 x D#define START 5,4,5,5 /初始迭代点0 i7 C9 K7 G% f7 Qconst int n=4; /变量个数6 9 q1 0 c3 h& H: gconst int nc=3; /不等式约束个数(也可根据题目需要增加等式约束个数修改约束构造函数)3 m6 O6 O0 J$ P, I: DL3 d9 O$ C# p p# k 9 yconst double epsilon1=1E-6;/一维搜索精度4 0 L(

3、 ?& N1 V0 y8 const double epsilon2=1E-4;/变尺度法收敛精度5 v* NZP; rconst double epsilon=1E-4;/惩罚函数法收敛精度, j/ t3 , a% p7 q: Adouble r=1; /惩罚因子# t X( A% s3 O( o: B, Udouble br=1; /变尺度法强行退出循环因子(可去除)3 d/ V) + C7 t/ d+ F4 N/ 9 cconst double zr=10; /递增系数% & |6 W! W F% B+ . udouble det=0.0014; /外惩罚函数法约束函数余量3 Q: U&

4、 I T+ U2 y B. O# Q1 9 _/(此det值仅适合当前A B C D，修改请将其恢复为0.0001或其他建议值)3 Y- J4 E4 6 h$ 9 v9 z; G/计算梯度# ! i, 3 f! o% mvoid comput_grad(double (*pf)(double *x),int n,double *point,double *grad); ( q; * I E4 A% . t& H d, f - q0 b4 $ t3 L) t4 X/无约束BFGS变尺度法（目标函数，维数，最小点）返回0值; W q1 wU* P% od* mdouble BFGS(double

5、(*pf)(double *x), int n, double *min_point);! l / s |) C$ a- h9 1 H ( e( K+ a; f Z0 s/ a# T /构造混合外点惩罚函数 返回惩罚函数值（构造点）9 v/ v* K( x# o4 V ) double funP(double *p); w5 L2 u C7 B S# e! M/返回目标函数值,并返回约束函数值存入ys （计算点，约束函数值存放处） / x7 w# X5 U: i+ F3 |double fun(double *x , double *ys );) . b# l9 V; a% c# x/ ; q

6、5 B# w/一维搜索函数 返回最优步长（目标函数，维数，初始点，方向）3 Z7 # N/ 5 gdouble line_search( double (*pf)(double *x),int n,double *start,double *direction); k. B7 Lh v- ; Y( W/配给与一维搜索函数 返回下一点的的惩罚函数值（计算点，计算方向，步长）: K- ?- v( L1 o- ddouble funPt(double *x,double *direction,double t);2 n) V: F9 G% W# ?/ F- Z9 pBvoid main()7 - J

7、. Y0 * C* d8 * C* z: n# h. y4 9 l v& C$ b+ U7 w: int i,j;% 4 N) t* y h F i double c=0;/表示外惩罚函数循环次数1 Q D: I6 x7 4 h) l1 + V: y, T! K2 double f_min,f_temp,f_punish;0 tT; E+ M2 l5 n9 y k7 mdouble *ys=new doublenc;O/ X% S3 d6 v! W. double x_x=0,f_ys=0;) |+ X5 X+ u4 e5 c; oR0 C double x_minn=START;9 u1 ,

8、 Y0 H6 I+ E7 U double x_tempn=START;2 e* n$ s L# Sf_temp=fun(x_temp,ys);& I8 t3 q; pfor(;)3 g3 F6 I# & r5 n / Y$ P7 m- ?. m% a2 5 y( Z- rc+;- z. 4 3 / M- hBFGS(funP,n,x_min);- P7 w9 # p9 i r3 Z# N/ n+ S5 f_min=fun(x_min,ys);! Z4 % k/ k W) a* Q$ efor(j=0;jn;j+)l! Z+ L! z% S x_x=x_x+pow(x_minj-x_tempj

9、),2);S9 d$ q$ r* W$ # i$ V W, pfor(i=0;inc;i+)L# q 8 d$ a f_ys+=fabs(ys);( ( - O) ?/ - yx_x=sqrt(x_x);: r: P; n2 z0 P/ M: + G2 s /参考值显示 ; r8 N2 h/ 8 U. R/ ncoutn外惩罚函数循环迭代第 c 轮endl. m$ a7 & F% Y3 目标函数: f_min=f_minendl( z: m3 C) E4 RP5 E* 自变量: x1=x_min0 x2=x_min1 x3=x_min2 x4=x_min3endl3 B1 P* e* n7 v

10、7 + D 约束函数: g1=ys0+det g2=ys1+det g3=ys2+detendl5 l0 z2 i6 ( f% K& l+ k% Q 判断收敛因子: x_x=x_xf_f=fabs(f_min-f_temp)det_f=fabs(f_min-f_temp)/f_min)f_ys=f_ysendl; 3 * d5 Z- P; ! x2 S0 d# D4 n/-|9 : 7 k# m+ h& m3 , E/收敛判断5 y0 V; x) c* iif(fabs(x_x)=epsilon)|(f_ys=epsilon); q1 |% a1 v T cout因为迭代xk与xk+1距离x_

11、x 或者 约束函数值f_ys 达到收敛精度 epsilon 而认为收敛endl;! o4 e# N5 w3 O break;$ c/ k0 F- v* zif(fabs(f_min)=1)&(fabs(f_min-f_temp)=epsilon)+ J% B! W1 yd: w4 _$ z. cout因为迭代xk与xk+1的函数值差距f_f达到收敛精度epsilon 而认为收敛 =1)&(fabs(f_min-f_temp)/f_min)=epsilon)2 z G% b7 d4 E: k0 U6 & O coutn 因为迭代xk与xk+1的函数值差距 det_f 达到收敛精度 epsilon

12、 而认为收敛 endl;: . T+ P( c7 n: V- O/ b7 Z4 r break;* |& g# w! k/ O, F9 q/-|+ P! I8 b! _) r=r*10;7 G- Q( l8 u, L1 G2 3 $ b/ t. Q4 e; pfor(i=0;if_temp)- I5 m! e. H& s% q) yf_min=f_temp;9 A; ; t3 l ofor(i=0;i=1E10)3 F$ t8 I) 9 t) u7 r0 _coutn!由于不可知的原因，导致结果发散!endl; l5 H/ F2 p& P: u3 zF) b2 B 请在源程序修改合适的精度值或

13、者其他nendl;& f. D: O1 7 c8 M, M$ 4 d% A5 ; * b coutn -外惩罚函数+BFGS变尺度法优化算法C+程序最终结果-endl0 g G t+ g2 z0 u( w0 u *也可根据实际需要在以上参考值中选择其他结果*endl;3 i6 Z1 j1 y, Z$ hi cout优化最小值:endl+ s1 & s8 Y( t/ F- |! . B 目标函数f=f_min 构造的惩罚函数值=f_punishendl! 2 n) j$ 2 F4 5 优化程序得到的最优点为:endl|4 B1 F3 c* A3 y1 |J3 j x1=x_min0 x2=x_m

14、in1 x3=x_min2 x4=x_min3endl# a3 a5 m1 B8 m8 P* 8 L约束条件情况:endl3 Z$ C# _; D1 e X0 r) J2 U* F3 约束1=ys0+det 约束2=ys1+det 约束3=ys2+detendl; 9 8 _$ t$ ?6 q/ n; ub2 m# C# l( ( W8 Z/ r/ I9 g t+ X/梯度计算模块2 cw5 Q$ o2 k E( k/参数：指向目标函数的指针，变量个数，求梯度的点，结果存放处. X+ a! % h( ; I# avoid comput_grad(double (*pf)(double *x),

15、1 a- A* g& & E7 S int n,/ R8 l6 jS) R, D4 7 v. Q8 w double *point,7 O: 2 R9 P; S) s7 0 |8 L$ 5 Y double *grad)& M9 a5 S( Y) l8 W7 L2 T) |8 _1 1 M 0 Q. o double h=1E-5;7 l5 n# T9 l2 | int i;/ M5 G! g2 Y/ a2 | B double *temp;) p/ : Y3 H+ ) G1 k temp = new doublen; X3 N5 _% p/ e j% m for(i=1;i=n;i+)/ |

16、- k9 & k9 E! V* L & H8 e7 g5 R+ K4 K8 L0 g8 6 u tempi-1=pointi-1; Q6 t! - I3 S& 9 M ) O7 d; Z( X5 d, D$ G2 J z$ _ for(i=1;i=n;i+)/ r h4 p0 h a3 B1 X9 e) E- Z$ y+ b# p) D tempi-1+=0.5*h;7 U$ J9 E( g 4 e- U1 c* D gradi-1=4*pf(temp)/(3*h); - f S X: 4 | tempi-1-=h;/ z B7 9 L) A. S+ s0 y* k- O2 gradi-1-=

17、4*pf(temp)/(3*h);) O) # z) D# Y l- w tempi-1+=(3*h/2);. Z2 S$ r% b$ hVB gradi-1-=(pf(temp)/(6*h);6 P* t% Q/ Z. tempi-1-=(2*h);4 v: $ | s* r2 s! P gradi-1+=(pf(temp)/(6*h);# q5 A, L8 Y3 O4 R tempi-1=pointi-1;4 ?4 Q7 X, b5 G& O) j | 7 B% t; Q1 l: m& g2 X delete temp;2 k# V7 + w: i* h: d& A- b- d& Q) e

18、+ L+ R8 T3 e3 t) JP+ j1 _ u0 u( G3 i/无约束变尺度法BFGS函数声明) F/ e Y: _3 u0 E0 u7 Y9 t3 i/参数：pf指向目标函数的指针，n变量个数，min_point接受初始点、存放结果) _% P5 0 A* W% 2 5 4 v d/返回0值$ e1 # a2 j Hdouble BFGS(. d9 Q$ / u! |6 D double (*pf)(double *x),3 k q/ p1 , n0 F: t+ u int n,3 X7 O9 j: 9 double *min_point0 wx4 H# P9 p& H0 X )*

19、 T, s# c! 5 Q& h/ G# O5 J# Q- u* I. L2 f. W int i,j;( V+ S/ z( j7 q, s int k=0;) q: jT* Y r6 _! e; R Ddouble e=epsilon2;) / y w( 6 G( S- L4 X/ Ndouble g_norm; /梯度的模3 I2 Z# U/ S k1 D . a a* W . i w0 j7 K7 X2 t+ A3 double *g0=new doublen; /梯度3 s1 + Q. u5 W% q double *g1=new doublen; O: U& v. a) i+ s/

20、V double *dg=new doublen;& j N9 j/ h4 N( Q& O double *p=new doublen; /搜索方向 =-g, x2 V7 z$ t: v0 hL8 D6 y! E double t; /一维搜索步长8 F% Y) 1 # K7 Dc _ Sdouble *x0=new doublen;4 Y; BU# qc s+ o% A. Edouble *x1=new doublen;! l) z* R) Aj- u% MH double *dx=new doublen;+ F 2 _# c) X9 Y2 J double *H=new double*n;

21、4 J& N% a$ , S6 r for (i=0; in; i+) H = new doublen;4 y. + J# E( B1 i7 h5 q double *gH=new doublen; /计算H变化量的中间变量-9 G1 d3 u8 r% N4 Y. & z double *Hg=new doublen;1 p; Q* ?+ E$ ! Y: 9 m L double num1;C9 1 L) c, F0 _6 v3 s2 double num2; /- |6 pm1 q/ n# R1 T$ j+ A. j4 ; ?: s9 e /初始化-o# u! Ge0 K* $ ?* s7

22、q for(i=0;in;i+)& A$ ) V! r* Z4 y/ F for(j=0;jn;j+); Q. D; r; + j 9 M6 U, u, j6 z- L 1 H if(i=j) Hj=1.0; / H0=I2 J5 x/ b |: H0 A$ I9 5 J% G else Hj=0.0;% F- _5 h9 h2 K7 D u/ U. / t; G* O# W9 e. M g1 d2 I/ P. n8 o8 c% U0 v for(i=0;in;i+)- 3 F: r. I. X9 X, S4 I1 |# D x0=min_point; : h6 c( k0 ?3 j) Q&

23、q comput_grad(pf,n,x0,g0);! _2 l3 ?) z/ ) B for(i=0;in;i+) p=-g0; O- z* g8 S8 E4 d! C) /-初始化结束 N. X% U7 I- r- x5 w do+ S( kr- & & ?I3 J1 Q 9 V0 D- L0 be/ Bt=line_search(pf,n,x0,p); : i: O/ j: q5 r! o4 Xfor(i=0;in;i+) , X. h3 & p* p; q* f x1=x0+t*p; 6 . r% k1 n3 z8 f* o comput_grad(pf,n,x1,g1); ( f;

24、u: N% |8 w# + /调试程序时参考输出8 Y4 T7 I) J4 z5 9 b/ coutk=kendlx00x01x02 x03 endlt=t; e=g_normendl;( t/ N; W1 f8 y2 B: x# c6 N1 p$ u A 5 I7 ?% u1 w0 K# |2 a for(i=0;in;i+)+ Y8 g5 H: o3 q9 : + B- , z$ t. u/ X0 a: U C% X9 Y. A dx=x1-x0;) W% |, / G5 E $ i dg=g1-g0;! N9 i5 n. N4 l$ p8 Vz + o h8 y; m u/ l8 x0

25、C: G /求H(k+1)的矩阵运算/% 6 K4 ( k: W- Z% u5 2 K /g*H,H*g ! b; U! # l for(i=0;in;i+) # H! F8 l- R- a0 X5 w; S$ x3 |7 G gH=0.0; /u* b7 K% g& D$ V, l Hg=0.0; /v- Z# Y7 g8 . q k0 l4 J6 E6 W) u- c! Q0 r4 g, O& w, P; E( C% i for(i=0;in;i+)8 P# n# a- & V: O5 Z 4 I7 + j3 G8 X) I. a! W3 E2 c0 s for(j=0;jn;j+)( 9

26、 ( F, v/ ) ) n, * p! C7 n 8 f, o6 ; % m a gH=gH+dgj*Hj;# $ l Y1 P9 N r7 R% Hg=Hg+Hj*dgj;8 y, - X a+ Q! 3 n! d! f! Y4 e& M y- W, D7 a ( Z/ w$ M/ j+ n0 K /num1,num21 G* A6 b I, b/ P E( s num1=0.0;. j. - $ q4 r- P+ & L num2=0.0;1 C/ a2 h7 v. z/ P for(i=0;in;i+). T9 k! E+ M/ g! p 4 R7 X9 f9 ?8 U! num1=n

27、um1+dx*dg;7 i8 ; Z- P+ p num2=num2+gH*dg;0 M3 |N% Q* k; N a- ug4 F # H/*if(k!=0&k%n=0) /由于在计算DFP变尺度矩阵的公式中，其分母会含有, t. M6 i3 t( S0 z3 V for(i=0;in;i+) /近似矩阵H，计算中容易引起数值不稳定，甚至可能8 W: L2 A1 F u: A for(j=0;jn;j+) /得到奇异矩阵，为了提高实际运算稳定性，通常将进行% M1 b8 I- p8 FE, / v if(i=j) Hj=1.0; /n次迭代作为一个循环，将变尺度矩阵重置单位阵( N& dYN

28、 $ d- J2 T m( 7 q elseHj=0.0; /左边代码作为DFP参考。* z# Y2 B3 r0 O. b; l | ! N J! . l$ & S; C: X- B8 Pelse*/ if(fabs(num1)!=0) /BFGS变尺度法修正公式可以改善DFP的计算稳定性问题- / Dk+ O7 z O- o% ( B ?, U( D8 O3 n6 f7 z, v0 l2 Z num2=1+num2/num1;* # A, H6 G) E3 q* Q# Hs for(i=0;in;i+) h/ | s, E3 BH( J for(j=0;jn;j+): |& T3 5 n9

29、i. P5 L. g Hj=Hj+(num2*dx*dxj-Hg*dxj-dx*gHj)/num1;% W9 p5 f j& N+ x, d1 i X) q+ t9 O6 jK* W$ J T3 q /Hk+1 计算完毕/. t5 I1 z# v$ M( I/ j / 计算P，迭代方向/ /% u5 I* n3 N% Z) z( J: V/ ffor(i=0;in;i+) p=0.0;) O1 k& T1 a7 B _. : E for(i=0;in;i+) Y) z) c3 : J( Z) . p4 k for(j=0;jn;j+)4 I! A y+ $ j& N6 n! b p=p-Hj*

30、g1j; ; _; h2 4 3 R; i8 R8 9 F& V+ l/迭代方向计算完毕/ F) m D+ W$ dk+;0 n7 O% |; Q, W/ for(i=0;in;i+)5 l3 |3 w0 a* . K; n. x* q) g( x0=x1;9 & H: J; S+ J Y/ j( T g0=g1;$ . 5 U* 1 S. Q7 V5 I9 / O- ; , y* g* k/ B# ?comput_grad(pf,n,x1,g1);* l$ S7 M3 t5 h8 n; l0 % g_norm=0.0; /计算判断收敛对象) 8 * D; S4 O* & for(i=0;i=

31、1e9)4 0 _% L, h5 E0 P coute);/ I1 d( / l7 + M; X; D( n for(i=0;in;i+) min_point=x1;f& o3 v7 k8 a! K; 2 b) delete g0; . q& # b+ I1 p5 p$ |, ?/ F5 Y+ Y delete g1;8 z- z3 s0 s% |- x I1 |: V delete dg;2 g( v$ S( s, V: M: X* u3 L) Z delete p; t% N s( f7 G; R5 V, . f& o delete x0;! l# x2 & e1 Rq : % M del

32、ete x1;2 i6 |) M: D) f4 x# w delete dx;t. q$ x S/ t; r2 g/ C for (i=0; in; i+) delete H;% t- M9 2 b7 R8 h: G2 H delete H;8 K1 Y: O R; R1 H/ I delete gH;% & F+ g7 r0 q, S: 7 Z* sq delete Hg;5 o2 U x; y- IS return 0; 4 R1 w, a1 f; D% fL4 m# 2 i& p/o; c; J& F9 Fdouble fun(double *x , double *ys )/返回目标函

33、数值, Q% , ?) C# C0 q- x/ w9 d6 d$ W?% o9 : _. 9 t* t* double f;, ? w ; d1 G- s! r( v8 u ys0=-x0*x0-x1*x1-x2*x2-x3*x3-x0+x1-x2+x3+A-det;+ r5 b5 t5 t& p ys1=-x0*x0-2*x1*x1-x2*x2-2*x3*x3+x0+x3+B-det;2 P K2 H j) j a/ Y0 Z ys2=-2*x0*x0-x1*x1-x2*x2-2*x0+x1+x3+C-det; W% U- k7 S; f. u& y. If=x0*x0+x1*x1+x2*x

34、2*2+x3*x3-5*x0-5*x1-21*x2+7*x3+D;6 r+ i# ?8 K! q. I* q F, V return f;1 E. b+ vU& n O , l q+ j# j l0 T) A G( Y, iX+ x3 Sdouble funP(double *p) /构造混合外点惩罚函数 返回惩罚函数值7 W/ : h! p k1 V2 S: k. s9 j3 N) j2 q2 ! & v$ P+ xint i;7 g& A9 _% w- _$ L3 i# j& p! O3 G double f;% E* y/ k/ f+ u( y! double *ys=new doubl

35、enc; L D) r. . 7 p) c& P9 P; i! B/ f f=fun(p,ys);$ * p# c X s0 IZ, i5 - ) j% G! N: t( & P5 G + S if(nc!=0)7 Y* p; + * w, Q ! u5 E. 2 N! c% jfor(i=0;in;i+)/ E3 M. D$ v8 Z$ z8 f! t! P6 if(ysf2)/向前搜索, _, M1 R) O! ; X/ b/ ) y5 t3=t2+h;f3=funPt(start,direction,t3); |4 j3 G, I- b0 Tu5 k2 Q8 J5 I while(1)3 |. |, A3 O3 m# I3 X+ Z if(f2=f3) a=t1,b=t3;break; I