大型連續(xù)Sylvester 方程外推的CSCS 迭代

劉仲云 張芳

(長(zhǎng)沙理工大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,長(zhǎng)沙,410114)

1 引言

大型連續(xù)Sylvester 方程[1]是最常用的線(xiàn)性矩陣方程之一.矩陣方程形如

這里A ∈Rn×n,B ∈Rm×m,E ∈Rn×m是給定的矩陣.眾所周知,當(dāng)系數(shù)矩陣A和?B沒(méi)有公共特征值時(shí),方程(1.1)有唯一解[2,3].

Sylvester 方程是設(shè)計(jì)Luenberger 觀測(cè)的經(jīng)典方法, 它廣泛應(yīng)用于信號(hào)處理、控制和系統(tǒng)理論[4–7],也常出現(xiàn)在計(jì)算不變量子空間的Riccati方程的線(xiàn)性和廣義特征值問(wèn)題中[8–10],也可用于設(shè)計(jì)常微分方程數(shù)值解的隱式龍格 庫(kù)塔積分公式和分塊多步公式中[11].矩形域上分離橢圓邊值問(wèn)題的有限差分離散產(chǎn)生的方程也可以寫(xiě)成Sylvester 方程[12,13],優(yōu)化理論、結(jié)構(gòu)動(dòng)力學(xué)、固體力學(xué)、地質(zhì)分子光學(xué)等領(lǐng)域都涉及到Sylvester 方程的求解[14].

求解連續(xù)Sylvester 方程(1.1)主要有兩種方法.第一種方法是將矩陣方程(1.1)轉(zhuǎn)換成一個(gè)線(xiàn)性方程Φx=c,這里Φ=Im ?A+BT ?In是矩陣A和BT的Kronecker 積的和,向量x和c是矩陣X和C的列向量,直接法和迭代法均可用于求解該線(xiàn)性方程. 第二種方法是對(duì)方程(1.1)進(jìn)行迭代.

當(dāng)系數(shù)矩陣階數(shù)比較大時(shí),線(xiàn)性方程Φx=c中的系數(shù)矩陣Φ 的階數(shù)會(huì)相當(dāng)大,一般會(huì)出現(xiàn)數(shù)據(jù)存儲(chǔ)和計(jì)算時(shí)間方面的困難.這說(shuō)明第一種方法主要用于中小維問(wèn)題.例如: 文獻(xiàn)[15]中提出的Bartel stewart 方法是基于使用特征值QR 算法將矩陣A和B化為實(shí)Schur 形式,然后使用直接法來(lái)求解幾個(gè)線(xiàn)性方程.文獻(xiàn)[16]提出的Hessenberg Schur 法是將矩陣A化為Hessenberg 形式,將矩陣B分解為擬三角Schur 形式,其速度比Bartels Stewart 方法快.這些方法都屬于直接法,并在MATLAB 中進(jìn)行了應(yīng)用.

對(duì)于大型稀疏Sylvester 方程,一般用迭代法求解. 比較常見(jiàn)的有古典迭代法,Smith 迭代法[17],ADI 迭代法[18]等. 最近比較流行的迭代法, 是利用系數(shù)矩陣A和B的分裂來(lái)建立類(lèi)似ADI 迭代格式的迭代法,例如,白中治在2011 年提出的HSS 迭代法[19].自該方法提出后,幾個(gè)更為有效的類(lèi)似的分裂迭代法相繼出現(xiàn),例如: 2013 年王湘提出的PSS 迭代法[20], 2014 年鄭青青提出的NSS 迭代法[21]等.這些方法都無(wú)條件收斂到方程的精確解.

文獻(xiàn)[1]考慮了(1.1)中系數(shù)矩陣A和B都是正定Toeplitz 矩陣的情況,提出了CSCS 迭代法.通過(guò)這種CSCS 迭代方法,將求解一般連續(xù)Sylvester 方程的問(wèn)題轉(zhuǎn)化為包含兩個(gè)循環(huán)矩陣和反循環(huán)矩陣的連續(xù)Sylvester 方程的求解問(wèn)題.通過(guò)這種轉(zhuǎn)化可以利用傅里葉算法快速求解.在對(duì)流擴(kuò)散方程的離散化中就出現(xiàn)了這種結(jié)構(gòu)的矩陣[13,22].

