求解一類隱式互補問題的加速模系矩陣分裂迭代法

殷俊鋒，丁 戩，李 蕊

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

1 背景介紹

互補問題廣泛地出現(xiàn)在科學(xué)和工程的許多應(yīng)用中，如彈性接觸、經(jīng)濟(jì)運輸、流體動力學(xué)的邊界問題、凸二次規(guī)劃以及反問題等［1］，因此受到越來越多研究者的關(guān)注和研究，并已經(jīng)取得了豐碩的成果。本文考慮如下一類隱式互補問題［2］：求向量u和r∈Rn滿足

式中：M∈Rn×n為給定矩陣；q=(q1，q2，…，qn) ∈Rn為給定實向量；f為Rn到它本身的一個映射。特別地，當(dāng)函數(shù)f= 0，問題（1）化為線性互補問題［3］。

針對線性互補問題，白中治［4］提出了模系矩陣分裂迭代法，并對系數(shù)矩陣為正定或H+‐矩陣下迭代法的收斂性進(jìn)行了分析。由于模系矩陣分裂迭代算法格式簡單且收斂速度較快，該方法吸引了眾多學(xué)者的關(guān)注并得到了進(jìn)一步的研究。如兩步模系矩陣分裂迭代法［5］、廣義模系矩陣分裂迭代法［6］和松弛模系矩陣分裂迭代法［7］等。為了更快速有效地求解大型稀疏線性互補問題，白中治等［8］通過引進(jìn)并行計算技術(shù)，構(gòu)造了模系矩陣同步多分裂迭代法，并建立了此算法在系數(shù)矩陣為H+‐矩陣時的收斂性分析。夏澤晨等［9］首次構(gòu)造了求解非線性互補問題的模系矩陣分裂迭代法，進(jìn)一步地還有兩步模系矩陣分裂迭代法［10-12］、松弛模系矩陣分裂迭代法［13］以及加速模系矩陣分裂迭代法［14-15］等也相繼被構(gòu)造。

李郴良等［16-17］構(gòu)造求解隱式互補問題（1）的模系矩陣分裂迭代法和模系矩陣同步多分裂迭代法，并建立了當(dāng)系數(shù)矩陣為正定或H+‐矩陣時迭代法的收斂性質(zhì)。曹陽等［18］對文獻(xiàn)［16］中的模系矩陣分裂迭代法的收斂理論給出了更加完整的證明。王艷等［19］采用兩步模系矩陣分裂迭代法對隱式互補問題（1）進(jìn)行求解，并且分析了系數(shù)矩陣為H+‐矩陣時算法的收斂性。曹陽等［20］同樣構(gòu)造了兩步模系矩陣分裂迭代法求解隱式互補問題（1）并且建立了迭代法在系數(shù)矩陣為正定或H+‐矩陣時迭代法的收斂性質(zhì)。

鄭寧等［21-22］通過對向量邊計算邊更替的思想，將線性互補問題轉(zhuǎn)換成了新的隱式不動點方程，提出了加速模系矩陣分裂迭代法求解線性互補問題，并從理論上分析了算法在系數(shù)矩陣為正定或H+‐矩陣時的收斂性質(zhì)。通過選擇合適的系數(shù)矩陣，該方法也可以得到文獻(xiàn)［4］中的傳統(tǒng)模系矩陣分裂迭代法。為了進(jìn)一步加快求解隱式互補問題的收斂速度，本文引入對向量邊計算邊更替的思想，提出了加速模系矩陣分裂迭代法。理論分析建立了新方法在系數(shù)矩陣為H+‐矩陣時的收斂性質(zhì)。數(shù)值實驗表明所提方法是有效的，并在迭代步數(shù)和時間上均優(yōu)于傳統(tǒng)的模系矩陣分裂迭代法。

2 加速模系矩陣分裂迭代法

先給出本文討論中使用的一些記號及基本概念。

給定兩個實矩陣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 =( mij)定義為

若對任意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‐相容分裂。

的等價不動點方程，其對于建立加速模系矩陣分裂迭代算法是至關(guān)重要的。

