王 艷， 殷俊鋒， 李 蕊,2

(1. 同濟(jì)大學(xué) 數(shù)學(xué)科學(xué)學(xué)院， 上海 200092； 2. 嘉興學(xué)院 數(shù)理與信息工程學(xué)院， 浙江 嘉興 314001)

給定大型稀疏矩陣M∈Rn×n和n維向量q∈Rn,考慮如下一類非線性互補(bǔ)問題:求向量z,w∈Rn滿足

z≥0,w∶=Mz+q+f(z)≥0,zTw=0

(1)

這里，f:Rn→Rn為給定的n元函數(shù)且滿足對角可微[1-2]，即fi，f的第i個分量僅僅是zi的可微函數(shù). 其中不等號是分量意義下的不等號,zT表示向量z的轉(zhuǎn)置,特別地,當(dāng)函數(shù)f為線性函數(shù)時,問題(1)化為線性互補(bǔ)問題.

互補(bǔ)問題在科學(xué)與工程等領(lǐng)域有著廣泛的應(yīng)用[3-4], 如變分不等式、彈性接觸問題、期權(quán)定價等自由邊界問題均可化為互補(bǔ)問題. 由于互補(bǔ)問題的重要性，引起了數(shù)學(xué)工作者的廣泛關(guān)注和高度重視. 近幾十年來， 其數(shù)值解的研究發(fā)展迅速，取得了豐碩的成果. Bai[5]提出了求解線性互補(bǔ)問題的模系矩陣分裂迭代法,并且為模系矩陣分裂迭代方法提供了一個基本框架. 由于其迭代格式簡單，在文獻(xiàn)[5]的基礎(chǔ)上，模系矩陣分裂迭代法得到了進(jìn)一步研究，比如，兩步模系矩陣分裂迭代法[6]、加速的模系矩陣分裂迭代法[7]和松弛模系矩陣分裂迭代方法[8]等. 隨著高性能并行計算機(jī)系統(tǒng)的快速發(fā)展,Bai等[9]利用并行計算的思想, 提出了模系同步多分裂迭代法、模系同步二級多分裂迭代方法[10]和兩步模系同步多分裂迭代方法[11].

針對非線性互補(bǔ)問題(1), 當(dāng)函數(shù)f對角可微時, Xia等在文獻(xiàn)[12]中首次將求解線性互補(bǔ)問題的模系矩陣分裂迭代法推廣至該非線性互補(bǔ)問題, 構(gòu)造了相應(yīng)的模系矩陣分裂迭代法, 并給出了系數(shù)矩陣為正定或H+-矩陣時的收斂理論. 隨后, Li等[13]采用加速模系矩陣分裂迭代法對該非線性互補(bǔ)問題進(jìn)行求解, 并且給出了矩陣M為H+-矩陣時的收斂性分析. 本文提出用松弛模系矩陣分裂迭代法求解該非線性互補(bǔ)問題, 討論該迭代法在系數(shù)矩陣M為H+-矩陣時的收斂條件， 并且給出了最優(yōu)參數(shù)的選取. 數(shù)值實驗表明, 對于大型稀疏非線性互補(bǔ)問題, 松弛模系矩陣分裂迭代法不僅是有效的, 而且收斂效果優(yōu)于模系矩陣分裂迭代法.

1 松弛模系矩陣分裂迭代法

算法1[12]令M=F-G為矩陣M的一個分裂，給定初始向量x(0)∈Rn, 對于k=1,2,…, 從下列方程組中解出x(k)：

為了提高計算效率，下面給出松弛模系矩陣分裂迭代方法.

對于任意給定的正對角矩陣Ω和Γ，令z=Ω(|x|+x),w=Γ(|x|-x)，MΩ=FΩ-GΩ為MΩ的一個分裂，x滿足如下不動點方程

(Γ+FΩ)x=GΩx+(Γ-MΩ)|x|-(q+f(z))

(2)

在不動點格式(2)對應(yīng)的迭代格式中，考慮對原向量和更新的向量做一個加權(quán)，通過引入一個非奇異的參數(shù)矩陣P∈Rn×n，建立如下松弛模系矩陣分裂迭代方法.

算法2 對于任意給定的正對角矩陣Ω和Γ，令MΩ=FΩ-GΩ為矩陣MΩ∈Rn×n的一個分裂. 給定初始向量x(0)∈Rn，對于k=1,2,…，通過求解如下系統(tǒng):

