求解BTTB系統(tǒng)的迭代算法

曹蓉

（汕頭職業(yè)技術(shù)學(xué)院 自然科學(xué)系，廣東 汕頭 515041）

求解BTTB系統(tǒng)的迭代算法

曹蓉

（汕頭職業(yè)技術(shù)學(xué)院 自然科學(xué)系，廣東 汕頭 515041）

BTTB矩陣在信號處理等工程問題中有著廣泛的應(yīng)用，因此，針對這種類型矩陣的特點(diǎn)，利用它們的結(jié)構(gòu)來設(shè)計一些數(shù)值穩(wěn)定的、收斂性能好的快速算法，具有極為重要的意義.文章討論了塊三角Toeplitz矩陣的一些性質(zhì)，給出了求解塊下三角Toeplitz矩陣逆的快速算法，并對其復(fù)雜性進(jìn)行了分析.利用這種求逆算法進(jìn)而給出了求解BTTB系統(tǒng)的塊Gauss-Seidel迭代算法和塊SOR迭代算法，并討論了其收斂性.數(shù)值實(shí)驗(yàn)得到驗(yàn)證.

BTTB；塊Gauss-Seidel迭代；塊SOR迭代

Toeplitz矩陣出現(xiàn)在數(shù)字信號處理領(lǐng)域的許多應(yīng)用中，是應(yīng)用最廣泛的特殊矩陣之一.這些問題很多都?xì)w結(jié)為對Toeplitz系統(tǒng)的求解.因此涌現(xiàn)了很多關(guān)于Toeplitz系統(tǒng)的快速求解，如[1-5].特殊的BTTB矩陣在高維信號、圖像復(fù)原、數(shù)值偏微分方程、折積型積分方程等有廣泛的應(yīng)用.因此，針對這類型矩陣的特點(diǎn)，利用它們的結(jié)構(gòu)來設(shè)計一些數(shù)值穩(wěn)定的、收斂性能好的快速算法具有極為重要的意義.現(xiàn)在已有一些人對塊Toeplitz系統(tǒng)和特殊的BTTB系統(tǒng)求其快速算法[6-9].因此本文在此基礎(chǔ)上利用B-FFT技術(shù)，討論BTTB系統(tǒng)的求解.對于這些線性方程組的求解主要分為直接法和迭代法.由于直接法求解的運(yùn)算量為O（N3），另外直接解法非常不穩(wěn)定，因此更多人選擇迭代法求解.對于迭代法，一個好的預(yù)條件子能使其運(yùn)算量大大降低，因此人們更專注于迭代法的研究.M.K.Ng在1997年給出了mn階BTTB方程組的帶狀預(yù)處理迭代法[10].H.Akaike和游兆永[11-12]已對塊Toeplitz矩陣求其逆給出算法.本文利用B-FFT技術(shù)給出塊三角BTTB矩陣的快速求逆算法，并將這個求逆算法應(yīng)用到塊Gauss-Seidel迭代算法和塊SOR迭代算法.

一個矩陣Am，n有如下結(jié)構(gòu)

則稱Am，n為一個m，n階的塊Toeplitz矩陣，其中Ai（1-n≤i≤n-1）為m階的任意矩陣，即Am，n主對角線上的分塊矩陣相同，且平行于主對角線線上的分塊矩陣也相同.若Ai（1-n≤i≤n-1）為Toeplitz矩陣，即則稱Am，n為 BTTB（Block-Toeplitz-Toeplitz-Block）矩陣.

對于求解線性方程組Am，nx=b（Am，n為（1）中的形式）的迭代格式一般可表示為

其中Am，n=M-N為矩陣Am，n的分裂.如果M非奇異，顯然，任給初始向量x0，由方程（2）生成的迭代序列{xk}收斂的充要條件是迭代矩陣H=M-1N的譜半徑ρ（H）＜1.下面首先討論塊下三角BTTB矩陣的逆.

1 塊下三角BTTB矩陣的逆

定義1[13]如果Am，n=（Ai）是一個mn階的Toeplitz矩陣，則Am，n=C（1∶mn，1∶mn），其中C是一個（2n-1）m階的塊循環(huán)矩陣且：

