周期三對角矩陣求逆的快速算法

唐 達

(上海電機學院 數理教學部， 上海 201306)

周期三對角矩陣求逆的快速算法

唐 達

(上海電機學院 數理教學部， 上海 201306)

利用t向量來求周期三對角矩陣之逆。求逆的運算量為2n2+O(n)乘除法及n2+O(n)加減法。該算法計算量小且計算精度高。若對t向量進行截斷、快速求逆，則求逆的計算量僅與n成正比。與現有快速算法相比，清除了電腦內存溢出的情況。文末列出了部分數值算例。

周期三對角矩陣； 逆矩陣； 溢出；t向量； 快速求逆

矩陣

(1)

稱為周期三對角矩陣。式中，n為矩陣的階；在許多科學與工程計算中，如在3次樣條周期邊界插值、某些拋物型方程差分解、常微分方程周期邊值問題的差分解、時間序列分析、圖論的特征值分析等領域中，常會使用周期三對角矩陣。因此，對這類矩陣的計算一直是學術界的熱點之一；周期三對角陣求逆[1-13]也不例外。文獻[1]中所述的解析逆，若用于計算逆陣的全部元素，則其計算量較大；文獻[2]中所述的求逆方法，其計算量為3n2+O(n)乘除法及2n2+O(n)加減法；文獻[3]中所述求逆方法，其計算量為3n2+O(n)乘除法及n2+O(n) 加減法；文獻[5]中提出了4種求解T-1的算法，其最小的計算量也與文獻[3]中的相同；文獻[6]中所述求逆的計算量為2n2+O(n)乘除法及n2+O(n)加減法，但此算法可能引起電腦計算溢出；文獻[7]中提出的算法計算量與文獻[6]中的相仿，但亦可能引起電腦計算上溢或下溢；文獻[8-9]中三對角矩陣求逆的算法同樣可能溢出。

本文利用t向量來求周期三對角矩陣之逆，其計算量為2n2+O(n)乘除法及n2+O(n)加減法，釆用了防止電腦計算溢出的措施。在上述諸多算法中，本文的算法計算量小，算法穩定，且具有較高的精度。若將t向量截斷，進行快速求逆，則求逆的計算量僅與n成正比。

1 周期三對角陣求逆

1.1方法概述

在式(1)中，若a1=cn=0，則周期三對角矩陣退化為三對角矩陣，令其為T1。現討論以下方程的解：

T1t=ej

(2)

式中，ej為n階單位陣的第j列；t為n維解向量。

為了求解方程(2)，可令

(3)

(4)

(5)

式中，B的第1列為t(1)的尾分量；B的第2列為t(2)的尾分量。B的行相應位于T1的第j及j+1行上。解二階方程(x1、x2為方程的未知數)：

(6)

為了求T的逆陣，令

T1t(3)=en，T1t(4)=e1

(7)

(8)

解出式(8)后，向量

(9)

就是T-1的第j列列向量。令j=1～n，即可求出T-1。

1.2計算量估計

求解式(9)的計算量最大，基本上為整個求逆的計算量。t(1)為t(3)的一段(或一部分)，它們間僅相差一個比例常數。令t(3)與t(1)對應分量之比為k1。計算式(9)時，t(3)可與t(1)對應分量合并計算：

(k1x1+x3)t(1)

同理，t(4)可與t(2)對應分量合并計算：

(k2x2+x4)t(2)

其中，k2為t(4)與t(2)對應分量之比。

上述2段的計算量共需n2+O(n)次乘法；t(3)剩下的一段乘以x1，t(4)剩下的一段乘以x2，共需n2乘法；兩者相加需n2+O(n)加減法；因此，計算T-1共需2n2+O(n)乘除法及n2+O(n)加減法。

1.3誤差分析

利用式(3)或式(4)計算t向量時，相當于解具有3條對角線的下三角或上三角方程，而解三角形方程組具有較高的相對精度[15-16],故利用式(3)或式(4)計算t向量具有較高的相對精度。求T逆的導出方程只有4階，其解不會產生太大的誤差，故利用本文方法求T-1具有較高精度。

2 防止內存溢出的方法

(10)