定理1［16］令M = F ?G 為矩陣M 的一個分裂，給定初始向量x(0)，對于隱式互補問題（1），下列結(jié)論成立：

是隱式互補問題（1）的解。

基于隱式不動點方程（2），李郴良等［16］構(gòu)造了如下求解隱式互補問題的模系矩陣分裂迭代法。

算法1 （模系矩陣分裂迭代法）

步驟1 定義U ={ u：u?f ( u)≥0，Mu+ q≥0 }，令M = F ?G為矩陣M的一個分裂。

步驟2 給定ε > 0，u(0)∈U，k = 0。

步驟3 計算解：

1）計算初始向量

2）通過下列迭代格式計算x(j+1，k)：

步驟4 如果u(k+1)滿足停止準(zhǔn)則，則迭代停止。否則，令k = k + 1且轉(zhuǎn)向步驟2。

為了加快求解隱式互補問題的傳統(tǒng)模系矩陣分裂迭代法的收斂速度，通過引入對向量邊計算邊更替的思想，可構(gòu)造如下加速模系矩陣分裂迭代法。

算法2 （加速模系矩陣分裂迭代法）

步驟1 定義U ={ u：u?f ( u)≥0，Mu+ q≥0 }，令M = F1?G1= F2?G2為 矩 陣M 的 兩 個分裂。

步驟2 給定ε > 0，u(0)∈U，k = 0。

步驟3 計算解u(k+1)：

1）計算初始向量

步驟4 如果u(k+1)滿足停止準(zhǔn)則，則迭代停止；否則，令k = k + 1且轉(zhuǎn)向步驟2。

3 收斂性分析

本節(jié)建立了系數(shù)矩陣為H+‐矩陣情況下加速模系矩陣分裂迭代法求解隱式互補問題（1）的收斂性分析。

引理1［23］令M ∈Rn×n是一個H+‐矩陣，則|M?1|≤M?1。

引理2［24］令M ∈Rn×n，則ρ( M )< 1當(dāng)且僅當(dāng)

由定理1知，若(u*，r*)為隱式互補問題（1）的解，則

定理2 令矩陣M ∈Rn×n為H+‐矩陣，M =F1?G1= F2?G2為矩陣M的兩個H‐相容分裂。假設(shè)Ω ∈Rn×n為一個正對角矩陣，γ是一個正常數(shù)，且f (u)是一個李普希茲連續(xù)函數(shù)，即存在非負(fù)矩陣N滿足|f ( u)?f ( v )|≤N|u?v|，?u，v ∈Rn，定義

如果Ω ≥diag (F2)且Ω + F1?| G2|是一個M‐矩陣及

因為M = F1?G1= F2?G2是矩陣M 的兩個H‐相容分裂，則

由M是H+‐矩陣可得F1為H+‐矩陣，且由引理1知

類似地，式（3）減去式（6），且兩邊同時加上絕對值得

或者等價的有

因為Ω + F1?| G2|為M‐矩陣，則Ω + F1也為M‐矩陣。 因此Ω + F1?| G2|是M‐分裂，且

進(jìn)一步知I ?(Ω + F1)?1| G2|為M‐矩陣，且其逆為非負(fù)矩陣，從而

類似可得

定義Z?= ?j1+1( I +|Ω?1M|+ N )+ ?3N，則式（7）可以表示為

顯然，如果ρ( Z?)< 1，迭代序列{ u(k)}∞k=0收斂到u*。事實上

由文獻(xiàn)［21］可知當(dāng)M = F1?G1= F2?G2為矩陣M 的兩個H‐相容分裂，且Ω + F1?| G2|是一 個M 陣 時，可 得ρ(?1)< 1。 由 引 理2 可 得= 0。故存在一個正整數(shù)j0，當(dāng)j ≥j0時有

于是，對于?j ≥j0有

類似地，可以建立下列加速的模系矩陣加速超松弛迭代法的收斂理論。

