一种构造多项式求解特征值问题的方法

李 宏，陈志宝

(华南理工大学数学系, 中国 广州 510640)

矩阵的特征值问题,无论在理论上还是在工程计算中,都是一个十分重要的问题.到目前为止人们已经对它进行了深入的讨论,得到了许多卓有成效的算法, 常用的主要有两类, 一类是正交相似变换的方法称为变换法如经典的Jacobi方法、QR方法[1-3]等,变换法是用来计算矩阵全部特征值的方法,它的理论依据是将A的特征值问题转化为容易求解的矩阵B的特征值问题,而转化过程是采用正交相似变换, QR方法是目前计算中小型矩阵全部特征值问题的最有效的方法之一，具有收敛快，算法稳定等特点,matlab中自带求特征值的指令就是依据这一原理编制而成；另一类是求解大型矩阵特征值问题的子空间迭代法[4-5]、对称矩阵的Lanczos方法[6]等.这两种方法是用来计算矩阵若干个端部特征值的方法，文[6]指出块Chebyshev-Lanczos方法是目前求解大型稀疏对称矩阵部分极端特征值的最有效方法.受Krylov子空间方法[6-7]思想的启发，本文介绍一种利用线性方程组的解构造多项式,通过求解该多项式的根从而获得矩阵全部或部分特征值并且利用该多项式做除法求得对应特征值的特征向量的方法.

1 相关定理及结果

由高等代数知识，容易得到如下两个结果.

定理A设A∈Rn×n,r0∈Rn×1,若能找到一个次数最小的m次非零多项式fm(x),使得fm(A)r0=0,则fm(x)=0的每个根一定都是矩阵A的特征值.

定理B设A∈Rn×n且其秩为m,r0为A的列向量,则存在非零m次多项式fm(x),使得fm(A)r0=0.

注记：A的特征向量的一种新求法：

2 算法

定理A表明：假如能够找到次数最小的非零多项式fm(x),使得fm(A)r0=0,那么矩阵A的全部或部分特征值可以通过求解fm(x)的根获得.定理B表明:这样的非零多项式总是存在的,且其次数不超过A的秩.于是得到了下述求解矩阵特征值和特征向量的算法.

算法步骤:

(1) 取A的某一列向量r0,计算Ar0,…,Anr0,形成矩阵K=(r0,Ar0,…,Anr0).

(2) 对矩阵K做LU分解.

(3) 若记U的秩为m,则求解线性方程组U(1:m,1:m)x=U(1:m,m+1).

(4) 得到多项式fm(z)的系数(见下面的注明),用代数的方法求该多项式的根作为A的特征值.

(5) 用注记中的方法求对应特征值的特征向量.

注上述算法中,U(1:m,1:m)表示U的前m行m列组成的矩阵，U(1:m,m+1)表示U的第m+1列中的前m个元素组成的列向量.若记步骤3中线性方程组的解为x=(a0,a1,…,am-1)T,则fm(x)=a0+a1x+…+am-1xm-1-xm.

计算量分析：用上述方法计算A的特征值需要计算Ar0,…,Anr0,这需要n次矩阵与向量的乘积运算(相当于求A2的计算量),然后求解一个齐次线性方程组的计算量,最后是代数的方法求多项式的根的计算量.特征向量的计算量从注记中可以看出其几乎依赖于Ar0,…,Anr0的计算和fm(x)的根求解,并不需要大量额外计算量.

3 数值算例

算例1利用上述算法计算下述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个根为λ1=-3.211 3,λ2=1.700 7,λ3=2.285 9,λ4=3.039 6,λ5=3.685 2.

用x-λ1对多项式fm(x)做除法得到f1(x)=-x4+10.711 3x3-41.897 9x2+70.798 7x-43.546 7,则对应λ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自带求特征值和特征向量指令得到的结果一样.

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

算例2考虑100阶魔方阵A=magic(100),求A的特征值.

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

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

fm(x)的3个根为λ1=50 005.000 0,λ2=-28 866.070 0,λ3=28 866.070 0,这与matlab指令求得的3个特征根λ1=0.500 050 000×105,λ2=-0.288 660 700×105,λ3=0.288 660 700×105也一样.

通过上面两个算例，验证了作者提出的利用线性方程组的解构造多项式从而求解矩阵特征值的方法的有效性和准确性,此方法对于低阶或者高阶低秩矩阵无疑是适用的，它运算量少并且易于编程实现.

致谢:作者对审稿人的建议和评论表示衷心感谢.

参考文献：

[1] 曹志浩. 数值线性代数[M].上海：复旦大学出版社, 1996.

[2] GOLUB G H， VAN LOAN G F. 矩阵计算[M].袁亚湘，译.北京：科学出版社, 2001.

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

[4] 用Chebyshev多项式加速的子空间迭代法[J].南京航空航天大学学报，2002，4(2)：197-199.

[5] 曲庆国. 对求解大型稀疏特征值问题的子空间迭代法的研究[D].南京：南京航空航天大学理学院，2005.

[6] 戴 华. 求解大规模矩阵问题的Krylov子空间方法[J].南京航空航天大学学报，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.