線性二階錐權(quán)互補(bǔ)問(wèn)題的非精確非單調(diào)光滑化牛頓法

遲曉妮,劉文麗,劉三陽(yáng),趙 敏

(1. 桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院,廣西密碼學(xué)與信息安全重點(diǎn)實(shí)驗(yàn)室,廣西 桂林 541004; 2. 桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院,廣西自動(dòng)檢測(cè)技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室,廣西 桂林 541004; 3. 西安電子科技大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,西安 710071)

0 引 言

作為互補(bǔ)問(wèn)題(CP)的推廣,權(quán)互補(bǔ)問(wèn)題(WCP)在大氣化學(xué)、多體動(dòng)力學(xué)、經(jīng)濟(jì)學(xué)等領(lǐng)域的一些均衡問(wèn)題中應(yīng)用廣泛[1]. 本文考慮線性二階錐權(quán)互補(bǔ)問(wèn)題(LWSOCCP),即給定M∈n×n,q∈n×n,權(quán)向量w∈K,找到一向量對(duì)(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范數(shù). 當(dāng)w=0時(shí),LWSOCCP(1)退化為線性二階錐互補(bǔ)問(wèn)題(LSOCCP):

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

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

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

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

1 預(yù)備知識(shí)

對(duì)任意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. 對(duì)任意x=(x1,x2)∈×n-1,定義箭形矩陣[10]為

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

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

(2)

w2?Kx2?w?Kx.

(3)

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

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

2 光滑函數(shù)及其性質(zhì)

對(duì)于LWSOCCP(1),構(gòu)造一個(gè)新的含參數(shù)光滑函數(shù)φτ(μ,x,s):+×n×n→n:

其中w∈K為給定權(quán)向量,參數(shù)τ∈[0,4). 當(dāng)w=0時(shí),φτ(μ,x,s)退化為L(zhǎng)SOCCP光滑函數(shù)

對(duì)任意τ∈[0,4),有

(5)

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

(6)

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

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

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

(7)

其中

(8)

(9)

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

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

故由文獻(xiàn)[12]中定理3.2易證結(jié)論成立.

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

Δμ=0,

MΔx-Δs=0,

(10)

(11)

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

整理式(12)有

從而

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

因此由引理1得

(15)

(16)

所以

(17)

又因?yàn)镸為半正定矩陣,故由式(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,結(jié)論成立. 證畢.

3 非精確非單調(diào)光滑化牛頓法

下面基于光滑函數(shù)(4),給出求解LWSOCCP(1)的非精確非單調(diào)光滑化牛頓法,且在適當(dāng)?shù)募僭O(shè)條件下證明算法適定.

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

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

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

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

設(shè)序列{υk}滿(mǎn)足υ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,轉(zhuǎn)步驟1).

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

1) {Ck}單調(diào)遞減;

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

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

證明：類(lèi)似文獻(xiàn)[13]中注3.1和文獻(xiàn)[14]中引理1可得結(jié)論.

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

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

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

(23)

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

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

(24)

對(duì)任意k≥0,由T(zk)的定義知,當(dāng)‖Gτ(zk)‖≤1時(shí),有

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

(25)

當(dāng)‖Gτ(zk)‖>1時(shí),有

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

(26)

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

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

(27)

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

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

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

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

4 收斂性分析

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

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

證明: 不失一般性,設(shè)當(dāng)k→∞時(shí),{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,…}. 故

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

(i) 若對(duì)任意k∈J有αk≥α*>0,其中α*為某固定常數(shù). 由式(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)矛盾. 故結(jié)論成立.

證明: 類(lèi)似文獻(xiàn)[4]中定理8可證.

5 數(shù)值算例

下面用MATLAB R2014在Intel(R) core(TM)i5-5200 CPU 2.7 GHz,4.0 GB內(nèi)存,Windows 7操作系統(tǒng)的計(jì)算機(jī)上編程,驗(yàn)證算法1的性能,且每個(gè)算例均隨機(jī)生成15個(gè)問(wèn)題進(jìn)行求解. 參數(shù)取值為σ=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. 隨機(jī)生成矩陣N∈n×n和向量q∈n,令矩陣M=NTN,s=Mx+q,以‖G(z)‖≤10-6為算法終止條件,(x0,s0)為初始點(diǎn),記AIT為平均迭代次數(shù),ACPU為平均CPU時(shí)間.

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

表1 算法1求解不同初始點(diǎn)和不同p值LWSOCCP的數(shù)值結(jié)果

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

表2 算法1和算法NS求解LWSOCCP的數(shù)值結(jié)果

綜合表1和表2可見(jiàn),算法1可有效求解LWSOCCP,且收斂速度快,迭代次數(shù)少,因而算法性能穩(wěn)定.