線性權互補問題的新全牛頓步可行內點算法

遲曉妮,張睿婕,劉三陽

(1.桂林電子科技大學數(shù)學與計算科學學院,廣西 桂林541004;2.桂林電子科技大學廣西密碼學與信息安全重點實驗室,廣西 桂林541004;3.桂林電子科技大學廣西自動檢測技術與儀器重點實驗室,廣西 桂林541004;4.西安電子科技大學數(shù)學與統(tǒng)計學院,陜西 西安710071)

1.引言

作為互補問題[1]的推廣,權互補問題[2]在科學和工程等領域有著廣泛的應用.Potra[2]證明了Fisher市場均衡問題可以建模為非線性互補問題[3],也可以建模為權互補問題,但后者在某些情況下能得到更高效的數(shù)值解方法[4].全牛頓步內點算法[5]是求解線性規(guī)劃的熱門算法.2003年,Darvay[6]基于連續(xù)可微函數(shù)提出求解線性規(guī)劃的新全內點算法,并給出了該算法的迭代復雜度.隨后,Darvay[7]等人又提出求解線性規(guī)劃的新全牛頓步內點算法,并分析了算法的多項式迭代復雜度.由于權向量非零,權互補問題的理論[8]和算法研究相比于互補問題更為復雜,因此關于權互補問題的算法尚不多見.Potra[9]設計了一種預估-校正內點算法[10]用于求解充分權互補問題.

本文將文[11]中P?(κ)線性互補問題全牛頓步可行內點算法的連續(xù)可微函數(shù)推廣到權互補問題,對中心路徑進行等價變換得到新搜索方向,給出求解Rn上線性權互補問題的可行內點算法.該算法采用全牛頓步,避免進行線搜索求解步長,節(jié)省了計算時間和內存.分析了算法的可行性,并證明其迭代復雜度與目前線性優(yōu)化最好的多項式時間復雜度相同.數(shù)值算例結果表明算法的有效性.

2.算法

考慮Rn上的線性權互補問題(LWCP)[1]：尋找向量對(x,y,s)∈Rn×Rm×Rn使得

其中A ∈Rm×n,b ∈Rm,c ∈Rn,ω ∈Rn++.定義LWCP(2.1)的嚴格可行域為}

令初始點(x0,y0,s0)∈F0,ω(t)=tx0s0+(1?t)ω,其中t ∈[0,1] 且t0=1.考慮(2.1)的擾動問題

這里e=(1,··· ,1)T.假定A為行滿秩矩陣,即R(A)=m.若內點條件(IPC)[5]成立,即(x,y,s)∈F0,則對任意參數(shù)t ∈(0,t0],方程組(2,2)有唯一解(x(t),y(t),s(t)).集合{(x(t),y(t),s(t))|t>0}稱為LWCP(2.1)的中心路徑.當t →0時ω(t)→ω,得LWCP(2.1)的最優(yōu)解.

考慮線性優(yōu)化問題內點算法中的連續(xù)可微函數(shù)[6]φ:R+→R+,并假設其反函數(shù)φ?1存在.用等價變換替換LWCP的擾動問題(2.2)中第三式,則新牛頓搜索方向(?x,?y,?s)應滿足方程組

定義

由(2.4)式可知,方程組(2.3)可化為

其中:=AV ?1X,V:=diag(υ),X:=diag(x),W(t)=diag(ω(t)).

將P?(κ)線性互補問題的全牛頓步可行內點算法中的函數(shù)[11]推廣到LWCP(2.1),則方程組(2.5)可化為

定義鄰近測度

令

引理2.1[5]若u與v正交,則

由方程組(2.6)知dTx ds=0.因此由(2.8)式得

由引理2.1和(2.7)式可知

引理2.2對任意υ ∈Rn,有

選取適當參數(shù)θ及任意初始點x0>0,s0>0,y0= 0,令ω(t0)=x0s0,其中t0= 1.顯然δ(x0,s0,ω(t0))=0.求解方程組(2.6)并結合(2.4)式得牛頓搜索方向(?x,?s,?y),其中t+=(1?θ)t,θ ∈(0,1).令新迭代點

滿足內點條件并且δ(x+,s+;ω(t+))<βt+.當∥xs ?ω∥≤ε時,算法終止.下面給出本文算法的具體步驟.

算法2.1求解LWCP的新全牛頓步可行內點算法

步1 選擇參數(shù)t0= 1,ε >0,β ∈(0,1),初始點(x0,y0,s0)∈F0,y0= 0且ω(t0)=x0s0.置k:=0.

步2 當∥xs ?ω∥≤ε,算法終止;否則轉步3.

步3 求解方程組(2.6)并結合(2.4)式得搜索方向(?x,?s,?y),令

步4 由下面的(4.1)式求得θ ∈(0,1),令

置k:=k+1.轉步2.

3.算法分析

引理3.1若δ(υ)<1,則(x+,y+,s+)∈F0.

