唐 達(dá)

(上海電機(jī)學(xué)院 數(shù)理教學(xué)部， 上海 201306)

周期三對(duì)角矩陣求逆的快速算法

唐 達(dá)

(上海電機(jī)學(xué)院 數(shù)理教學(xué)部， 上海 201306)

利用t向量來(lái)求周期三對(duì)角矩陣之逆。求逆的運(yùn)算量為2n2+O(n)乘除法及n2+O(n)加減法。該算法計(jì)算量小且計(jì)算精度高。若對(duì)t向量進(jìn)行截?cái)?、快速求逆，則求逆的計(jì)算量?jī)H與n成正比。與現(xiàn)有快速算法相比，清除了電腦內(nèi)存溢出的情況。文末列出了部分?jǐn)?shù)值算例。

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

矩陣

(1)

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

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

1 周期三對(duì)角陣求逆

1.1方法概述

在式(1)中，若a1=cn=0，則周期三對(duì)角矩陣退化為三對(duì)角矩陣，令其為T(mén)1?，F(xiàn)討論以下方程的解：

T1t=ej

(2)

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

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

(3)

(4)

(5)

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

(6)

為了求T的逆陣，令

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

(7)

(8)

解出式(8)后，向量

(9)

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

1.2計(jì)算量估計(jì)

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

(k1x1+x3)t(1)

同理，t(4)可與t(2)對(duì)應(yīng)分量合并計(jì)算：

(k2x2+x4)t(2)

其中，k2為t(4)與t(2)對(duì)應(yīng)分量之比。

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

1.3誤差分析

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

2 防止內(nèi)存溢出的方法

(10)

(11)

令

(12)

(13)

(14)

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

(15)

(16)

(17)

從而可求出

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

3 快速求逆

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

4 算法歩驟

周期三對(duì)角矩陣求逆算法如下：

(1) 輸入矩陣T；

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

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

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

(5) 列出式(15)求逆的導(dǎo)出方程的系數(shù)陣，其右端項(xiàng)為(0,1,0,0)T及(0,0,1,0)T；

(6) 計(jì)算導(dǎo)出方程之解x1,x2,x3,x4；

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

(8) nextj；

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

(10) 輸出T-1。

5 數(shù)值算例

本文算例用PC機(jī)CPU為AMD 2.59GB，內(nèi)存2GB。用Matlab 7.1語(yǔ)言編程計(jì)算。每一算例均列出本文算法、文獻(xiàn)[3]與文獻(xiàn)[7]中算法、Matlab內(nèi)部函數(shù)Inv求逆算法的比較。若令

εi,j=TT-1-I

(18)

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

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

(19)

各例的計(jì)算時(shí)間及算法誤差如表1所示。

表1 算例的計(jì)算時(shí)間及精度

注： *有2例溢出

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

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

例3Toeplitz周期三對(duì)角陣，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周期三對(duì)角陣，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函數(shù)計(jì)算T-1時(shí)，T應(yīng)是稀疏矩陣；其求逆的時(shí)間與文獻(xiàn)[3]中有較大差別，原因如下： ① 文獻(xiàn)[3]中是采用滿陣求逆的；② 由于所用Matlab版本不同，故計(jì)算的時(shí)間差別較大。另外，利用Matlab語(yǔ)言編制的程序是解釋程序；若計(jì)算量相同，利用Matlab語(yǔ)言編制的程序運(yùn)行要比其內(nèi)部函數(shù)運(yùn)行慢得多。

6 結(jié) 語(yǔ)

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

[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] 余承依，陳耀輝，趙立群.求三對(duì)角和周期三對(duì)角矩陣逆陣的-種新算法[J].長(zhǎng)江大學(xué)學(xué)報(bào)： 自然科學(xué)版,2010,7(1)： 126-128.

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

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

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

[6] 袁志杰.有關(guān)特殊矩陣的計(jì)算問(wèn)題及性質(zhì)[D].西安： 西北工業(yè)大學(xué)，2005.

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

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

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

[10] 車 毅,徐 仲,雷小娜.分塊周期三對(duì)角矩陣逆矩陣的新算法[J].紡織高?；A(chǔ)科學(xué)學(xué)報(bào),2011,24(1)： 15-20.

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

[12] 黃卓紅,劉建州.嚴(yán)格對(duì)角占優(yōu)周期三對(duì)角矩陣逆元素的上下界估計(jì)[J].工程數(shù)學(xué)學(xué)報(bào),2008,25(4)： 703-707.

[13] 袁志杰,徐 仲.嚴(yán)格對(duì)角占優(yōu)周期三對(duì)角矩陣逆元素的上界估計(jì)[J].工程數(shù)學(xué)學(xué)報(bào),2004,21(8)： 68-72.

[14] 唐 達(dá).三對(duì)角矩陣計(jì)算[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào),1997,19(2)： 97-104.

[15] Willkinson J H.代數(shù)特征值問(wèn)題[M].石鐘慈,鄧健新,譯.北京： 料學(xué)出版社,2001： 259-261.

[16] 唐 珍.舍入誤差分析引論[M].上海： 上海科學(xué)技術(shù)出版社,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] 雷光耀.強(qiáng)主元稀疏陣的近似求逆[J].科學(xué)通報(bào),1991,36(8)： 572-574.

[20] 季 青，王能超.解循環(huán)三對(duì)角線性方程組的追趕法[J].小型微型計(jì)算機(jī)系統(tǒng),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

唐 達(dá)(1935-)，男，工程師，主要研究方向?yàn)楣こ虜?shù)值計(jì)算，E-mail: tangdayouxiang@sina.cn