李 劍,張 娜,李瑞霞

(陜西科技大學(xué) 數(shù)學(xué)與數(shù)據(jù)科學(xué)學(xué)院, 陜西 西安 710021)

0 引言

考慮以下3×3塊大型稀疏線性系統(tǒng)的快速求解:

(1)

式中:矩陣A∈Rn×n是對(duì)稱正定的,矩陣B∈Rm×n和C∈Rp×m都是行滿秩的,u∈Rn、v∈Rm、w∈Rp是給定的向量。這些假設(shè)保證了線性系統(tǒng)(1)的解的存在唯一性[1]。 通過以下兩種形式的劃分:

(2)

和

(3)

可以發(fā)現(xiàn),3×3塊線性系統(tǒng)(2)可以歸結(jié)為經(jīng)典鞍點(diǎn)問題[1]。而由于線性系統(tǒng)(3)中的(2,2)塊是反對(duì)稱矩陣,而(1,2)塊是秩虧的,使其不同于傳統(tǒng)鞍點(diǎn)問題。對(duì)于具有如此特殊分塊的3×3塊線性系統(tǒng)(1),由于其形似鞍點(diǎn)或雙鞍點(diǎn)的結(jié)構(gòu),故也稱為雙鞍點(diǎn)問題。它廣泛來源于工程應(yīng)用和科學(xué)計(jì)算中的大多數(shù)問題,如約束二次規(guī)劃[2]、橢圓偏微分方程的混合有限元方法[3]、計(jì)算流體動(dòng)力學(xué)[4-5]以及最小二乘問題[6]等。

近年來,為了提高Krylov子空間方法求解傳統(tǒng)鞍點(diǎn)系統(tǒng)的收斂速度,已經(jīng)有很多預(yù)處理子被學(xué)者們陸續(xù)提出,如HSS預(yù)處理子及其變形[7-8]、移位分裂預(yù)處理子[9]、結(jié)構(gòu)預(yù)處理子[10]等。但是考慮到結(jié)構(gòu)上存在的差異,傳統(tǒng)求解鞍點(diǎn)問題的預(yù)處理子并不能直接應(yīng)用于雙鞍點(diǎn)問題,因此尋找新的有效預(yù)處理形式是非常有必要的。

最近,針對(duì)雙鞍點(diǎn)系統(tǒng)(1)的快速求解,HUANG等[11]提出了兩個(gè)塊對(duì)角(Block-Diagonal,簡稱BD)預(yù)處理子:

文獻(xiàn)[12]分析了相應(yīng)預(yù)處理矩陣的譜性質(zhì),表明這些預(yù)處理子會(huì)顯著加快GMRES方法的收斂速度。但當(dāng)內(nèi)部系統(tǒng)被不精確求解時(shí),計(jì)算時(shí)間會(huì)急劇增加。ABDOLMALEKI等[13]提出以下形式的塊對(duì)角預(yù)處理子:

與HUANG等[11]提出的塊對(duì)角(BD)預(yù)處理子相比,該預(yù)處理子的主要優(yōu)點(diǎn)是不涉及Schur補(bǔ)矩陣,所以計(jì)算過程易于實(shí)施。作者從3×3分塊矩陣的本征系統(tǒng)得到了一個(gè)三次方程,給出了相應(yīng)預(yù)處理矩陣的特征值的分布范圍。

此外,CAO[14]基于矩陣分裂,得到了一類移位分裂(Shift-Splitting,簡稱SS)預(yù)處理子:

式中:α>0是給定的常數(shù),I是具有適當(dāng)維數(shù)的單位矩陣。作者證明了相應(yīng)迭代法的無條件收斂性,并給出了相應(yīng)預(yù)處理矩陣的特征值分布范圍。與文獻(xiàn)[11]提出的塊對(duì)角預(yù)處理子相比,移位分裂預(yù)處理子在計(jì)算效率上表現(xiàn)出良好的性能。