推論1 設(shè)矩陣M 為H+‐矩陣，令M=D?L?U=D?B，這里D、?L和?U分別是矩陣M的對角部分，嚴(yán)格下三角部分及嚴(yán)格上三角部分。 假設(shè)Ω為一個正對角矩陣滿足Ω ≥D，γ 是一個正常數(shù)，若為M‐矩陣，松弛系數(shù)α 和β 滿足下列條件：

其中λ同定理2中定義。

則任意給定初始向量u(0)∈U，算法2生成的迭代序列{ u(k)}∞k=0收斂到隱式互補問題（1）的唯一解u*。

用式（10）減去式（9），可得

由條件可知Ω ≥D且0≤β ≤α，在式（11）兩邊同時加上絕對值可得

式（12）等價于

與定理2類似可得

由文獻(xiàn)［21］可知，當(dāng)α、β滿足條件（8）時，有ρ(Δ1)<1。由引理2 可得jl→im∞Δj1+1= 0。故存在一個正整數(shù)j0，當(dāng)j≥j0時

4 數(shù)值實驗

本節(jié)將給出一些數(shù)值算例來驗證加速模系矩陣分裂迭代法的有效性。數(shù)值實驗將從迭代步數(shù)（IT）、所消耗時間（CPU）和殘差范數(shù)（RES）3 個方面來反映算法的優(yōu)劣。其中殘差范數(shù)定義為

所有計算均在一臺CPU 為1. 8 GHz 和內(nèi)存為8 GB 的機(jī)器上運行，編程語言為MATLAB。在本實驗中，所有的初始向量取為：u(0)=(0，0，…，0，0，…)T∈Rn. 并且內(nèi)迭代的停止準(zhǔn)則為‖x(j+1，k)?x(j，k)‖≤10?4，外 迭 代 的 停 止 準(zhǔn) 則 為RES(u(k))≤10?6。α為加速模系矩陣逐步超松弛迭代法與模系矩陣逐步超松弛迭代法的迭代系數(shù)，最優(yōu)參數(shù)為使得這兩種迭代法迭代步數(shù)與時間最小的參數(shù)。

例1 設(shè)m為給定的正整數(shù)，n=m2。在式（1）中取M=M?+μI，q=?Mz?，其中

為塊三對角矩陣，其中S=tridiag（?1. 5，4，?0. 5）∈Rm×m為三對角矩陣，取z?=（1，2，1，2，…，1，2，…）T∈Rn和f（u）=（arctan（u1），…，arctan（un））。

例2 設(shè)m為給定的正整數(shù)，n=m2。在式（1）中取M=M?+μI，q=?Mz?其中

為塊三對角矩陣，其中S=tridiag（?1，4，?1）∈Rm×m為三對角矩陣，取z?=（1，2，1，2，…，1，2，

表1列出了例1和例2中使用的6種測試方法及其縮寫符號。表2 列出了μ= 2，α以0. 2 為間隔從1. 0 到2. 2 變化以及矩陣維數(shù)增加時，MSOR 及AMSOR算法迭代步數(shù)與迭代時間的數(shù)值結(jié)果。

表1 測試方法及其縮寫Tab. 1 Abbreviations of testing methods

從表2 中可以看到當(dāng)μ= 2，矩陣維數(shù)為400、1 600、3 600 和6 400 時，從最小化迭代時間考慮，MSOR 相對應(yīng)的最優(yōu)參數(shù)分別為2. 0、2. 0、2. 0 和1. 8，然而AMSOR 的最優(yōu)參數(shù)分別為1. 8、1. 6、1. 6和1. 4。

表3和表4中列出了當(dāng)μ分別為2和4，矩陣維數(shù)增加時在最優(yōu)參數(shù)下6 種方法的迭代步數(shù)、迭代時間以及殘差范數(shù)的比較結(jié)果。

表3 例1中最優(yōu)參數(shù)下6種方法的數(shù)值結(jié)果（μ=2）Tab. 3 Numerical results for six methods with optimal parameters for example 1（μ=2）

表2 例1中當(dāng)α變化時MSOR和AMSOR的迭代步數(shù)與時間（μ=2）Tab. 2 IT and CPU for MSOR and AMSOR with the change in α for example 1（μ=2）

