吳思婷，鮑 亮

(華東理工大學(xué)數(shù)學(xué)學(xué)院,上海 200237)

1 引言

大型稀疏線性方程組在許多科學(xué)與工程領(lǐng)域都有著非常重要的作用，比如在最優(yōu)化問題、流體力學(xué)和電磁學(xué)等領(lǐng)域都需要對此類問題進(jìn)行求解[1-7]。因此，如何快速有效地求解大型稀疏線性方程組已經(jīng)成為當(dāng)下非常重要的課題之一。

本文考慮如式(1)所示的大型稀疏正定線性方程組的求解：

Ax=b

(1)

其中,A∈Cn×n為非Hermitian正定矩陣，x∈Cn為未知向量，b∈Cn為已知向量。當(dāng)方程組(1)的系數(shù)矩陣規(guī)模較小時(shí)，使用直接法更加方便。而對于較大規(guī)模的稀疏非Hermitian正定線性方程組，研究人員一般采用迭代法進(jìn)行求解。

在求解某些特定問題時(shí)，求解系數(shù)矩陣αI+S的線性方程組可能會出現(xiàn)一些問題，甚至可能會同求解原始的線性方程組Ax=b一樣困難。而選擇恰當(dāng)?shù)木仃嘖會使得矩陣αI+K+S比矩陣αI+S更對角占優(yōu)，因此求解系數(shù)矩陣αI+K+S的線性方程組會比求解系數(shù)矩陣αI+S的線性方程組更加方便?；谏鲜霭l(fā)現(xiàn)，本文將系數(shù)矩陣進(jìn)行廣義的HS分裂，再進(jìn)行非對稱的二步迭代。這可以將二步迭代中求解系數(shù)矩陣αI+S的線性方程組簡化為求解系數(shù)矩陣αI+K+S的線性方程組，大大加快了迭代方法的收斂速度。

本文將系數(shù)矩陣進(jìn)行廣義的Hermitian和反Hermitian分裂，再通過引入新的變量同時(shí)結(jié)合外推的技術(shù)，給出一種外推的廣義Hermitian和反Hermitian迭代方法EGHSS(Extrapolated and Generalized HSS)。還理論分析了本文方法的收斂性，給出了該方法收斂的充要條件，并將該方法迭代矩陣的譜半徑與GHSS方法迭代矩陣的譜半徑和EHSS方法迭代矩陣的譜半徑進(jìn)行比較。數(shù)值實(shí)驗(yàn)結(jié)果表明，在處理某些實(shí)際問題中，EGHSS迭代方法比GHSS迭代方法和EHSS迭代方法更有效，且在合適的參數(shù)值下新方法的收斂效率可以大大提高。

2 EGHSS迭代方法

首先簡要回顧EHSS迭代方法[13]。

方法1EHSS迭代方法。

(2)

其中,α>0和ω>0均為給定的常數(shù)。

事實(shí)上，可以將式(2)改寫成如式(3)～式(5)所示的等價(jià)形式：

x(k+1)=M1(α,ω)x(k)+N1(α,ω)b,

k=0,1,2,…

(3)

[HS-(1-ω)α(H+S)+α2I]

(4)

(5)

這里的M1(α,ω)是EHSS迭代方法的迭代矩陣。

考慮到上述迭代方法中求解系數(shù)矩陣為αI+S的線性方程組仍有一定的難度，為了將其簡化為求解系數(shù)矩陣αI+K+S的線性方程組，本文將HS分裂中的Hermitian矩陣H作如式(6)所示的分裂：

H=G+K

(6)

其中，G和K為Hermitian半正定矩陣，從而得到系數(shù)矩陣A∈Cn×n有如式(7)所示廣義的HS分裂：

A=G+K+S

(7)

其中，G和K為Hermitian半正定矩陣，S為反Hermitian矩陣。一種常見的矩陣分裂格式如式(8)所示：

(8)