總的來說,雖然塊對(duì)角預(yù)處理子M不涉及Schur補(bǔ)矩陣的計(jì)算,但與系數(shù)矩陣A的近似程度不高。而移位分裂預(yù)處理子PSS在計(jì)算過程中雖然與系數(shù)矩陣A的近似程度較高,但又涉及Schur補(bǔ)矩陣的求解。因此,構(gòu)造不涉及Schur補(bǔ)矩陣的計(jì)算且與系數(shù)矩陣的近似程較高的預(yù)處理子,對(duì)于快速求解線性系統(tǒng)(1)是非常有意義的。

本文基于塊對(duì)角預(yù)處理子M的構(gòu)造特點(diǎn),并結(jié)合移位分裂的思想,提出一種新的不涉及Schur補(bǔ)計(jì)算的參數(shù)化塊移位分裂預(yù)處理子,用以求解線性方程組(1)。該預(yù)處理子具有易于實(shí)現(xiàn)且計(jì)算效率高等特點(diǎn)。此外,給出預(yù)處理矩陣特征值的上下界,并通過數(shù)值算例驗(yàn)證新預(yù)處理子的有效性。

1 一類新的參數(shù)化塊移位分裂預(yù)處理子

將系統(tǒng)(1)中的系數(shù)矩陣A0進(jìn)行如下分裂:

式中:α、β>0為迭代參數(shù)。該分裂可導(dǎo)出如下形式的塊移位分裂預(yù)處理子:

(4)

由此,預(yù)處理子P的算法實(shí)現(xiàn)過程如下:

算法1 求解Ps=r

1.求解As1=r1;

2.求解(αI+βBBT+CTC/α)s2=r2+CTr3/α;

3.計(jì)算s3=(r3-Cs2)/α。

算法1表明,利用預(yù)處理子P加速Krylov子空間方法來求解線性系統(tǒng)(1)時(shí),只涉及以矩陣A和αI+βBBT+(1/α)CTC為系數(shù)矩陣的子線性系統(tǒng)的求解,不涉及顯式Schur補(bǔ)的運(yùn)算,可期望達(dá)到較高的計(jì)算效率。一般來說,預(yù)處理矩陣的譜分布越聚集,Krylov子空間方法的收斂速度就越快,下面給出預(yù)處理矩陣的譜分析。

引理1[15]設(shè)

是一個(gè)復(fù)系數(shù)一元多項(xiàng)式,λ是p(z)的任意根,那么|λ|滿足以下不等式(柯西上下界):

|λ|≤max{|a0|,1+|a1|,…,1+|an-1|}。

定理1 設(shè)A∈Rn×n是對(duì)稱正定矩陣,B∈Rm×n和C∈Rp×m都是行滿秩矩陣,λ為預(yù)處理矩陣P-1A的特征值,[u*,v*,w*]*為相應(yīng)的特征向量,則預(yù)處理矩陣P-1A有n個(gè)特征值為1,剩余的特征值分布范圍如下:

(5)

式中:

λmin(BBT)≤b≤λmax(BBT),

λmin(CTC)≤c≤λmax(CTC)。

證明考慮廣義特征值方程P-1Az=λz,即

式中:λ是P-1A的特征值,z=[u*,v*,w*]*是相應(yīng)的特征向量,通過計(jì)算得到

(λ-1)Au-BTv=0,

(6)

Bu+λ(αI+βBBT)v+(1-λ)CTw=0,

(7)

(λ-1)Cv+αλw=0。

(8)

由于預(yù)處理矩陣P是非奇異的,故λ≠0。首先令u=0,則BTv=0,由B是行滿秩矩陣可知,v=0。又因α≠0,由式(7)可得w=0,這與[u*,v*,w*]*是特征向量相矛盾。接著令u≠0,v=0,則由式(6)和式(8)可以推斷出

(λ-1)Au=0,αλw=0,

這表明λ=1,其相應(yīng)的特征向量的形式為[u*,0,0]*且u∈null(B),則預(yù)處理矩陣P-1A有n個(gè)特征值為1。故令u≠0,v≠0,由式(6)可知λ≠1。而A是非奇異的,由式(6)得

由式(8)得

代入式(7)得

(9)

等式(9)兩邊同時(shí)乘αλ(λ-1)v*/(v*v),可得

進(jìn)一步計(jì)算可得到以下關(guān)于λ的一個(gè)三次方程:

λ3(α2+αβb+c)-λ2(α2+αβb+3c)+λ(αa+3c)-c=0,

