快速投影Hessian矩阵算法

文章编号: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 文献标识码:A

Q uick projection method with H essian matrix

T 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 linear

conver gence

1 投影Hessian 矩阵算法的缺点

考察等式约束非线性规划问题: m in x R

n

f (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 )-∃m

i=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.3

Jul.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 O

x k

k

=-

g k-A k k

c k

,

x k+1 k+1=

x k

k

+

x k

k

,∀

W k-A k A T k O p k

k

=-

g k

c k

,

x k+1=x k+x k,∀

Q T k W k-A k

A T k O

Q k Q T k

x k

k+1

=-Q T k

g k

c k

,

其中,x k=y k p y

k +z k p z

k

,

Q T k=(y k,z k)O

O I

,∀

y T k W k y k y T k W k z k-R k

z T k W k y k z T k W k z k O

R T k O O

p y

k

p z

k

k+1

=-

y T k g k

z T k g k

c k

,(3)

式(3)是m+n个未知数的线性方程组,由此可知

R k p y

k =-c k∀得出p y

k

,(4a)

z T k W k y k p y

k +z T k W k z k p z

k

=

-z T k g k∀得出p z

k

,(4b)

x k=y k p y

k +z k p z

k

,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 k

R k O

p y

k

p z

k

=-

z T k g k

c 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)

R0

O

.

Step2 解方程组

R k p y

k

=-c k(7)

得出p y

k

,解方程组

B k p z

k

=-z T k g k(8)

得出p z

k

,置x k+1=x k+y k p y

k

+z k p z

k

.

Step3 计算f k+1,g k+1,c k+1,A k+1,分解

A k+1=(y k+1 z k+1)

R k+1

O

,

计算

%

z T k+1

c k+1%2

.

Step4 如果

%

z T k+1

c 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)=

1

2

x2(1)-∀x(1)x(2)+

1

2

x2(2)-

(x(1)-∀)3

3∀

,∀为常数,

约束c(x)=1

2-x(2)

-1=0.

(9)

由于约束相当于x(2)=1,不难看出本例的最优

解为x*=(∀,1)T.

∋

19

∋

第29卷 第3期 汤大林:快速投影H essian矩阵算法