唐 達

（上海電機學院 數(shù)理教學部，上海201306）

在許多科學與工程計算中，常需要計算矩陣的特征值。而對于對稱矩陣，一般是將其約化為三對角矩陣后再求其特征值；另外，三對角矩陣的特征問題也常作為原始問題出現(xiàn)。因此，三對角矩陣的特征問題長期來為許多學者所關注，如美籍華人數(shù)學家顧明關于三對角矩陣分－治算法的研究成果［1］在SIAM第六屆應用線性代數(shù)會議上獲最佳論文。

當今，計算三對角矩陣特征問題的方法很多，有QL或 QR 算法［2］372－373［3－4］、二（多）分法［2］391－397［5－9］、分－治算法［2］391－397［1，7，10－11］、同 倫 法［7，12］以 及 其 他 的一些迭代算法［7，13］。其中，二（多）分法、分－治算法、同倫法等都能用于并行計算。

一般來講，計算一個n階三對角矩陣的n個特征對，其計算量總不小于O（n2）；而本文利用一類三對角矩陣特征對的局限性質(zhì)，來計算大型非對稱矩陣，其計算量僅為O（n）。本文的算法也非常適合并行計算，其并行效率較高。

1 三對角矩陣t向量的衰減性質(zhì)

設A為n階實三對角矩陣，t＝（t1，t2，…為n維實向量。若

式（1）中右端向量僅第n個分量w不為零，則稱t為A所對應的t向量[14-15]?？梢宰C明，對角占優(yōu)三對角矩陣（i＝2，3，…，n）。也就是說，t向量之分量的絕對值是隨i的減小而衰減的。實際上，除對角占優(yōu)矩陣外，有很大一類三對角矩陣具有這種衰減性質(zhì)。

考慮元素為正態(tài)分布N（0，1）上的隨機三對角矩陣A，其t向量的平均增長率為

設n＝120，進行104次的計算，即樣本容量為104，由此繪出的的頻率密度直方圖[16]如圖1所示。由 MATLAB函數(shù) mean及var可算得此時樣本的均值為1.69，方差為0平.0均3 11 6／1。.6由9此的可速知率，衰在減這。種情況下將以

圖1 頻率密度直方圖Fig.1 Frequency histogram

若隨機矩陣A的元素為（0，1）上均勻分布的隨機數(shù)，則用與上述相同方法可算得的均值為1.40，方差為0.016 5，故此時t向量旳衰減率比上述要低。

以上所述，t向量是從矩陣A的左上方計算到右下方。若t向量是從矩陣A的右下方計算到左上方，即式（1）右端向量的第1個分量變?yōu)閣、第n個分量變?yōu)榱?，由于A是隨機矩陣，故此時的均值及方差均不變。

2 一類三對角矩陣特征對的局限性質(zhì)

仍以元素服從N（0，1）分布的隨機三對角矩陣A為例，設n＝1 000。若任意截取A中一段n1＝200階的子矩陣A1，則稱A為主矩陣、A1為子矩陣。圖2所示為A1的2個特征向量之分量的絕對值或模（取對數(shù)值）隨分量下標i變化的曲線。特征向量之分量的絕對值或模也同t向量一樣具有相同的衰減率。因此，由圖2可見，特征向量曲線在其最大值處向兩旁逐漸減小。其中曲線1在子矩陣A1的兩端已衰減，接近于零，故曲線1所對應的A1的特征對也可近似地看做主矩陣的特征對，這就是特征對的局限性質(zhì)。此時，要計算主矩陣A的特征對，可以由計算子矩陣A1的特征對取得；由于截斷，曲線2在子段的左端特征向量的分量并不趨于零，故該曲線對應的特征對不能看做主矩陣的特征對。由此可見，子矩陣A1的特征對只有一部分是屬于主矩陣A的。數(shù)值實驗表明，此時A1中約有82%的特征值是屬于A的。子矩陣的階數(shù)越大，百分比越大；反之，越小。

圖2 特征向量曲線Fig.2 Curve of eigenvector

除上述呈正態(tài)分布的隨機矩陣外，服從均勻分布、t分布等的隨機三對角矩陣及Wilkinson測試矩陣等其特征對也都具局限性質(zhì)。在特征對具局限性質(zhì)的三對角矩陣中，隨機矩陣占大多數(shù)，而當今隨機矩陣無論在理論上還是實用上都具有一定的價值[17-18]。

3 分段算法

本文將特征對具有局限性質(zhì)的主矩陣A分成互相有重疊的k段，以形成子矩陣A1，A2，…，Ak；再從這些子矩陣的特征對中分離出屬于主矩陣A的特征對，由此求得A的全部特征對。

