林 婷，柯藝芬，張 振，馬昌鳳

(1.福建師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，福建，福州 350117；2.中國(guó)科學(xué)院大學(xué)計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100049)

考慮如下形式的非線(xiàn)性互補(bǔ)問(wèn)題(NCP)：求x∈Rn，滿(mǎn)足

x≥0,f(x) ≥0,xTf(x)=0,

(1)

非線(xiàn)性互補(bǔ)問(wèn)題在運(yùn)籌學(xué)[1-2]、工程設(shè)計(jì)[3-4]和經(jīng)濟(jì)均衡[4-5]等方面有許多重要的應(yīng)用，并發(fā)展了許多求解非線(xiàn)性互補(bǔ)問(wèn)題的數(shù)值方法.

近年來(lái)，光滑函數(shù)法在求解NCP問(wèn)題上引起了廣泛的關(guān)注[6]，主要原因是其優(yōu)越的數(shù)值性能．光滑函數(shù)法的思想是使用一個(gè)光滑函數(shù)將非線(xiàn)性互補(bǔ)問(wèn)題(1)重新表述為光滑方程組，在迭代步中使用牛頓法近似求解光滑方程組.

在光滑函數(shù)中,Fischer-Burmeister (FB)函數(shù)受到廣泛的關(guān)注.FB函數(shù)的定義如下

其中(a，b)∈R2.

然而，由于在牛頓法中求解雅可比矩陣及其逆矩陣需要大量的計(jì)算工作，因此求解光滑方程組代價(jià)是昂貴的．為了減少計(jì)算量，無(wú)導(dǎo)數(shù)法以及梯度法是目前比較有效的方法.

Jiang[7]利用平方FB函數(shù)轉(zhuǎn)化的等價(jià)非線(xiàn)性方程組，在非線(xiàn)性映射方向可微且強(qiáng)單調(diào)情況下建立了一種無(wú)導(dǎo)數(shù)下降法；Mangasarian 等[8]通過(guò)最小化隱式拉格朗日函數(shù)，給出了強(qiáng)單調(diào)NCP的無(wú)導(dǎo)數(shù)下降方法，并建立了該方法的全局收斂性；Ma 等[9-10]提出了求解NCP的光滑Broyden-like 方法，該方法基于光滑逼近的FB 函數(shù)，并利用了無(wú)導(dǎo)數(shù)線(xiàn)搜索規(guī)則，證明了該算法在適當(dāng)?shù)臈l件下具有全局和超線(xiàn)性收斂性.

最近，Yu 等[11]提出了一類(lèi)求解凸約束的非線(xiàn)性單調(diào)方程的改進(jìn)譜梯度投影方法，其步長(zhǎng)由長(zhǎng)Barzilli-Borwein (BB)步長(zhǎng)和短BB步長(zhǎng)的凸組合決定，該算法具有良好的數(shù)值性能.

基于上述工作的啟發(fā)，本文提出了一類(lèi)求解非線(xiàn)性互補(bǔ)問(wèn)題的新修正譜梯度投影方法．首先，將問(wèn)題(1)重新表述為一個(gè)基于FB 函數(shù)的單調(diào)且利普希茲連續(xù)的非線(xiàn)性方程組．然后，結(jié)合線(xiàn)搜索方案[12]和多元譜梯度投影方法[13]以及改進(jìn)譜梯度投影方法[11]，對(duì)得到的單調(diào)非光滑系統(tǒng)提出了一種新修正譜梯度投影方法，并對(duì)所提出的方法建立了全局收斂性．通過(guò)數(shù)值算例，表明了所提出算法的有效性.

本文符號(hào)說(shuō)明如下:‖·‖記為歐幾里得范數(shù),〈·,·〉記為歐幾里得內(nèi)積,(·)T記為向量或矩陣的轉(zhuǎn)置,[n]∶={1，2，…，n}.

1 等價(jià)非線(xiàn)性方程組的建立

為了便于建立和分析本文的主要結(jié)果，先給出如下定義.

(x-y)T(f(x)-f(y))≥0.

(xi-yi)(fi(x)-fi(y)) ≥0,?i∈[n].

‖F(xiàn)(x)-F(y)‖ ≤L‖x-y‖.

(b) 問(wèn)題(1)中的fi(i∈[n]) ，存在一個(gè)正常數(shù)L滿(mǎn)足

