一種構(gòu)造多項(xiàng)式求解特征值問題的方法

李 宏，陳志寶

(華南理工大學(xué)數(shù)學(xué)系, 中國 廣州 510640)

矩陣的特征值問題,無論在理論上還是在工程計(jì)算中,都是一個(gè)十分重要的問題.到目前為止人們已經(jīng)對它進(jìn)行了深入的討論,得到了許多卓有成效的算法, 常用的主要有兩類, 一類是正交相似變換的方法稱為變換法如經(jīng)典的Jacobi方法、QR方法[1-3]等,變換法是用來計(jì)算矩陣全部特征值的方法,它的理論依據(jù)是將A的特征值問題轉(zhuǎn)化為容易求解的矩陣B的特征值問題,而轉(zhuǎn)化過程是采用正交相似變換, QR方法是目前計(jì)算中小型矩陣全部特征值問題的最有效的方法之一，具有收斂快，算法穩(wěn)定等特點(diǎn),matlab中自帶求特征值的指令就是依據(jù)這一原理編制而成；另一類是求解大型矩陣特征值問題的子空間迭代法[4-5]、對稱矩陣的Lanczos方法[6]等.這兩種方法是用來計(jì)算矩陣若干個(gè)端部特征值的方法，文[6]指出塊Chebyshev-Lanczos方法是目前求解大型稀疏對稱矩陣部分極端特征值的最有效方法.受Krylov子空間方法[6-7]思想的啟發(fā)，本文介紹一種利用線性方程組的解構(gòu)造多項(xiàng)式,通過求解該多項(xiàng)式的根從而獲得矩陣全部或部分特征值并且利用該多項(xiàng)式做除法求得對應(yīng)特征值的特征向量的方法.

1 相關(guān)定理及結(jié)果

由高等代數(shù)知識，容易得到如下兩個(gè)結(jié)果.

定理A設(shè)A∈Rn×n,r0∈Rn×1,若能找到一個(gè)次數(shù)最小的m次非零多項(xiàng)式fm(x),使得fm(A)r0=0,則fm(x)=0的每個(gè)根一定都是矩陣A的特征值.

定理B設(shè)A∈Rn×n且其秩為m,r0為A的列向量,則存在非零m次多項(xiàng)式fm(x),使得fm(A)r0=0.

注記：A的特征向量的一種新求法：

2 算法

定理A表明：假如能夠找到次數(shù)最小的非零多項(xiàng)式fm(x),使得fm(A)r0=0,那么矩陣A的全部或部分特征值可以通過求解fm(x)的根獲得.定理B表明:這樣的非零多項(xiàng)式總是存在的,且其次數(shù)不超過A的秩.于是得到了下述求解矩陣特征值和特征向量的算法.

算法步驟:

(1) 取A的某一列向量r0,計(jì)算Ar0,…,Anr0,形成矩陣K=(r0,Ar0,…,Anr0).

(2) 對矩陣K做LU分解.

(3) 若記U的秩為m,則求解線性方程組U(1:m,1:m)x=U(1:m,m+1).

(4) 得到多項(xiàng)式fm(z)的系數(shù)(見下面的注明),用代數(shù)的方法求該多項(xiàng)式的根作為A的特征值.

(5) 用注記中的方法求對應(yīng)特征值的特征向量.

注上述算法中,U(1:m,1:m)表示U的前m行m列組成的矩陣，U(1:m,m+1)表示U的第m+1列中的前m個(gè)元素組成的列向量.若記步驟3中線性方程組的解為x=(a0,a1,…,am-1)T,則fm(x)=a0+a1x+…+am-1xm-1-xm.

計(jì)算量分析：用上述方法計(jì)算A的特征值需要計(jì)算Ar0,…,Anr0,這需要n次矩陣與向量的乘積運(yùn)算(相當(dāng)于求A2的計(jì)算量),然后求解一個(gè)齊次線性方程組的計(jì)算量,最后是代數(shù)的方法求多項(xiàng)式的根的計(jì)算量.特征向量的計(jì)算量從注記中可以看出其幾乎依賴于Ar0,…,Anr0的計(jì)算和fm(x)的根求解,并不需要大量額外計(jì)算量.