從子矩陣的特征對中分離出n個屬于主矩陣的特征對，而其計算量又要小于O（n2），是一件較困難的事。所用的軟件功能不同，分離的方法也不盡一樣。本文采用Matlab軟件進行分離。

在圖2中，可以把特征向量之分量絕對值或模最大處的下標i形象地看成其對應的特征值的位置。因此，分離特征值就是把每個子矩陣（子段）兩端的特征值剔除，留下中間屬于主矩陣的特征值。

假定將主矩陣A分成A1、A2、A33段有互相重疊部分的子矩陣。每段長度（階數(shù)）為1.5m，其中，m為分段長參數(shù)，重疊部分為0.5m。再將主矩陣分成有互相重疊部分的B1、B22段，每段長度（階數(shù)）為2.5m。主矩陣的長度（階數(shù)）為3.5m，如圖3所示。m的值不能太小，否則，分離出的特征值誤差太大。如對于服從N（0，1）分布的隨機三對角矩陣，m不應小于120。

圖3 矩陣分段Fig.3 Piecewise matrix

先求出子矩陣A1、A2、B1的全部特征值，然后找出B1的特征值與A1、A2中數(shù)值最接近的2.5m個特征值，令這些特征值的集合為Ω1。Ω1中屬于A1的特征值也應屬于主矩陣A，記此集合為Λ1，則Ω1中屬于A2的特征值的集合Λ2＝Ω1-Λ1。再求出子矩陣A2、A3與B2的全部特征值，找出B2的特征值與A2、A3中數(shù)值最接近的2.5m個特征值，記此集合為Ω2。Ω2中屬于A3的特征值也應屬于主矩陣A，記此集合為Λ4。Ω2中屬于A2的特征值的集合Λ3＝Ω2-Λ4。再求Λ2與Λ3的交集Ψ1＝Λ2∩Λ3。Ψ1中的特征值也應屬于主矩陣A。因此，集合Ψ＝Λ1＋Ψ1＋Λ4中的特征值就是主矩陣A中的3.5m個特征值。

Ψ中的特征值所對應的子矩陣的特征向量并不能作為主矩陣A的特征向量，因為特征向量衰減得較慢，若以此作為主矩陣的特征向量將帶來較大誤差。解決的辦法是將每段1.5m的長度再向兩端延長一些，在延長的子矩陣中求特征向量，這就可以近似地看作主矩陣的特征向量。

4 分段算法程序

假定將主矩陣A分段成k個子矩陣A1，A2，…，Ak，則有如上文所述特征值的計算層次，如圖4所示。

圖4 計算層次Fig.4 Computational hierarchy

先計算Ωi，將Ωi中屬于2個子矩陣的特征值分成Λ2i-1、Λ2i2個集合，求Λ2i、Λ2i＋1的交集

Ψi＝Λ2i∩Λ2i＋1，i＝1，2，…，k-2

最后，集合

中的特征值就是主矩陣A的全部特征值。

算法程序是用MATLAB語言編寫的，其中相關的函數(shù)參見文獻[19] 。程序步驟如下：

（1）輸入矩陣A、分段長參數(shù)m、分段數(shù)k、矩陣A的階數(shù)n＝（k＋0.5）m以及為了精確計算特征向量而將子段端部延長的數(shù)值Δm；

（2）將A以1.5m的長度（階數(shù)）分段成A1、A2、…、Akk個子矩陣，再以2.5m的長度分段成B1、B2、…、Bk-1k-1個子矩陣，各矩陣間的相互位置如圖3所示；

（3）計算A1、A2及B1的特征值，利用dsearchn函數(shù)求出集合Ω1，利用find函數(shù)求出Λ1及Λ2；

（4）將子矩陣A1的右端延長Δm，以集合Λ1中的特征值用逆迭代法求其特征向量；

（5）Fori＝2∶k－1；

（6）計算Ai、Ai＋1、Bi的特征值，求出集合Ωi、Λ2i－1及Λ2i；

（7）利用intersect函數(shù)求出Λ2i－2、Λ2i－1的交集Ψi－1；

（8）將子矩陣Ai的兩端延長Δm，以集合Ψi－1中的特征值用逆迭代法求其特征向量；

（9）end

（10）計算集合Λ2k－2；

（11）將子矩陣Ak的左端延長Δm，以集合Λ2k－2中的特征值用逆迭代法求其特征向量；

（12）輸出用式（3）計算得到主矩陣A的n個特征值，以上算得的特征向量即為A的全部特征向量。

5 數(shù)值算例

本文算例使用的PC機CPU為AMD 2.59GB，內(nèi)存2GB。用Matlab 7.12.0語言編程計算，共計算3例。

例1 元素服從N（0，1）分布的隨機三對角矩陣，為非對稱矩陣，其特征對可能是復數(shù)。其中，m＝120，Δm＝0.45m。

