線性二階錐權互補問題的非精確非單調光滑化牛頓法

遲曉妮,劉文麗,劉三陽,趙 敏

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

0 引 言

作為互補問題(CP)的推廣,權互補問題(WCP)在大氣化學、多體動力學、經濟學等領域的一些均衡問題中應用廣泛[1]. 本文考慮線性二階錐權互補問題(LWSOCCP),即給定M∈n×n,q∈n×n,權向量w∈K,找到一向量對(x,s)∈n×n,使得

x∈K,s=Mx+q∈K,x°s=w,

(1)

這里K為n維二階錐,即

K∶={x=(x1,x2)∈×n-1: ‖x2‖≤x1},

其中‖·‖表示Euclid范數. 當w=0時,LWSOCCP(1)退化為線性二階錐互補問題(LSOCCP):

x∈K,s=Mx+q∈K,x°s=0.

由于WCP中存在非零權向量,其理論和算法比CP更復雜,因此目前關于WCP的研究文獻報道較少,且大多僅限于n上. Potra[1]最早提出了WCP的概念,并研究了求解WCP的兩種內點算法；Tang[2]給出了求解線性WCP的非單調光滑化牛頓法；Chi等[3]研究了Euclid Jordan代數上水平線性WCP解的存在性與唯一性.

在優化問題求解[4]中,光滑化牛頓法是一種有效的算法[5-6],其主要思想為先將原問題轉化為與之等價的方程組,然后用光滑化算法對該方程組進行求解,從而得到原問題的解. 與內點算法[7-8]不同,光滑化牛頓法不要求初始點和中間迭代點嚴格可行[9]. 此外,光滑化牛頓算法每次迭代都需精確求解方程組,但當初始點離解較遠時,精確線搜索通常效率不高,因此本文采用非精確線搜索,從而減少計算工作量.

本文結合Fischer-Burmeister函數和Natural-Residual函數[2],給出二階錐上權互補問題的一個新的含參數光滑函數,并基于該函數,提出新的求解LWSOCCP的非精確非單調光滑化牛頓法. 算法采用非單調線搜索技術,提高了找到全局最優解的可能性及收斂速度,數值算例結果表明算法有效.

1 預備知識

對任意x=(x1,x2)∈×n-1,s=(s1,s2)∈×n-1,定義K上的Jordan積[10]為

x°s=(xTs,x1s2+s1x2),

其單位元e∶=(1,0,…,0)T∈n,且x2∶=x°x. 對任意x=(x1,x2)∈×n-1,定義箭形矩陣[10]為

其中I為(n-1)階單位矩陣. 易證x°s=Lxs,Lx+s=Lx+Ls.Lx是正定(半正定)矩陣當且僅當x∈intK(x∈K).

引理1[10]對任意x,s∈n,w?K0,有

(2)

w2?Kx2?w?Kx.

(3)

若將式(2)和式(3)中“?”替換為“±”,結論也成立.

引理2[11]對任意a,b,u,v∈n,若a?K0,b?K0,a°b?K0滿足〈u,v〉≥0和a°u+b°v=0,則u=v=0.

2 光滑函數及其性質

對于LWSOCCP(1),構造一個新的含參數光滑函數φτ(μ,x,s):+×n×n→n:

其中w∈K為給定權向量,參數τ∈[0,4). 當w=0時,φτ(μ,x,s)退化為LSOCCP光滑函數

對任意τ∈[0,4),有

(5)

下面基于含參數光滑函數(4),將LWSOCCP(1)轉化為含參數光滑方程組. 令z∶=(μ,x,s)∈+×n×n,定義函數

(6)

由式(1)和式(4)～(6)知,(x*,s*)是LWSOCCP(1)的解,且μ*=0 ?Gτ(z*)=0.

定理1令z∶=(μ,x,s)∈+×n×n,τ∈[0,4),則下列結論成立:

1)Gτ(z)全局Lipschitz連續,處處強半光滑,且在點z∶=(μ,x,s)∈++×n×n處連續可微,其Jacobi矩陣為

