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

(1.福建師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，福建，福州 350117；2.中國科學(xué)院大學(xué)計算地球動力學(xué)重點實驗室，北京 100049)

考慮如下形式的非線性互補問題(NCP)：求x∈Rn，滿足

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

(1)

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

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

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

其中(a，b)∈R2.

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

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

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

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

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

1 等價非線性方程組的建立

為了便于建立和分析本文的主要結(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) 問題(1)中的fi(i∈[n]) ，存在一個正常數(shù)L滿足

|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)的.

證明 對于所有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)滿足假設(shè)1，那么由式(2)定義的F(x)是利普希茲連續(xù)的.

證明 對于任意的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)滿足假設(shè)1，則有

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

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

因此，有

證畢.

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

F(x)=0.

(3)

2 算法

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

(4)

(5)

進一步，結(jié)合線搜索方案[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 計算搜索方向

dk=-αkF(xk).

(6)

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

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

(7)

步4 更新迭代

(8)

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

(9)

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

下面引理說明算法1適定.

引理1在假設(shè)1 下，對于所有的k≥1 ，存在一個非負數(shù)mk滿足算法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)不滿足任意的非負整數(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 的一個解并且

那么有

‖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*‖ } 是一個遞減序列．對于所有k≥1，有‖xk-x*‖ ≤ ‖x0-x*‖，意味著 {xk} ?{x∈Rn：‖x-x*‖ ≤ ‖x0-x*‖ }．因此， {xk} 是一個有界序列．因為F(x) 是連續(xù)的，{ ‖F(xiàn)(xk) ‖} 是有界的.

(16)

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

由式(15)，有

(17)

由式(8)，有

定理4在假設(shè)1下，設(shè)序列 {xk} 由算法1 生成，且ε=0， 則序列 {xk} 收斂到某個點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ù)值實驗

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

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

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

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

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

其中τ=10-8，且

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

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

下面給出具體實例.

例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的取值會影響NMGPM算法的迭代次數(shù)，對算法的數(shù)值性能產(chǎn)生影響，因此需要選擇適合的λk.

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

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

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

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

5 結(jié)論

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