總體極小向后擾動(dòng)塊方法求解多右端對(duì)稱(chēng)線性方程組

朱景福，李欣，孟亞輝

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

在理論研究和工程應(yīng)用中，比如控制理論、工程振動(dòng)反問(wèn)題、結(jié)構(gòu)力學(xué)等許多領(lǐng)域的應(yīng)用中經(jīng)常遇到求解多右端對(duì)稱(chēng)項(xiàng)線性方程組

其中A?Rn×n是非奇異對(duì)稱(chēng)矩陣，X，B?Rn×S。

塊krylov子空間方法[1]是求解多右端線性方程組的常用方法。塊方法的基本思想是，在求解多右端線性方程組時(shí)，不是對(duì)每一個(gè)方程都單獨(dú)應(yīng)用迭代法的方式，而是采用同時(shí)對(duì)所有方程應(yīng)用迭代，這樣的計(jì)算方法會(huì)更有效。塊方法一經(jīng)產(chǎn)生就引起了數(shù)學(xué)家們的很大興趣。在上世紀(jì)八十年代Underwood.R和Saad Y[2-3]就提出了塊lanczos方法。1980年O’Leary在共軛梯度法（CG）的基礎(chǔ)上提出了的塊CG方法（BCG）[4]，九十年代初Sadkan.M和Vital.B提出求解對(duì)稱(chēng)正定方程組的塊Davidson方法[5]，進(jìn)一步論證了塊方法的有效性。1996年Simoncini V和Gallopoulos E[6-7]提出了塊GMRES方法（BGMRES）及其變形。1997年Freund R W和Malhotra M[8]提出求解多右端的非Hermitian線性方程組的塊QMR方法（BQMR），同年Chan T F和Wan W L[9]對(duì)求解多右端線性方程組的投影方法做了較全面的分析。1999年顧桂定等提出了求解多右端非對(duì)稱(chēng)線性方程組的塊EN方法（BEN）[10-11]。2000年戴華[12]教授研究了求解多右邊大型稀疏對(duì)稱(chēng)線性方程組的塊Lanczos算法和塊Minres算法，并建立了兩種算法之間的關(guān)系。2008年鐘寶江等[13]提出了一種更簡(jiǎn)單的塊GMRES方法，該方法避免了塊上Hessenberg矩陣的因式分解。因此，編程要簡(jiǎn)單得多，所需的工作也少得多。同年張建華等[14]提出求解具有多個(gè)右端項(xiàng)線性方程組的總體CGS算法。Taherian A和Toutounian F[15]提出了求解多右端非對(duì)稱(chēng)線性方程組的塊GPBi-CG方法。一般地，在討論算法的收斂性和穩(wěn)定性時(shí)，Krylov子空間方法通常用殘量范數(shù)作為判斷算法終止的條件，即極小殘量法[16-18]。通常，若近似解是精確的，殘量范數(shù)是小的，但是反過(guò)來(lái)不一定。為克服殘量范數(shù)作為終止條件的不足，數(shù)學(xué)工作者們考慮采用系數(shù)矩陣擾動(dòng)的方法求解線性方程組，1995-1997年Kasenally[19-20]在用GMRES方法解非對(duì)稱(chēng)線性方程組時(shí)，考慮求滿足擾動(dòng)方程的近似解使向后擾動(dòng)范數(shù)‖△A‖F(xiàn)極小化，給出了求解大型非對(duì)稱(chēng)線性方程組的廣義極小向后擾動(dòng)方法。1998年曹志浩[22]推廣了廣義極小向后擾動(dòng)方法，把總體向后擾動(dòng)方程的精確解看做所求方程的精確解，給出求解大型非對(duì)稱(chēng)線性方程組的總體GMBACK方法[21]。孫繼廣在2001年出版了《矩陣擾動(dòng)分析》，系統(tǒng)地研究了矩陣特征值擾動(dòng)分析和求解方程組的擾動(dòng)分析問(wèn)題。李欣等[23-24]提出了求解對(duì)稱(chēng)線性方程組的總體最小擾動(dòng)方法和求解對(duì)稱(chēng)多右端位移方程組的種子投影方法。Zhang Lei-Hong等[25]提出向后攝動(dòng)分析和基于殘差的誤差界的線性相關(guān)特征值問(wèn)題。