(3)

得到x(k)，令z(k)=Ω(|x(k)|+x(k)).

系統(tǒng) (3) 等價于

(Γ-MΩ)|x(k-1)|-q-f(Ω(|x(k-1)|+x(k-1)))]

(4)

算法2提供了松弛模系矩陣分裂迭代方法的一個基本框架. 特別地，取

當(dāng)α,β分別取α=β,α=β=1和α=1,β=0時，稱為松弛模系逐步超松弛迭代格式(RMSOR)、 松弛模系高斯-賽德爾迭代格式(RMGS)和松弛模系雅克比迭代格式(RMJ).

特別地，令P=εI, 系統(tǒng) (3) 等價為如下形式：

(5)

這里，ε∈R{0}, 并且ε隨著k的變化而變化.

2 收斂性分析

H+-矩陣在數(shù)學(xué)物理問題、控制論、電力系統(tǒng)理論、經(jīng)濟(jì)數(shù)學(xué)以及彈性力學(xué)等眾多領(lǐng)域中都有廣泛應(yīng)用，例如經(jīng)濟(jì)價值模型、 反網(wǎng)絡(luò)分析模型、經(jīng)濟(jì)學(xué)中的投入產(chǎn)出增長模型和概率統(tǒng)計中的Markov鏈等問題.

給定兩個實矩陣M=(mij),N=(nij)∈Rm×n, 如果對任意的1≤i≤m,1≤j≤n都有mij≥nij(mij>nij), 則記M≥N(M>N). 本文用記號|M|=(|mij|)∈Rm×n表示M的絕對值.

令M∈Rn×n為一個實n×n矩陣，它的比較矩陣〈M〉=(〈m〉ij)定義為

若對任意i≠j, 有mij≤0, 則稱M為Z-矩陣. 若M為非奇異的Z-矩陣且M-1≥O,O為零矩陣, 則稱M為M-矩陣. 如果〈M〉為M-矩陣, 則稱M為H-矩陣. 對于H-矩陣M, 有M非奇異, 且|M-1|≤〈M〉-1. 特別地, 對角元為正的H-矩陣稱為H+-矩陣.

給定M∈Rn×n, 設(shè)M=F-G, 如果F非奇異, 則稱M=F-G為矩陣M的一個分裂; 如果F是非奇異的M-陣并且G≥0, 則稱分裂為M-分裂；當(dāng)〈F〉-|G|是一個M-矩陣，則稱M=F-G是一個H-分裂；如果〈M〉=〈F〉-|G|, 則稱分裂為H-相容分裂.

引理1[14]如果矩陣M∈Rn×n是嚴(yán)格對角占優(yōu)的，對于任意矩陣N∈Rn×n，

引理2[15]矩陣M是一個對角元為正的Z-陣,M是一個M-陣當(dāng)且僅當(dāng)存在一個正對角矩陣D, 使得MD是對角元為正的嚴(yán)格對角占優(yōu)矩陣.

引理3[8]令M∈Rn×n是一個H+-矩陣，MΩ=FΩ-GΩ為矩陣MΩ的H-分裂.則存在一個正對角矩陣D使得(〈FΩ〉-|GΩ|)D和(Γ+〈FΩ〉)D是兩個嚴(yán)格對角占優(yōu)矩陣，且有

[(〈MΩ〉+〈FΩ〉-|GΩ|)De]i>0,i=1,2,…,n

下面,考慮算法2的收斂性. 假設(shè)對角可微函數(shù)f滿足

(Γ-MΩ)|x*|-q-f(Ω(|x*|+x*))]

(6)

由于f是對角可微的，根據(jù)均值定理，有

f(z(k))-f(z*)=f(Ω(|x(k)|+x(k)))-

f(Ω(|x*|+x*))=

J(k)Ω(|x(k)|-|x*|+x(k)-x*)

x(k)-x*=(I-P)(x(k-1)-x*)+

(Γ-MΩ)(|x(k-1)|-|x*|)-

(f(Ω(|x(k-1)|+x(k-1)))-

f(Ω(|x*|+x*)))]=

(I-P)(x(k-1)-x*)+

(Γ-MΩ)(|x(k-1)|-|x*|)-

