幂法求解矩阵主特征值的加速方法 毕业论文.doc

共20页 河南理工大学数学与信息科学学院本科毕业论文 第20页幂法求解矩阵主特征值的加速方法 摘要：本论文主要研究的是幂法求解矩阵的主特征值和特征向量。物理、力学和工程技术中有许多需要我们求矩阵的按模最大的特征值（及称为主特征值）和特征向量。幂法是计算一个矩阵的模最大特征值和对应的特征向量的一种迭代方法。它最大的优点是方法简单，适合于大型稀疏矩阵的主特征值，但是收敛速度非常慢。所以我们要用加速的方法来加速收敛，加速方法包括原点平移加速、rayleigh商加速和aitken加速算法。关键词：幂法；原点平移加速；rayleigh商加速；aitken加速算法1 引言我们来介绍矩阵特征值和特征向量的计算方法，大家知道求一个矩阵的特征值的问题实质上是求一个多项式的根的问题，而数学上已经证明5阶以上的多项式的根一般不能用有限次运算求得。因此，矩阵特征值的计算方法本质上都是迭代，而对于大型的稀疏矩阵就需要用幂法求解最简单。但是由于收敛速度非常的慢所以我们需要用加速的方法加快收敛速度而本次论文也是针对加速问题来通过对几种方法的研究讨论。并且通过算法的实现来说明那种加速算法收敛得快，哪个计算量比较节省。其实本文主要讨论的问题是一个应用中常见的一类数值计算问题。2 加速算法的背景2.1幂法幂法是计算一个矩阵的模最大特征值和对应的特征向量的一种迭代方法。它适用于大型稀疏矩阵。为了说明其基本思想我们先假设是可对角化的，即有如下分解其中非奇异，再假定现任取一向量由于的列向量构成的一组基，故可表示为这里这样，我们有由此可知 这表明，当而且充分大时，向量这是的一个很好的近似特征向量。然而，实际计算时，这是行不通的，其原因有二：一是我们事先并不知道的特征值；二是对充分大的计算的工作量太大。所以找出一种工作量较小的方法，而幂法求解收敛速度很慢所以我们还要找出一种加快速度的算法。迭代格式： 是的模最大分量，其中是任意给定的初始向量，通常要求定理 2.1.1 设有个线性无关的特征向量，主特征值满足则对任何非零初始向量按下面构造的向量序列则有(1)(2) (注：此定理证明参阅文献5)例1 计算矩阵的主特征值。用幂法求解矩阵a的计算结果如下表17.812600027.623501237.2003614246.0001200680256.0000006893266.0000004232276.0000003831由此求得主特征值2.2幂法的应用物理，力学和工程技术中的很多问题在数学上都归结为求矩阵特征值问题。例如，振动问题（大型桥梁或建筑物的振动，机械振动，电磁振荡等），物理学中某些临界值得确定，这些问题都归结为求矩阵的特征值的数学问题。而幂法求解实际应用矩阵特征值是十分有效方法之一，但是收敛速度太慢，所以在实际应用中它所需要的时间非常的长，而且计算过程中所消耗的时间造成了实际问题的完成进度。因而我们需要通过用加速算法来加快收敛速度，让实际问题提前或者按时完成。为了加快幂法求解矩阵主特征值的收敛速度，让幂法更有效广泛的运用在实际应用生活中，我们现在就来认识几种加速方法，如原点平移法、rayleigh商加速、aitken加速算法、一种改进的aitken加速算法和一种新的改进的aitken加速算法并且对他们进行比较，看哪种加速方法收敛得快，哪种计算量比较节省等。下面我们就来说说这几种加速方法。3 常见的几种加速算法3.1原点平移法定理 3.1.1 设,个互不相同的特征值满足并且模最大特征值是半单的（即的几何重数等于它的代数重数）。如果初始向量在的特征子空间上的投影不为零，则定理（2.1.2）产生的向量序列收敛到的一个特征向量，而且由定理（2.1.2）产生的数值序列收敛到。（注：此定理证明参阅1）由定理（2.1.1）可知幂法的收敛速度主要取决于的大小。在定理（2.1.1）的条件下，这个数总是小于1的，它越小收敛也就越快。当它接近于1时收敛是很慢的。所以为了加快幂法的收敛速度，通常用位移的方法，即应用幂法于上。如果适当选取可使之模最大特征值与其他特征值之模的距离更大，就可起到加速的目的。首先我们引进矩阵其中为选择参数。设的特征值为则的相应特征值为而且，的特征向量相同。如果需要计算的主特征值，就要适当选择，使仍然是的主特征值，且使对应用幂法，使得在计算的主特征值的过程中得到加速。这种方法通常称为原点平移法。对于的特征值的某种分布，它是十分有效的。对于参数的选择依赖于对矩阵特征值分布的大致了解。通常可以用gerschgorin（盖尔）圆盘定理得到矩阵的特征值分布情况。定理3.1.2（gerschgorin圆盘定理）设为阶实矩阵，则（1） 的每一个特征值必定属于下述个闭圆盘（称为gerschgorin圆）的并集之中；（2） 由矩阵的所有gerschgorin圆组成的连通部分中任取一个，如果由个gerschgorin圆构成，则在这个连通部分中有且仅有的个特征值（gerschgorin圆相重时重复计数，特征值相同时也重复计数）。求得矩阵的主特征值后，可得矩阵的主特征值同时得到对应的特征向量的近似。（注：此定理证明参阅文献1）事实上，如果对于矩阵的特征值能够分离得很清楚，就可以利用原点平移法求得矩阵的所有特征值及其相应的特征向量。但需要说明的是，虽然常常能够选择有利的值，使幂法得到加速，但设计一个自动选择适当参数的过程是困难的。下面考虑的特征值是实数时，怎样选择使幂法计算得到加速。设的特征值满足 （3.1.1）则不管如何，的主特征值为或当我们希望计算及时，首先应选择使且使收敛速度的比值显然，当即时为最小，这时收敛速度的比值为当的特征值满足（3.1.1）且能初步估计时，我们就能确定的近似值。当希望计算时，应选择使得应用幂法计算得到加速。3.2 rayleigh商加速定义 设为阶实对称矩阵，对于任一非零向量，称为对应于向量的rayleigh商。定理 3.2.1 设为对称矩阵（其特征值次序记为）则（1） (对任何非零)；（2）（3）（注：此定理证明参阅文献5）由定理3.2.1知，对称矩阵的及可用rayleigh商的极值来表示。下面我们将把rayleigh商应用到用幂法计算实对称矩阵的主特征值的加速收敛上来。定理 3.2.2 设为对称矩阵，特征值满足对应的特征向量满足应用幂法（公式（2.1.2）计算的主特征值，则规范化向量的rayleigh商给出的（注：此定理证明参阅5）3.3 aitken加速算法aitken加速算法在诸多方面得到了广泛的应用，尤其是从对幂法加速的诸多方法中突颖而出。但是在实际应用中，由于幂法针对的大多是大型矩阵，而计算速度要求较快，精度要求较高，传统的aitken方法越来越不能满足需要。3.3.1 aitken加速算法设序列线性收敛到极限，而且对所有有如果存在实数，且满足：则定义为的序列收敛到，且比快，而且算法实现:把上述方法应用于幂法迭代格式中的序列即可做到对幂法的加速。具体算法伪代码如下:(1) 输入初始向量误差限,最大迭代次数.(2) 置(3) 求整数,使(4) 计算置(5) 计算(6) 若输入停机；否则转（7）.(7) 若置转（3）；否则停止.例2 用此加速算法计算矩阵的主特征值。计算结果如下表17.812526.266666666636.062546.0153846153116.0000037368126.0000002335136.0000000583146.0000000364156.0000000145166.0000000091由此a的主特征值3.3.2改进的aitken加速算法文献6给出一种改进的aitken 算法，具体算法伪代码如下:(1) 输入初始向量误差限,最大迭代次数.(2) 置(3) 求整数,使(4) 计算置(5) 计算(6) 若输入停机；否则转（7）.(7) 若置转（3）；否则停止.例3 用此改进的加速算法计算矩阵的主特征值。计算结果如下表格：16.410256410224.399999999934.166666666644.0769230769365.9999964976376.0000000700385.9999999562396.0000000730405.9999999945415.9999999993由此得到结果经过比较我发现改进的加速算法还不如不改进的。所以我重新阅读了文献6发现证明中有一步假设是不合适的，而定理的结论是依赖于这个假设的，因此，我们认为，这样的算法在实际应用中不具有应用价值。具体的解释说明请关注附录。3.3.3新的改进的aitken加速算法由于改进的加速没有达到我们想要的结果。所以我又参阅文献7了解到一种解非线性方程的新算法，经过反复与本文对比发现此新算法可以推及应用到求矩阵的主特征值上，具体计算步骤如下：(1) 输入初始向量误差限,最大迭代次数.(2) 置(3) 求整数,使(4) 计算置(5) 置(6) 计算置(7) 计算(8) 若输入停机；否则转（7）.(9) 若置转（3）；否则停止.例4.用此改进的加速算法计算矩阵的主特征值。计算结果如下表格：16.009554140126.000429664036.0000246277106.0000000009116.0000000001125.9999999999136由此得到结果新的加速算法既可以加快收敛速度而且使结果更接近于真实值。4 结论通过用vc+将几种加速算法实现，并作了一个对比得出以下结论：(1)原点平移加速算法是一个矩阵变换方法。这种变换容易计算，又不破坏矩阵的稀疏性而且收敛速度比较快，但是参数的选择依赖于对的特征值分布的大致了解，而且设计一个自动选择适当参数的过程比较困难，所以在参数的选择过程中会遇到很多麻烦和困难造成花费很多时间。(2)rayleigh商加速算法是系统特征值问题中一个非常重要的概念。该方法适用于求解最大特征值和相应的特征向量，一般情况下收敛速度是二阶，矩阵对称正定时到达三阶。rayleigh商加速算法的收敛速度与反幂法相同，而且不用选取特定的迭代初值是求解这类问题更出色的一种方法。(3)aitken加速算法比原点平移加速收敛得更快，而且它的计算量也不是太大。但是有时aitken加速方法可能会失败，如当迭代向量起伏较大的时候初值和真实值有较大距离时aitken加速就可能失败。（4）改进的aitken加速算法由于没有达到加速的目的，所以我们现在还不能用此改进的加速算法。我们对于此改进的加速算法还需完善。（5）新的改进的aitken加速算法既大大减少了计算量又加快了收敛速度和让计算结果更接近于真实值。总之这几种方法都可以达到使幂法加速求解矩阵主特征值的目的，而其中最好的是新的改进的加速算法，并且我认为这种方法用于实际问题和其它算法领域的加速当中，相信会有更好的效果。而本文介绍的改进的aitken加速还需要进一步完善。附录：程序（1）#include#includeint vector_max_component(const int n,const double *v) double max=fabs(v0); for(int i=0,index=0;in;i+)if(maxfabs(vi)index=i;max=fabs(vi); return index;double vector_inner_product(const int n,const double *u,const double *v) double value=0; for(int i=0;in;i+)value+=ui*vi; return value;double vector_infty_norm(const int n,const double *v)double value=fabs(v0);for(int i=1;in;i+)if(valuefabs(vi)value=fabs(vi); return value;double *matrix_multiply_vector(const int& m,const int& n,double *a,double *v) double *vector; vector=new double m; for(int i=0;im;i+)vectori=vector_inner_product(n,ai,v); return vector;double *square_matrix_multiply_vector(const int& n,double *a,double *v) return matrix_multiply_vector(n,n,a,v);double *matrix_principal_eigen(const int n,double *a) double *v=new doublen,*eig=new double*2;/先定义两个向量，分别存储乘幂前后的迭代向量; eig0=new double1;/存放主特征值; eig1=new doublen;/存放主特征向量; eig00=0; for(int i=0;in;i+)eig1i=1+i;/用（1,2,.,n)作为初始试探向量v； double old_lambda; int index;/定义一个指标变量，来提取乘幂后的最大分量指标； int step=0; do step+; old_lambda=eig00; /double norm=vector_1_norm(n,eig1);/规范化，防止特征值大的时候计算机存储有溢出； double norm=vector_infty_norm(n,eig1);/规范化，防止特征值大的时候计算机存储有溢出； for(i=0;in;i+)vi=eig1i/norm;/更新； eig1=square_matrix_multiply_vector(n,a,v);/乘幂； index=vector_max_component(n,eig1);/提取最大分量指标； eig00=eig1index/vindex;/计算这次乘幂后最大分量与乘幂前该分量的比值； /cout第step+步，分量vindex = vindex, vvindex = vvindex,主特征值lambda的近似值为：lambdaendl; /while(fabs(vvindex)1e-15&step1e-8&step10000);/x循环终止条件为二选一； cout经过step步乘幂方法计算,得到主特征值近似值为:eig00endl; return eig;void main()cout.precision(16);int i,n;cout请输入矩阵规模n:n;double *a=new double*n;cout请输入 n x n 的矩阵元素：endl;for(i=0;in;i+)ai=new doublen;for(int j=0;jaij;double *eig=matrix_principal_eigen(n,a);cout矩阵主特征值为：eig00endl;cout矩阵主特征向量为：endl;for(i=0;in;i+)couteig1iendl;程序（2）#include#includeint vector_max_component(const int n,const double *v) double max=fabs(v0); for(int i=0,index=0;in;i+)if(maxfabs(vi)index=i;max=fabs(vi); return index;double vector_inner_product(const int n,const double *u,const double *v) double value=0; for(int i=0;in;i+)value+=ui*vi; return value;double vector_infty_norm(const int n,const double *v)double value=fabs(v0);for(int i=1;in;i+)if(valuefabs(vi)value=fabs(vi); return value;double *matrix_multiply_vector(const int& m,const int& n,double *a,double *v) double *vector; vector=new double m; for(int i=0;im;i+)vectori=vector_inner_product(n,ai,v); return vector;double *square_matrix_multiply_vector(const int& n,double *a,double *v) return matrix_multiply_vector(n,n,a,v);double *matrix_principal_eigen(const int n,double *a) double *v=new doublen,*eig=new double*2;/先定义两个向量，分别存储乘幂前后的迭代向量; eig0=new double1;/存放主特征值; eig1=new doublen;/存放主特征向量; eig00=0; for(int i=0;in;i+)eig1i=1;/用（1,2,.,n)作为初始试探向量v； double a,a0=0,a1=0,a2,lmd,lmd0=1; int index;/定义一个指标变量，来提取乘幂后的最大分量指标； for(i=0;i100;i+)index=vector_max_component(n,eig1);a=fabs(eig1index);for(int j=0;jn;j+)eig1j/=a;eig1=square_matrix_multiply_vector(n,a,eig1);index=vector_max_component(n,eig1);a2=fabs(eig1index);/lmd=a0-pow(a1-a0,2.0)/(2*a2-3*a1+a0);/new-aitken method;lmd=a0-pow(a1-a0,2.0)/(a2-2*a1+a0);/aitken method;eig00=lmd;cout经过i步计算，求得主特征值为：lmdendl;if(fabs(lmd-lmd0)1e-8)break;a0=a1;a1=a2;lmd0=lmd; return eig;void main()cout.precision(16);int i,n;cout请输入矩阵规模n:n;double *a=new double*n;cout请输入 n x n 的矩阵元素：endl;for(i=0;in;i+)ai=new doublen;for(int j=0;jaij;double *eig=matrix_principal_eigen(n,a);cout矩阵主特征值为：eig00endl;cout矩阵主特征向量为：endl;for(i=0;in;i+)couteig1iendl;程序（3）#include#includeint vector_max_component(const int n,const double *v) double max=fabs(v0); for(int i=0,index=0;in;i+)if(maxfabs(vi)index=i;max=fabs(vi); return index;double vector_inner_product(const int n,const double *u,const double *v) double value=0; for(int i=0;in;i+)value+=ui*vi; return value;double vector_infty_norm(const int n,const double *v)double value=fabs(v0);for(int i=1;in;i+)if(valuefabs(vi)value=fabs(vi); return value;double *matrix_multiply_vector(const int& m,const int& n,double *a,double *v) double *vector; vector=new double m; for(int i=0;im;i+)vectori=vector_inner_product(n,ai,v); return vector;double *square_matrix_multiply_vector(const int& n,double *a,double *v) return matrix_multiply_vector(n,n,a,v);double *matrix_principal_eigen(const int n,double *a) double *v=new doublen,*eig=new double*2;/先定义两个向量，分别存储乘幂前后的迭代向量; eig0=new double1;/存放主特征值; eig1=new doublen;/存放主特征向量; eig00=0; for(int i=0;in;i+)eig1i=1;/用（1,2,.,n)作为初始试探向量v； double a0=0,a1=0,a2,a3,lmd,lmd0=1; int index;/定义一个指标变量，来提取乘幂后的最大分量指标；double e=1e-14; for(i=0;i100;i+)index=vector_max_component(n,eig1);a0=eig1index;couta0endl;for(int j=0;jn;j+)eig1j/=a0;eig1=square_matrix_multiply_vector(n,a,eig1);index=vector_max_component(n,eig1);a1=eig1index;couta1endl;a2=(a0+a1)/2.0;couta2endl;for(j=0;jn;j+)eig1j/=a1;eig1=square_matrix_multiply_vector(n,a,eig1);index=vector_max_component(n,eig1);a3=eig1index;for(j=0;jn;j+)eig1j/=a3;couta3endl;/lmd=a0-pow(a1-a0,2.0)/(a3-2*a1+a0);/aitken method;lmd=a0-(a1-a0)*(a2-a0)/(a3+a0-a1-a2);/my method;eig00=lmd;cout经过i步计算，求得主特征值为：lmdendl;if(fabs(lmd-lmd0)e)break;lmd0=lmd; return eig;void main()cout.precision(16);int i,n;cout请输入矩阵规模n:n;double *a=new double*n;cout请输入 n x n 的矩阵元素：endl;for(i=0;in;i+)ai=new doublen;for(int j=0;jaij;double *eig=matrix_principal_eigen(n,a);cout矩阵主特征值为：eig00endl;cout矩阵主特征向量为：endl;for(i=0;in;i+)couteig1iendl;(四)改进的aitken加速算法的证明：我们只需证明以上证明用到了的假设，并且其结论是以此为基础成立的，而这个假设我们认为是不合适的（易知这样的条件甚至会导致原来简单迭代序列的发散），这样情况下，讨论迭代的加速显然是不恰当的。因此，我们为了将此算法改善，结合参考文献7中关于非线性方程迭代方法的加速技巧，给出了3.3.3算法。致谢：本论文是在牛老师的悉心指导下完成的，牛老师对学术的严谨和精益求精的工作作风给我留下了深刻的印象。从选题后的题目分析到开题报告，从写作提纲，再到毕业论文的编写、修改，每一步都有牛老师的细心指导和认真解析，在此我表示衷心的感谢。四年大学生活即将结束，回顾几年的历程，老师们给了我们很多指导和帮助。他们严谨的治学，优良的作风和敬业的态度，为我们树立了为人师表的典范。在此，我对所有数信学院的老师表示感谢，祝您们身体健康，工作顺利！参考文献1徐树方，高立，张平文.数值线性代数m.北京大学出版社.2000.2曹志浩.矩阵特征值问题m.上海科学技术出版社.1980.3王萼方，石生明.高等代数(第三版)m.高等教育出版社.2003.4华东师范大学师范系.数学分析(第三版)m.高等教育出版社.2001.5李庆扬，王能