(11)

令

(12)

(13)

(14)

用s(1)、s(2)代替t(1)、t(2)后，求T逆的導出方程的系數陣為

(15)

(16)

(17)

從而可求出

如果矩陣T1可均，則根據文獻，用本文算法亦可方便地求出T-1，且其計算量比不可約的矩陣小。

3 快速求逆

絕大多數常見的三對角矩陣，其t(1)向量之分量的絕對值是按行序增大的方向遞增的；則按行序減小的方向，其絕對值就是衰減的(對t(2)向量也有類似結論)；因此，T-1的元素也具有衰減的性質。文獻[17-19]中從不同角度也闡述了逆矩陣元素的衰減性質。在利用式(9)計算T-1的元素時，當t向量之分量衰減到非常小(如相鄰兩分量的絕對值均小于10-18)時，就對t向量進行截斷，將其以后的分量視為零而不予計算，這樣做不會影響求逆的精度。因此，求逆的計算量就由與n2成正比減少為與n成正比。此時，T-1是一稀疏矩陣，其非零元的數目不是n2，而是與n成正比。文獻[20]中也提出用截斷方法來快速求解周期三對角線性方程組，但對截斷的條件要求比較苛刻，其實只要t向量具有衰減性，就可根據情況截斷。

4 算法歩驟

周期三對角矩陣求逆算法如下：

(1) 輸入矩陣T；

(2) 根據式(11)計算s(1)及s(2)；

(3) 根據式(16)計算ui及vi，i=1，2,…,n；

(4) forj=2 step 2 ton-1；

(5) 列出式(15)求逆的導出方程的系數陣，其右端項為(0,1,0,0)T及(0,0,1,0)T；

(6) 計算導出方程之解x1,x2,x3,x4；

(7) 按照式(9)及1.2節所述方法計算T-1的第j及j+1列，其中t向量的各分量根據式(13)、(14)遞推求得。在遞推時，若檢測到相鄰兩個t分量之值均小于10-18，則進行截斷，結束遞推；否則，對t(1)、t(3)按行序減小方向一直遞推到T的第1行，對t(2)、t(4)按行序増大方向一直遞推到T的第n行；

(8) nextj；

(9) 計算T-1的第1列及n列，其導出方程只有二階，算法與上相仿；

(10) 輸出T-1。

5 數值算例

本文算例用PC機CPU為AMD 2.59GB，內存2GB。用Matlab 7.1語言編程計算。每一算例均列出本文算法、文獻[3]與文獻[7]中算法、Matlab內部函數Inv求逆算法的比較。若令

εi,j=TT-1-I

(18)

式中,I為n階單位陣；ε為算法誤差。則

ε=max|εi,j|，i,j=1,2,…,n

(19)

各例的計算時間及算法誤差如表1所示。

表1 算例的計算時間及精度

注： *有2例溢出

例1T的元素為正態分布N(0,1)中之隨機數。共計算10次，計算誤差及時間t為10次的平均值。

例2式(1)中的T為Toeplitz周期三對角陣，ai=1,bi==6(i=1,2,…,n),ci=5(i=1,2,…,n-1),cn=1。本例取自文獻[3]中。

例3Toeplitz周期三對角陣，a1=0.5,cn=0.5,ai=1(i=2，3,…,n)，bi=2.01(i=1,2,…,n),ci=1(i=1,2,…,n-1)。

例4Toeplitz周期三對角陣，a1=1,cn=1,ai=2(i=2,3,…,n),bi=2(i=1,2,…,n),ci=-3(i=1,2,…,n-1)。

在表1中，利用MATLAB Inv函數計算T-1時，T應是稀疏矩陣；其求逆的時間與文獻[3]中有較大差別，原因如下： ① 文獻[3]中是采用滿陣求逆的；② 由于所用Matlab版本不同，故計算的時間差別較大。另外，利用Matlab語言編制的程序是解釋程序；若計算量相同，利用Matlab語言編制的程序運行要比其內部函數運行慢得多。

6 結 語

在當今諸多求T-1的算法中，本文的算法計算時間較短，特別對于大型矩陣是這樣；計算精度也較高；同時，消除了電腦計算溢出的弊端。因此，本文的算法是一種較好的求T-1的實用算法。