表4 例1中最優(yōu)參數(shù)下6種方法的數(shù)值結(jié)果（μ=4）Tab. 4 Numerical results for six methods with optimal parameters for example 1（μ=4）

從表3 和圖1 中可以看到所有的方法都快速收斂，并且所有迭代法的迭代時間隨著矩陣維數(shù)的增加而增加。測試的6種算法中AMSOR 需要最少的迭代步數(shù)與時間，并且所有的加速模系矩陣分裂迭代法與相對應(yīng)的模系矩陣分裂迭代法相比迭代步數(shù)與迭代時間都更少。

圖1 例1兩種算法的殘差隨迭代步數(shù)變化的比較Fig. 1 Residual comparison of two algorithms for example 1

當(dāng)μ= 4時系數(shù)矩陣變得更加對角占優(yōu)。從表4可以看出，在相同的矩陣維數(shù)下，6 種方法的迭代步數(shù)與時間少于表3 中的結(jié)果，并且加速模系矩陣分裂迭代法仍優(yōu)于相對應(yīng)的模系矩陣分裂迭代法。

表5 列出了μ= 2 時，α以0. 2 為間隔從1. 0 到2. 2 變化以及矩陣維數(shù)增加時，MSOR 及AMSOR算法迭代步數(shù)與迭代時間的數(shù)值結(jié)果。

從表5 中可以看到，當(dāng)μ= 2，矩陣維數(shù)為400、1 600、3 600 和6 400 時，MSOR 相對應(yīng)的最優(yōu)參數(shù)分別為1. 8、1. 8、1. 6和1. 6，而AMSOR的最優(yōu)參數(shù)分別為1. 6、1. 6、1. 8和1. 8。

表6和表7中列出了當(dāng)μ分別為2和4，矩陣維數(shù)增大時在最優(yōu)參數(shù)下6 種方法的迭代步數(shù)、迭代時間以及殘差范數(shù)的比較結(jié)果。

表5 例2中當(dāng)α變化時MSOR和AMSOR的迭代步數(shù)與時間（μ=2）Tab. 5 IT and CPU for MSOR and AMSOR with the change in α for example 2（μ=2）

表7 例2中最優(yōu)參數(shù)下6種方法的數(shù)值結(jié)果（μ=4）Tab. 7 Numerical results for six methods with optimal parameters for example 2（μ=4）

表6 例2中最優(yōu)參數(shù)下6種方法的數(shù)值結(jié)果（μ=2）Tab. 6 Numerical results for six methods with optimal parameters for example 2（μ=2）

表6 和表7 中的數(shù)值結(jié)果進(jìn)一步驗證了從表3和表4得到的結(jié)論，并且6種方法對于對稱矩陣的迭代步數(shù)與時間略多于非對稱矩陣。測試的6種算法中AMSOR 算法的計算效率優(yōu)于其他5 種方法，并且所有的加速模系矩陣分裂迭代法仍優(yōu)于相對應(yīng)的模系矩陣分裂迭代法。

5 結(jié)論

本文構(gòu)造了求解一類隱式互補問題的加速模系矩陣分裂迭代法。理論分析建立了系數(shù)矩陣為H+‐矩陣時的收斂性。數(shù)值結(jié)果表明，加速模系矩陣分裂迭代法在迭代步數(shù)和迭代時間上均優(yōu)于傳統(tǒng)的模系矩陣分裂迭代法。

然而對于隱式互補問題的求解，無論是模系矩陣分裂迭代法還是加速的模系矩陣分裂迭代法，內(nèi)迭代的收斂速度在整個迭代過程中都有著十分重要的影響。合適的內(nèi)迭代停止準(zhǔn)則或迭代步數(shù)可以加速求解隱式互補問題的計算效率。目前，對于內(nèi)迭代最優(yōu)的停止準(zhǔn)則以及停止的迭代步數(shù)還沒有很好的結(jié)論，并且內(nèi)迭代與外迭代之間如何平衡也是一個十分有挑戰(zhàn)性的課題，這需要筆者進(jìn)一步的研究。