上述這些情況激發(fā)了我們對(duì)連續(xù)Sylvester 方程進(jìn)一步的研究興趣. 本文提出一個(gè)外推的CSCS(ECSCS)迭代.文章的結(jié)構(gòu)如下: 第二節(jié)介紹基本定義,第三節(jié)介紹CSCS 迭代,第四節(jié)給出ECSCS 迭代法,并分析其收斂性,第五節(jié)通過(guò)數(shù)值實(shí)驗(yàn)驗(yàn)證該方法的有效性.

2 基本定義

對(duì)于任意的矩陣K ∈Rn×n,K?,ρ(K)和λ(K)分別表示它的共軛轉(zhuǎn)置矩陣、譜半徑和特征值.Ki,j是矩陣K的第i行j列元素.若x ∈R,則x的實(shí)部記為Re(x),虛部記為Im(x).如果矩陣(K+K?)是正定的,則矩陣K也是正定的.該正定條件也等價(jià)于: 存在一個(gè)非零向量z ∈R,使得Re(z?Kz)>0,即對(duì)于矩陣K的任意特征值λ,有Re(λ)>0[23].

任意一個(gè)Toeplitz 矩陣T都有一個(gè)循環(huán)反循環(huán)分裂,即:T=CT+ST,這里CT是循環(huán)矩陣,ST是反循環(huán)矩陣,定義如下:

眾所周知,循環(huán)矩陣C和反循環(huán)矩陣S可由傅里葉矩陣F和酉矩陣?F=FD進(jìn)行如下對(duì)角化:

3 CSCS 迭代

若方程(1.1)中的系數(shù)矩陣A ∈Rn×n,B ∈Rm×m是正定Toeplitz 矩陣,則矩陣A,B可以進(jìn)行如下循環(huán)反循環(huán)分裂:

給定常數(shù)α,β及適當(dāng)維數(shù)的單位矩陣I,則有

由此可知,連續(xù)Sylvester 方程(1.1)可以等價(jià)地寫(xiě)成不動(dòng)點(diǎn)矩陣方程

文獻(xiàn)[1]中提出的求解(1.1)的CSCS 迭代定義如下.

CSCS 迭代給定一個(gè)初始矩陣X(0)∈Rn×m,使用以下迭代格式計(jì)算X(k+1)∈Rn×m(k=0,1,2,···):

對(duì)于上面的CSCS 迭代,文獻(xiàn)[1]證明了如下收斂性定理.

因而,對(duì)任意γ>0,有ρ(T(γ))≤σ(γ)<1.

4 ECSCS 迭代

為了進(jìn)一步加快CSCS 迭代(3.1)的收斂速度,本節(jié)給出ECSCS 迭代.

ECSCS 迭代給定一個(gè)初始矩陣X(0)∈Rn×m,使用以下迭代格式計(jì)算X(k+1)∈Rn×m(k=0,1,2,···):

用矩陣向量形式,ECSCS 迭代可以等價(jià)為

進(jìn)而變成一步迭代格式

關(guān)于ECSCS 迭代,有如下的收斂性定理:

定理2在定理1 的假設(shè)下,(4.1)中產(chǎn)生的ECSCS 迭代序列X(k)收斂到(1.1)的精確解X?,其中ω ∈(0,).

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

例1 考慮連續(xù)Sylvester 方程(1.1),其中矩陣M,N ∈Rn×n都是三對(duì)角矩陣:M= tridiag(?1,2,?1),N= tridiag(?0.5,0,0.5). 數(shù)值結(jié)果見(jiàn)表1和表2.

表1 SSOR,HSS,CSCS 和ECSCS 的IT 與CPU

表2 SSOR,HSS,CSCS 和ECSCS 的IT 與CPU

例2 考慮Sylvester 方程(1.1),其中矩陣

數(shù)值結(jié)果見(jiàn)表3.

表3 SSOR,HSS,CSCS 和ECSCS 的IT 與CPU

從例1 和例2 中可以看到: CSCS 迭代在迭代次數(shù)和計(jì)算效率上都優(yōu)于HSS 迭代和SSOR 迭代,而ECSCS 比CSCS 收斂的效果更好,并且系數(shù)矩陣的階數(shù)越大,效果越明顯,因此外推的CSCS迭代很好地提高了CSCS 迭代法的收斂速度.