|fi(x)-fi(y)|≤L|xi-yi|, ?x,y∈Rn.

設(shè)

F(x)∶=x+f(x)-g(x),

(2)

定理1若f(x) 是一致P0-映射，則式(2)定義的F(x) 是單調(diào)的.

證明 對(duì)于所有x，y∈Rn，i∈[n]，考慮以下兩種情況.

若xi=yi=0，顯然有(xi-yi)(fi(x)-Fi(y))=0.

若xi≠0 或yi≠0，則有

可得

和

根據(jù)假設(shè)條件，有

和

(xi-yi)(fi(x)-fi(y)).

因此，有

(xi-yi)(fi(x)-Fi(y))=(xi-yi)[(xi+fi(x)-gi(x))-(yi+fi(y)-gi(y))]=

(xi-yi)2+(xi-yi)(fi(x)-fi(y))-(xi-yi)(gi(x)-gi(y))=

由上述證明，有

證畢．

定理2若fi(x)滿(mǎn)足假設(shè)1，那么由式(2)定義的F(x)是利普希茲連續(xù)的.

證明 對(duì)于任意的x，y∈Rn,i∈[n],令u=(xi，fi(x))T,v=(yi,fi(y))T,e= (1,1)T，則有

φFB(xi，fi(x))=φFB(u)=‖u‖-eTu,

和

φFB(yi，fi(y))=φFB(v)=‖v‖-eTv.

由于fi(x)滿(mǎn)足假設(shè)1，則有

|Fi(x)-Fi(y) |=|φFB(u)-φFB(v) |=

|‖u‖-eTu-‖v‖+eTv|≤

因此，有

證畢.

根據(jù)定理1 和定理2，在假設(shè)1 的條件下，求解問(wèn)題(1) 等價(jià)于求解一個(gè)基于FB 函數(shù)的單調(diào)且利普希茲連續(xù)的非線(xiàn)性方程組

F(x)=0.

(3)

2 算法

采用許多有效的算法來(lái)解決由問(wèn)題(1) 轉(zhuǎn)化的等價(jià)的非線(xiàn)性方程組(3)，如牛頓法、無(wú)導(dǎo)數(shù)法和梯度法．由于在牛頓法中求解雅可比矩陣及其逆矩陣需要大量的計(jì)算工作，因此求解光滑方程組代價(jià)是昂貴的．為了減少計(jì)算量，無(wú)導(dǎo)數(shù)法和梯度法是目前比較有效的方法．在本節(jié)中，提出了一種新修正譜梯度投影方法(New Modified Gradient Projection Method，簡(jiǎn)稱(chēng)NMGPM)來(lái)解決問(wèn)題(1)．該方法結(jié)合了線(xiàn)搜索技術(shù)[12]和多元譜梯度投影方法[13]以及改進(jìn)譜梯度投影方法[11].

(4)

(5)

進(jìn)一步，結(jié)合線(xiàn)搜索方案[12]和多元譜梯度投影方法[13]，得出如下算法.

算法1(NMGPM 算法)

步0 給定ε>0,x0∈Rn，設(shè)σ,β∈(0,1),α0=1,r∈[0,1],置k∶=0.

步1 若‖F(xiàn)(xk)‖≤ε，停止，并接受xk為近似解.

步2 計(jì)算搜索方向

dk=-αkF(xk).

(6)

步3 計(jì)算zk=xk+θkdk，其中θk∶=βmk,mk是最小的非負(fù)整數(shù)使得

-F(xk+θkdk)Tdk≥σγkθk‖dk‖2,

(7)

步4 更新迭代

(8)

步5 選取λk∈[0,1]，更新譜向量

(9)

步6 置k:=k+1，轉(zhuǎn)步1.

下面引理說(shuō)明算法1適定.

引理1在假設(shè)1 下，對(duì)于所有的k≥1 ，存在一個(gè)非負(fù)數(shù)mk滿(mǎn)足算法1.

(10)

(11)

(L2+Lr+(L+r)r)〈sk-1，sk-1〉=(L+r)2〈sk-1，sk-1〉.

由αk的定義式(9)，可得

(12)

下面用反證法證明結(jié)論．假設(shè)存在k0≥1 使得式(7)不滿(mǎn)足任意的非負(fù)整數(shù)m，即

