朱景福，李欣，郭東山

(廣東石油化工學(xué)院 理學(xué)院，廣東 茂名 525000)

本文研究多右端對(duì)稱線性方程組AX=B，其中A∈Rn×n為非奇異對(duì)稱矩陣；X,B∈Rn×s的數(shù)值求解方法。

多右端線性方程組在許多領(lǐng)域都有著重要的應(yīng)用背景，如控制論、結(jié)構(gòu)力學(xué)、電磁場等領(lǐng)域，都涉及多右端線性方程組的求解。多右端線性方程組的求解研究主要有精確解的研究和數(shù)值解的研究。在數(shù)值解法中常見的迭代法有種子投影方法[1,2]和塊迭代法。塊迭代法主要有塊CG方法(BCG)[3]、塊GMRES方法(BGMRES)[4]及其變形、塊QMR方法(BQMR)[5]、塊EN方法(BEN)[6]、塊Lanczos方法[7]等，這些塊方法都是建立在塊krylov方法[8](krylov子空間方法[9]的推廣)的基礎(chǔ)之上的。在利用krylov子空間方法求解線性方程組AX=b(即AX=B的右端項(xiàng)矩陣B只有1列)時(shí)，通常用殘量范數(shù)作為判斷算法終止的條件。若近似解的精確度比較高，則殘量范數(shù)較小，反之不一定。為了克服用殘量作為算法終止條件的不足，采用極小向后擾動(dòng)范數(shù)代替殘量作為算法終止的條件，Kasenally[10]給出了求解大型非對(duì)稱線性方程組的廣義極小向后擾動(dòng)方法。Cao[11]推廣了廣義極小向后擾動(dòng)方法，給出了求解大型非對(duì)稱線性方程組的總體GMBACK方法。

本文研究多右端對(duì)稱線性方程組AX=B的右端項(xiàng)矩陣B的列數(shù)s>1時(shí)的求解方法。新方法在采用塊Lanczos方法的過程中結(jié)合極小向后擾動(dòng)技術(shù)，討論AX=B的左端矩陣A的向后擾動(dòng)形式和極小向后擾動(dòng)范數(shù)的求法，將極小向后擾動(dòng)范數(shù)作為迭代終止的判定條件，提出求解多右端對(duì)稱線性方程組的一種新算法，并對(duì)新算法做殘量分析，通過數(shù)值對(duì)比實(shí)驗(yàn)驗(yàn)證新方法的有效性。

1 預(yù)備知識(shí)

首先給出本文所用到的符號(hào)約定如下：本文涉及的矩陣A均為n階實(shí)對(duì)稱矩陣，λ1,λ2,…,λn為A的特征值，x1,x2,…,xn為相應(yīng)的標(biāo)準(zhǔn)正交特征向量，并且λ1≥λ2≥…≥λn。En表示n階單位矩陣，在不產(chǎn)生混淆的情況下用E表示?！珹‖F(xiàn)、‖A‖2分別表示矩陣A的Frobenius范數(shù)和2-范數(shù)。xT和XT分別表示向量x和矩陣X的轉(zhuǎn)置。νec(A)表示矩陣A的列拉直(列展開)。A?B表示A與B的Kronecker積(直積，張量積)。A+表示矩陣A的Moore-Penrose廣義逆。

引理1 設(shè)A∈Cm×n，則‖A‖F(xiàn)=‖νec(A)‖2。

引理2設(shè)A∈Cm×n，B∈Cn×p，C∈Cp×q則νec(ABC)=(CT?A)νec(B)[12]。

引理4設(shè)A∈Cm×n，B∈Cp×q，C∈Cm×q，矩陣方程AXB=C有解的充要條件為AA+CB+B=C，在有解的情況下，AXB=C的通解為X=A+CB++Y-A+AYBB+

(1)

其中Y∈Cn×p是任意的；如果AXB=C無解，則AXB=C的最小二乘解仍為式(1)的形式，并且AXB=C的極小最小二乘解為X=A+CB+[12]。