定理1[13]一個塊Toeplitz矩陣乘一個塊向量可轉(zhuǎn)化為一個塊循環(huán)矩陣乘這個塊向量.

即若Y=CX，C為（3）中的塊循環(huán)矩陣，X，Y分別為mn維的塊向量，X（i），Y（i）表示塊向量X，Y的第i個m維向量，并且X=（n+1∶2n-1）=0，0為（n-1）m維的零向量，則.

推論1[13]塊循環(huán)矩陣Cm，n可塊對角化，即：

其中F(n)=Fn?Im，F(xiàn)n是n階離散Fourier變換（Dis?crete Fourier Transform即DFT）矩陣，滿足?表示Kronecker積，Λ表示一個mn階的對角陣，即Λ=Blkdiag（F（n）Cm，n（∶，1∶m）），這里Cm，n（∶，1∶m）是塊矩陣Cm，n的前m列（即第1列塊），Blkdiag是由括號里的塊向量生成的塊對角陣.由參考文獻(xiàn)[14]知，若vec(Y)=(Fn?Im)vec(X)，則能得到

因此形如z=Cm，nx的這種一個塊循環(huán)矩陣與一個向量相乘可以使用B-FFT技術(shù)，由公式（4）知，

引理1[1]表示階數(shù)為mn階的塊下三角To?eplitz矩陣的集合：

由此可得出，A-1的遞歸方程：

根據(jù)引理1、推論1和推論2得到塊下三角BTTB矩陣逆的快速算法.從推論2中公式（6）知，若計算塊下三角BTTB矩陣AM的逆，當(dāng)知道了時，則，只需要計算BN.由引理1的（c）知，是一個塊下三角Toeplitz矩陣，因此，只需求出BN的第1列塊，它的第1列塊等于積的第1列塊，而與CN都是塊Toeplitz矩陣，能嵌入到塊循環(huán)矩陣.另外注意到，BTTB矩陣Am，n的每一小塊也是To?eplitz矩陣.

下面給出塊下三角BTTB矩陣的逆的算法.不妨設(shè)n=2q，min=2l，（0≤l＜q）.對于n非2的冪，見后面的分析.

算法1 塊下三角BTTB矩陣的逆的快速算法.

算法1中的函數(shù)hti（·）指的是Doubling算法[15]，cirsys（·）指的是利用FFT技術(shù)計算循環(huán)矩陣乘以一個向量，strass（·）表示兩個矩陣相乘的Strassen算法，blkcirsys（·）與bttbcirsys（·）均利用B-FFT技術(shù)計算塊循環(huán)矩陣乘以塊向量.

當(dāng)n非2的冪時（如2q-1＜n＜2q），計算塊下三角BTTB矩陣Am，N（N=2q）的逆，它的第1列塊是Am，n的第1列塊和m階的零塊矩陣，即（A0，A1，…，An-1，O，…，O）.引理1的（d）保證了Am，N的第1列塊的前n個與Am，n的第1列塊相同.

因此這三個算法的復(fù)雜度近似為

先討論塊循環(huán)矩陣乘以塊向量的復(fù)雜度.由推論1知，計算y=Cm，nx，x=（x1，x2，…，xn）T，xi是個m維向量，可以使用B-FFT計算，即

因此一個mn階的塊循環(huán)矩陣乘以塊向量使用B-FFT和基2-FFT后的復(fù)雜度為

此時的矩陣大小為（n-1）m×（n-1）m，因此計算的第1列塊的復(fù)雜度為

2 BTTB矩陣的兩種迭代算法

a）塊Gauss-Seidel迭代算法

給定方程組Am，nx=b（Am，n為（1）中的形式），令A(yù)m，n=D-L-U，其中

如果在方程組（2）中取，M=D-L，N=U那么就得到如下求解Am，nx=b的塊Gauss-Seidel迭代格式：

可以看到D-L是塊下三角BTTB矩陣，故可以利用算法1求其逆，因此得到塊Gauss-Seidel算法如下：

算法2 塊Gauss-Seidel迭代算法

1）給定Am，n，x0，b，n=0，精度e；

2）利用算法1求D-L逆H；

3）使用B-FFT計算H*b；

4）Fork=1，2，…

5）使用基2-FFT和兩次B-FFT計算HUxk；

6）xk+1=HUxk+Hb n=n+1；

