黎健玲,王培培

(广西大学数学与信息科学学院,广西南宁　530004)



二次半定规划一个原始对偶路径跟踪算法*

黎健玲**,王培培

(广西大学数学与信息科学学院,广西南宁530004)

(College of Mathematics and Information Science,Guangxi University,Nanning,Guangxi,530004,China)

摘要:本文提出求解二次半定规划的一个基于H..K..M方向的原始对偶路径跟踪算法.文中首先导出确定H..K..M方向的线性方程组,并证明该搜索方向的存在唯一性;然后给出算法的具体步骤,并证明算法产生的迭代点列落在中心路径的某个邻域内.最后采用Matlab(R2011b)数学软件编程对算法进行数值试验.数值结果表明算法是有效的.

关键词：二次半定规划原始对偶算法路径跟踪中心路径

0　引言

【研究意义】(线性)半定规划(简记SDP)既是定义在半正定矩阵锥上的规划,也是线性规划、凸二次规划、二阶锥规划等的一种推广,在控制论、组合优化等诸多领域有着广泛的应用.而二次半定规划是半定规划的拓展，它在最优控制、证券、金融风险分析等领域中都有着广泛的应用.【前人研究进展】原始对偶内点算法是求解半定规划的有效方法[1-10].由于确定搜索方向的方法不同,因此导出多种形式的原始对偶内点法.目前常用的搜索方向有以下3种:(i)AHO搜索方向[2];(ii)H..K..M搜索方向[3-5];(iii)Nesterov-Todd搜索方向(简称NT方向)[6-7].基于SDP的原始对偶内点算法,文献[11]提出一个求解二次半定规划的势减少(potential reduction)算法;文献[12]结合Dikin型方向和牛顿中心步提出一个求解二次半定规划的预估校正算法;文献[13]提出一个求解二次半定规划的投影收缩算法;文献[14]将基于NT方向的原始对偶内点算法推广到二次半定规划.【本研究切入点】本文考虑如下二次半定规划问题(简记为QSDP):

aTHsvec(X)+C·X

(0.1)

s.t.Ai·X =bi,i=1,…,m;

XT=X0,

1　二次半定规划问题的对偶规划

利用凸规划问题的Wolfe对偶理论,可得QSDP(0.1)的对偶规划(简记为DSDP)为

(1.1)

分别记问题(0.1),(1.1)的可行集为FP,FD,即

FP={X|Ai·X =bi,X0},

记FP和FD的严格内部为F°P和F°D.

易知DSDP(1.1)也是一个凸二次半定规划问题.对于任意的(X,y,Z)∈FD,若X∈FP,则对偶间隙为

(1.2)

2　H..K..M搜索方向

假设：

(A1)设矩阵组{Ai,i=1,…,m}线性无关.

(A2)Slater条件成立,即存在X≻0,Z≻0和y∈Rm使得X∈FP,(X,y,Z)∈FD.

由(1.2)式以及文献[13]的定理3.5知,在假设(A2)下,X∈FP是QSDP(0.1)的最优解当且仅当存在y∈Rm,Z∈Sn,使得

(2.1a)

XZ=0.

(2.1b)

点对(X,Z)如果满足X∈FP,(X,y,Z)∈FD且

XZ=σμI,

Ai·X =bi,i=1,…,m,

(2.2a)

(2.2b)

XZ=σμI.

(2.2c)

设当前迭代点为(X,y,Z)，且X∈F°P,(y,Z)∈F°D,用一步Newton法求解非线性方程组(2.2)得到的系统如下:

Ai·ΔX=0,i=1,…,m,

(2.3a)

(2.3b)

ΔXZ+XΔZ=σμI-XZ-ΔXΔZ.

(2.3c)

忽略(2.3c)式中的ΔXΔZ,得到

ΔXZ+XΔZ=σμI-XZ.

(2.4)

注意到X,Z不可交换,即XZ≠ZX,而内点法的一个关键是产生的矩阵ΔX,ΔZ要满足对称性.文献[9]给出对称化矩阵M的一般通式:

Ai·ΔX=0,i=1,…,m,

(2.5a)

(2.5b)

(2.5c)

GT=(svec(A1),…,svec(Am)),

(2.6a)

(2.6b)

(2.6c)

(2.6d)

(2.7)

证明已知

rank(FHTH+E)=rank(F-1(FHTH+E))=rank(HTH+F-1E),

(2.8)

往证F-1E是对称正定矩阵.由对称Kronecker积的性质(文献[1]附录)知

(Z⊗sX-1)]=X-1⊗sZ,