選取方程組AX=B的初始矩陣X0，并由QR分解得到X0=Q1T1。

塊Lanczos方法(BLanczos方法)過程如下[7]：

可得

(2)

2 極小向后擾動(dòng)塊方法

2.1 向后擾動(dòng)分析

本文把方程組AX=B的近似解看作向后擾動(dòng)方程組(A-Δ)X=B的精確解，其中Δ稱為AX=B的向后擾動(dòng)。下面討論在用塊Lanczos方法解方程組AX=B時(shí)，考慮把極小化向后擾動(dòng)范數(shù)作為算法終止的判定條件，以克服Lanczos過程產(chǎn)生的特征向量失去正交性的不足，構(gòu)造求解式AX=B的新算法。

證明：設(shè)Xk為方程(A-Δ)X=B的解，則AXk-B=ΔXk所以-Rk=ΔXk

(3)

(4)

證畢。

2.2 極小向后擾動(dòng)塊Lanczos方法(BMINBACK方法)

本文以極小向后擾動(dòng)范數(shù)式(4)作為新算法迭代終止的判定條件，給出求解方程組AX=B的極小向后擾動(dòng)塊Lanczos方法(BMINBACK方法)，算法過程如下：

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

本文的數(shù)值實(shí)驗(yàn)均在主頻1.8 GHz的酷睿i7CPU、內(nèi)存8 GB、Win10企業(yè)版操作系統(tǒng)的筆記本電腦上采用Matlab R2014a編程實(shí)現(xiàn)。

例1 已知對(duì)稱矩陣

圖1 BMINBACK算法殘量隨迭代次數(shù)的變化

由圖1可知，采用BMINBACK方法求解方程組AX=B計(jì)算的殘量達(dá)到小數(shù)點(diǎn)后5位，驗(yàn)證了極小向后擾動(dòng)塊方法是有效的。

例2 已知矩陣A、多右端項(xiàng)B和初始估計(jì)X0(同例1)，當(dāng)n=100，k=8時(shí)，多右端項(xiàng)B的列數(shù)s變化時(shí)，BLanczos方法和BMINBACK方法的殘量結(jié)果見表1。

由表1可知，在求解方程組AX=B時(shí)，BMINBACK方法計(jì)算速度接近BLanczos方法的計(jì)算速度，BMINBACK方法的殘量略高于BLanczos方法的殘量，但BMINBACK方法的執(zhí)行過程中，采用了極小向后擾動(dòng)范數(shù)作為循環(huán)迭代的條件，克服了BLanczos方法采用殘量范數(shù)作為迭代條件而容易引起算法中斷的不足，進(jìn)一步驗(yàn)證了新方法是有效性的。

表1 BLanczos方法與BMINBACK方法殘量sBMINBACK方法CPU/s殘量范數(shù)BLanczos方法CPU/s殘量范數(shù)20.08761.81×10-50.01603.73×10-640.08832.35×10-70.00171.68×10-860.09559.45×10-50.00252.24×10-1080.13226.31×10-50.00271.69×10-7

4 結(jié)論

本文提出了求解多右端對(duì)稱線性方程組的新方法——極小向后擾動(dòng)塊方法(BMINBACK方法)，并對(duì)新算法做了殘量分析。數(shù)值實(shí)驗(yàn)結(jié)果表明，新算法求解方程組的殘量范數(shù)達(dá)到小數(shù)點(diǎn)后5位，在新算法的執(zhí)行過程中，采用極小向后擾動(dòng)范數(shù)作為循環(huán)迭代的條件，克服了BLanczos方法以殘量范數(shù)作為循環(huán)迭代條件而容易引起算法中斷的不足。而新方法的殘量不是特別小，在滿足一定精度的要求下，理論和數(shù)值實(shí)驗(yàn)說明了新方法的有效性。