-〈F(xk0+βmdk0) ，dk0〉<σγk0βm‖dk0‖2.

-〈F(xk0) ，dk0〉 ≤0.

(13)

然而，由dk的定義式(6)和式(12)，有

(14)

這與式(13)矛盾.證畢.

3 全局收斂性

引理2[14]設(shè)F(x) 是單調(diào)的，且x，z∈Rn使得F(z)T(x-z)>0．設(shè)x*為F(x)=0 的一個(gè)解并且

那么有

‖x+-x*‖2≤ ‖x-x*‖2-‖x+-x‖2.

定理3在假設(shè)1下，設(shè)序列 {xk} 由算法1 生成，且ε=0，則有

證明 根據(jù)假設(shè)1，f(x) 是一致P0-映射，所以F(x) 是單調(diào)的．由算法1 中的式(7)，有

F(zk)T(xk-zk)=-βmkF(zk)Tdk≥σγkβmk‖dk‖2>0.

由引理2 有

‖xk+1-x*‖2≤‖xk-x*‖2-‖xk+1-xk‖2,

(15)

其中x*是F(x)=0 的解．因此，{ ‖xk-x*‖ } 是一個(gè)遞減序列．對(duì)于所有k≥1，有‖xk-x*‖ ≤ ‖x0-x*‖，意味著 {xk} ?{x∈Rn：‖x-x*‖ ≤ ‖x0-x*‖ }．因此， {xk} 是一個(gè)有界序列．因?yàn)镕(x) 是連續(xù)的，{ ‖F(xiàn)(xk) ‖} 是有界的.

(16)

由式(16)， {dk} 也是有界的且 {zk} 也是有界的．因此，存在常數(shù)M>0 使得對(duì)于任意的k，有 ‖F(xiàn)(zk) ‖ ≤M．

由式(15)，有

(17)

由式(8)，有

定理4在假設(shè)1下，設(shè)序列 {xk} 由算法1 生成，且ε=0， 則序列 {xk} 收斂到某個(gè)點(diǎn)x*使得F(x*)=0.

證明 由定理3，有

(18)

由算法1和式(12)，有

(19)

由式(16)， {dk} 是有界的．下面分兩種情況證明.

當(dāng)

由式(19)，有

由式(18)，可得

由算法1 ，有

-〈F(xk+βmk-1dk) ，dk〉<σγkβmk-1‖dk‖2.

(20)

(21)

4 數(shù)值實(shí)驗(yàn)

為了測(cè)試所提出的NMGPM方法的數(shù)值性能，本節(jié)給出了一些初步的數(shù)值結(jié)果.

首先討論凸組合系數(shù)λk的選擇，并為數(shù)值實(shí)驗(yàn)做了準(zhǔn)備．最簡(jiǎn)單的方案是在所有迭代中固定λk，即在[0,1] 內(nèi)選擇一個(gè)常量作為λk的值．對(duì)于所有的k，分別設(shè)置λk=0.1，0.3，0.618和0.9．實(shí)驗(yàn)中則選擇實(shí)驗(yàn)效果最好的λk作為最后的結(jié)果.

將算法1 與多元譜梯度投影法(簡(jiǎn)稱(chēng)MSGP)[13]和改進(jìn)的梯度投影法(簡(jiǎn)稱(chēng)MGP)[15]進(jìn)行了比較.在MATLAB R2018a軟件上進(jìn)行數(shù)值實(shí)驗(yàn)．若 ‖F(xiàn)(xk)‖ ≤10-5，則終止運(yùn)行.

算法參數(shù)設(shè)置如下：

(1) 對(duì)于MSGP算法，設(shè)置ρ=0.5，σ=0.001，=10-10，r=0.01,

其中τ=10-8，且

(3)對(duì)于NMGPM算法，設(shè)置α0=1,β=0.618,σ=0.01,r=0.001．

所有的實(shí)驗(yàn)起始點(diǎn)為x0=rand(n，1)，實(shí)驗(yàn)結(jié)果為5 次實(shí)驗(yàn)的平均值.

下面給出具體實(shí)例.

例1函數(shù)f(x)的組成部分描述如下

fi(x)=ln(xi)-sin(|xi-2|)，i∈[n].

例2函數(shù)f(x)的組成部分描述如下