7）當(dāng)‖xk+1-xk‖2＜e時，停止；

8）輸出xk+1，n.

b）塊SOR迭代算法

將線性方程組Am，nx=b兩邊兩邊乘以ω變?yōu)棣谹m，nx=ωb，且在ωAm，n=M-N中取M=D-ωL，N=[ωU+（1-ω）D]，那么得到如下求解ωAm，nx=ωb的塊SOR的迭代格式

其中（0＜ω＜2），由于D-ωL是塊下三角BTTB矩陣，所以可以由算法1求得其逆.由此給出塊SOR迭代算法如下：

算法3 塊SOR迭代算法

1）給定Am，n，x0，b，n=0，精度e；

2）利用算法1求D-ωL逆H；

3）使用B-FFT計算H*b；

4）Fork=1，2，…

7）當(dāng)‖xk+1-xk‖2＜e時，停止；

8）輸出xk+1，n.

3 迭代法的收斂性

定義2[16]一個矩陣A=(aij) ∈Rn×n，稱為：

（1）非負(fù)矩陣，如果aij≥0，i，j=1，2，…，n；

（2）Z型矩陣，如果aij≤0，i≠j，i，j=1，2，…，n；

（3）M矩陣，如果A是非奇異的Z矩陣且A-1≥0；

（4）H矩陣，如果其比較矩陣〈A〉是非奇異的且它是M矩陣，這里〈A〉=（αij），其中

定義3[17]若A,B∈Rm×n且滿足A-B≤O，則稱A小于等于B，記為A≤B.

定義4[17]若A=(aij)∈Rm×n，則定義A的絕對值為每個元素的絕對值，記為|A|=(|αij|).

注意：不能與“方陣的行列式”混淆.

定義2 一個n階矩陣A的分裂A=M-N稱為：

（1）弱正則分裂，如果M-1≥0，且M-1N≥0；

（2）P-正則分裂，如果MT+N正定.

引理2 如果A=(aij)∈Rn×n對稱正定，且A=M-N，P-正則分裂，則ρ（M-1N）＜1；如果A為M矩陣，且A=M-N，弱正則分裂，則ρ（M-1N）＜1.

引理3 令A(yù),B∈Rn×n，

（1）若A是一個M-矩陣，B∈Zn×n（全體n階Z型矩陣構(gòu)成的集合），且A≤B，則B也是一個M-矩陣；

（2）若A是一個H-矩陣，則；

（3）若|A|≤B，則ρ（A）≤ρ（B），ρ（A）表示A的譜半徑.

定理2 如果（1）式中Am，n為H-矩陣或?qū)ΨQ正定，則對任給初始向量x0，則方程（7）生成的迭代序列{xk}收斂.

證明 根據(jù)引理2只要驗(yàn)證塊Gauss-Seidel分裂滿足收斂性的條件即可.

當(dāng)Am，n為對稱正定時，D也對稱正定.又（D-L）T+U=D，即Am，n=（D-L）-U為P-正則分裂，故滿足收斂條件.

當(dāng)Am，n為H-矩陣時，〈Am，n〉為M-矩陣.由引理 3（1）知，〈D-L〉為M矩陣，|U|為非負(fù)矩陣，于是有〈D-L〉-1≥0，.即〈Am，n〉=〈D-L〉-|U|為弱正則分裂，于是

又ρ(B)≤ρ(|B|)和ρ（B）≤ρ（C），如果C≥B≥0以及矩陣不等式有

定理3 如果式（1）中Am，n為H-矩陣或?qū)ΨQ正定，則任給初始向量x0，由方程（8）生成的迭代序列{xk}收斂.

證明 根據(jù)引理2只要驗(yàn)證塊SOR分裂滿足收斂性的條件就行.

當(dāng)Am，n為對稱正定時，又0＜ω＜2，所以（2-ω）Am，n對稱正定，容易知道塊對角矩陣（2-ω）D也對稱正定.又

即ωAm，n=（D-ωL）-[ωU+（1-ω）D]為P-正則分裂，故滿足收斂條件.

當(dāng)Am，n為H-矩陣時，ωAm，n也為H-矩陣，〈ωAm，n〉為M-矩陣，〈D-ωL〉為M-矩陣，|ωU+(1 -ω)D|為非負(fù)矩陣，

