幂法求解矩阵主特征值的加速方法傅鹏河南理工大学数学与信息科学学院信息与计算科学专业 2009 级 1 班摘要：本论文主要研究的是幂法求解矩阵的主特征值和特征向量。物理、力学和工程技术中有许多需要我们求矩阵的按模最大的特征值（及称为主特征值）和特征向量。幂法是计算一个矩阵的模最大特征值和对应的特征向量的一种迭代方法。它最大的优点是方法简单，适合于大型稀疏矩阵的主特征值，但是收敛速度非常慢。所以我们要用加速的方法来加速收敛，加速方法包括原点平移加速、Rayleigh 商加速和 Aitken 加速算法。关键词：幂法；原点平移加速；Rayleigh 商加速；Aitken 加速算法1 引言我们来介绍矩阵特征值和特征向量的计算方法，大家知道求一个矩阵的特征值的问题实质上是求一个多项式的根的问题，而数学上已经证明 5 阶以上的多项式的根一般不能用有限次运算求得。因此，矩阵特征值的计算方法本质上都是迭代，而对于大型的稀疏矩阵就需要用幂法求解最简单。但是由于收敛速度非常的慢所以我们需要用加速的方法加快收敛速度而本次论文也是针对加速问题来通过对几种方法的研究讨论。并且通过算法的实现来说明那种加速算法收敛得快，哪个计算量比较节省。其实本文主要讨论的问题是一个应用中常见的一类数值计算问题。2 加速算法的背景2.1 幂法幂法是计算一个矩阵的模最大特征值和对应的特征向量的一种迭代方法。它适用于大型稀疏矩阵。为了说明其基本思想我们先假设是可对角化的，即有如下分解n nACA1AXX其中112,n n ndiagXxxC LL非奇异，再假定12.nL现任取一向量由于的列向量构成的一组基，故可表示为0.nuCXnC0u这里这样，我们有01 122,nnuxxxL.iC0 2111 1 21n kk jj jn k jjj jknjk jj jA uA xxxx 由此可知 0 1 1 1lim.kkkA ux这表明，当而且充分大时，向量这是的一个很好的近似特征向量。然而，实10k01kkkA uuA际计算时，这是行不通的，其原因有二：一是我们事先并不知道的特征值；二是对充分大的A1计算的工作量太大。所以找出一种工作量较小的方法，而幂法求解收敛速度很慢所以我们还kkA 要找出一种加快速度的算法。迭代格式： 1,kkyAu是的模最大分量， ,kk kjjky,kkkyu其中是任意给定的初始向量，通常要求0nuC01.u定理 2.1.1 设有个线性无关的特征向量，主特征值满足则对任n nARn112,nL何非零初始向量按下面构造的向量序列0010 ,vu ,:kkuv则有0010, 1,2,max(),kkkkk k kvuvAu kvvu L(1) 11lim,maxkkxux(2) 1lim.kk (注：此定理证明参阅文献5)例 1 计算矩阵的主特征值。41405130102A 用幂法求解矩阵 A 的计算结果如下表Kmax()kv17.27.37.LL246.256.266.276.由此求得主特征值6.0000003831.2.2 幂法的应用物理，力学和工程技术中的很多问题在数学上都归结为求矩阵特征值问题。例如，振动问题（大型桥梁或建筑物的振动，机械振动，电磁振荡等） ，物理学中某些临界值得确定，这些问题都归结为求矩阵的特征值的数学问题。而幂法求解实际应用矩阵特征值是十分有效方法之一，但是收敛速度太慢，所以在实际应用中它所需要的时间非常的长，而且计算过程中所消耗的时间造成了实际问题的完成进度。因而我们需要通过用加速算法来加快收敛速度，让实际问题提前或者按时完成。为了加快幂法求解矩阵主特征值的收敛速度，让幂法更有效广泛的运用在实际应用生活中，我们现在就来认识几种加速方法，如原点平移法、Rayleigh 商加速、Aitken 加速算法、一种改进的Aitken 加速算法和一种新的改进的 Aitken 加速算法并且对他们进行比较，看哪种加速方法收敛得快，哪种计算量比较节省等。下面我们就来说说这几种加速方法。3 常见的几种加速算法3.1 原点平移法定理 3.1.1 设,个互不相同的特征值满足并且模最大特征值是n nACp12,nL1半单的（即的几何重数等于它的代数重数） 。如果初始向量在的特征子空间上的投影不为零，10u1则定理（2.1.2）产生的向量序列收敛到的一个特征向量，而且由定理（2.1.2）产生的 ku11x数值序列收敛到。 k1（注：此定理证明参阅1）由定理（2.1.1）可知幂法的收敛速度主要取决于的大小。在定理（2.1.1）的条件下，21 这个数总是小于 1 的，它越小收敛也就越快。当它接近于 1 时收敛是很慢的。所以为了加快幂法的收敛速度，通常用位移的方法，即应用幂法于上。如果适当选取可使之模最大特AIuAI征值与其他特征值之模的距离更大，就可起到加速的目的。首先我们引进矩阵BApI其中为选择参数。设的特征值为则的相应特征值为pA12,n LB而且，的特征向量相同。12,npppLAB如果需要计算的主特征值，就要适当选择，使仍然是的主特征值，且使A1p1pB2211.p p 对应用幂法，使得在计算的主特征值的过程中得到加速。这种方法通常称为原点平移BB1p法。对于的特征值的某种分布，它是十分有效的。A对于参数的选择依赖于对矩阵特征值分布的大致了解。通常可以用 Gerschgorin（盖尔）pA圆盘定理得到矩阵的特征值分布情况。A定理 3.1.2（Gerschgorin 圆盘定理）设为阶实矩阵，则An（1）的每一个特征值必定属于下述个闭圆盘（称为 Gerschgorin 圆）An的并集之中；1,1,2,niiiij jj iarainL（2）由矩阵的所有 Gerschgorin 圆组成的连通部分中任取一个，如果由个 Gerschgorin 圆Ak构成，则在这个连通部分中有且仅有的个特征值（Gerschgorin 圆相重时重复计数，Ak特征值相同时也重复计数）。求得矩阵的主特征值后，可得矩阵的主特征值同时得到对应的特征向B11pA11, p量的近似。1x（注：此定理证明参阅文献1）事实上，如果对于矩阵的特征值能够分离得很清楚，就可以利用原点平移法求得矩阵的所有特征值及其相应的特征向量。但需要说明的是，虽然常常能够选择有利的值，使幂法得到加速，但p设计一个自动选择适当参数的过程是困难的。下面考虑的特征值是实数时，怎样选择使幂pAp法计算得到加速。1设的特征值满足A（3.1.1）121,nnL则不管如何，的主特征值为或当我们希望计算及时，首先应选择pBApI1p.np11x使p1,npp且使收敛速度的比值1p211max,minnpp pp显然，当即时为最小，这时收敛速度的比值为211,npp pp 1 2npp221112.2nnnpp pp 当的特征值满足（3.1.1）且能初步估计时，我们就能确定的近似值。A2,n p当希望计算时，应选择n11,2npp使得应用幂法计算得到加速。n3.2 Rayleigh 商加速定义 设为阶实对称矩阵，对于任一非零向量，称Anx , ,Ax xR xx x为对应于向量的 Rayleigh 商。x定理 3.2.1 设为对称矩阵（其特征值次序记为）则n nAR121nnL（1） (对任何非零)； 1, ,nAx x x xnxR（2） 1,0,max;,nx RxAx x x x （3） ,0,min;,nnx RxAx x x x （注：此定理证明参阅文献5）由定理 3.2.1 知，对称矩阵的及可用 Rayleigh 商的极值来表示。下面我们将把A1nRayleigh 商应用到用幂法计算实对称矩阵的主特征值的加速收敛上来。A定理 3.2.2 设为对称矩阵，特征值满足n nAR121,nnL对应的特征向量满足应用幂法（公式（2.1.2） ）计算的主特征值，则规范化向,ijijx xA1量的 Rayleigh 商给出的ku1 22 1 1,.,kkkkkAu u u u（注：此定理证明参阅5）3.3 Aitken 加速算法Aitken 加速算法在诸多方面得到了广泛的应用，尤其是从对幂法加速的诸多方法中突颖而出。但是在实际应用中，由于幂法针对的大多是大型矩阵，而计算速度要求较快，精度要求较高，传统的 Aitken 方法越来越不能满足需要。3.3.1 Aitken 加速算法设序列线性收敛到极限，而且对所有有如果存在实数，且 0nnPp0,n 0.nppA满足：1,A 则定义为1lim.nnnppApp 2 1212nn nn nnnppqpppp的序列收敛到，且比快，而且 0nnqp 0nnplim0.nnnpp pq算法实现:把上述方法应用于幂法迭代格式中的序列即可做到对幂法的加速。具体算 1kkx法伪代码如下:(1) 输入初始向量误差限,最大迭代次数. ,ijAa12,nxx xxLN(2) 置001,0,1.ka(3) 求整数,使r11max,.riri nxxxa (4) 计算置,xyxAya2.rxa(5) 计算.2a0122 1 0aaaaan (6) 若输入停机；否则转（7）.,0,x(7) 若.kN置转（3）；否则停止.10210,1aa aakk 例 2 用此加速算法计算矩阵的主特征值。41405130102A 计算结果如下表Kmax()kv17.812526.36.062546.LL116.126.136.146.156.166.由此 A 的主特征值6.0000000091.3.3.2 改进的 Aitken 加速算法文献6给出一种改进的 Aitken 算法，具体算法伪代码如下:(1) 输入初始向量误差限,最大迭代次数. ,ijAa12,nxx xxLN(2) 置0101,0,1.kaa(3) 求整数,使r 1max,.riri nxxxa (4) 计算置,xyxAya2.rxa(5) 计算2 10 0 210a.23aa aaa(6) 若输入停机；否则转（7）.,0,x(7) 若.kN置转（3）；否则停止.10210,1aa aakk 例 3 用此改进的加速算法计算矩阵的主特征值。41405130102A 计算结果如下表格：Kmax()kv16.24.34.44.LL365.376.385.396.405.415.由此得到结果5.9999999993.经过比较我发现改进的加速算法还不如不改进的。所以我重新阅读了文献6发现证明中有一步假设是不合适的，而定理的结论是依赖于这个假设的，因此，我们认为，这样的算法在实际应用中不具有应用价值。具体的解释说明请关注附录。3.3.3 新的改进的 Aitken