J(k-1)Ω(|x(k-1)|-|x*|+x(k-1)-x*)]

(7)

|x(k)-x*|=

(8)

(9)

(10)

進(jìn)而，有

{[Γ+〈FΩ〉-(Γ+DFΩ)|P-I|-

|BFΩ||P-I|-|BGΩ|P-|BMΩ|P]De}i.

(11)

情況1：如果(Pe)i≥1，由式(11)可得到

|BGΩ|P-|BMΩ|P]De}i≥{[Γ+〈FΩ〉-

(Γ+DFΩ)(P-I)-|DGΩ|P-(Γ-DMΩ)P-

|BFΩ|(P-I)-|BGΩ|P-|BMΩ|P]De}i=

{[2Γ+2DFΩ-(2Γ+|FΩ|+|GΩ|-

〈MΩ〉)P]De}i>0

(12)

情況2：如果0<(Pe)i<1，有

|BGΩ|P-|BMΩ|P]De}i≥{[Γ+〈FΩ〉-

(Γ+DFΩ)(I-P)-|DGΩ|P-(Γ-DMΩ)P-

|BFΩ|(I-P)-|BGΩ|P-|BMΩ|P]De}i=

{[(〈MΩ〉+|FΩ|-|GΩ|)P-

2|BFΩ|]De}i>0.

(13)

結(jié)合式(12)和式(13)，得到不等式(10).證畢.

基于引理4，建立如下收斂性定理.

ai<ε

其中，D由引理3給出，則

i=1,2,…,n

(14)

由算法2生成的迭代序列收斂到非線性互補(bǔ)問題(1)的解x*.

[(2Γ+2DFΩ)De]i-[(2Γ+|FΩ|+|GΩ|-

〈MΩ〉)De]i=[(〈MΩ〉+〈FΩ〉-

|GΩ|)De]i>0

以及

[(〈MΩ〉+|FΩ|-|GΩ|)De]i-2(|BFΩ|De)i=

[(〈MΩ〉+〈FΩ〉-|GΩ|)De]i>0

然后，有0

注1 定理1給出了松弛參數(shù)ε的一個區(qū)域，可以保證算法2的收斂. 在文獻(xiàn)[11]中，給出求解線性互補(bǔ)問題的理論最優(yōu)參數(shù)，

由于x*未知，在具體實現(xiàn)時，x*用x(k-1/2)代替，并且用k-1代替k. 因此，松弛參數(shù)ε的一個合理的選擇為

(15)

在第k步，松弛參數(shù)ε的選取方法如下：

(16)

這里,a,b由定理1給出.

3 數(shù)值實驗

本節(jié)通過數(shù)值實驗驗證幾種松弛模系矩陣分裂迭代法求解非線性互補(bǔ)問題(1)的有效性.

在數(shù)值實驗中, 初始向量取為x(0)=(1,1,…,1)T∈Rn,Γ=2DM，對于MSOR方法和RMSOR方法，松弛參數(shù)取α=1.2， 停止準(zhǔn)則設(shè)定為

所有計算均在一臺CPU為2.40 GHz和內(nèi)存為4.00 GB的機(jī)器上運行, 編程語言為MATLAB.

為塊三對角矩陣, 其中S=tridiag(-1,4,-1)∈Rm×m為三對角矩陣,

表1給出了松弛參數(shù)ε的選取策略.對于例1，式 (14) 中取D=I. 基于定理1，通過簡單地計算，對于RMJ，RMGS，RMSOR方法, 策略2中的(a,b)為

(0,1.200 0),(0.333 3,1.200 0),

(0.428 6,1.133 3)

(17)

策略3中(a,b)為

(0,1.333 3),(0,1.333 3),(0,1.259 3)

(18)

表2分別列出當(dāng)n=1 600, 松弛參數(shù)ε在0.5～1.5之間變化時幾種迭代法的迭代步數(shù)(記作IT)和迭代時間(記作CPU).

從表2中可以看出，當(dāng)松弛參數(shù)ε<1時，算法2需要的迭代步數(shù)和計算時間比算法1多，當(dāng)ε增大時，對于RMJ方法，迭代步數(shù)和計算時間是不斷下降的，而RMGS方法和RMSOR方法的迭代步數(shù)和計算時間先減少后增大. 當(dāng)ε不在收斂區(qū)域時，迭代方法不收斂.

