用離散化方法證明半定規劃的拉格朗日強對偶定理

羅 丹,羅洪林

(重慶師范大學 數學科學學院， 重慶 401331)

1 背景

(P) mincTx

s.t.F(x)≥0

半定規劃也稱為帶有半正定錐約束的線性規劃,其退化情形包括線性規劃和凸二次規劃。半定規劃廣泛地存在于系統與控制理論[1]、金融工程[2]、量子化學[3]、信號處理[4]等諸多領域。 半定規劃的多項式內點算法為求解組合優化領域的某些中小規模的NP難問題(如著名的旅行商問題[5]和最大割問題[6])提供了有效的解決途徑。 自20世紀90年代初期至今,半定規劃一直是優化領域的一個熱門的研究課題。

半定規劃的對偶理論在半定規劃的理論研究和算法設計中都扮演著十分重要的角色。半定規劃的對偶理論研究除了常見的拉格朗日對偶和共軛對偶理論外,具有代表性的包括Ramana等提出的廣義拉格朗日slater對偶理論[7]和極小錐對偶理論[8]。 本文主要從算法的角度重新考察半定規劃的拉格朗日對偶理論。

原始半定規劃問題(P)的拉格朗日對偶模型為:

(D) max -TrF0Z

s.t.TrFiZ=ci,i=1,…,m

Z≥0

記val(D)為(D)的最優目標函數值， 稱(D)可行是指其可行域非空,即：存在Z≥0使得TrFiZ=ci,i=1,…,m; 稱(D)嚴格可行是指存在Z>0使得TrFiZ=ci,i=1,…,m成立。

內點算法是求解中小規模的半定規劃的十分有效的求解算法。內點算法利用半定規劃的原始對偶模型的1階最優性條件(KKT條件)建立可行或不可行中心路徑的方式提出了各種可行內點算法和不可行內點算法。 最近,楊洋等[9]利用不可行中心路徑給出了一種寬領域不可行內點算法。 內點算法由于需要儲存和分解牛頓矩陣,這將占用大量的內存空間且十分耗時，在2002年,Helmberg[10]指出:在一般的計算機處理系統上通過合理的等待時間,用內點算法能有效計算的半定規劃問題,其矩陣變量n≤200,約束的個數m≤3 000。 截至2013年,Huang等[11]指出,目前的內點算法一般可以有效求解的中小規模的半定規劃的矩陣變量n≤1 000,約束的個數m≤10 000。

如何為大規模的半定規劃問題設計有效的求解算法是一個有待解決的問題。 為了擺脫牛頓矩陣的儲存和分解,受Yang[12]的啟發,本文將從強對偶定理的一種新的證明方法來闡述半定規劃的一種新的求解算法。

本文的第1部分利用離散化方法將半定規劃問題近似地轉換為一個線性規劃問題并給出了2個問題的最優解的近似程度的刻畫。本文的第2部分給出了強對偶定理的一個新的證明方法以及半定規劃的一種新的離散化求解算法及其收斂性證明。 針對本文提出的半定規劃的一種新的理論算法,在第3部分提出了該算法數值實現時可能面臨的幾個待解決的問題。 附錄部分呈現了半定規劃的強對偶定理的經典證明過程以方便讀者比較它與離散化證明方法之間的異同。

2 半定規劃的離散化方法

本節主要介紹對偶半定規劃(D)的離散化方法以及(D)的最優解與離散化問題的最優解之間的誤差刻畫。為了半定規劃的離散化方法的描述完整性,首先引入Yang[12]關于原始半定規劃問題(P)的離散化過程:

第1步利用{x∈Rm|F(x)≥0}與{x∈Rm|yTF(x)y≥0,?y∈Rn}的等價性,將(P)等價地轉化為如下形式的優化問題:

mincTx

s.t.yTF(x)y≥0

y∈Rn

(1)

第2步將集合{y∈Rn|yTF(x)y≥0}正則化為{y∈Rn|yTF(x)y≥0,yTy=1},則問題(1)可等價地表示為如下的線性半無限規劃問題:

mincTx

s.t.aT(y)x+b(y)≥0

y∈Y

(2)

其中Y={y∈Rn|yTy=1},aT(y)=(yTF1y,yTF2y,…,yTFmy),bT(y)=yTF0y

記

問題(2)關于κ∈R的下水平集定義為:

第3步選取有限網格Yd?Y離散化問題(2),得到一個線性規劃問題:

mincTx

s.t.aT(y)x+b(y)≥0

y∈Yd

(3)

網格Yd與集合Y的Hausdorff距離定義為

記

通過上述離散化方法可將(P)近似地轉換為一列線性規劃問題進行近似求解,求解的精度取決于網格Yd的取法,具體參見文獻[12]中的引理2。

眾所周知,不論是半定規劃還是線性規劃,研究其對偶理論的主要目的之一就是化解問題求解的難度,即:相較于原問題,當對偶問題的求解更加容易時,在強對偶定理的理論支撐下,可以通過求解其對偶問題而獲得原問題的解。基于此,利用Dattorro[13]中的:

(4)

給出對偶半定規劃問題(D)的離散化過程:

第1步利用等式(4)將對偶半定規劃問題(D)等價地轉化成如下形式的優化問題:

(5)

第2步通過bj≥0,j=1,2,…的適當選取(仍然記為bj≥0,j=1,2,…),將向量zj∈Rn,j=1,2,…單位化,則式(5)可等價地轉換為:

(6)

記

i=1,…,m,zj∈Z′}

問題(6)關于λ∈R的下水平集定義為:

第3步選取有限網格Zk?Z′,不妨設Zk含有k個向量z1,z2,…,zk,則與之對應的k個非負實數為b1,b2,…,bk,現通過該網格離散化問題(6)得到一個線性規劃問題:

(7)

網格Zk與集合Z的Hausdorff距離定義為:

記

i=1,…,m,z∈Zk}

為了利用離散化方法證明半定規劃的強對偶定理,需要建立以下幾個結論：

定理1 如果(P)嚴格可行,(D)可行,那么對任意給定的λ∈R,水平集

有界。

證明假設存在一個λ0∈R使得水平集

無界,則存在非負數列{bj}?L≥(BD(Z),F0,λ0)滿足:

(8)

(9)

且

zj∈Z′={z∈Rn|||z||=1}=Y

對任意給定的y∈Y,則存在序列{vj}滿足

(10)

(11)

由級數收斂的必要條件和式(9)可得:

再結合式(10)以及式(11)可推得:

yTF0y=0,?y∈Y

(12)

(13)

(14)

(15)

由級數收斂的必要條件和(ref{eq7})可得:

再結合式(14)與(15)可推得:

yTFiy=0,?y∈Y,i=1,2,…,m

(16)

故,由式(12)和(16)可得:

?y∈Y

由于(P)與式(2)等價,不難發現上式與半定規劃問題(P)的嚴格可行性假設矛盾。

由于(D)的嚴格可行性必蘊含其可行性,則由定理1可得如下推論：

推論1 問題(P)和(D)都嚴格可行,則對λ∈R,水平集

有界。

命題1 若問題(P)和(D)都嚴格可行,則對κ∈R,水平集

有界。

證明由Yang[12]中的引理1和半定規劃(P)嚴格可行性蘊含問題的可行性得結論成立。

如下2個命題揭示了:用離散化方法可以求得半定規劃問題(P)和(D)滿足任意精度的近似解。

證明因為(P)與式(2)等價,(D)與式(6)等價,由定理1及文獻[14]中的定理4.4可以證得。

類似地,可以利用推論2和命題1以及文獻[14]中的定理4.4可以證明下面這個命題：

命題3 若問題(P)和(D)都嚴格可行,那么,

3 強對偶定理的離散化證明方法及其在算法設計中的應用

當半定規劃問題(P)和(D)滿足一定的條件時,若能推出對偶間隙:

val(P)-val(D)=0

則稱為半定規劃的強對偶定理。 關于半定規劃的強對偶定理的經典證明方法請參見附錄。2005年,Yang[12]通過原始半定規劃(P)的離散化方法證明了在原始半定規劃(P)可行和對偶半定規劃(D)嚴格可行的假設條件下的強對偶定理。 值得指出的是,可以利用Yang[12]關于原始半定規劃(P)的離散化方法為(P)的求解設計一種新的求解算法:利用文獻[14-15]中關于半無限規劃的離散化算法近似地求解(P)。

下面利用對偶半定規劃問題(D)的離散化方法重新給出半定規劃的另外2個強對偶定理的證明,進而為對偶半定規劃問題(D)設計一種新的求解算法。

3.1 強對偶定理的證明

定理2 如果(P)嚴格可行,(D)可行,則val(P)=val(D),且(D)的最優解Z*最優可達。