基于上述系數(shù)矩陣的廣義HS分裂，下面介紹EGHSS迭代方法。

方法2EGHSS迭代方法。

設(shè)A∈Cn×n為非Hermitian正定矩陣，給定一個(gè)初始向量x(0)∈Cn，對于k=0,1,2,…, 直到迭代序列{x(k)}收斂，計(jì)算如式(9)所示：

(9)

其中,α>0和ω>0均為給定的常數(shù)。

當(dāng)ω=0時(shí)，EGHSS迭代方法就成為GHSS迭代方法；當(dāng)K=0時(shí)，EGHSS迭代方法就成為EHSS迭代方法；當(dāng)ω=0且K=0時(shí)，EGHSS迭代方法就成為HSS迭代方法。

觀察EHSS方法的迭代格式，不難發(fā)現(xiàn)式(9)等價(jià)式(10)～式(12)所示的矩陣-向量形式：

x(k+1)=M(α,ω)x(k)+N(α,ω)b,k=0,1,2,…

(10)

M(α,ω)=(αI+K+S)-1[K+S-(1-ω)αI+

(2-ω)α(αI+G)-1(αI-K-S)]=

(αI+K+S)-1(αI+G)-1{(αI+G)[K+

S-(1-ω)αI]+(2-ω)α(αI-K-S)}=

(αI+K+S)-1(αI+G)-1{G(K+S)-

(1-ω)α(G+K+S)+α2I}

(11)

(12)

這里的M(α,ω)為EGHSS迭代方法的迭代矩陣。

當(dāng)然EGHSS迭代方法也可由系數(shù)矩陣A經(jīng)過如式(13)～式(15)所示的分裂得到：

A=B(α,ω)-C(α,ω)

(13)

(14)

(1-ω)α(G+K+S)+α2I]=

(15)

3 收斂性分析

本節(jié)討論了EGHSS迭代方法的收斂性，并給出了迭代方法收斂的充要條件。

引理1[12]設(shè)A∈Cn×n為正定矩陣，令A(yù)=G+K+S,G和K為Hermitian半正定矩陣，S為反Hermitian矩陣，M(α)為GHSS迭代方法的迭代矩陣。如果G或K為正定矩陣，且α>0，則迭代矩陣M(α)的譜半徑ρ(M(α))滿足式(16)：

(16)

即GHSS迭代方法收斂到線性方程組式(1)的精確解x*∈Cn。

從上述EGHSS方法的迭代格式推導(dǎo)過程不難發(fā)現(xiàn)，EGHSS方法的迭代矩陣與GHSS方法的迭代矩陣存在一定的關(guān)系。為了證明EGHSS迭代方法的收斂性，本節(jié)根據(jù)EGHSS迭代方法的迭代矩陣M(α,ω)的形式，進(jìn)一步分析了EGHSS迭代方法的迭代矩陣M(α,ω)與GHSS迭代方法的迭代矩陣M(α)之間的關(guān)系。

定理1設(shè)A∈Cn×n為非Hermitian正定矩陣，令A(yù)=G+K+S,G和K為Hermitian半正定矩陣，S為反Hermitian矩陣，對任意的α>0和ω≥0，GHSS迭代方法的迭代矩陣如式(17)所示：

(17)

那么EGHSS迭代方法的迭代矩陣M(α,ω)滿足式(18)：

(18)

證明

ωI+(2-ω)M(α)=

ωI+(2-ω)(αI+K+S)-1(αI-G)

(αI+G)-1(αI-K-S)=

(αI+K+S)-1(αI+G)-1[ω(αI+G)

(αI+K+S)+(2-ω)(αI+G)(αI-G)(αI+G)-1

(αI-K-S)]=(αI+K+S)-1(αI+G)-1

[ω(αI+G)(αI+K+S)+

(2-ω)(αI-G)(αI+G)(αI+G)-1

(αI-K-S)]=(αI+K+S)-1(αI+G)-1