此外，ε≥1時三種方法比ε<1時收斂得快，這表明在公式 (15) 中，精確解更靠近x(k-1/2). 松弛參數(shù)的三種選取策略數(shù)值實驗顯示在表2的最后三行. 從表2最后三行可以看出，基于三種選取策略的算法2明顯優(yōu)于算法1，且策略1是最優(yōu)策略.

表1 3種策略描述Tab.1 Description of three strategies

表2 例1的3種迭代法的迭代時間和迭代步數(shù)隨參數(shù)ε的變化Tab.2 IT and CPU for three methods with the change of ε for Example 1

當(dāng)矩陣維數(shù)不斷增大時，數(shù)值實驗結(jié)果顯示在表3中. 從表3中可以得出，無論是計算時間還是迭代步數(shù)，基于2種策略的算法2明顯優(yōu)于算法1.

表3 矩陣維數(shù)變化時例1兩種算法比較Tab.3 Comparison of two algorithms with the change of n for Example 1

例2 設(shè)m為給定的正整數(shù),n=m2. 在式(1)中取

表4 例2的三種迭代法的迭代時間和迭代步數(shù)隨參數(shù)ε的變化Tab.4 IT and CPU for three methods with the change of ε for Example 2

為塊三對角矩陣,S=tridiag(-1,4,-1)∈Rm×m為三對角矩陣,

對于例2，式(14)中取D=diag[(〈FΩ〉-|GΩ|)-1e]. 基于定理1，通過簡單的計算，對于RMJ，RMGS，RMSOR方法, 策略2的中(a,b)為

(0,1.006 2),(0.930 8,1.006 2)

(18)

(0.748 9,1.062 5)

(19)

策略3中的(a,b)為

(0,1.294 8),(0,1.294 8),(0,1.199 1)

(20)

表4分別列出當(dāng)n=1 600, 松弛參數(shù)ε在0.5～1.5之間變化時幾種迭代法的迭代步數(shù)和迭代時間.從表4中可以看出，當(dāng)松弛參數(shù)ε<1時，算法2需要的迭代步數(shù)和計算時間比算法1多，當(dāng)ε增大時，對于RMJ方法，迭代步數(shù)和計算時間是不斷下降的，而RMGS方法和RMSOR方法的迭代步數(shù)和計算時間先減少后增大. 算法2的實際收斂區(qū)域比式(18)、式(19)要大，因此，收斂區(qū)域還可以進(jìn)行改進(jìn).

松弛參數(shù)的三種選取策略數(shù)值實驗顯示在表4的最后三行. 從表4最后三行可以看出，基于三種選取策略的算法2明顯優(yōu)于算法1，且策略1是最優(yōu)策略.

當(dāng)矩陣維數(shù)不斷增大時，數(shù)值實驗結(jié)果顯示在表5中. 從表5可以得出，無論是計算時間還是迭代步數(shù)，基于三種策略的算法2明顯優(yōu)于算法1.

根據(jù)實驗結(jié)果，得出策略1是最優(yōu)策略，在圖1和圖2中，分別畫出基于策略1的算法2和算法1隨著迭代步數(shù)變化的殘差下降曲線. 從兩張圖中可以看出,算法2收斂速度明顯快于算法1，表明提出的新方法是有效的且明顯優(yōu)于模系矩陣分裂迭代方法.

表5 矩陣維數(shù)變化時例2兩種算法比較Tab.5 Comparison of two algorithms with the change of n for Example 2

圖1 例1兩種算法的殘差比較 Fig.1 Residual comparison of two algorithms for Example 1

圖2 例2兩種算法的殘差比較 Fig.2 Residual comparison of two algorithms for Example 2

4 結(jié)論

本文構(gòu)造了求解一類非線性互補(bǔ)問題的松弛模系矩陣分裂迭代法,理論上分析了當(dāng)矩陣M為H+-矩陣時的收斂性，并且給出了最優(yōu)參數(shù)的選取方法. 數(shù)值實驗進(jìn)一步驗證了松弛模系矩陣分裂迭代法的有效性. 數(shù)值結(jié)果表明,松弛模系矩陣分裂迭代法無論在迭代步數(shù)還是迭代時間上均優(yōu)于模系矩陣分裂迭代法.