黎健玲,王培培

(廣西大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,廣西南寧　530004)

?

二次半定規(guī)劃一個(gè)原始對(duì)偶路徑跟蹤算法*

黎健玲**,王培培

(廣西大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,廣西南寧530004)

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

摘要:本文提出求解二次半定規(guī)劃的一個(gè)基于H..K..M方向的原始對(duì)偶路徑跟蹤算法.文中首先導(dǎo)出確定H..K..M方向的線性方程組,并證明該搜索方向的存在唯一性;然后給出算法的具體步驟,并證明算法產(chǎn)生的迭代點(diǎn)列落在中心路徑的某個(gè)鄰域內(nèi).最后采用Matlab(R2011b)數(shù)學(xué)軟件編程對(duì)算法進(jìn)行數(shù)值試驗(yàn).數(shù)值結(jié)果表明算法是有效的.

關(guān)鍵詞：二次半定規(guī)劃原始對(duì)偶算法路徑跟蹤中心路徑

0　引言

【研究意義】(線性)半定規(guī)劃(簡(jiǎn)記SDP)既是定義在半正定矩陣錐上的規(guī)劃,也是線性規(guī)劃、凸二次規(guī)劃、二階錐規(guī)劃等的一種推廣,在控制論、組合優(yōu)化等諸多領(lǐng)域有著廣泛的應(yīng)用.而二次半定規(guī)劃是半定規(guī)劃的拓展，它在最優(yōu)控制、證券、金融風(fēng)險(xiǎn)分析等領(lǐng)域中都有著廣泛的應(yīng)用.【前人研究進(jìn)展】原始對(duì)偶內(nèi)點(diǎn)算法是求解半定規(guī)劃的有效方法[1-10].由于確定搜索方向的方法不同,因此導(dǎo)出多種形式的原始對(duì)偶內(nèi)點(diǎn)法.目前常用的搜索方向有以下3種:(i)AHO搜索方向[2];(ii)H..K..M搜索方向[3-5];(iii)Nesterov-Todd搜索方向(簡(jiǎn)稱NT方向)[6-7].基于SDP的原始對(duì)偶內(nèi)點(diǎn)算法,文獻(xiàn)[11]提出一個(gè)求解二次半定規(guī)劃的勢(shì)減少(potential reduction)算法;文獻(xiàn)[12]結(jié)合Dikin型方向和牛頓中心步提出一個(gè)求解二次半定規(guī)劃的預(yù)估校正算法;文獻(xiàn)[13]提出一個(gè)求解二次半定規(guī)劃的投影收縮算法;文獻(xiàn)[14]將基于NT方向的原始對(duì)偶內(nèi)點(diǎn)算法推廣到二次半定規(guī)劃.【本研究切入點(diǎn)】本文考慮如下二次半定規(guī)劃問題(簡(jiǎn)記為QSDP):

aTHsvec(X)+C·X

(0.1)

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

XT=X0,

1　二次半定規(guī)劃問題的對(duì)偶規(guī)劃

利用凸規(guī)劃問題的Wolfe對(duì)偶理論,可得QSDP(0.1)的對(duì)偶規(guī)劃(簡(jiǎn)記為DSDP)為

(1.1)

分別記問題(0.1),(1.1)的可行集為FP,FD,即

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

記FP和FD的嚴(yán)格內(nèi)部為F°P和F°D.

易知DSDP(1.1)也是一個(gè)凸二次半定規(guī)劃問題.對(duì)于任意的(X,y,Z)∈FD,若X∈FP,則對(duì)偶間隙為

(1.2)

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

假設(shè)：

(A1)設(shè)矩陣組{Ai,i=1,…,m}線性無關(guān).

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

由(1.2)式以及文獻(xiàn)[13]的定理3.5知,在假設(shè)(A2)下,X∈FP是QSDP(0.1)的最優(yōu)解當(dāng)且僅當(dāng)存在y∈Rm,Z∈Sn,使得

(2.1a)

XZ=0.

(2.1b)

點(diǎn)對(duì)(X,Z)如果滿足X∈FP,(X,y,Z)∈FD且

XZ=σμI,

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

(2.2a)

(2.2b)

XZ=σμI.

(2.2c)

設(shè)當(dāng)前迭代點(diǎn)為(X,y,Z)，且X∈F°P,(y,Z)∈F°D,用一步Newton法求解非線性方程組(2.2)得到的系統(tǒng)如下:

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,而內(nèi)點(diǎn)法的一個(gè)關(guān)鍵是產(chǎn)生的矩陣ΔX,ΔZ要滿足對(duì)稱性.文獻(xiàn)[9]給出對(duì)稱化矩陣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是對(duì)稱正定矩陣.由對(duì)稱Kronecker積的性質(zhì)(文獻(xiàn)[1]附錄)知

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