[ω(α2I+αG+αK+αS+GK+GS)+(2-ω)

(αI-G)(αI-K-S)]=(αI+K+S)-1(αI+G)-1

[2GK+2GS+2α2I-2(1-ω)α(G+K+S)]

從而可以得到式(19)：

[G(K+S)-(1-ω)α(G+K+S)+α2I]=

(19)

□

由引理1和定理1可以得出EGHSS迭代方法收斂的充要條件。

定理2[14,15]設(shè)A∈Cn×n為非Hermitian正定矩陣，令A(yù)=G+K+S,G和K為Hermitian半正定矩陣，S為反Hermitian矩陣，如果0≤ω<2，那么對任意的α>0，EGHSS方法的迭代格式式(9)收斂。

證明由定理1知EGHSS迭代方法的迭代矩陣M(α,ω)與GHSS迭代方法的迭代矩陣M(α)之間存在如式(20)所示的關(guān)系：

(20)

設(shè)λ和μ分別為GHSS迭代方法的迭代矩陣M(α)和EGHSS迭代方法的迭代矩陣M(α,ω)的特征值，則有式(21):

(21)

若λ為實(shí)數(shù)，則式(22)成立:

(22)

若λ為復(fù)數(shù)，設(shè)λ=a+bi,i為虛數(shù)單位，則式(23)成立:

(23)

綜上所述,式(24)成立：

(24)

又由引理1知?α>0，ρ(M(α))<1,所以當(dāng)0≤ω<2時(shí)，即0<2-ω≤2時(shí)，有式(25):

(25)

因此,EGHSS迭代方法收斂。

□

由定理2可知，EGHSS迭代方法是條件收斂的。由定理2的證明過程可知，EGHSS方法的迭代矩陣M(α,ω)的特征值與GHSS方法的迭代矩陣M(α)的特征值存在著一定的聯(lián)系。接下來進(jìn)一步分析EGHSS方法的迭代矩陣M(α,ω)的譜半徑與GHSS方法的迭代矩陣M(α)譜半徑之間的關(guān)系。

定理3設(shè)A∈Cn×n為非Hermitian正定矩陣，令A(yù)=G+K+S,G和K為Hermitian半正定矩陣，S為反Hermitian矩陣，λ=a+bi為GHSS方法的迭代矩陣M(α)的特征值，則有：

證明根據(jù)定理2可得式(26):

(26)

(27)

(28)

此時(shí)|μ|<|λ|，于是當(dāng)a2+b2<1時(shí)，有0≤ω<2，則式(29)成立:

ρ(M(α,ω))≤ρ(M(α))<1

(29)

ρ(M(α))≤ρ(M(α,ω))

(30)

ρ(M(α))≤ρ(M(α,ω))<1

(31)

□

根據(jù)定理3可知，通過選擇合適的參數(shù)，EGHSS方法的迭代矩陣M(α,ω)的譜半徑會小于GHSS方法的迭代矩陣M(α)的譜半徑。下一節(jié)將通過具體的數(shù)值實(shí)驗(yàn)對EGHSS迭代方法、GHSS迭代方法和EHSS迭代方法的收斂速度和迭代次數(shù)進(jìn)行比較，并詳細(xì)分析矩陣規(guī)模、α和ω對EGHSS迭代方法收斂速度的影響。

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

本節(jié)對EGHSS迭代方法進(jìn)行了數(shù)值實(shí)驗(yàn)，并將其結(jié)果與GHSS迭代方法和EHSS迭代方法進(jìn)行比較。本次數(shù)值實(shí)驗(yàn)是在英特爾四核處理器(2.40 GHz，8 GB RAM)環(huán)境下應(yīng)用Matlab編程語言實(shí)現(xiàn)的。

例1考慮文獻(xiàn)[16]中的數(shù)值例子,如式(32)所示：

-u″+qu′=f

(32)

邊界條件為齊次邊界條件。根據(jù)中心差分方法對式(32)進(jìn)行離散，可以得到如式(33)所示的系數(shù)矩陣：