(7)

其中

(8)

(9)

2) 若M為半正定矩陣,則Gτ(z)在任意點z∶=(μ,x,s)∈++×n×n處可逆.

證明: 1) 因為

故由文獻[12]中定理3.2易證結論成立.

2) 令Δz∶=(Δμ,Δx,Δs)∈×n×n是Gτ(z)零空間中任一向量,則只需證Δz=0. 由式(7)知

Δμ=0,

MΔx-Δs=0,

(10)

(11)

將式(11)兩端左乘Lc,得

整理式(12)有

從而

對任意μ>0,由c的定義和w∈K,有

因此由引理1得

(15)

(16)

所以

(17)

又因為M為半正定矩陣,故由式(10)知〈Δx,Δs〉=〈Δx,MΔx〉≥0,從而有

〈Δx+μΔs,μΔx+Δs〉=μ(‖Δx‖2+‖Δs‖2)+(1+μ2)〈Δx,Δs〉≥0.

(18)

由式(14)～(18)和引理2有Δx+μΔs=0,μΔx+Δs=0,故由式(10)得Δx=0,Δs=0,結論成立. 證畢.

3 非精確非單調光滑化牛頓法

下面基于光滑函數(4),給出求解LWSOCCP(1)的非精確非單調光滑化牛頓法,且在適當的假設條件下證明算法適定.

定義函數Hτ:+×n×n→+為Hτ(z)∶=‖Gτ(z)‖2.

算法1LWSOCCP的非精確非單調光滑化牛頓法(INS).

選擇常數σ,δ∈(0,1),p∈[0,1],μ0>0. 取ρ∈(0,1)和ξ∶=‖Gτ(z0)‖+1,使得μ0ρξ<1. 令z0∶=(μ0,x0,s0)∈×n×n為任意初始點,

Q0=1,C0∶=Hτ(z0),T(z0)=ρmin{1,‖Gτ(z0)‖1+p}.

設序列{υk}滿足υk∈[0,1-μ0ρξ]. 選取ηmin和ηmax,使得0≤ηmin≤ηmax<1. 置k∶=0.

步驟1) 若Hτ(zk)=0,則算法停止;

步驟2) 解下列方程組得搜索方向Δzk∶=(Δμk,Δxk,Δsk)∈×n×n:

(19)

其中T(zk)=min{ρ,ρ‖Gτ(zk)‖1+p,T(zk-1)},‖r(zk)‖≤υk‖Gτ(zk)‖;

步驟3) 令αk為1,δ,δ2,…中使得下式成立的最大值:

Hτ(zk+αkΔzk)≤[1-σ(1-μ0ρξ-υk)αk]Ck;

(20)

步驟4) 令ηk∈[ηmin,ηmax],zk+1∶=zk+αkzk,且

Qk+1∶=ηkQk+1,

(21)

(22)

置k∶=k+1,轉步驟1).

定理2令{Ck},{T(zk)},{zk}是算法1生成的迭代序列,則:

1) {Ck}單調遞減;

2) 對任意k≥0,有Hτ(zk)≤Ck且{T(zk)}單調遞減;

3) 對任意k≥0,有μk>0且μ0T(zk)≤μk.

證明：類似文獻[13]中注3.1和文獻[14]中引理1可得結論.

定理3若M為半正定矩陣,則算法1適定.

下證算法1中步驟3)適定. 對任意α∈(0,1],定義

L(α)=Gτ(zk+αΔzk)-Gτ(zk)-αGτ(zk)Δzk.

(23)

由定理1中1)知Gτ(·)在任意點zk=(μk,xk,sk)∈++×n×n處連續可微,則

‖L(α)‖=ο(α).

(24)

對任意k≥0,由T(zk)的定義知,當‖Gτ(zk)‖≤1時,有

T(zk)≤ρ‖Gτ(zk)‖1+p≤ρ‖Gτ(zk)‖;

(25)

當‖Gτ(zk)‖>1時,有

T(zk)=ρ≤ρ‖Gτ(zk)‖.

(26)