(2.9)

因X-1和Z是對(duì)稱正定陣,故X-1?sZ,即F-1E是對(duì)稱正定陣.又因?yàn)镠TH是對(duì)稱半正定陣,因此HTH+F-1E是對(duì)稱正定陣,結(jié)合(2.8)式即知矩陣FHTH+E是可逆的.

證明因?yàn)?/p>

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

(HTH+F-1E)-1,

引理2.3假設(shè)(A1)和(A2)成立,X∈F°P,(X,y,Z)∈F°D,則線性方程組(2.7)有唯一解.

證明首先證明線性方程組(2.7)的系數(shù)矩陣可逆,即證明方程組(2.7)對(duì)應(yīng)的齊次線性方程組

(2.10)

只有平凡解.記

B=

利用塊高斯消去法對(duì)(2.10)式進(jìn)行化簡(jiǎn),得到

(2.11)

于是有

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

(2.12)

由引理2.2知(FHTH+E)-1F是對(duì)稱正定陣,又因?yàn)镚為行滿秩,所以G(FHTH+E)-1FGT也是對(duì)稱正定陣,于是由(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假設(shè)(A1)和(A2)成立,設(shè)X∈F°P,(X,y,Z)∈F°D,(ΔX,Δy,ΔZ)是線性方程組(2.7)對(duì)于某個(gè)給定的矩陣W(∈Rn×n)的解,則以下結(jié)論成立:

①ΔZ·ΔX≥0;

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

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

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

證明①根據(jù)(2.5b)式和(2.5a)式,并結(jié)合矩陣內(nèi)積的性質(zhì)可得

最后一個(gè)不等式成立是因?yàn)镠TH是對(duì)稱半正定矩陣.

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

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

③根據(jù)矩陣內(nèi)積的線性性質(zhì)以及矩陣跡的性質(zhì),再結(jié)合結(jié)論(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),

即結(jié)論(3)成立.

3　H..K..M搜索方向的計(jì)算

H..K..M搜索方向(ΔX,Δy,ΔZ)可以通過求解線性方程組(2.7)得到,而方程組(2.7)包含m+n(n+1)個(gè)線性方程,通過塊高斯消去法將其化簡(jiǎn)為如下的方程

(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是對(duì)稱正定陣.事實(shí)上,因?yàn)?/p>

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

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

Z是對(duì)稱正定陣知X-1?sZ也是對(duì)稱正定陣,而HTH是對(duì)稱半正定陣,因此(HTH+(X-1?sZ))-1是對(duì)稱正定陣,故M是對(duì)稱正定陣.

于是可利用Cholesky分解求解線性方程組(3.2)得到Δy,進(jìn)而由(3.1b)式和(3.1c)式可求出

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

(3.3a)

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

(3.3b)

4　短步原始對(duì)偶路徑跟蹤算法

這一節(jié)我們將給出基于H..K..M方向的短步原始對(duì)偶路徑跟蹤算法.

算法的具體步驟如下:

算法A

當(dāng)μk>2-ημ0時(shí)執(zhí)行以下步驟:

步驟1記X=Xk,(y,Z)=(yk,Zk),μ=μk.

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

步驟4選擇α≥0,使得

算法產(chǎn)生的迭代點(diǎn)列將落在中心路徑的如下鄰域內(nèi):

(4.1)

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

(4.2a)

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

(4.2b)

(4.2c)

則有

(4.3)

其中I為n階單位陣.

證明假設(shè)α≥0是給定的.根據(jù)引理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)

于是,根據(jù)(4.2c)式,(4.4)式,(4.5)式得到

最后一個(gè)等式成立是根據(jù)(2.5c)式.

‖A+B‖F(xiàn)≤‖A‖F(xiàn)+‖B‖F(xiàn),‖AB‖F(xiàn)≤‖A‖F(xiàn)‖B‖F(xiàn),‖AB‖F(xiàn)≤‖A‖F(xiàn)‖B‖.

(4.6)

(4.7)

證明根據(jù)(2.5c)式得

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

所以

以下定理說明本文的短步路徑跟蹤算法產(chǎn)生的迭代點(diǎn)列將落在中心路徑鄰域NF(γ)內(nèi).

(4.8)

證明因?yàn)?/p>

(4.9)

從而

(4.10)

(4.11)

根據(jù)(4.2b)式,(4.2c)式和(4.3)式得

再根據(jù)(4.6)式,有

(4.12)

由矩陣Frobenious范數(shù)的定義知

即