證令α ∈[0,1],記

由(2.4),(2.8)式和方程組(2.6)得

因為ω(t)>0,(1?α)υ2≥0,所以由(3.1)式知,若

則x(α)s(α)>0.由(2.7)和(2.9)式,及δ(υ)<1可知

又由(2.7)式和引理2.2得

當δ(υ)<1時,由(3.3)式得

因此若δ(υ)<1時,則由(3.1),(3.2)和(3.4)式可得x(α)s(α)>0.

由于x(0)=x>0,s(0)=s>0且x(α),s(α)與α呈線性關系,所以對α ∈[0,1]有x(α),s(α)>0,相應地,x(1),s(1)>0.證畢.

引理3.2令,則

證由(3.1)式和引理2.2得

令f(η)=η2+η ?η3.由f′(η)=2η+1?3η2=?(3η+1)(η?1)知,f(η)在(0,1)為單調遞增函數(shù).故由f(η)≤f(1)=1,η ∈(0,1)可得

因此,由(2.10),(3.5)和(3.6)式,得

引理3.3令若δ(υ)<1,則

證由(2.6),(2.8)和(3.1)式得

因為δ(υ)<1,所以由引理2.2可得

由(2.7),(2.8),(2.9),(3.7)和(3.8)式知,當δ(υ)<1時有

引理3.4令,則

證因為ω(t+)=ω(t)+θt(ω ?x0s0),所以

故由引理3.3得

引理3.5令則

證由引理3.2和引理3.4得

引理3.6給定常數(shù)β ∈(0,1),令若δ(υ)≤βt,則δ(υ+)<βt+.

證由引理3.5得

因為(3.11)式右端函數(shù)關于δ(υ)單調增加,所以若δ(υ)≤βt,則

不妨設

其中t+=(1?θ)t,θ ∈(0,1)化簡(3.13)式得

即

因為β ∈(0,1),所以故由式(3.14)得

4.復雜度分析

選取參數(shù)β ∈(0,1),t ∈[0,1],令

因為t0= 1,則ω(t0)=x0s0,故δ(υ0)= 0<βt0.又t+= (1?θ)t,由引理3.1和引理3.6可知下一迭代點(x+,y+,s+)嚴格可行且δ(υ+)<βt+.

引理4.1設參數(shù)θ,K按(4.1)式選取.若δ(υ)≤βt,則

證由引理3.3得

因為β ∈(0,1),t ∈[0,1]所以若δ(υ)≤βt,則

引理4.2設參數(shù)θ,K按(4.1)式選取.若LWCP(2.1)存在最優(yōu)解,則算法2.1至多經過

次迭代得到LWCP(2.1)的ε-近似解.

證因為∥ω(t)∥∞≤max{max(x0s0),max(ω)},所以由引理4.1可得

5.數(shù)值算例

為檢驗算法2.1的有效性,在Intel(R)Core i5 CPU2.3GHz 8.0內存,IOS操作系統(tǒng)的計算機上運用MATLABR 2016b編程進行數(shù)值實驗.

在算法2.1中令參數(shù)t0= 1,ε= 10?5.隨機生成5個不同規(guī)模的LWCP,且每種規(guī)模產生10個問題,分別取β= 0.7,β= 0.8和β= 0.9進行求解.隨機生成行滿秩矩陣A ∈Rm×n.選取任意起始點(x0,y0,s0)∈F0,權向量ω >0,按照(4.1)式選擇參數(shù)θ,K.算法的終止準則為∥xs ?ω∥≤ε,記Gap=∥xs ?ω∥.如表5.1,5.2所示,在算法2.1求解相同規(guī)模的LWCP時K值上界與迭代次數(shù)和運行時間呈正相關關系.在K值上界確定的情況下,m和n對運行時間和迭代次數(shù)也有較大影響;從數(shù)值試驗結果可知,在t從1減少到0的過程中,θ單調遞增,且β的取值越趨向1,算法2.1求解同一LWCP所需的迭代次數(shù)就越少.表中數(shù)據(jù)均為求解不同規(guī)模的LWCP10次結果的平均值.

表5.1 K ≤2時求解不同規(guī)模和β值的LWCP的數(shù)值結果

表5.2 K ≤0.5時求解不同規(guī)模和β值的LWCP的數(shù)值結果

例5.1考慮R6上的LWCP(2.1),其中

取初始點x0= (2,1,11,5,7,3)T,y0= (0,0,0,0)T,s0= (7,8,3,5,6,4)T,t0= 1,β=0.9.用算法2.1求解例5.1,得到最優(yōu)解x?= (2.100,0.738,10.986,4.836,7.105,3.092)T,y?=?(0.141,0.109,0.035,0.026)T,s?=(7.618,9.486,2.731,3.722,6.334,4.851)T.圖5.1 可知,隨著t的減小,Gap逐漸減小并趨于0,且每步迭代鄰近測度都小于βt.

圖5.1 迭代過程中鄰近測度及Gap的值