曹蓉

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

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

曹蓉

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

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

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

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

一個矩陣Am，n有如下結構

則稱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矩陣乘一個塊向量可轉化為一個塊循環(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是由括號里的塊向量生成的塊對角陣.由參考文獻[14]知，若vec(Y)=(Fn?Im)vec(X)，則能得到

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

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

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

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

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

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

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

當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列塊相同.

因此這三個算法的復雜度近似為

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

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

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

2 BTTB矩陣的兩種迭代算法

a）塊Gauss-Seidel迭代算法

給定方程組Am，nx=b（Am，n為（1）中的形式），令Am，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）當‖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）當‖xk+1-xk‖2＜e時，停止；

8）輸出xk+1，n.

3 迭代法的收斂性

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

（1）非負矩陣，如果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,B∈Rn×n，

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

[2]曹蓉,童細心.求解Toeplitz系統(tǒng)循環(huán)和反循環(huán)分裂算法[J].云南民族大學學報:自然科學版,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].西南民族大學學報:自然科學版,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ù)學理論與應用,2012,32(1):1-5.

[9]趙敏.塊Toeplitz矩陣的一種快速Q(mào)R分解及算法實現(xiàn)[J].長江大學學報:自然科學版,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三角線性方程組求解的復雜性[J].數(shù)學研究與評論,1989,9（1）:101-106.

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

[14]張凱院,徐仲.數(shù)值代數(shù)[M].2版.北京:科學出版社,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].北京:清華大學出版社,2004:330-354.

責任編輯：畢和平

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è)技術學院科研基金資助項目（SZK2014Y35）