將上式代入(4.12)式,再結(jié)合(4.11)式,(4.9)式和(4.8)式,得

(4.13)

根據(jù)(4.10)式,(4.9)式和定理4.1的條件,得

(4.14)

5　數(shù)值試驗(yàn)

采用Matlab(R2011b)數(shù)學(xué)軟件編程對(duì)算法A進(jìn)行數(shù)值試驗(yàn),程序運(yùn)行環(huán)境為Windows7(64bite),Intel(R)Core(TM)i3-2330M CPU@2.20GHZ,RAM:4G.

考慮如下兩類測(cè)試問題:(1)隨機(jī)二次半定規(guī)劃問題;矩陣Ai,Hj,C都是隨機(jī)生成的對(duì)陣矩陣,初始點(diǎn)(X,y,Z)的選取需在FP,FD內(nèi).記該類測(cè)試問題為RP.(2)最近相關(guān)矩陣Nearest correlation matrix(簡(jiǎn)記為NCM)問題[16-17].

s.t.diag(X)=e,

其中G是對(duì)稱正定陣,對(duì)角線上的元素為1,對(duì)角線外的元素在-1到1之間.e是所有元素都為1的n維向量.

NCM問題具有廣泛的應(yīng)用,例如在市場(chǎng)營(yíng)銷以及經(jīng)濟(jì)學(xué)方面.但是遺憾的是,由于缺乏數(shù)據(jù)等信息，不能得到一個(gè)完整的矩陣G,即矩陣G的一些元素是未知的.于是通過求解NCM問題得到一個(gè)有效的、且與G最近的相關(guān)矩陣.在數(shù)值測(cè)試中,取G是隨機(jī)矩陣且G∶=(1-α)B+α E[18],其中α∈(0,1),E是元素屬于[-1,1]的隨機(jī)對(duì)稱陣,B是具有指定特征值的測(cè)試矩陣,在NCM問題的數(shù)值測(cè)試中取n=10,20,30.

利用對(duì)稱Kronecker積的矩陣形式

將矩陣的對(duì)稱Kronecker積轉(zhuǎn)化為矩陣的Kronecker積計(jì)算,其中Q的取法詳見文獻(xiàn)[1]附錄.

在數(shù)值測(cè)試中,取算法A中的參數(shù)δ=0.3,γ=0.3,終止準(zhǔn)則是Zk· Xk≤10-6.

測(cè)試的數(shù)值結(jié)果見表1(所有結(jié)果都是每個(gè)維數(shù)測(cè)試10次,然后取平均值得到的),其中表中符號(hào)含義如下:prob,測(cè)試問題的編號(hào);n,矩陣X的維數(shù);m,不等式約束的個(gè)數(shù);Itr,算法迭代次數(shù);Nf,目標(biāo)函數(shù)值的計(jì)算次數(shù);time(s),算法迭代所需CPU時(shí)間(s).

表1數(shù)值結(jié)果

Table 1Numerical results

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

上述數(shù)值結(jié)果表明本文提出的基于H..K..M方向的路徑跟蹤算法是有效的，能在合理的時(shí)間內(nèi)求解中等規(guī)模的NCM問題.為進(jìn)一步提高本文算法的數(shù)值效果，我們將在后續(xù)的研究中深入探討步長(zhǎng)α的更有效的確定方法，即線搜索技術(shù).

參考文獻(xiàn):

[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]關(guān)秀翠,刁在筠.二次半定規(guī)劃問題及其投影收縮算法[J].全國(guó)高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào),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]徐鳳敏,徐成賢.求解二次半定規(guī)劃的原對(duì)偶內(nèi)點(diǎn)算法[J].工程數(shù)學(xué)學(xué)報(bào),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.

(責(zé)任編輯：米慧芝)

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

作者簡(jiǎn)介：黎健玲(1965-),女,博士,教授,主要從事最優(yōu)化理論與方法研究,E-mail:jianlingli@126.com。

中圖分類號(hào):C934

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1005-9164(2016)05-0396-08

修回日期：2016-09-25

*國(guó)家自然科學(xué)基金項(xiàng)目(No.11561005)和廣西自然科學(xué)基金項(xiàng)目(2016GXNSFAA380248,2014GXSFFA118001)資助。

**通信作者。

廣西科學(xué)Guangxi Sciences 2016,23(5):396～403

網(wǎng)絡(luò)優(yōu)先數(shù)字出版時(shí)間：2016-11-21【DOI】10.13656/j.cnki.gxkx.20161121.007

網(wǎng)絡(luò)優(yōu)先數(shù)字出版地址：http://www.cnki.net/kcms/detail/45.1206.G3.20161121.1520.014.html