A=tridiag(-1+qh/2,2,-1-qh/2)∈RN×N

(33)

其中，tridiag(·,·,·)表示三對角矩陣。

例2考慮區(qū)域?yàn)棣?[0,1]×[0,1]×[0,1]的一類三維對流擴(kuò)散方程[8],如式(34)所示：

-(uxx+uyy+uzz)+q(ux+uy+uz)=

f(x,y,z)

(34)

其中,q=1000為給定常數(shù)，且區(qū)域邊界滿足Dirichlet邊界條件。下面采用向前差分離散方程得到如式(35)所示的系數(shù)矩陣：

A=Tx?I?I+I?Ty?I+

I?I?Tz∈RN3×N3

(35)

4.1 收斂速度分析

表1～表3分別給出了例1中GHSS迭代方法、EHSS迭代方法和EGHSS迭代方法的收斂迭代次數(shù)和收斂時(shí)間，并且給出了3種迭代方法近似最優(yōu)變量的取值。其中N的取值表示的矩陣規(guī)模分別為256×256階，512×512階,1024×1024階，2048×2048階。從表1～表3可以發(fā)現(xiàn)，當(dāng)數(shù)值例子相同且矩陣規(guī)模相同時(shí)，EGHSS迭代方法的譜半徑、迭代次數(shù)和迭代時(shí)間普遍都小于GHSS迭代方法的和EHSS迭代方法的。且隨著qh和N的不斷增大，EGHSS迭代方法的迭代矩陣譜半徑、迭代次數(shù)和迭代時(shí)間都遠(yuǎn)遠(yuǎn)小于GHSS迭代方法的和EHSS迭代方法的。因此，不難從表1～表3中發(fā)現(xiàn)，在合適參數(shù)下EGHSS迭代方法會相比GHSS迭代方法和EHSS迭代方法有更好的效果。

表4給出了例2中GHSS迭代方法、EHSS迭代方法和EGHSS迭代方法的收斂迭代次數(shù)和收斂時(shí)間，并且給出了3種迭代方法近似最優(yōu)變量的取值。其中N的取值表示的矩陣規(guī)模為1728×1728階和4096×4096階。從表4可以發(fā)現(xiàn)，當(dāng)矩陣規(guī)模相同時(shí)，EGHSS迭代方法的迭代速度比GHSS迭代方法的和EHSS迭代方法的更快。甚至當(dāng)?shù)螖?shù)相近時(shí)，EGHSS迭代方法的迭代時(shí)間明顯小于GHSS迭代方法的。因此不難發(fā)現(xiàn)，在選擇合適參數(shù)下EGHSS迭代方法會相比GHSS迭代方法和EHSS迭代方法有更好的效果。

Table 1 IT and CPU of GHSS iterative method for example 1

Table 2 IT and CPU of EHSS iterative method for example 1

Table 3 IT and CPU of EGHSS iterative method for example 1

Table 4 IT and CPU of GHSS,EGHSS and EHSS iterative method for example 2

圖1展示了例1取qh=10時(shí)，系數(shù)矩陣規(guī)模分別為256×256階和512×512階的情況下，EGHSS迭代方法、GHSS迭代方法和EHSS迭代方法的殘量隨著迭代步數(shù)的變化趨勢。顯然EGHSS迭代方法的殘量下降曲線位于GHSS迭代方法和EHSS迭代方法的殘量下降曲線之下。這說明對上述qh=10的數(shù)值例子，EGHSS迭代方法相較于GHSS迭代方法和EHSS迭代方法迭代次數(shù)更少，迭代時(shí)間更短。并且系數(shù)矩陣規(guī)模越大，EGHSS迭代方法相較于其他2種方法的優(yōu)勢越明顯。因此，在處理某些實(shí)際問題中，EGHSS迭代方法比GHSS迭代方法和EHSS迭代方法具有迭代次數(shù)更少、迭代時(shí)間更短等優(yōu)點(diǎn)。