3 數(shù)值算例

算例1利用上述算法計(jì)算下述A的特征值和特征向量:

取r0=A(:,1)=(1/2,1,-1,3/2,-3/2)T,算得

U的秩是5，解得x=(-139.84,183.81,-63.75,-7.5,7.5)T,從而

fm(x)=-139.84+183.81x-63.75x2-7.5x3+7.5x4-xm,

fm(x)的5個(gè)根為λ1=-3.211 3,λ2=1.700 7,λ3=2.285 9,λ4=3.039 6,λ5=3.685 2.

用x-λ1對多項(xiàng)式fm(x)做除法得到f1(x)=-x4+10.711 3x3-41.897 9x2+70.798 7x-43.546 7,則對應(yīng)λ1的特征向量f1(A)r0=K(1∶5,1∶5)(-43.546 7,70.798 9,-41.897 9,10.711 3,-1.000)T=(0.555 7,-0.463 9,0.387 7,-0.401 2,0.405 9)T(已單位化處理,下同),同理解得

f2(A)r0=(0.679 3,0.686 4,-0.200 2,-0.138 3,-0.090 9)T,

f3(A)r0=(-0.222 6,0.434 4,0.858 6,-0.101 5,-0.119 2)T,

f4(A)r0=(0.331 1,-0.057 7,0.240 9,0.899 7,0.140 0)T,

f5(A)r0=(-0.265 8,0.348 7,-0.120 1,0.013 8,0.890 6)T.

這與表1 matlab自帶求特征值和特征向量指令得到的結(jié)果一樣.

表1 matlab指令求得的特征值和特征向量

算例2考慮100階魔方陣A=magic(100),求A的特征值.

由于rank(A)=3,根據(jù)定理B,只需要形成K=(r0,Ar0,A2r0,A3r0)，求得

fm(x)=-4.166 7×1014+8.332 5×108x+5.000 5×105x2-x3，

fm(x)的3個(gè)根為λ1=50 005.000 0,λ2=-28 866.070 0,λ3=28 866.070 0,這與matlab指令求得的3個(gè)特征根λ1=0.500 050 000×105,λ2=-0.288 660 700×105,λ3=0.288 660 700×105也一樣.

通過上面兩個(gè)算例，驗(yàn)證了作者提出的利用線性方程組的解構(gòu)造多項(xiàng)式從而求解矩陣特征值的方法的有效性和準(zhǔn)確性,此方法對于低階或者高階低秩矩陣無疑是適用的，它運(yùn)算量少并且易于編程實(shí)現(xiàn).

致謝:作者對審稿人的建議和評論表示衷心感謝.

參考文獻(xiàn)：

[1] 曹志浩. 數(shù)值線性代數(shù)[M].上海：復(fù)旦大學(xué)出版社, 1996.

[2] GOLUB G H， VAN LOAN G F. 矩陣計(jì)算[M].袁亞湘，譯.北京：科學(xué)出版社, 2001.

[3] WILKINSON J H. The algebraic eigenvalue problem[M]. Oxford：Oxford University Press, 1965.

[4] 用Chebyshev多項(xiàng)式加速的子空間迭代法[J].南京航空航天大學(xué)學(xué)報(bào)，2002，4(2)：197-199.

[5] 曲慶國. 對求解大型稀疏特征值問題的子空間迭代法的研究[D].南京：南京航空航天大學(xué)理學(xué)院，2005.

[6] 戴 華. 求解大規(guī)模矩陣問題的Krylov子空間方法[J].南京航空航天大學(xué)學(xué)報(bào)，2001，4(2)：139-145.

[7] SIMONCINI V, SZYLD D B. Recent computational developments in Krylov subspace methods for linear systems[J]. Numerical Linear Algebra with Applications, 2007, 14：1-59.