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

遲曉妮,張睿婕,劉三陽

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

1.引言

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

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

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,則對任意參數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)的最優解.

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

定義

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

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

將P?(κ)線性互補問題的全牛頓步可行內點算法中的函數[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,有

選取適當參數θ及任意初始點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 選擇參數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)為單調遞增函數.故由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給定常數β ∈(0,1),令若δ(υ)≤βt,則δ(υ+)<βt+.

證由引理3.5得

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

不妨設

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

即

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

4.復雜度分析

選取參數β ∈(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設參數θ,K按(4.1)式選取.若δ(υ)≤βt,則

證由引理3.3得

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

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

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

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

5.數值算例

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

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

表5.1 K ≤2時求解不同規模和β值的LWCP的數值結果

表5.2 K ≤0.5時求解不同規模和β值的LWCP的數值結果

例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,得到最優解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的值