王 艷， 殷俊鋒， 李 蕊,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ù)且滿足對(duì)角可微[1-2]，即fi，f的第i個(gè)分量?jī)H僅是zi的可微函數(shù). 其中不等號(hào)是分量意義下的不等號(hào),zT表示向量z的轉(zhuǎn)置,特別地,當(dāng)函數(shù)f為線性函數(shù)時(shí),問題(1)化為線性互補(bǔ)問題.

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

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

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

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

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

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

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

(2)

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

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

(3)

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

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

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

(4)

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

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

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

(5)

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

2 收斂性分析

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

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

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

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

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

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

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

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

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

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

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

(6)

由于f是對(duì)角可微的，根據(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ù)ε的一個(gè)區(qū)域，可以保證算法2的收斂. 在文獻(xiàn)[11]中，給出求解線性互補(bǔ)問題的理論最優(yōu)參數(shù)，

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

(15)

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

(16)

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

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

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

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

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

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

表1給出了松弛參數(shù)ε的選取策略.對(duì)于例1，式 (14) 中取D=I. 基于定理1，通過簡(jiǎn)單地計(jì)算，對(duì)于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í)幾種迭代法的迭代步數(shù)(記作IT)和迭代時(shí)間(記作CPU).

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

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

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

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

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

表3 矩陣維數(shù)變化時(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ù)隨參數(shù)ε的變化Tab.4 IT and CPU for three methods with the change of ε for Example 2

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

對(duì)于例2，式(14)中取D=diag[(〈FΩ〉-|GΩ|)-1e]. 基于定理1，通過簡(jiǎn)單的計(jì)算，對(duì)于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í)幾種迭代法的迭代步數(shù)和迭代時(shí)間.從表4中可以看出，當(dāng)松弛參數(shù)ε<1時(shí)，算法2需要的迭代步數(shù)和計(jì)算時(shí)間比算法1多，當(dāng)ε增大時(shí)，對(duì)于RMJ方法，迭代步數(shù)和計(jì)算時(shí)間是不斷下降的，而RMGS方法和RMSOR方法的迭代步數(shù)和計(jì)算時(shí)間先減少后增大. 算法2的實(shí)際收斂區(qū)域比式(18)、式(19)要大，因此，收斂區(qū)域還可以進(jìn)行改進(jìn).

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

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

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

表5 矩陣維數(shù)變化時(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+-矩陣時(shí)的收斂性，并且給出了最優(yōu)參數(shù)的選取方法. 數(shù)值實(shí)驗(yàn)進(jìn)一步驗(yàn)證了松弛模系矩陣分裂迭代法的有效性. 數(shù)值結(jié)果表明,松弛模系矩陣分裂迭代法無(wú)論在迭代步數(shù)還是迭代時(shí)間上均優(yōu)于模系矩陣分裂迭代法.