[1] EI-Shehawey M A,EI-Shreef Gh A,AI-Henawy A Sh.Analytical inversion of general periodic tridiagonal Matrices[J].Journal of Mathematical Analysis and Application,2008,345(1): 123-134.

[2] 余承依，陳耀輝，趙立群.求三對角和周期三對角矩陣逆陣的-種新算法[J].長江大學學報： 自然科學版,2010,7(1)： 126-128.

[3] 余承依,陳耀輝.周期三對角矩陣逆的一種新算法[J].數學的實踐與認識,2010,40(22)： 192-198.

[4] 陳耀輝,趙立群,余承依.周期三對角矩陣行列式和逆矩陣新算法： 英文[J].數學研究,2011,44(1)： 27-35.

[5] 趙立群.一些稀疏矩陣的逆和行列式的計算[D].漳州： 漳州師范學院，2011.

[6] 袁志杰.有關特殊矩陣的計算問題及性質[D].西安： 西北工業大學，2005.

[7] 杜永恩,陸 全,徐 仲.求分塊周期三對角矩陣逆矩陣的新算法[J].計算機工程與應用，2012，48(17)： 41-43.

[8] 冉瑞生,黃庭祝.三對角矩陣的逆[J].哈爾濱工業大學學報,2006,38(5)： 815-817.

[9] 冉瑞生,黃庭祝,劉興平,等.三對角矩陣求逆的算法[J].應用數學和力學,2009,30(2)： 238-244.

[10] 車 毅,徐 仲,雷小娜.分塊周期三對角矩陣逆矩陣的新算法[J].紡織高校基礎科學學報,2011,24(1)： 15-20.

[11] 陳 芳,徐 仲.求分塊三對角矩陣和分塊周期三對角矩陣逆陣的快速算法[J].數學的實踐與認識,2005,32(12)： 183-187.

[12] 黃卓紅,劉建州.嚴格對角占優周期三對角矩陣逆元素的上下界估計[J].工程數學學報,2008,25(4)： 703-707.

[13] 袁志杰,徐 仲.嚴格對角占優周期三對角矩陣逆元素的上界估計[J].工程數學學報,2004,21(8)： 68-72.

[14] 唐 達.三對角矩陣計算[J].高等學校計算數學學報,1997,19(2)： 97-104.

[15] Willkinson J H.代數特征值問題[M].石鐘慈,鄧健新,譯.北京： 料學出版社,2001： 259-261.

[16] 唐 珍.舍入誤差分析引論[M].上海： 上海科學技術出版社,1987： 207-209.

[17] Nabben R.Decay rates of the inverse of nonsymmetric tridiagonal and band matrix[J].Journal on Matrix Analysis and Applications,1999,20(3)： 820-837.

[18] Meurant G.A review on the inverse of symmetric tridiagonal and block tridiagonal matrices[J].Journal on Matrix Analysis and Applications,1992, 13(3)： 707-728.

[19] 雷光耀.強主元稀疏陣的近似求逆[J].科學通報,1991,36(8)： 572-574.

[20] 季 青，王能超.解循環三對角線性方程組的追趕法[J].小型微型計算機系統,2002,23(11)： 1393-1395.

Fast Inversion of Periodic Tridiagonal Matrices

TANGDa

(Department of Mathematics and Physics, Shanghai Dianji University, Shanghai 201306, China)

In this paper, inversion of periodic tridiagonal matrices is discussed using thetvector. Computational complexity of the inversion is 2n2+O(n) of multiplication and division andn2+O(n) of addition and subtraction. The algorithm has a low computational cost as well as high precision. The computational cost is further reduced to become proportional tonif thetvector is truncated and fast inversed. In contrast to the existing fast algorithms, the out-of-memory errors no longer occur. Numerical examples are presented.

periodic tridiagonal matrix; inverse matrix; overflow;t-vector; fast inversion

2095-0020(2013)05 -0300-05

O 241.6

A

2013-08-19

唐 達(1935-)，男，工程師，主要研究方向為工程數值計算，E-mail: tangdayouxiang@sina.cn