因此由式(25)～(26)得

T(zk)≤ρ‖Gτ(zk)‖.

(27)

由定理2中2)和式(6)知,對任意k≥0,有

故eμk≤ξ. 從而由式(19),(23),(24),(27)知,對任意α∈(0,1],有

Hτ(zk+αΔzk)≤[1-σ(1-μ0ρξ-υk)α]Ck.

故算法1中步驟3)適定. 從而算法1適定.

4 收斂性分析

下面證明算法1的全局收斂性和局部超線性收斂性.

定理4設M為半正定矩陣,{zk}是算法1生成的迭代序列,則{zk=(μk,xk,sk)}的任一聚點z*∶=(μ*,x*,s*)都是Gτ(z)=0的解.

證明: 不失一般性,設當k→∞時,{zk=(μk,xk,sk)}收斂到z*∶=(μ*,x*,s*). 由定理2中1)知{Ck}為收斂序列,即

(28)

由定理2中2)知

0≤‖Gτ(zk)‖2=Hτ(zk)≤Ck≤Ck-1≤C0,

(29)

所以序列{Hτ(zk)}有界,因此必存在收斂子列{Hτ(zk)}k∈J,其中J?{0,1,2,…,k,…}. 故

反證法. 設Gτ(z*)>0,則C*>0,考慮下列兩種情形:

(i) 若對任意k∈J有αk≥α*>0,其中α*為某固定常數. 由式(20)～(22)得

由{Ck}有界知

(31)

又由式(21)得

(32)

(33)

從而由定理2中2)有

(34)

由定理2中3)知

的解. 由式(29),(33)得Hτ(z*)≤C*≤Hτ(z*),因此

Hτ(z*)=C*.

(35)

從而

2Gτ(z*)TGτ(z*)Δz*≥-σ(1-μ0ρξ-υ*)C*.

(36)

又由式(27)和式(35)知

T(z*)‖Gτ(z*)‖≤ρHτ(z*)=ρC*,

(37)

于是由式(36)～(37)有

即σ(1-μ0ρξ-υ*)≥2(1-μ0ρξ-υ*),與σ∈(0,1)矛盾. 故結論成立.

證明: 類似文獻[4]中定理8可證.

5 數值算例

下面用MATLAB R2014在Intel(R) core(TM)i5-5200 CPU 2.7 GHz,4.0 GB內存,Windows 7操作系統的計算機上編程,驗證算法1的性能,且每個算例均隨機生成15個問題進行求解. 參數取值為σ=0.04,μ0=0.01,δ=0.5,ρ=0.02,υk=2-k,τ=1. 令0=(0,0,…,0)T,e=(1,0,…,0)T,1=(1,…,1)T. 隨機生成矩陣N∈n×n和向量q∈n,令矩陣M=NTN,s=Mx+q,以‖G(z)‖≤10-6為算法終止條件,(x0,s0)為初始點,記AIT為平均迭代次數,ACPU為平均CPU時間.

例1給定權向量w=(1,0.1×1),令ηk=0.1,n=100. 表1列出了算法1求解不同初始點和不同p值LWSOCCP的數值結果. 由表1可見,不同初始點對算法所需的AIT和ACPU影響明顯,而p值對算法影響較小.

表1 算法1求解不同初始點和不同p值LWSOCCP的數值結果

例2給定權向量w=e,令p=1,初始點(x0,s0)=(e,e),問題規模從200維到600維. 用算法1求解不同ηk值的LWSOCCP,并比較算法1(INS)與采用精確線搜索的光滑化算法NS(nonmontone smoothing)[13]求解LWSOCCP的數值結果,結果列于表2. 由表2可見,算法1采用非單調線搜索(ηk≠0)比采用單調線搜索(ηk=0)所需的ACPU和AIT少,且在相同條件下,算法1比算法NS求解LWSOCCP的效果更好.

表2 算法1和算法NS求解LWSOCCP的數值結果

綜合表1和表2可見,算法1可有效求解LWSOCCP,且收斂速度快,迭代次數少,因而算法性能穩定.