(2.9)

因X-1和Z是对称正定阵,故X-1⊗sZ,即F-1E是对称正定阵.又因为HTH是对称半正定阵,因此HTH+F-1E是对称正定阵,结合(2.8)式即知矩阵FHTH+E是可逆的.

证明因为

(FHTH+E)-1F=(F-1(FHTH+E))-1=

(HTH+F-1E)-1,

引理2.3假设(A1)和(A2)成立,X∈F°P,(X,y,Z)∈F°D,则线性方程组(2.7)有唯一解.

证明首先证明线性方程组(2.7)的系数矩阵可逆,即证明方程组(2.7)对应的齐次线性方程组

(2.10)

只有平凡解.记

B=

利用块高斯消去法对(2.10)式进行化简,得到

(2.11)

于是有

(G(FHTH+E)-1FGT)Δy=0.

(2.12)

由引理2.2知(FHTH+E)-1F是对称正定阵,又因为G为行满秩,所以G(FHTH+E)-1FGT也是对称正定阵,于是由(2.12)式知Δy=0.

将Δy=0代入(2.11)式,即得

svec(ΔX)=(FHTH+E)-1FGTΔy=0,

svec(ΔZ)=-F-1E svec(ΔX)=0,

所以(2.10)式只有平凡解,从而(2.7)式的解存在并且唯一.

引理2.4假设(A1)和(A2)成立,设X∈F°P,(X,y,Z)∈F°D,(ΔX,Δy,ΔZ)是线性方程组(2.7)对于某个给定的矩阵W(∈Rn×n)的解,则以下结论成立:

①ΔZ·ΔX≥0;

②X·ΔZ+Z·ΔX=Tr(W);

(X+αΔX)·(Z+αΔZ)=(1-α+α σ)(X·

Z)+α2(ΔX·ΔZ).

证明①根据(2.5b)式和(2.5a)式,并结合矩阵内积的性质可得

最后一个不等式成立是因为HTH是对称半正定矩阵.

2(X·ΔZ+Z·ΔX),

即X·ΔZ+Z·ΔX=Tr(W).

③根据矩阵内积的线性性质以及矩阵迹的性质,再结合结论(2)得到

(X+αΔX)·(Z+αΔZ)=X·Z+α σμn-αTr(XZ)+α2(ΔX·ΔZ)=X·Z+α σμn-

α(X· Z)+α2(ΔX·ΔZ)=(1-α+α σ)(X·Z)+α2(ΔX·ΔZ),

即结论(3)成立.

3　H..K..M搜索方向的计算

H..K..M搜索方向(ΔX,Δy,ΔZ)可以通过求解线性方程组(2.7)得到,而方程组(2.7)包含m+n(n+1)个线性方程,通过块高斯消去法将其化简为如下的方程

(G(FHTH+E)-1FGT)Δy=-G(FHTH+E)-1rc,

(3.1a)

GTΔy-(F-1E+HTH)svec(ΔX)=-F-1rc,

(3.1b)

Esvec(ΔX)+Fsvec(ΔZ)=rc,

(3.1c)

记M=G(FHTH+E)-1FGT,h=-G(FHTH+E)-1rc,则(3.1a)式可改写为

MΔy=h.

(3.2)

可证M是对称正定阵.事实上,因为

M=G(F-1(FHTH+E))-1GT=G(HTH+F-1E)-1GT=G(HTH+(X-1⊗sZ))-1GT,

上述最后一个等式成立是由(2.9)式得到.由X-1,

Z是对称正定阵知X-1⊗sZ也是对称正定阵,而HTH是对称半正定阵,因此(HTH+(X-1⊗sZ))-1是对称正定阵,故M是对称正定阵.

于是可利用Cholesky分解求解线性方程组(3.2)得到Δy,进而由(3.1b)式和(3.1c)式可求出

ΔX=smat((FHTH+E)-1(FGTΔy+rc)),

(3.3a)

ΔZ=smat(HTHsvec(ΔX)-GTΔy).

(3.3b)

4　短步原始对偶路径跟踪算法

这一节我们将给出基于H..K..M方向的短步原始对偶路径跟踪算法.

算法的具体步骤如下:

算法A

当μk>2-ημ0时执行以下步骤:

步骤1记X=Xk,(y,Z)=(yk,Zk),μ=μk.

步骤3求解线性方程组(2.7)得到H..K..M搜索方向(ΔX,Δy,ΔZ).

步骤4选择α≥0,使得

算法产生的迭代点列将落在中心路径的如下邻域内:

(4.1)

