熊 露，張 婷，李勝坤

(成都信息工程大學(xué) 應(yīng)用數(shù)學(xué)學(xué)院，四川 成都 610225)

0 引言

本文討論一般的大型稀疏Stein矩陣方程

X+AXB=C

(1)

的數(shù)值解法，其中系數(shù)矩陣A∈n×n,Β∈s×s,C∈n×s是給定的矩陣，X∈n×s是待求的未知矩陣．方程(1)存在唯一解的充分必要條件為λi(A)λj(B)≠-1,其中λi(A),λj(B)分別是矩陣A,B的特征值[1]．

Stein矩陣方程是一類特殊的矩陣方程，在離散時間系統(tǒng)、概率、統(tǒng)計、光譜分析，神經(jīng)網(wǎng)絡(luò)和圖象復(fù)原等領(lǐng)域有著廣泛應(yīng)用[2-9]．同時Stein矩陣方程也是控制系統(tǒng)設(shè)計的計算工具[10]．因此，求解Stein矩陣方程顯得尤為重要．

對于低階稠密的Stein矩陣方程，通常采用Kronecker積將其轉(zhuǎn)化為線性方程組后直接求解．由于實際問題產(chǎn)生的Stein矩陣方程通常都是高階稀疏的，因此迭代法[11-13]具有更大的優(yōu)勢．對Stein矩陣方程而言，系數(shù)矩陣間階數(shù)的不同對迭代法的構(gòu)造有著重要的影響．

近年來，一些求解不同類型Stein矩陣方程的迭代方法相繼被提出．當(dāng)矩陣B的階數(shù)遠(yuǎn)遠(yuǎn)小于矩陣A的階數(shù)時，文獻(xiàn)[14]利用塊Arnoldi正交化過程和非對稱塊Lanczos雙正交化過程將Stein方程降階后直接求解，提出了塊Arnoldi型Stein方法和非對稱塊Lanczos型Stein方法．文獻(xiàn)[15]則采用推廣的塊Arnoldi型方法降階求解．對低秩Stein矩陣方程，即矩陣C=EFT的秩很小時，文獻(xiàn)[16-17]利用全局Arnoldi正交化過程對Stein矩陣方程降階求解，提出了Stein全局Arnoldi算法，得到其低秩解．塊Arnoldi型正交化方法[18]，推廣的塊Arnoldi型正交化方法[15]，塊Hessenberg型方法[19]也相繼被提出．對于一般的Stein矩陣方程，文獻(xiàn)[20]提出了Smith型迭代法求解Stein矩陣方程，但是這種方法對系數(shù)矩陣有一定的要求，即當(dāng)ρ(Α)ρ(B)<1時方法才收斂，其中ρ(Α)和ρ(B)分別表示矩陣A和B的譜半徑．文獻(xiàn)[21]利用全局Arnoldi正交化過程提出了全局FOM(the global full orthogonalization method)和全局GMRES算法(the global generalized minimal residual)．文獻(xiàn)[22]利用矩陣Krylov子空間的位移不變性，提出了位移型全局FOM和位移型全局GMRES算法，降低了全局FOM和全局GMRES算法的計算量．

Hessenberg過程[23]是另外一種構(gòu)造Krylov子空間基的方法，通常比Arnoldi正交化過程需要更少的計算量．緊接著，塊Hessenberg過程[19,24]和全局Hessenberg過程[25-26]相繼被提出．本文基于全局Heseenberg過程，提出了位移型全局Hess方法和位移型全局CMRH方法求解Stein矩陣方程．?dāng)?shù)值結(jié)果表明新方法相對于位移型全局Arnoldi型方法而言占有一定的優(yōu)勢．

對給定的兩個實矩陣X,Y∈m×n,定義Frobenius內(nèi)積為〈X,Y〉F=tr(XTY),相應(yīng)的F范數(shù)用‖·‖F(xiàn)表示．如果〈X,Y〉F=0,就稱X,Y正交．另外，X?Y=[xijY]∈mm×nn表示矩陣X和Y的Kronecker積．