例2 元素服從N（0，1）中均勻分布的隨機三對角矩陣，為對稱矩陣。其中，m＝160，Δm＝0.45m。

例3 Wilkinson測試矩陣，其非對角元素均為1，對角元素依次為n／2＋1－i（i＝1，2，…，n）。其中，m＝100，Δm＝0.45m。

計算結果列于表1中。每個子矩陣的特征值用Matlab eig函數(shù)計算，而特征向量則用逆迭代法計算得到。例1與例2每例均重復計算10次，計算時間為10次的平均值；例3只計算1次。表中，ε為本文算法分段算得的特征值與Matlab eig函數(shù)算得的特征值之差的絕對值或模；e為本文算法算得的特征向量與Matlab eig函數(shù)算得的特征向量其對應的各分量之差的絕對值或模的最大值。

表1 算例的計算精度及計算時間Tab.1 Computation precision and consumed time for the examples

由表可見，計算精度是滿意的，而計算時間僅供參考。由于同一個問題用Matlab內(nèi)部函數(shù)計算要比用Matlab語言編程計算的速度約快300%～1 000%。但從表中可見，本文算法的計算用時與矩陣階數(shù)n大致呈線性關系，而Matlab是用滿矩陣求特征對的，故計算用時約與n3成比例。一般來講，用QR迭代法或二分法計算三對角矩陣的特征值，其計算量為O（n2）。

6 結 語

本文所述分段快速算法適合于大型非對稱三對角矩陣，對n個特征對的計算復雜性僅為O（n），這是目前文獻中不曾見的，只要矩陣的特征向量具局限性質(zhì)，就能采用本文算法。從理論分析及算例均可見，計算具有較高精度。本文算法也非常適合并行計算。

[1] Gu Ming，Eisenstat S C.A divide-and-conquer algorithm for the symmetric tridiagonal eigenproblem[J] .SIAM Journal on Matrix Analysis and Application，1995，16（1）：172-191.

[2] Golub G H，Van Loan C F.矩陣計算[M] .袁亞湘譯.3版.北京：人民郵電出版社，2011：372-373，391-397.

[3] 何光渝，高永利.Visual Fortran常用數(shù)值算法集[M] .北京：科學出版社，2002：357-363.

[4] 何渝.計算機常用數(shù)值算法與程序[M] .北京：人民郵電出版社，2003：137-140.

[5] Wilkinson J H.代數(shù)特征值問題[M] .石鐘慈，鄧健新，譯.北京：科學出版社，2001：310-318.

[6] 曹志浩，張玉德，李瑞遐.矩陣計算和方程求根[M] .2版.北京：高等教育出版社，1984：159-177.

[7] 鄧健新.解實對稱矩陣特征值問題的并行算法[J] .數(shù)值計算與計算機應用，1997，18（4）：246-258.

[8] 劉艷紅，呂全義.Jacobi矩陣特征值的并行算法[J] .紡織高?；A科學學報，2011，24（1）：21-25.

[9] Roopamala T D，Katti S K.Eigenvalue of tridiagonal matrix using Strum sequence and Gerschgorin theorem[J] .International Journal on Computer Science and Engineering，2011，3（12）：3722-3727.

[10] 羅曉廣，李曉梅.求解對稱三對角矩陣特征值的一種新的分而治之算法[J] .數(shù)值計算與計算機應用，1997，18（1）：74-80.

[11] Coakley Ed S，Rokhlin V.A fast divide-and-conquer algorithm for computing the spectra of real symmetric tridiagonal matrices[J] .Applied and Computational Harmonic Analysis，2013，34（3）：379-414.

[12] Brockman P，Carson T，Cheng Yun，et al.Homotopy method for the eigenvalues of symmetric tridiagonal matrices[J] .Journal of Computational and Applied Mathematics，2013，237（1）：644-653.

[13] Matsekh A M.The Godunov-inverse iteration：A fast and accurate solution to the symmetric tridiagonal eigenvalue problem[J] .Applied Numerical Mathematics，2005，54（2）：208-221.

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

[15] 唐達.周期三對角矩陣求逆的快速算法[J] .上海電機學院學報，2013，16（5）：300-304.

[16] 何正風.MATLAB概率與數(shù)理統(tǒng)計分析[M] .2版.北京：機械工業(yè)出版社，2012：87.

[17] 曾杏元.有關隨機矩陣領域最新研究動態(tài)與進展的綜述報告[J] .數(shù)學理論與應用，2011，31（3）：7-19.

[18] 王磊，鄭寶玉，崔景伍.隨機矩陣理論與無線通訊[J] .南京郵電大學學報：自然科學版，2010，30（3）：90-96.

[19] 徐東艷，孟曉剛.MATLAB函數(shù)庫查詢辭典[M] .北京：中國鐵道出版社，2006：161，216，321.