考慮在求解多右端對(duì)稱(chēng)線性方程組時(shí)，采用塊Lanczos過(guò)程與總體極小向后擾動(dòng)相結(jié)合的方法，分析論證總體向后擾動(dòng)的格式及總體極小向后擾動(dòng)范數(shù)的求法，建立求解多右端對(duì)稱(chēng)線性方程組的總體極小向后擾動(dòng)塊方法。本文的結(jié)構(gòu)安排如下：首先引言和預(yù)備知識(shí)；在求解多右端對(duì)稱(chēng)線性方程組的塊Lanczos過(guò)程中，做總體向后擾動(dòng)分析，提出用總體極小向后擾動(dòng)范數(shù)作為判斷算法終止的準(zhǔn)則，給出總體極小向后擾動(dòng)塊lanczos方法；最后通過(guò)數(shù)值實(shí)驗(yàn)驗(yàn)證新算法的可行性和有效性。

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

下面先引入一些所采用的記號(hào)、基本概念和引理。

A+表示矩陣A的Moore-Penrose廣義逆。

‖A‖F(xiàn)、‖A‖2分別表示矩陣A的Frobenius范數(shù)和2-范數(shù)。

xT和XT分別表示向量x和矩陣X的轉(zhuǎn)置。

νec（A）表示矩陣A的列拉直（列展開(kāi)）。

A?B表示A與B的Kronecker積（直積，張量積）。

上面定義的記號(hào)均見(jiàn)參考文獻(xiàn)[26]。

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

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

證明由矩陣的Frobenius范數(shù)定義和2-范數(shù)定義結(jié)合矩陣的列拉直定義易得。

引理2.3[26]設(shè)A∈Cm×n，B∈Cp×q，C∈Cm×q，矩陣方程AXB=C有解的充分必要條件為AA+CB+B=C，在有解的情況下AXB=C的通解為

其中Y∈Cn×p是任意的，A+為矩陣A的Moore-Penrose廣義逆；如果AXB=C無(wú)解，則AXB=C的最小二乘解仍為（2.1），并且AXB=C的極小最小二乘解為X=A+CB+。

引理2.4[26]設(shè)A是任意的m×n矩陣，若A的奇異值分解為

其中U∈Rm×m和V∈Rn×n，Σ=diag（σ1，σ2，…σr）（σ1，σ2，…σr是A的正奇異值），則的Moore-Penrose廣義逆A+

證明由Moore-Penrose廣義逆的定義直接驗(yàn)證可得。

下面回顧塊Lanczos過(guò)程。首先選取多右端項(xiàng)對(duì)稱(chēng)線性方程組AX=B（1.1）的初始矩陣X0，并由QR分解得到X0=Q1T1。

算法2.1塊Lanczos過(guò)程[2-3，26]

（1）選?。?.1）的初始矩陣X0，X0=Q1T1；令Q0=0，計(jì)算M1=Q1TAQ1；

（2）for j=1:k

（2.1）計(jì)算Pj+1=AQj-Qj-1TTj-QjMj；

（2.2）將Pj+1進(jìn)行QR分解Pj+1=Qj+1Tj+1，其中Qj+1是n×s規(guī)范正交矩陣，Tj+1是s×s的上三角矩陣；

（2.3）計(jì)算Mj+1=QTj+1AQj+1；

End for

上式兩邊同時(shí)左乘以分塊矩陣（Q1Q2…Qk），

因?yàn)镼k是規(guī)范正交矩陣，則有

由矩陣的分塊乘法易得

用通式表示為三項(xiàng)遞推公式

變形為

若記Qj+1Tj+1=Pj+1，即對(duì)Pj+1做QR分解可得Qj+1和Tj+1，再計(jì)算Mj+1=QTj+1AQj+1，即可得Mj+1，這樣可構(gòu)造出矩陣及。

在（2.2）式中的三項(xiàng)遞推公式可用如下矩陣形式來(lái)表示

則（2.3）式可以表示為

設(shè)x0（i）（i=1，…，s）是有s個(gè)右端項(xiàng)的對(duì)稱(chēng)線性方程組（1.1）的初始估計(jì)，記X0=[x0（1），x0（2），…，x0（s）]，則R0=B-AX0=[r0（1），…，r0（s）]。