4.2 參數(shù)α和ω的靈敏度分析

通常來說，迭代矩陣的譜半徑越小，迭代方法的收斂速度越快。因此，本節(jié)主要討論各個(gè)變量對EGHSS方法的迭代矩陣譜半徑的影響。首先分析單一變量α對EGHSS迭代方法、GHSS迭代方法和EHSS迭代方法的迭代矩陣譜半徑的影響；再分析2個(gè)變量α和ω共同對EGHSS迭代方法的譜半徑的影響。本節(jié)參數(shù)靈敏度分析僅考慮例1中的參數(shù)α和ω對EHSS迭代方法譜半徑的影響。

圖2a表示例1選取矩陣規(guī)模為256×256階，且qh=10和ω=0.6時(shí)，EGHSS迭代方法、GHSS迭代方法和EHSS迭代方法中參數(shù)α與迭代矩陣譜半徑的關(guān)系。圖2b表示例1選取矩陣規(guī)模為512×512階，且qh=10,ω=0.5時(shí)，EGHSS迭代方法、GHSS迭代方法和EHSS迭代方法中參數(shù)α與迭代矩陣譜半徑的關(guān)系。觀察圖2可以發(fā)現(xiàn)，3種方法迭代矩陣的譜半徑均隨著α增大先減小后增大。除此之外,EGHSS迭代方法的最小譜半徑明顯小于GHSS迭代方法和EHSS迭代方法的最小譜半徑(都小于1)，這說明了EGHSS迭代方法相比GHSS迭代方法和EHSS迭代方法有更好的迭代效果。

圖3a表示例1選取系數(shù)矩陣規(guī)模為256×256階且qh=10時(shí)，EGHSS迭代方法中的參數(shù)α和ω與迭代矩陣譜半徑的關(guān)系。圖3b表示例1選取系數(shù)矩陣規(guī)模為512×512階且qh=10時(shí)，EGHSS迭代方法中的參數(shù)α和ω與迭代矩陣譜半徑的關(guān)系。觀察圖3可以發(fā)現(xiàn)，當(dāng)α和ω處于一個(gè)適當(dāng)?shù)娜≈捣秶鷷r(shí)，EGHSS迭代方法的譜半徑小于1，符合定理2的收斂性證明。當(dāng)系數(shù)矩陣規(guī)模為256×256階且qh=10時(shí)，α=1.6和ω=0.6時(shí),EGHSS迭代方法的譜半徑近似取得最小值(小于1)。當(dāng)系數(shù)矩陣規(guī)模為512×512階且qh=10時(shí)，α=1.1和ω=0.5時(shí),EGHSS迭代方法的譜半徑近似取得最小值(小于1)。因此,圖3a和圖3b分別選擇ω=0.6和ω=0.5是合理的。

5 結(jié)束語

本文提出了一種求解大型稀疏非Hermitian正定線性方程組的EGHSS迭代方法，給出了該方法收斂的充要條件，并進(jìn)行了數(shù)值實(shí)驗(yàn)。在與GHSS迭代方法和EHSS迭代方法的譜半徑、迭代步數(shù)及迭代時(shí)間的比較過程中不難發(fā)現(xiàn)，EGHSS迭代方法在處理某些問題時(shí)相較于其他2種迭代方法具有一定的優(yōu)越性。除此之外，本文還分析了單一參數(shù)α對EGHSS方法迭代矩陣譜半徑的影響，以及2個(gè)參數(shù)α和ω共同對EGHSS方法迭代矩陣譜半徑的影響。當(dāng)然，理論上的最優(yōu)參數(shù)α和ω的確定還有待進(jìn)一步的研究?？偟膩碚f，EGHSS迭代方法是一種較GHSS迭代方法和EHSS迭代方法更有競爭力的方法，特別是選取了合適的參數(shù)α和ω后，EGHSS方法的收斂速度大大提高。