文章编号:1671 1114(2009)03 0018 04快速投影Hessian 矩阵算法收稿日期:2008 03 10基金项目:天津市高校发展基金项目(20060402)作 者:汤大林(1965 ),男,高级工程师,主要从事数学建模及应用方面的研究.汤大林(天津理工大学理学院,天津300191)摘 要:分析了求解等式约束非线性规划问题的投影H essian 矩阵算法,找出了算法两步Q 超线性收敛的原因,并用BY RD 的例子说明此算法的收敛效果较差,即甚至不是线性收敛;对算法进行了合理的改进,并用改进后的算法求解BY RD 问题,得到了满意的收敛效果,即Q 超线性收敛.借助数值试验验证了改进算法的快速收敛性.关键词:等式约束非线性规划;投影H essian 矩阵算法;超线性收敛中图分类号:O 221.2 文献标识码:AQ uick projection method with H essian matrixT AN G Dalin(School of Science,T ianjin University of Techn ology,T ian jin 300191,China)Abstract:T he project ion method wit h H essian mat rix used to so lve nonlinear prog ramming w ith equality co nstr aint is analyzed and the r easo n w hy the method is superlinear conver gent by two steps is found o ut.Its bad converg ent effect at linear ity is illuminated by BYRD's example.T he method is impr ov ed and quickly super linear conver gence o f the impr ov ed metho d is illuminated using BY RD's ex ample.T he quickly co nv erg ent effect o f t he impro ved method is verified by a numerical experiment.Key words:nonlinear pr og ramming w ith equality constra int ;project ion method wit h H essian mat rix ;super linearconver gence1 投影Hessian 矩阵算法的缺点考察等式约束非线性规划问题: m in x Rnf (x ),约束c(x)=0,(1)其中,目标函数f (x ):R n !R,约束c(x):R n !R m 是二次可微函数,且m ∀n,即m 个等式约束.为叙述方便,引入如下记号:x =(x (1),x (2),#,x (n)),c(x)=(c (1)(x ),c (2)(x ),#,c (m)(x ))T, g(x)= f (x )=( f x (1), f x (2),#, f x (n))T,A(x)= c(x)=c (1)x (1)c (2)x (1)# c (m) x (1) c (1) x (2) c (2) x (2)# c (m) x (2) ! c (1) x (n)c (2) x (n)#c (m) x (n ),L (x, )=f (x )-∃mi=1(i)c (i)(x ),其中, (i)为拉格朗日乘子,i =1,2,#,m.将A(x)QR 分解为A(x)=(y(x),z(x))(R(x)O ),其中y (x),z(x)均为n 阶正交矩阵,R(x)为m 阶上三角矩阵.V ol.29N o.3Jul.2009第29卷 第3期2009年7月 天津师范大学学报(自然科学版)Jour nal of T ianjin N orma l U niver sity (N atural Science Edit ion)定理1[1] 1)x*是问题(1)的解的一阶必要条件是存在拉格朗日乘子 ,满足q(x, )= x, L(x, )=(g(x)-A(x)c(x))=0,(2)记W(x, )= 2x L(x, )= 2f(x)-∃m i=1 (i) c2(i)(x).2)x*是问题(1)的解的二阶必要条件是投影H essian阵Z T(x)W(x, )Z(x)半正定.3)x*是问题(1)的解的二阶充分条件是式(2)成立,且Z T(x)W(x, )Z(x)正定.假设x*是问题(1)的局部解,且A*=A(x*)列满秩.为求解问题(1),T.F.Colemen等设计了投影H essian矩阵修正算法[2].下面由该算法的简化推导过程说明该算法的缺点.由定理1的1)知,求解问题(1)即求式(2)的零点,可用牛顿法.为了叙述方便,在点(x k, k)处,记g k=g(x k),c k=c(x k),W k=W(x k, k),其他类似.W k-A k A T k Ox kk=-g k-A k kc k,x k+1 k+1=x kk+x kk,∀W k-A k A T k O p kk=-g kc k,x k+1=x k+x k,∀Q T k W k-A kA T k OQ k Q T kx kk+1=-Q T kg kc k,其中,x k=y k p yk +z k p zk,Q T k=(y k,z k)OO I,∀y T k W k y k y T k W k z k-R kz T k W k y k z T k W k z k OR T k O Op ykp zkk+1=-y T k g kz T k g kc k,(3)式(3)是m+n个未知数的线性方程组,由此可知R k p yk =-c k∀得出p yk,(4a)z T k W k y k p yk +z T k W k z k p zk=-z T k g k∀得出p zk,(4b)x k=y k p yk +z k p zk,x k+1=x k+x k,(4c)R k k+1=y T k g k+y T k W k x k,(4d)式(4d)为 k+1的二阶估计,也可由其一阶估计即最小二乘估计给出,即由m in%A k+1 -g k+1%2给出,k+1=R-1k+1y T k+1g k+1,(5)将式(5)代入式(4d)得z T k W k y k z T k W k z kR k Op ykp zk=-z T k g kc k,(6)由定理1的二阶必要条件知z T k W k z k为半正定阵,其计算量很大.若用一个正定矩阵B k近似z T k W k z k似乎是可行的方法,但是这样又会导致z T k W k y k无法处理,解决该问题的最简单的办法就是直接略去z T k W k y k,如此便得到投影H essian矩阵算法.算法1 投影H essian矩阵算法.Step1 1)选择精度!,置k=0.2)置初始点x0,B0,其中B0为正定近似投影阵.3)分别计算f0,g0,c0,A0,分解A0=(y0 z0)R0O.Step2 解方程组R k p yk=-c k(7)得出p yk,解方程组B k p zk=-z T k g k(8)得出p zk,置x k+1=x k+y k p yk+z k p zk.Step3 计算f k+1,g k+1,c k+1,A k+1,分解A k+1=(y k+1 z k+1)R k+1O,计算%z T k+1c k+1%2.Step4 如果%z T k+1c k+1%2<!,则结束.否则,用DFT或BFGS修正公式计算B k+1,并置k=k+1,返回Step2.本算法被证明[3 4]在适当的假设条件下两步Q 超线性收敛于问题(1)的一个局部极小值点x*,即%x k+1-x*%2%x k-x*%2!0,k!&.关于此算法能否一步Q 超线性收敛,下面通过BYRD的一个例子给出了否定的回答.例min f(x)=12x2(1)-∀x(1)x(2)+12x2(2)-(x(1)-∀)33∀,∀为常数,约束c(x)=12-x(2)-1=0.(9)由于约束相当于x(2)=1,不难看出本例的最优解为x*=(∀,1)T.∋19∋第29卷 第3期 汤大林:快速投影H essian矩阵算法现用算法1求解,取x0=(∀,1+#)T,B0为投影H essian矩阵,于是x1=(∀+∀#,1+#2)T,x2=(∀,1+#4)T,x3=(∀+∀#4,1+#8)T,#,当k为偶数时,x k=(∀,1+#2k)T,x k+1=(∀+∀#2k,1+#2k+1)T,易知%x k-x*%2=#2k,%x k+1-x*%2=#2k∀2+#2k,%x k+1-x*%2%x k-x*%2=∀2+#2k,(10)由式(10)知,只有当|#|∀1时算法才收敛,但只是局部收敛.另外当∀(0时,算法即使收敛也不是Q 超线性收敛,当∀>1时甚至不是线性收敛.2 算法的改进从上例看出算法1的收敛效果不能令人满意,这是因为算法1在式(6)中简单地去掉了z T k W k y k,这样势必导致不精确.重新考虑z T k W k y k,在式(6)中将其移到等式右边,并做近似,则式(6)变为R k p yk=-c k,z T k W k z k-p zk =-z T k g k-z T k W k y k p yk)-z T k x L(x k+y k p yk , k),(11)如此处理比直接从式(6)中简单地去掉z T k W k y k提高了一阶精度,从而得到算法1的改进算法.算法2 改进的投影H essian矩阵算法Step1 同算法1.Step2 解方程组R k p yk =-c k,得出p yk,置x∗k=x k+y k p yk.解方程组B k p zk =-z T k x L(x∗k, k),得出p zk,置x k+1=x k+y k p yk +z k p zk.Step3 同算法1.Step4 同算法1.下面用算法2来求解BYRD的例子,说明算法2的优点.易知A k=( 01(2-x(2)k)2)=0110(1(2-x(2)k)20),y k=01,z k=1.取B k=2fx2(1)k=1-2(x(1)k-∀)∀,g k=x(1)k-∀x(2)k-(x(1)k-∀)2∀-∀x(1)k+x(2)k.第一次迭代:x0=(∀,1+#)T,B0=1,R0=1(1-#)2,p y=-(R T0)-1c0=-(1-#)2(11-#-1)=#(#-1),x∗0=(∀,1+#)T+(0,1)T#(#-1)=(∀,1+#2)T,g∗0=(-∀#2,1+#2-∀2)T,p z=-B-10(1,0)(-∀#2,1+#2-∀2)T=∀#2,x1=x∗0+z0p z=(∀+∀#2,1+#2)T.第二次迭代:B1=1-2#2,R1=1(1-#2)2,c1=#21-#2,p y1=-(R T1)-1c1=-(1-#2)#2,x∗1=(∀+∀#2,1+#2)T+(0,1)T#2(#2-1)=(∀+∀#2,1+#4)T,g∗1=∀#2-2∀#4-∀(1+#2)+1+#4,p z1=-B-11z1g∗1=-∀#2,x2=x∗1+z1p z1=(∀,1+#4)T.第三、四次迭代:x3=(∀+∀#8,1+#8)T,x4=(∀,1+#16)T,#.由归纳法知:x k=(∀,1+#2k)T, k为偶数,(∀+∀#2k,1+#2k)T,k为奇数.故%x k-x*%=#2k, k为偶数,#2k∀+1,k为奇数.%x k+1-x*%2%x k-x*%2=#2k∀2+1, k为偶数,#2k∀2+1, k为奇数,k!&0.(12)式(12)表明,无论∀取何值,只要|#|<1,算法2即为一步Q 超线性收敛.3 数值试验利用BYRD的例子,应用数学软件M athematica6.0编程,分两种情形对两种算法进行数值对比试验.∋20∋天津师范大学学报(自然科学版)2009年7月情形1:取迭代初值x0=(0,1)T,∀=0.2+0.2+3,收敛精度为%x k-x*%∀max(10-9,10-8%x k%)且3f(x k)∀10-9.试验结果见表1.情形2:取迭代初值x0=(∀,1)T,∀=0.2+0.2+3,收敛精度为%x k-x*%∀max(10-9,10-8%x k%)且3f(x k)∀10-9.试验结果见表2.两种情形的对比试验结果均表明:1)同等条件下,算法2比算法1收敛速度快;2)不论∀<1或>1,算法2均收敛,而算法1只在∀<1时收敛.表1 情形1∀算法2(x0=(0,1)T)算法1(x0=(0,1)T)迭代次数x*最优值迭代次数x* 最优值0.22x(1)!0,x(2)!10.513333 5x(1)!0.07779,x(2)!0.01556 0.00594670.42x(1)!0,x(2)!10.553333 5x(1)!0.16473,x(2)!0.06589 0.02224930.62x(1)!0,x(2)!10.620000 5x(1)!0.27502,x(2)!0.16502 0.04327120.82x(1)!0,x(2)!10.713333 5x(1)!0.44287,x(2)!0.35429 0.05428281.02x(1)!0,x(2)!10.833333 23x(1)!1.00000,x(2)!1.00000 01.22x(1)!0,x(2)!10.980000 2x(1)!2.71726,10103,x(2)!8.09698,1018-5.5730120,10309 1.42x(1)!0,x(2)!11.153330 101x(1)!2.28077,10102,x(2)!1.37916,1020-2.8248700,10306 1.62x(1)!0,x(2)!11.353330 101x(1)!1.69299,10103,x(2)!9.95408,1019-1.0109318,103091.82x(1)!0,x(2)!11.580000 101x(1)!1.88499,10103,x(2)!1.10902,1020-1.2403130,103092.02x(1)!0,x(2)!11.833330 101x(1)!1.30476,10103,x(2)!-2.22676,1020-3.7019995,10308 2.22x(1)!0,x(2)!12.113330 101x(1)!1.48876,10103,x(2)!2.12286,1019-4.9995037,10308 2.42x(1)!0,x(2)!12.420000 101x(1)!1.81431,10103,x(2)!1.23266,1020-8.2947708,10308 2.62x(1)!0,x(2)!12.753330 101x(1)!5.28182,10103,x(2)!3.32167,1019-1.8891057,103102.82x(1)!0,x(2)!13.113330 101x(1)!2.77578,10103,x(2)!4.03269,1019-2.5460911,103093.02x(1)!0,x(2)!13.500000 101x(1)!3.20881,10103,x(2)!5.90223,1019-3.6710456,10309表2 情形2∀算法2(x0=(∀,1)T)算法1(x0=(∀,1)T)迭代次数x*最优值迭代次数x* 最优值0.22x(1)!0.2,x(2)!1 0.48 6x(1)!0.07779,x(2)!0.01556 0.00594670.42x(1)!0.4,x(2)!1 0.42 6x(1)!0.16473,x(2)!0.06589 0.02224930.62x(1)!0.6,x(2)!1 0.32 6x(1)!0.27502,x(2)!0.16502 0.04327120.82x(1)!0.8,x(2)!1 0.18 5x(1)!0.44287,x(2)!0.35429 0.05428281.02x(1)!1.0,x(2)!1 0 2x(1)!1.00000,x(2)!1.00000 01.22x(1)!1.2,x(2)!1-0.22 2x(1)!4.99756,10103,x(2)!-3.29194,1018-3.4671348,10310 1.42x(1)!1.4,x(2)!1-0.48 101x(1)!3.81463,10103,x(2)!-2.10787,1018-1.3216242,10310 1.62x(1)!1.6,x(2)!1-0.78 101x(1)!3.58878,10103,x(2)!3.54846,1019-9.6293836,103091.82x(1)!1.8,x(2)!1-1.12 101x(1)!3.64849,10103,x(2)!6.18958,1019-8.9938505,103092.02x(1)!2.0,x(2)!1-1.50 101x(1)!4.09807,10103,x(2)!3.84194,1019-1.1470649,10310 2.22x(1)!2.2,x(2)!1-1.92 101x(1)!3.28328,10103,x(2)!1.57167,1020-5.3626603,10309 2.42x(1)!2.4,x(2)!1-2.38 101x(1)!3.62176,10103,x(2)!1.47605,1020-6.5982074,10309 2.62x(1)!2.6,x(2)!1-2.88 101x(1)!5.14778,10103,x(2)!1.68313,1019-1.7488987,103102.82x(1)!2.8,x(2)!1-3.42 101x(1)!6.50098,10103,x(2)!1.70756,1019-3.2708260,103103.02x(1)!3.0,x(2)!1-4.00 101x(1)!2.19182,10104,x(2)!5.10395,1019-1.1699615,10312参考文献:[1] Broyden C G.T he con vergen ce of a clas s of double rank m inimiz ation algorith ms[J].J In st M ath Appl,1970,6:76 90.[2] C oleman T F,Conn A R.On the local convergence of a qu asiN ew ton m ethod for the n onlinear programm ing problem[J].S IAM J Num er Anal,1984,21:769 775.[3] Fletcher R.Practical methods of optimization[M].NewYork:John Wiley&Sons Inc,1980.[4] Shannon D F.Conditioning of quasi New ton methods for functionmini mization[J].M ath ematics of Com puting,1970,24: 647 656.(责任编校 马新光)∋21∋第29卷 第3期 汤大林:快速投影H essian矩阵算法。