(X(α),Z(α),y(α))≡(X,Z,y)+α(ΔX,ΔZ,Δy),

(4.2a)

μ(α)≡(X(α)· Z(α))/n,

(4.2b)

(4.2c)

则有

(4.3)

其中I为n阶单位阵.

证明假设α≥0是给定的.根据引理2.4(3)有

X(α)·Z(α)=(X+αΔX)·(Z+αΔZ)=(1-α+α σ)(X·Z)+α2(ΔX·ΔZ),

所以由(4.2b)即知

因此,

(1-α)(XZ-μI)+α(XZ-σμI)+α(XΔZ+

(4.4)

类似可以得到

Z(α)X(α)-μ(α)I=(1-α)(ZX-μI)+

(4.5)

于是,根据(4.2c)式,(4.4)式,(4.5)式得到

最后一个等式成立是根据(2.5c)式.

‖A+B‖F≤‖A‖F+‖B‖F,‖AB‖F≤‖A‖F‖B‖F,‖AB‖F≤‖A‖F‖B‖.

(4.6)

(4.7)

证明根据(2.5c)式得

由引理2.4(1)知ΔZ·ΔX≥0,再根据(4.6)式以及引理4.2的已知条件得到

所以

以下定理说明本文的短步路径跟踪算法产生的迭代点列将落在中心路径邻域NF(γ)内.

(4.8)

证明因为

(4.9)

从而

(4.10)

(4.11)

根据(4.2b)式,(4.2c)式和(4.3)式得

再根据(4.6)式,有

(4.12)

由矩阵Frobenious范数的定义知

即

将上式代入(4.12)式,再结合(4.11)式,(4.9)式和(4.8)式,得

(4.13)

根据(4.10)式,(4.9)式和定理4.1的条件,得

(4.14)

5　数值试验

采用Matlab(R2011b)数学软件编程对算法A进行数值试验,程序运行环境为Windows7(64bite),Intel(R)Core(TM)i3-2330M CPU@2.20GHZ,RAM:4G.

考虑如下两类测试问题:(1)随机二次半定规划问题;矩阵Ai,Hj,C都是随机生成的对阵矩阵,初始点(X,y,Z)的选取需在FP,FD内.记该类测试问题为RP.(2)最近相关矩阵Nearest correlation matrix(简记为NCM)问题[16-17].

s.t.diag(X)=e,

其中G是对称正定阵,对角线上的元素为1,对角线外的元素在-1到1之间.e是所有元素都为1的n维向量.

NCM问题具有广泛的应用,例如在市场营销以及经济学方面.但是遗憾的是,由于缺乏数据等信息，不能得到一个完整的矩阵G,即矩阵G的一些元素是未知的.于是通过求解NCM问题得到一个有效的、且与G最近的相关矩阵.在数值测试中,取G是随机矩阵且G∶=(1-α)B+α E[18],其中α∈(0,1),E是元素属于[-1,1]的随机对称阵,B是具有指定特征值的测试矩阵,在NCM问题的数值测试中取n=10,20,30.

利用对称Kronecker积的矩阵形式

将矩阵的对称Kronecker积转化为矩阵的Kronecker积计算,其中Q的取法详见文献[1]附录.

在数值测试中,取算法A中的参数δ=0.3,γ=0.3,终止准则是Zk· Xk≤10-6.

测试的数值结果见表1(所有结果都是每个维数测试10次,然后取平均值得到的),其中表中符号含义如下:prob,测试问题的编号;n,矩阵X的维数;m,不等式约束的个数;Itr,算法迭代次数;Nf,目标函数值的计算次数;time(s),算法迭代所需CPU时间(s).

表1数值结果

Table 1Numerical results

probnmItrNftime(s)RP501261261.009955e+00601461462.871151e+00NCM10101631631.842785e+0120202442441.767764e+0230303073071.529803e+03

上述数值结果表明本文提出的基于H..K..M方向的路径跟踪算法是有效的，能在合理的时间内求解中等规模的NCM问题.为进一步提高本文算法的数值效果，我们将在后续的研究中深入探讨步长α的更有效的确定方法，即线搜索技术.

参考文献:

[1]TODD M J,TOH K C,TÜTÜNCÜ R H.On the nesterov-todd direction in SDP[J].SIAM Journal on Optimization,1998,8(3):769-796.

[2]ALIZADEH F,HAEBERLY J P,OVERTON M L.Primal-dual interior-point methods for semidefinite programming:Convergence rates,stability and numerical results[J].SIAM Journal on Optimization,1998,8:746-768.

[3]KOJIMA M,SHINDOH S,HARA S.Interior-point me- thods for the monotone semidefinite linear complementarity problem in symmetric matrices[J].SIAM Journal on Optimization,1997,7(1):86-125.

[4]HELMBERG C,RENDL F,VANDERBEI R J,et al.An interior-point methods for semidefinite programming[J].SIAM Journal on Optimization,1996,6(2):342-361.

[5]MONTEIRO R D C.Primal-dual path-following algorithms for semidefinite programming[J].SIAM Journal on Optimization,1997,7(3):663-678.

[6]NESTEROV Y E,TODD M J.Primal-dual interiorpoin- t methods for self-scaled cones[J].SIAM Journal on Optimization,1998,8(2):324-364.

[7]NESTEROV Y E,TODD M J.Self-scaled barriers and interior-point methods for convex programming[J].Mathematics of Operations Research,1997,22(1):1-42.

[8]STURM J F,ZHANG S Z.Symmetric primal-dual path-following algorithms for semidefinite programming[J].Applied Numerical Mathematics,1999,29(3):301-315.

[9]ZHANG Y.On extending some primal-dual interior- point algorithms from linear programming to semidefinite programming[J].SIAM Journal on Optimization,1998,8(2):365-386.

[10]NESTEROV Y,NEMIROVSKII A S.Interior-Point Polynomial Methods in Convex Programming[M].Philadelphia PA:SIAM,1994.

[11]NIE J W,YUAN Y X.A potential reduction algorithm for an extended SDP problem[J].Science in China Series A Mathematics,2000,43(1):35-46.

[12]NIE J W,YUAN Y X.A Predictor-corrector algorith- m for QSDP combining dikin-type and newton centering steps[J].Annals of Operations Research,2001,103(1/2/3/4):115-133.

[13]关秀翠,刁在筠.二次半定规划问题及其投影收缩算法[J].全国高等学校计算数学学报,2002,24(2):97-108. GUAN X C,SIAO Z Y.The quadratic semi-definite programming problem and its projection and contraction algorithm[J].Numerical Mathematics a Journal of Chinese Universities,2002,24(2):97-108.

[14]徐凤敏,徐成贤.求解二次半定规划的原对偶内点算法[J].工程数学学报,2006,23(4):590-598. XU F M,XU C X.Primal-dual algorithm for quadratic semi-definite programming[J].Chinese Journal of Engineering Mathematics,2006,23(4):590-598.

[15]HORN R A,JOHNSON C R.Matrix Analysis[M].New York:Cambridge University Press,1985.

[16]HIGHAM N J.Computing the nearest correlation matrixa problem from finance[J].IMA Journal of Numerical Analysis,2002,22(3):329-343.

[17]TOH K C.An inexact primal-dual path following algorithm for convex quadratic SDP[J].Mathematical Programming,2008,112(1):221-254.

[18]ZHAO X Y.A Semismooth Newton-CG Augmented Lagrangian Method for Large Scale Linear and Convex Quadratic SDPs[D].Singapore:National University of Singapore,2009.

(责任编辑：米慧芝)

A Primal-dual Path-following Algorithm for Quadratic Semi-definite Programming

LI Jianling,WANG Peipei

Key words：quadratic semi-definite programming,primal-dual,algorithm,path-following,central path

Abstract：A primal-dual path-following algorithm based on H..K..M direction for quadratic semi-definite programming problems(QSDP) is proposed.Firstly,the system of linear equations yielding the H..K..M direction are derived,and the existence and uniqueness of the search direction are shown;Secondly,the algorithm is described in detail.We show that the iterates generated by the algorithm can fall into some neighborhood of the central path under some mild conditions.Finally,a preliminary numerical experiment is performed for the algorithm by using Matlab (R2011b) mathematical software,and the numerical results show that the proposed algorithm is effective.

收稿日期：2016-08-05

作者简介：黎健玲(1965-),女,博士,教授,主要从事最优化理论与方法研究,E-mail:jianlingli@126.com。

中图分类号:C934

文献标识码：A

文章编号：1005-9164(2016)05-0396-08

修回日期：2016-09-25

*国家自然科学基金项目(No.11561005)和广西自然科学基金项目(2016GXNSFAA380248,2014GXSFFA118001)资助。

**通信作者。

广西科学Guangxi Sciences 2016,23(5):396～403

网络优先数字出版时间：2016-11-21【DOI】10.13656/j.cnki.gxkx.20161121.007

网络优先数字出版地址：http://www.cnki.net/kcms/detail/45.1206.G3.20161121.1520.014.html