為方便敘述，定義算子

M:X→AXB,

即方程(1)等價轉(zhuǎn)化為

X+M(X)=C.

(2)

1 全局Hessenberg過程

本節(jié)介紹全局Hessenberg過程構(gòu)造矩陣Krylov子空間的基的過程．

給定矩陣V∈n×s和算子M,定義矩陣Krylov子空間

Km(M,V)=span{V,M(V),…,Mm-1(V)},

其中，Mi(V)=M(Mi-1(V)),i=2,…,m-1.

全局Hessenberg過程構(gòu)造基{V1,V2,…,Vm}滿足的正交條件為：

Vk+1⊥FY1,Y2,…,Yk,k=1,2,…,m,

即

算法1 全局Hessenberg過程[24-25]1.令θ=F=(V)1,1,V1=V/θ;2.Forj=1,2,…,m do3.W=M(Vj);4.Fori=1,2,…,j do5.hi,j=F/F=(W)i,1;6.W=W-hi,jVi;7.Endfor8.hj+1,j=F=(W)j+1,1;若hj+1,j=0,則停機;9.Vj+1=W/hj+1,j10.Endfor

假設(shè)算法1在第m步之前不中斷，則算法1生成的V1,V2,…,Vm構(gòu)成Krylov子空間Km(M,V)=span{V,M(V),…,Mm-1(V)}的一組非正交基．

若Zm∈Km(M,V),則Zm可以表示為

其中，

算法2 帶主元策略的全局Hessenberg過程[24-25]1.定義i0,j0,滿足(V)i0,j0=max{|vi,j|}1≤j≤s1≤i≤n;2.θ=(V)i0,j0,V1=V/θ;p1,1=i0,p1,2=j0;3.Forj=1,2,…,m do4.W=M(Vj);5.Fori=1,2,…,j do6.hi,j=(W)pi,1,pi,1;7.W=W-hi,jVi;8.Endfor9.定義i0,j0,滿足(W)i0,j0=max{|wi,j|}1≤j≤s1≤i≤n;10.hj+1,j=(W)i0,j0;11.Vj+1=W/hj+1,j;12.pj+1,1=i0,pj+1,2=j0;13.Endfor

不難看出，全局Hessenberg過程(算法2)構(gòu)造搜索空間

Km(M,V)=span{V,M(V),…,Mm-1(V)}的基V1,V2,…,Vm滿足：

因此，下列關(guān)系式成立．

定理1令

(3)

(4)

其中，em=(0,…,0,1)T∈m.

2 位移型全局Hess算法

本節(jié)給出求解Stein矩陣方程(2)的位移型全局Hess算法．給定初始近似值X0∈n×s,R0=C-X0-M(X0)為相應(yīng)的殘差．為了能夠充分利用初始值的相關(guān)信息，第m步迭代解Xm在仿射空間X0+Km(M,R0)中尋找，即Xm∈X0+Km(M,R0),使得殘差Rm滿足

Rm⊥FY1,Y2,…,Ym.

利用全局Hessenberg過程(算法2)，可構(gòu)造矩陣Krylov子空間Km(M,R0)的一組基V1,V2,…,Vm,從而

相應(yīng)殘差為

根據(jù)正交條件Rm⊥FY1,Y2,…,Yk,k=1,2,…,m,可得

(Im+Hm)ym=θe1,

從而

在實際計算時，隨著迭代步數(shù)m的增加，算法的計算量和存儲量會增加．因此，采用固定步數(shù)(m步)的重啟型，相應(yīng)的位移型全局Hess算法(簡記為SGl-Hess(m))描述如下：

算法3 SGl-Hess(m)1.選擇初始近似矩陣X0和精度ε;2.計算R0=C-X0-M(X0);3.運行算法2,得到一組基V1,V2,…,Vm和Hm;4.計算(Im+Hm)ym=θe1,得到y(tǒng)m;5.計算Xm=X0+Vm(ym Is),Rm=C-Xm-M(Xm);6.若‖Rm‖F(xiàn)<ε,停止迭代,輸出Xm,否則繼續(xù);7.令X0=Xm,R0=Rm,返回步驟3.