(10)

式中:

(11)

對(duì)于任意具有實(shí)特征值的方陣Φ,其最小和最大特征值分別用λmin(Φ)和λmax(Φ)表示,由式(11)以及文獻(xiàn)[16]可得a、b、c的上下界:

λmin(BBT)≤b≤λmax(BBT),

λmin(CTC)≤c≤λmax(CTC)。

為了得到點(diǎn)(1,0)附近特征值的聚類性質(zhì),令μ=λ-1,式(10)可改寫為以下形式:

根據(jù)引理1,可以得到μ的柯西下界和上界,進(jìn)而得到λ的柯西下界和上界(式(5))。證畢。

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

本節(jié)通過數(shù)值算例比較新提出的預(yù)處理子與已有預(yù)處理子的有效性。所有的數(shù)值實(shí)驗(yàn)都是在裝有Intel(R) Core(TM) i7-8700 CPU @3.20 GHz 3.19 GHz 8.00 GB的電腦上使用Matlab R2019b進(jìn)行計(jì)算。給出數(shù)值例子來分析新的塊移位分裂預(yù)處理子的性能,z(0)的初始值設(shè)置為零向量,并在最大迭代次數(shù)kmax=1000時(shí)停止;或滿足相對(duì)殘差限時(shí)停止,

例1考慮3×3塊鞍點(diǎn)問題(1),其中

B=(I?FF?I)∈Rp2×2p2,

C=E?F∈Rp2×p2,

且

E=diag(1,p+1,…,p2-p+1),

?表示Kronecker張量積,h=1/(p+1)為網(wǎng)格尺寸,總維數(shù)為4p2,右端向量b=(1,1,…,1)T∈R4p2。

數(shù)值實(shí)驗(yàn)是通過設(shè)置6個(gè)遞增的剖分p=16、32、64、128、256、512生成的。表1給出了6種預(yù)處理子加速GMRES的數(shù)值結(jié)果,其中IT表示迭代次數(shù)、CPU表示計(jì)算時(shí)間、RES表示相對(duì)殘差,符號(hào)“—”表示超過最大迭代次數(shù)或者超出存儲(chǔ)需求。

表1 不同預(yù)處理子的數(shù)值結(jié)果Tab. 1 Numerical results of different preconditioners

對(duì)于預(yù)處理子M,設(shè)置α=0.001和β=1。對(duì)于預(yù)處理子PSS,設(shè)置α=0.01。而設(shè)置α=0.000 01和β=0.001作為本文提出的預(yù)處理子P的實(shí)驗(yàn)最佳選擇。從表1可以看出,在迭代次數(shù)方面,新提出的預(yù)處理子P優(yōu)于PBD2預(yù)處理子和M預(yù)處理子;在計(jì)算時(shí)間方面,新提出的預(yù)處理子P優(yōu)于其他預(yù)處理子。特別地,在剖分p=256和p=512時(shí),計(jì)算時(shí)間仍然有較優(yōu)的結(jié)果,計(jì)算效率高于目前已有的預(yù)處理子。

圖1分別繪制了當(dāng)p=16時(shí)系數(shù)矩陣A0在PBD2、P2、P3、M、PSS、P作用下的預(yù)處理矩陣的特征值分布?？梢钥闯?新提出的預(yù)處理矩陣P的特征值具有較好的聚集性。

圖1 不同預(yù)處理子的特征值分布圖(p=16)Fig. 1 Eigenvalue distribution of different preconditioners (p=16)

3 結(jié)束語

針對(duì)一類3×3塊鞍點(diǎn)問題,提出了一種新的不涉及Schur補(bǔ)的參數(shù)化塊移位分裂預(yù)處理子。與現(xiàn)有的預(yù)處理子相比,該預(yù)處理子易于實(shí)現(xiàn)且具有計(jì)算上的優(yōu)勢(shì)。文中討論了相應(yīng)預(yù)處理矩陣非零特征值的分布范圍及特征值1的代數(shù)重?cái)?shù)。數(shù)值結(jié)果表明,所提出的預(yù)處理子是非常有效的,特別是選取了合適的參數(shù)值之后,可以很大程度上提高Krylov子空間方法的收斂速度。