fi(x)=exi-1，i∈[n].

例3函數(shù)f(x)的組成部分描述如下

例4函數(shù)f(x)的組成部分描述如下

例5函數(shù)f(x)的組成部分描述如下

fi(x)=2xi-sin(|xi|)，i∈[n].

例6函數(shù)f(x)的組成部分描述如下

fi(x)=xi-sin(|xi-1|)，i∈[n].

例7函數(shù)f(x)的組成部分描述如下

表1給出了例2-7的數(shù)值結(jié)果.

表1 例2-7的數(shù)值結(jié)果Tab.1 The numerical results for examples 2-7

圖1 例1中λk與IT關(guān)系圖Fig.1 The relationship between λk and IT for example 1

圖1給出了例1的數(shù)值結(jié)果.λk的取值會(huì)影響NMGPM算法的迭代次數(shù)，對(duì)算法的數(shù)值性能產(chǎn)生影響，因此需要選擇適合的λk.

數(shù)值結(jié)果表明，對(duì)于例2-7，所提出的NMGPM方法與MSGP方法相比，在CPU 時(shí)間上有明顯的改進(jìn)，因?yàn)镸SGP方法中需要對(duì)譜梯度αk的每個(gè)元素進(jìn)行賦值，需要花費(fèi)一定的時(shí)間.而NMGPM方法直接對(duì)αk進(jìn)行整體的賦值，花費(fèi)的時(shí)間更少.3 種方法中MGP方法使用的時(shí)間最少，原因是NMGPM方法有MGP方法中沒(méi)有的線(xiàn)搜索技術(shù)，在線(xiàn)搜索步中尋找mk的過(guò)程需要花費(fèi)一定的時(shí)間.MSGP方法在階數(shù)為10 000時(shí)沒(méi)有結(jié)果，原因是算法運(yùn)行過(guò)程中對(duì)譜梯度αk的每個(gè)元素進(jìn)行賦值，階數(shù)越高花費(fèi)的時(shí)間越多，在階數(shù)為10 000時(shí)需要花費(fèi)超過(guò)1 000 s的CPU時(shí)間.

在迭代次數(shù)方面，所提出的NMGPM方法的次數(shù)最少，與MSGP方法比較相差較小，與MGP比較有明顯的減少．并且NMGPM方法的迭代次數(shù)在實(shí)驗(yàn)階數(shù)10 倍增長(zhǎng)的情況下沒(méi)有明顯的增長(zhǎng),原因是NMGPM方法使用一種新的線(xiàn)搜索方法，找到了3 種算法中最優(yōu)的下降方向,因此NMGPM方法的迭代次數(shù)最小.

在誤差方面，3 種方法都達(dá)到了終止條件中的要求，無(wú)明顯差別.

總體來(lái)說(shuō)，NMGPM 方法在CPU 時(shí)間和迭代次數(shù)上都有良好的表現(xiàn)，尤其在迭代次數(shù)上有明顯優(yōu)勢(shì)，該算法對(duì)于求解非線(xiàn)性互補(bǔ)問(wèn)題是有效的.

5 結(jié)論

本文將非線(xiàn)性互補(bǔ)問(wèn)題重新表述為非線(xiàn)性方程組，并研究了一定假設(shè)條件下，得到的非線(xiàn)性方程組中映射的單調(diào)性與利普希茲連續(xù)性．在此基礎(chǔ)上，提出了一種新修正的譜梯度投影方法來(lái)求解NCP，并建立了該方法的全局收斂性．與原來(lái)的多元譜梯度投影方法相比，所提出的NMGPM 算法中譜向量的更新方案成本更低，線(xiàn)搜索規(guī)則提高了數(shù)值性能．與改進(jìn)的梯度投影法相比，所提出的NMGPM 算法具有更少的迭代步數(shù)，新的線(xiàn)搜索規(guī)則發(fā)揮了很大作用．初步的數(shù)值結(jié)果表明，該方法優(yōu)于現(xiàn)有的一些方法，對(duì)于求解非線(xiàn)性互補(bǔ)問(wèn)題是有效的．然而，所提出的NMGPM 算法中λk的選擇還有待進(jìn)一步改進(jìn)，若能找到最優(yōu)的λk，相信算法的性能會(huì)有更好的表現(xiàn).