4 結(jié)論

本文主要利用B-FFT給出塊下三角BTTB矩陣快速求逆算法，利用其算法給出BTTB系統(tǒng)的兩種迭代算法，分析其收斂性.不足之處，由于BTTB矩陣本身的特殊性，是否可將其收斂的條件進(jìn)行改進(jìn)！

[1]徐仲,張凱院,陸全.Toeplitz矩陣類的快速算法[M].西安:西安工業(yè)大學(xué)出版社,1999:38-45.

[2]曹蓉,童細(xì)心.求解Toeplitz系統(tǒng)循環(huán)和反循環(huán)分裂算法[J].云南民族大學(xué)學(xué)報:自然科學(xué)版,2012,21（5）:356-360.

[3]Ng M K.Iterative Methods for Toeplitz Systems[M].Ox?ford:Oxford University Press,2004.

[4]Chan R H,Ng M K.Conjugate gradient methods for To?eplitz systems[J].SIAM Rev,1996,38(3):427-482.

[5]畢永青.三對角對稱Toeplitz矩陣的解析逆陣[J].西南民族大學(xué)學(xué)報:自然科學(xué)版,2003,29（4）:390-393.

[6]Shi Yongjie,Pi Xuebo.New preconditioners for systems of linear equations with Toeplitz structure[J].Calcolo,2014,51(1):31-55.

[7]Lin furong,Wang chixi.BTTB preconditioners for BTTB systems[J].Numerical Algorithms,2012,60(1):153-167.

[8]馮月華,劉成志,劉仲云.塊Toeplitz方程組的快速塊Gauss-Seidel迭代算法[J].數(shù)學(xué)理論與應(yīng)用,2012,32(1):1-5.

[9]趙敏.塊Toeplitz矩陣的一種快速Q(mào)R分解及算法實(shí)現(xiàn)[J].長江大學(xué)學(xué)報:自然科學(xué)版,2007,4(2):4-5.

[10]Ng M K.Band precondtioners for Block-Toeplitz-To?eplitz-Block Systems[J].Linear algebra and its applications,1997,259:307-327.

[11]Akaike H.Block Toeplitz matrix inversion[M].SIAM J.Ap?pl.Math,1973,2(4):234-241.

[12]游兆永,路浩.塊Toeplitz三角陣求逆及塊Toeplitz三角線性方程組求解的復(fù)雜性[J].數(shù)學(xué)研究與評論,1989,9（1）:101-106.

[13]金小慶.Toeplitz系統(tǒng)預(yù)處理方法[M].北京:高等教育出版社,2010:59-72.

[14]張凱院,徐仲.數(shù)值代數(shù)[M].2版.北京:科學(xué)出版社,2006:36-45.

[15]Morf M.Doubling algorithms for Toeplitz and related equa?tions[J].In:IEEE Conference on Acoustics Speech and Sig?nal Processing,1980:954-959.

[16]Commges D,Monsion M.Fast Inversion of Triangular To?eplitz Matrices[J].IEEE Transactions on automatic control,1984,29(3):250-251.

[17]方保镕,周繼東,李醫(yī)民.矩陣論[M].北京:清華大學(xué)出版社,2004:330-354.

責(zé)任編輯：畢和平

The Iterative Algorithm for Block-Toeplitz-Toeplitz-Block Systems

CAO Rong
（Natural Sciences，Shantou Vocational and Technical College，Shantou515041，China）

BTTB matrices have a wide range of engineering applications such as in signal processing.In view of the charac?teristics of this type of matrix,it is very significant that we design some fast algorithms with numerical stability and the good property of convergence.Firstly,we discussed some properties of block triangular Toeplitz matrices,then we presented fast algorithms for computing the inverse of such a class of matrices and also analyzed the complexity of this algorithm.Using the inverse algorithm,we gave Block-Gauss-Seidel iteration algorithm and Block-SOR iteration algorithm for solving the BTTB system.

BTTB；Block Gauss-Seidel iteration；Block SOR iteration

O 241.6

A

1674-4942(2015)02-0134-05

2015-03-01

汕頭職業(yè)技術(shù)學(xué)院科研基金資助項目（SZK2014Y35）