證明定理1可知,對任意給定的λ∈R,水平集L≥(BD(Z),F0,λ)是有界閉的,則L≥(BD(Z),F0,λ)是緊集且當λ充分大時非空。 因此,(D)的最優解可以在可行域中取到。

現定義網格Zk={z1,z2,…,zk},則式(7)可以等價地表示成如下形式的線性規劃問題:

(17)

maxcTx

s.t.aT(yj)x+b(yj)≥0

yj∈Yk,j=1,…,k

(18)

由半定規劃的弱對偶定理有:

val(P)≥val(D)

(19)

(20)

因此,對于充分小的dZk,val(Zk)有界。故由線性規劃的強對偶定理,有

val(Zk)=val(xk)

(21)

由于(P)的可行解包含式(18)的可行解,則有

val(P)≤val(xk)

(22)

結合式(20)、(21)和(22),令dYk→0可得:

val(P)≤val(D)

(23)

綜合式(19)和(23)得: val(P)=val(D)

定理3 如果(P)和(D)都嚴格可行,則 val(P)=val(D),且(P)的最優解x*和(D)的最優解Z*都最優可達。

證明由命題1和推論1可知,對任意給定的λ∈R,水平集L≥(XP(Y),c,κ)和L≥(BD(Z),F0,λ)是有界閉的,則L≥(XP(Y),c,κ),L≥(BD(Z),F0,λ)是緊集且分別當κ、λ充分大時非空。 因此,(P)和(D)的最優解可以在其各自的可行域中取到。

現定義網格Zk={z1,z2,…,zk},則式(7)可以等價地表示成如下形式的線性規劃問題:

(24)

maxcTx

s.t.aT(yj)x+b(yj)≥0

yj∈Yk,j=1,…,k

(25)

(26)

(27)

因此,對于充分小的dYk, val(xk)有界；對于充分小的dZk, val(Zk)有界。 故由線性規劃的強對偶定理,有

val(Zk)=val(xk)

(28)

結合式(26)、(17)和(28),令dYk→0,dZk→0可得:

val(P)=val(D)

(29)

3.2 半定規劃的離散化算法

求解一般的中小規模的半定規劃問題時可以用CVX[16]中的SDPT3或SeDuMi求得滿足任意精度的近似解。 不論是SDPT3還是SeDuMi,其求解原理都是利用原始對偶內點算法,需要儲存和分解牛頓矩陣,這將占用大量的內存空間且十分耗時。

根據文獻[14]對半無限規劃問題的離散化算法的描述,當(P)嚴格可行,(D)可行,針對對偶半定規劃問題(D),設計如下算法：

算法1 首先將半定規劃對偶問題(D)轉換為等價的線性半無限規劃問題(16)，然后,給定一個步長向量h∈Rm,其中hj>0,j=1,…,m,并且固定一個z0∈Rm,定義網格

Gh={z|(z-z0)j=αjhj,αj∈N,j=1,…,m}

并且

Zh=Z∩Gh

具體步驟如下:

① 設hi+1=(1/ni)hi,ni∈N,ni≥2。

④ 如果i>i0(預先規定步數)則停止,否則繼續第(i+1)步。

注：1) 該離散化方法在求解對偶半定規劃問題(D)時,無需儲存和分解牛頓矩陣;

下面簡單地給出該算法的收斂性證明。

4 結論和展望

利用離散化方法,本文給出了強對偶定理的一種新的證明方法,并利用該方法設計了一種半定規劃的求解算法。 本文僅從理論上給出了半定規劃的求解算法及其收斂性分析,值得指出的是，雖然該算法擺脫了傳統的內點算法對于牛頓矩陣的儲存和分解,但是在算法步驟中涉及隨機變量,這可能會導致數值實驗結果的不穩定性。若讓算法中的隨機變量滿足一定的概率分布,對于滿足均勻分布或標準正態分布等是否可以使得數值實驗結果更加穩定， 如何改進該算法以適用于有效地求解大規模的半定規劃問題， 這些是接下來準備研究的問題。

參考文獻：

[1] YAO D,ZHANG S Z,ZHOU X Y.Control theory stochastic linear-quadratic control via semidefinite programming [J].Society for Industrial and Applied Mathematics,2001,40:1-23.

[2] DALAKOURAS G V,KWON R H,PARDALOS P M.Semidefinite programming approaches for bounding asian option prices [J].Computational Methods in Financial Engineering,2008,103-116.