3 位移型全局CMRH算法

本節(jié)構(gòu)造位移型全局CMRH算法求解Stein矩陣方程(2)．同樣的，給定初始近似值X0,R0為相應(yīng)的殘差．第m步迭代解

其中，

這里，ym采用極小化‖Rm‖F(xiàn)來計算，即

類似于位移型全局Hess算法，位移型全局CMRH算法(簡記為SGl-CMRH(m))描述如下：

算法4 SGl-CMRH(m)1.選擇初始近似矩陣X0和精度ε;2.計算R0=C-X0-M(X0);3.運行算法2,得到一組基V1,V2,…,Vm和Hm+1,m;4.計算miny∈Rm‖θ e1-I～my-Hm+1,my‖2,得到y(tǒng)m;5.計算Xm=X0+Vm(ym Is),Rm=C-Xm-M(Xm);6.若‖Rm‖F(xiàn)<ε,停止迭代,輸出Xm,否則繼續(xù);7.令X0=Xm,R0=Rm,返回步驟3.

下面給出位移型全局CMRH算法的殘差范數(shù)與位移型全局GMRES算法的殘差范數(shù)之間的一個關(guān)系．

其中

由范數(shù)的性質(zhì)可得

由最小范數(shù)條件可得

因此，

4 數(shù)值算例

下面給出兩個數(shù)值算例說明位移型全局Hess算法(SGl-Hess(m))和位移型全局CMRH算法(SGl-CMRH(m))的有效性，并且與位移型全局FOM(SGl-FOM(m))和位移型全局GMRES方法(SGl-GMRES(m))[22]在迭代步數(shù)和計算時間等方面進(jìn)行比較．所有數(shù)值計算都是用Matlab R2012a編程．在所有的數(shù)值算例中，選擇的初始近似矩陣X0為零矩陣，停止迭代的條件為‖Rm‖F(xiàn)<10-8或者重啟次數(shù)不超過200次．

表1 例1的數(shù)值結(jié)果

圖1 例1的收斂曲線

通過表1可以看出，在給定的相同停止迭代條件下，SGl-Hess(m)和SGl-CMRH(m)方法在迭代步數(shù)上略少于SGl-FOM(m)和SGl-GMRES(m)方法，但是在CPU運行時間方面的優(yōu)勢非常明顯．原因在于全局Hessenberg過程比全局Arnoldi正交化過程的計算量小，在迭代步數(shù)相差不多的情況下，運行時間更少．另外通過圖1可以看出，SGl-Hess(m)方法的殘差光滑性比其他三種方法稍微差點．

矩陣C的選取滿足精確解X=I∈n×s.在這里，取

固定步數(shù)分別取m=5和m=10，計算結(jié)果見表2和圖2．

表2 例2的數(shù)值結(jié)果

圖2 例2的收斂曲線

通過表2可以看出，SGl-Hess(5)和SGl-CMRH(5)方法在重啟次數(shù)和CPU運行時間方面都優(yōu)于SGl-FOM(5)和SGl-GMRES(5)方法；SGl-Hess(10)，SGl-CMRH(10)，SGl-FOM(10)和SGl-GMRES(10)方法的迭代次數(shù)相當(dāng)，但CPU運行時間方面SGl-Hess(10)和SGl-CMRH(10)方法明顯優(yōu)于SGl-FOM(10)和SGl-GMRES(10)方法．圖2表明SGl-Hess(m)方法的殘差出現(xiàn)一些微小的震蕩．

5 總結(jié)

基于全局Hessenberg過程，利用矩陣Krylov子空間的位移不變性，構(gòu)造了新的算法，提出了位移型全局Hess算法和位移型全局CMRH算法求解一般的大型Stein矩陣方程，給出了相應(yīng)的殘差估計．?dāng)?shù)值例子表明，新方法相對于SGl-FOM(m)和SGl-GMRES(m)方法而言，在運行時間和迭代步數(shù)上占有一定的優(yōu)勢．