可令根據(jù)塊Lanczos方法產(chǎn)生的方程組（1.1）的近似解有如下形式：

用矩陣表示為

其中Yk=[yk（1），yk（2），…，yk（s）]，由Yk=EkR1求出相應(yīng)的Yk。殘量矩陣表示為

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

2.1 總體極小向后擾動(dòng)分析

本節(jié)討論如何建立求解多右端對(duì)稱(chēng)線性方程組（1.1）總體極小向后擾動(dòng)的塊方法，如何構(gòu)造總體向后擾動(dòng)的格式，以及總體極小向后擾動(dòng)范數(shù)的求法，在新方法執(zhí)行塊Lanczos的過(guò)程中如何與總體極小向后擾動(dòng)相結(jié)合。

首先把多右端對(duì)稱(chēng)線性方程組AX=B（1.1）的近似解看作擾動(dòng)方程

的精確解。這里我們把[△，δ]稱(chēng)為總體向后擾動(dòng)。這種求解方程組（1.1）近似解的方法，稱(chēng)作擾動(dòng)方法，多年來(lái)擾動(dòng)方法得到廣大數(shù)學(xué)工作者的重視，見(jiàn)參考文獻(xiàn)[19，20，21，23]。

下面討論把塊Lanczos過(guò)程與總體極小向后擾動(dòng)方法相結(jié)合求解（2.6）式中的Yk，即求解滿足min{‖△，δ‖F(xiàn)｜（A-△）X=B+δ}的Yk值。

首先以定理的形式給出總體向后擾動(dòng)的格式，以及總體極小向后擾動(dòng)范數(shù)的求法。

定理3.1若Xk∈Rn×S是方程組（3.1）的精確解，則總體向后擾動(dòng)[△，δ]=-Rk（XkTXk+E）-1（XkTE），且總體極小向后擾動(dòng)范數(shù)為

這里νec（*）表示矩陣*的列拉直?！?‖F(xiàn)表示矩陣*的Frobenius范數(shù)，‖*‖2表示矩陣*的2范數(shù)。具體定義見(jiàn)參考文獻(xiàn)[26]。

證明：設(shè)Xk∈Rn×S是（3.1）的精確解，代入擾動(dòng)方程（3.1），則

若記Rk=B-AXk表示方程組（1.1）的殘量，則上式可化為

由矩陣乘法，上面的（3.2）式可以看作

由引理2.3方法解此矩陣方程，得總體向后擾動(dòng)

下面考慮求解總體向后擾動(dòng)[△，δ]的范數(shù)的極小值，由引理2.4得

由Frobenius范數(shù)和2-范數(shù)的定義及引理2.1和引理2.2

化簡(jiǎn)為

則

即

證畢。

2.2 總體極小向后擾動(dòng)塊方法

根據(jù)3.1節(jié)的分析，本文在求解方程組（1.1）的塊Lanczos過(guò)程中，把（3.5）式的總體極小向后擾動(dòng)范數(shù)作為迭代終止的條件，在滿足一定精度要求時(shí)，給出求解多右端對(duì)稱(chēng)線性方程組的如下算法——總體極小向后擾動(dòng)塊Lanczos方法（TBMINBACK）。

算法3.1總體極小向后擾動(dòng)塊Lanczos方法（TBMINBACK）

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

數(shù)值實(shí)驗(yàn)均在個(gè)人計(jì)算機(jī)上實(shí)現(xiàn)，具體配置：CPU：intel（R）Core（TM）i7-8550u，主頻：1.8GHz，內(nèi)存：8GB，系統(tǒng)：Win10企業(yè)版，軟件：Matlab R2014a

由于Matlab軟件的優(yōu)越性，在算法的執(zhí)行過(guò)程中，比如QR分解、求矩陣的范數(shù)等，都可以直接調(diào)用現(xiàn)成的命令得到，這對(duì)算法的執(zhí)行帶來(lái)了更便利的條件。

下面利用算法2.1塊Lanczos方法（記為BLanczos）和算法3.1（TBMINBACK）兩種方法求解方程組（1.1），把兩種方法的計(jì)算速度和殘量做對(duì)比，比較兩種方法的有效性。

