李 宏,雷秀仁,陈志宝

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

构造多项式求解矩阵特征值问题的方法

李 宏,雷秀仁*,陈志宝

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

提出了2种构造多项式,通过求该多项式的根,获得矩阵部分特征值的方法,数值算例比较了这2种方法,表明了这2种方法对于求解较小规模矩阵特征值的有效性和准确性.

特征值问题; 最小多项式; 残差多项式; QR分解; 多项式求根

矩阵的特征值问题无论在理论上还是在工程计算中都非常重要.前人进行了深入的讨论,得到了许多卓有成效的算法,除了经典的Jacobi方法、QR方法[1-2]外,还有近几十年来发展起来的子空间迭代法、Krylov方法[3-5]等等,在本质上属于变换法或迭代法.本文则利用线性方程组的解构造向量r0的最小多项式,通过求解该多项式的根从而获得矩阵的部分特征值的方法,另外,当矩阵满秩时,利用GMRES算法中的残差多项式可以得到向量r0的最小多项式的一个非常好的近似,数值算例表明了这2种方法对于求解较小规模矩阵特征值是非常有效和准确的方法.

1 多项式的构造

1.1通过求解线性方程组构造多项式

为求得向量r0的最小多项式，文献[1]给出了一种方法，称为Krylov方法，具体地说是找到最小的正整数m，满足r0,Ar0,…,Amr0线性相关，然后求解

[r0,Ar0,…,Amr0](x0,x1,…,xm-1,1)T=0,

(1)

最后得到r0的最小多项式

f(z)=x0+x1z+…+xm-1zm-1+zm.

(2)

当矩阵A的阶数较大时，向量组r0,Ar0,…,Amr0的线性相关性开始变得模糊，并且Krylov矩阵[r0,Ar0,…,Anr0]不可避免的变得病态，即不容易确定r0的最小多项式次数和精确求解线性方程组(1)，因此我们试图归一化Krylov矩阵的列，以期降低其病态程度，然后对其做QR分解，以R的秩作为r0的最小多项式的次数并求解

(3)

最后得到r0的最小多项式

(4)

注记：A的特征向量的一种求法:记i为多项式f(z)的根,用z-i对f(z)做带余除法(显然这里余数为0)得到fi(z)=f(z)/(z-i),记fi(z)=bm-1×zm-1+bm-2zm-2+…+b0,则fi(A)r0=(r0,Ar0,…,Am-1r0)(b0,b1,…,bm-1)T,显然fi(A)r0即为A的对应特征值i的特征向量.

1.2直接利用GMRES残差多项式

证明由向量的线性相关性可知,存在最小的正整数m(m≤n),满足r0,Ar0,…,Amr0线性相关,从而存在一组非零系数a0,a1,…,am,使得(a0+a1A+…+amAm)r0=0.若a0非零，则结论显然成立;若a0为零，则r0的极小多项式变为a1x+…+amxm,显然该多项式有零根,即A有零特征值,这与条件A为满秩矛盾，故结论成立.

(5)

(6)

其中

gm-1(x)=a0+a1x+…+am-1xm-1,

(7)

注意到GMRES中的Arnoldi迭代

(8)

(9)

Pm(x)=1-a0x-a1x2-…-am-1xm.

(10)

2 算法

利用1.1节的方法可以得到下述求解矩阵特征值和特征向量的算法1.

算法1

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

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

(4)得到多项式f(z)的系数,用代数的方法求该多项式的根作为A的特征值.

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

注明:上述算法中,R(1∶m,1∶m)表示R的前m行m列组成的矩阵,R(1∶m,m+1)表示R的第m+1列中的前m个元素组成的列向量,若步骤3中解得x=(x0,x1,…,xm-1)T.则

利用1.2节的GMRES残差多项式，可以得到下述算法2.

算法2

(2)迭代，对k=1,2,…,m,计算

iv)qk+1=rk/hk+1,k,

v) 按照式(9)计算Cm的第k列.

3 数值试验

考察20阶的对称矩阵

分别用Krylov方法[1]、算法1和算法2求A的特征值,其结果如表1所示.

表1 不同算法求得的特征值比较Table 1 Comparison of eigenvalue obtained by different algorithm

表1中Matlab指令求得的是精确特征值，从表中看出，Krylov方法只能求得18个特征值且与精确特征值存在较大误差，算法1是在Krylov方法的基础上改进而来，能够求出A的所有特征值且与精确特征值之间误差明显减小,算法2是利用GMRES方法中的残差多项式作为向量r0的最小多项式，通过求解残差多项式的根作为矩阵的特征值，用该方法求得的特征值与精确特征值之间的差别很小，较前2种方法更为精确.但是,算法1和算法2虽然能够用来计算矩阵部分特征值，但受到较为严格的限制，对于算法1，正如文献[1]提到的，只有当特征值的分布满足某种条件时才可能应用到较大规模的矩阵中，否则，当矩阵的阶较高时，计算得到的特征值舍入误差的影响而与实际特征值相差很大，对于算法2，矩阵必须为满秩且受舍入误差的影响，步骤2的迭代次数不能过大，即m不能取太大，否则步骤4的残差检验将难以满足，这些因素直接制约了矩阵的阶数不能太大，除非能够找到一个很好的初始向量r0使得步骤2的迭代次数较少且残差检验也能够满足.

[1] WILKINSON J H.The algebraic eigenvalue problem[M].Oxford：Oxford University Press,1965:365-381.

[2] 曹志浩.数值线性代数[M].上海：复旦大学出版社，1996：175-215.

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

[4] 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.

[5] 戴华.求解大规模矩阵问题的Krylov子空间方法[J].南京航空航天大学学报，2001，33(2)：139-145.

DAI Hua.Krylov subspace methods for solving large scale matrix problems[J].Journal of Nanjing University of Aeronautics & Astronautic,2001，33(2)：139-145.

[6] 全忠，向淑晃.基于GMRES的多项式预处理广义极小残差法[J].计算数学,2006，28(4):365-376.

QUAN Zhong,XIANG Shuhuang.A GMRES based polynomial preconditioning algorithm[J].Mathematica Numerica Sinica,2006，28(4):365-376.

Keywords： eigenvalue problem; minimal polynomial; residual polynomial; QR decomposition; polynomial roots

【责任编辑 庄晓琼】

METHODSFORSOLVINGEIGENVALUEPROBLEMSBYUSINGCONSTRUCTEDPOLYNOMIAL

LI Hong,LEI Xiuren*,CHEN Zhibao

(Department of Mathematics, South China University of Technology, Guangzhou 510640, China)

Two methods to construct a polynomial are first introduced,and then the roots of the polynomial are calculated in order to obtain some of the eigenvalues.These two methods are compared by numerical experiments which show that the methods are effective and accurate to solve the eigenvalues of smaller matrix.

2010-04-29

*通讯作者，maxlei@scut.edu.cn

1000-5463(2011)02-0043-03

O241.6

A