[3] DING Y,GE D,WOLKOWICZ H.On equivalence of Semidefinite relaxations for quadratic matrix programming [J].Mathematics of Operations Research,2011,36:8-104.

[4] BISWAS P,TOH K C,YE Y.A distributed Semidefinite approach for large-scale noisy anchor-free graph realization with applications to molecular conformation [J].Society for Industrial and Applied Mathematics,2008:1251-1277.

[5] KLERK E D,PASECHNIK D V,SOTIROV R.On semidefinite programming relaxations of the traveling salesman problem [J].Society for Industrial and Applied Mathematics,2008,19:1559-1573.

[6] GRIPPO L,PALAGI L,PICCIALLI V.An unconstrained minimization method for solving low-rank SDP relaxations of the maxcut problem [J].Mathematical Programming,2011,126:119-146.

[7] RAMANA M V.An exact duality theory for semidefinite programming and its complexity implications [J].Mathematical Programming,1997,77:129-162.

[8] Ramana M V,Tuncxelz L,Wolkowicz H.Strong duality for semidefinite programming [J].Society for Industrial and Applied Mathematics,1997,7(3):641-662.

[9] 楊洋,羅洪林,羅慧林.關于半定規劃的一種寬鄰域不可行內點算法的注記 [J].運籌學學報,2016,20(2):79-87.

[10] HELMBERG C.Invited Review,Semidefinite programming [J].European Journal of Operational Research,2002,137:461-482.

[11] HUANG A Q,XU C X.A trust region method for solving semidefinite programs [J].Computational Optimization and Applications,2013,55:49-71.

[12] YANG Q Z.A new proof of the strong duality theorem for semidefinite programming [J].it Mathematical analysis and applications,2005,303:622-626.

[13] DATTORRO J.Convex optimization and euclidean distance geometry[Z].Meboo Publishing USA,2005.

[14] HETTICH R,KORTANEK K O.Semi-infinite programming:Theory,methods,and application [J].Society for Industrial and Applied Mathematics,1993,35:380-429.

[15] HETTICH R.An implementation of a discretization method for semi-infinite programming [J].Math Programming,1986,34:354-361.

[16] GRANT M,BOYD S.The CVX Users’ Guide (Release 2.0 (Beta))[EB/OL].[2018-01-01].http://cvxr.com/cvx/(221) March 2013.

[17] ALIZADEH F.Interior point methods in semidefinite programming with applications to combinatorial optimization [J].Society for Industrial and Applied Mathematics,1995,5:13-51.

附錄A

定理4[17]令

z1=inf{C·X:Ai·X=bi,X≥0}

z2=sup{bTy:C-∑yiAi≥0}

如果存在一個m維向量y,使得ATy>0,那么z2=z1。

Mat(ATy1)≤0且bTy1>0

意味著對偶問題是無界的。由于任意的對偶可行解y能增加一個任意大的正乘子在y1前,得到另一個可行解有巨大的目標函數值。 因此,z2=z1=+∞。 可以假設z1和z2都是有限的， 假設z2

C·X=z2

AvecX=b

X≥0

是不可行的。 因此,通過推廣的Farkas引理2.3,存在一個標量y0和m維向量y,使得

且z2y0+bTy<0

(19)

其中vecAi是A的第i行。 對y0,下列之一成立。

1) 若y0=0,式(19)等價于

Mat(ATy)≤0 且bTy>0

通過推廣的Farkas引理得AvecX=b且X≥0是不可行的， 因此z1=+∞。

2) 若y0=0,則式(19)除以y0得

C-Mat(AT(-y/y0))≥0且

z2-bT(y/y0)<0

意味著z2不是對偶問題的最優值。

3) 若y0<0,則式(19)除以-y0得

-C+Mat(AT(-y/y0))≥0且

-z2+bT(-y/y0)<0

事實上,由于是嚴格不等式,則存在ε>0,使得

-C+Mat(AT(-y/y0))≥0且

-z2+bT(-y/y0)<-ε

同樣地,通過z2的最優性,存在一個y*,使得

C-Mat(ATy*)≥0且z2-bTy*<ε

最后2個式子相加得

Mat(AT(-y/y0-y*))≥0且

bT(-y/y0-y*)<0

再次通過推廣的Farkas引理得原問題不可行且z1=∞，與假設z2