例4.1已知對(duì)稱(chēng)矩陣

表1 兩種方法殘量對(duì)比Table 1 Residual comparison of two methods

由表1可見(jiàn)，在求解多右端對(duì)稱(chēng)方程組（1.1）時(shí)，TBMINBACK方法和BLanczos方法的計(jì)算速度都達(dá)到小數(shù)點(diǎn)后的第三位，但兩者差距很小，這是因?yàn)門(mén)BMINBACK方法是在BLanczos的基礎(chǔ)上改進(jìn)的，兩者都采用了塊Lanczos過(guò)程，新方法的運(yùn)算步驟相對(duì)增加，自然運(yùn)算速度會(huì)有相應(yīng)的增加。但是，在目前計(jì)算機(jī)性能日新月異的改進(jìn)中，如此小的速度差異可以忽略不計(jì)。更重要的是，顯然TBMINBACK方法計(jì)算的殘量范數(shù)更小。

例4.2已知對(duì)稱(chēng)矩陣

解：選取方程組（1.1）的初始估計(jì)X0=當(dāng)矩陣A的規(guī)模n=100，塊Lanczos過(guò)程循環(huán)次數(shù)k=8時(shí)，隨著右端項(xiàng)B的列數(shù)S變化，用算法2.1（BLanczos）和算法3.1（TBMINBACK）兩種方法求解方程組（1.1），實(shí)驗(yàn)結(jié)果見(jiàn)表2。

表2 兩種方法殘量對(duì)比Table 2 Residual comparison of two methods

由表2進(jìn)一步得出TBMINBACK方法的殘量范數(shù)比BLanczos方法的更小。特別是，當(dāng)多右端項(xiàng)的列數(shù)S比較大時(shí)，TBMINBACK方法的殘量更加小。說(shuō)明新方法對(duì)于多右端項(xiàng)更大型的對(duì)稱(chēng)線性方程組的求解更有效。

4 結(jié)論

提出了求解多右端對(duì)稱(chēng)線性方程組的新方法——總體極小向后擾動(dòng)塊Lanczos方法（TBMINBACK）。新算法在塊Lanczos方法的基礎(chǔ)上做了改進(jìn)，把總體極小向后擾動(dòng)范數(shù)作為迭代的終止條件，建立了總體極小向后擾動(dòng)塊Lanczos方法（TBMINBACK），論文給出了總體向后擾動(dòng)的格式和總體極小向后擾動(dòng)范數(shù)的求法，并對(duì)新算法做了殘量分析，從理論證明了新方法的可行性。

數(shù)值實(shí)驗(yàn)結(jié)果表明總體極小向后擾動(dòng)塊Lanczos方法（TBMINBACK）在求解方程組（1.1）時(shí)，殘量范數(shù)達(dá)到小數(shù)點(diǎn)后10位且明顯優(yōu)于塊Lanczos方法（BLanczos）。數(shù)值實(shí)驗(yàn)進(jìn)一步驗(yàn)證了新方法的可行性、有效性和優(yōu)越性。特別是，當(dāng)多右端項(xiàng)的列數(shù)S比較大時(shí)，TBMINBACK方法的殘量更加小。說(shuō)明新方法對(duì)于多右端項(xiàng)更大型的對(duì)稱(chēng)線性方程組的求解更有效。

需要說(shuō)明的是在數(shù)值實(shí)驗(yàn)環(huán)節(jié)，本文選取的多右端對(duì)稱(chēng)線性方程組（1.1）的系數(shù)矩陣A都是比較大型的，同時(shí)也考慮了系數(shù)矩陣A的非零元素比例的情況（稠密或者稀疏）。通過(guò)多組實(shí)驗(yàn)發(fā)現(xiàn)，在系數(shù)矩陣A比較稀疏的情況下，總體極小向后擾動(dòng)塊Lanczos方法（TBMINBACK）更有優(yōu)勢(shì)。稀疏矩陣方程組的求解問(wèn)題，在大型科學(xué)工程計(jì)算、圖像處理等研究領(lǐng)域經(jīng)常會(huì)遇到，說(shuō)明新算法具有一定的實(shí)踐意義。