数值分析矩阵特征值问题计算

数值分析矩阵特征值问题计算第一页，共八十四页，2022年，8月28日8.1引言工程技术中有多种振动问题，如桥梁或建筑物的振动，机械零件、飞机机翼的振动，及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题.下面先复习一些矩阵的特征值和特征向量的基础知识.第二页，共八十四页，2022年，8月28日定义1⑴已知n阶矩阵A=(aij)，则称为A的特征多项式.一般有n个根(实的或复的，复根按重数计算)称为A的特征值.用λ(A)表示A的所有特征值的集合.

A的特征方程第三页，共八十四页，2022年，8月28日⑵设λ为A的特征值，相应的齐次方程组

注：当A为实矩阵时，

(λ)=0为实系数n次代数方程，其复根是共轭成对出现.的非零解x称为矩阵A的对应于λ的特征向量.

例1求A的特征值及特征向量，其中第四页，共八十四页，2022年，8月28日

解矩阵A的特征方程为求得矩阵A的特征值为：对应于各特征值矩阵A的特征向量分别为：第五页，共八十四页，2022年，8月28日

定理1设λ为A∈Rn×n的特征值,且Ax=λx(x0)，则有⑵λ-p为A-pI的特征值，即(A-pI)x=(λ-p)x；⑴cλ为的cA特征值(c≠0为常数)；下面叙述有关特征值的一些结论：⑶λk为Ak的特征值，即Akx=λkx；⑷设A为非奇异矩阵，那么λ≠0,且λ-1为A-1的特征值，即A-1x=λ-1x.第六页，共八十四页，2022年，8月28日

定理2设λi(i=1,2,,n)为n阶矩阵A=(aij)的特征值，则有⑴称为A的迹；⑵定理3设A∈Rn×n，则有

定理4设A为分块上三角矩阵，即其中每个对角块Aii均为方阵，则第七页，共八十四页，2022年，8月28日

定理5设A与B为相似矩阵（即存在非奇异矩阵P使B=P-1AP），则定理5说明，一个矩阵A经过相似变换，其特征值不变.一个亏损矩阵是一个没有足够特征向量的矩阵，亏损矩阵在理论上和计算上都存在困难.⑴A与B有相同的特征值；⑵如果y是B的特征向量，则Py是A的特征向量.

定义2如果实矩阵A有一个重数为k的特征值λ,且对应于λ的A的线性无关的特征向量个数<k，则A称为亏损矩阵.第八页，共八十四页，2022年，8月28日

定理6⑴A∈Rn×n可对角化，即存在非奇异矩阵P使的充分必要条件是A具有n个线性无关的特征向量.⑵如果A∈Rn×n有m个(m≤n)不同的特征值λ1,λ2,,λm，则对应的特征向量x1,x2,,xm线性无关.第九页，共八十四页，2022年，8月28日

定理7(对称矩阵的正交约化)设A∈Rn×n为对称矩阵，则⑶存在一个正交矩阵P使的且λ1,λ2,,λn为A的特征值，而P＝(u1,u2,,un)列向量uj为A的对应于λj

的单位特征向量.⑴A的特征值均为实数；⑵A有n个线性无关的特征向量；第十页，共八十四页，2022年，8月28日

定义3

设n阶矩阵A=(aij)，令下面讨论矩阵特征值界的估计.⑴；⑵集合称为复平面上以aii为圆心，以ri为半径的n阶矩阵A的n个Gerschgorin圆盘.第十一页，共八十四页，2022年，8月28日

定理8(Gerschgorin圆盘定理)特别地，如果A的一个圆盘Di是与其它圆盘分离(即孤立圆盘)，则Di中精确地包含A的一个特征值.⑴设n阶矩阵A＝(aij)，则A的每一个特征值必属于下面某个圆盘之中⑵如果A有m个圆盘组成一个连通的并集S，且S与余下n-m个圆盘是分离的，则S内恰包含A的m个特征值.或者说A的特征值都在n个圆盘的并集中.第十二页，共八十四页，2022年，8月28日证明只就⑴给出证明.设λ为A的特征值，即Ax=λx，其中x=(x1,x2,,xn)T0.或记，考虑Ax=λx的第k个方程，即于是即第十三页，共八十四页，2022年，8月28日这说明，A的每一个特征值必位于A的一个圆盘中，并且相应的特征值λ一定位于第k个圆盘中(其中k是对应特征向量x绝对值最大的分量的下标).利用相似矩阵性质，有时可以获得A的特征值进一步的估计，即适当选取非奇异对角阵并做相似变换.适当选取可使某些圆盘半径及连通性发生变化.第十四页，共八十四页，2022年，8月28日

例2估计矩阵A的特征值范围，其中

解矩阵A的3个圆盘为由定理8，可知A的3个特征值位于3个圆盘的并集中，由于D1是孤立圆盘，所以D1内恰好包含A的一个特征值λ1(为实特征值)，即A的其它两个特征值λ2,λ3包含在D2,D3的并集中.第十五页，共八十四页，2022年，8月28日现在取对角阵做相似变换矩阵A1的3个圆盘为第十六页，共八十四页，2022年，8月28日显然，3个圆盘都是孤立圆盘，所以，每一个圆盘都包含A的一个特征值(为实特征值)，且有估计第十七页，共八十四页，2022年，8月28日当A为实矩阵，如果限制用正交相似变换，由于A有复的特征值，A不能用正交相似变换约化为上三角阵.用正交相似变换能约化到什么程度呢？定理9(Schur定理)设A∈Rn×n，则存在酉矩阵U使其中rii(i=1,2,,n)为A的特征值.下面给出理论上有关通过酉相似变换及正交变换可以约化一般矩阵A到什么程度的问题.第十八页，共八十四页，2022年，8月28日其中Rii(i=1,2,,m)为一阶或二阶方阵，且每个一阶Rii是A的实特征值，每个二阶对角块Rii的两个特征值是A的两个共轭复特征值.定理10(实Schur分解)设A∈Rn×n，则存在正交矩阵Q使第十九页，共八十四页，2022年，8月28日定义4设A∈Rn×n为对称矩阵，对于任一非零向量x，称我们转向实Schur型的实际计算.为对应于向量x的瑞利(Rayleigh)商.定理11设A∈Rn×n为对称矩阵(其特征值次序记为λ1≥λ2≥≥λn)，则1.(对任何非零x∈Rn)；2.；3..第二十页，共八十四页，2022年，8月28日

证明只证1，关于2,3自己作练习.由于A为实对称矩阵，可将λ1,λ2,,λn

对应的特征向量x1,x2,,xn

正交规范化，则有(xi,xj)=δij，设x0为Rn中任一向量，则有于是从而1成立.结论1说明瑞利商必位于λn和λ1之间.第二十一页，共八十四页，2022年，8月28日关于计算矩阵A的特征值问题，当n＝2,3时，我们还可按行列式展开的办法求(λ)=0的根.但当n较大时，如果按展开行列式的办法，首先求出(λ)的系数，再求(λ)的根，工作量就非常大，用这种办法求矩阵的特征值是不切实际的，由此需要研究求A的特征值及特征向量的数值解法.本章将介绍一些计算机上常用的两类方法，一类是幂法及反幂法（迭代法），另一类是正交相似变换的方法（变换法）.第二十二页，共八十四页，2022年，8月28日幂法与反幂法都是求实矩阵的特征值和特征向量的向量迭代法，所不同的是幂法是计算矩阵的主特征值(矩阵按模最大的特征值称为主特征值，其模就是该矩阵的谱半径)和相应特征向量的一种向量迭代法，而反幂法则是计算非奇异(可逆)矩阵按模最小的特征值和相应特征向量的一种向量迭代法.下面分别介绍幂法与反幂法.8.2幂法及反幂法第二十三页，共八十四页，2022年，8月28日现讨论求λ1及x1的方法.设实矩阵A=(aij)有一个完全的特征向量组，即A有n个线性无关的特征向量，设矩阵A的特征值为λ1,λ2,,λn,相应的特征向量为x1,x2,,xn.已知A的主特征值λ1是实根，且满足条件

8.2.1幂法（又称乘幂法）第二十四页，共八十四页，2022年，8月28日幂法的基本思想是:任取非零的初始向量v0

,由矩阵A构造一向量序列{vk}称为迭代向量，由假设，v0可唯一表示为第二十五页，共八十四页，2022年，8月28日于是其中由假设故从而为λ1的特征向量.第二十六页，共八十四页，2022年，8月28日所以当k充分大时，有即为矩阵A的对应特征值1的一个近似特征向量.用(vk)i

表示vk的第i个分量，则当k充分大时，有即为A的主特征值1的近似值.由于这种由已知非零向量v0及矩阵A的乘幂Ak构造向量序列{vk}以计算A的主特征值1(2.7)及相应特征向量(2.5)的方法就称为幂法.第二十七页，共八十四页，2022年，8月28日迭代公式实质上是由矩阵A的乘幂Ak与非零向量v0相乘来构造向量序列{vk}={Akv0}，从而计算主特征值λ1及其对应的特征向量，这就是幂法的思想.的收敛速度由比值来确定，r越小收敛越快，但当r≈1时收敛可能很慢.第二十八页，共八十四页，2022年，8月28日定理12设A∈Rn×n有n个线性无关的特征向量，主特征值λ1满足条件|λ1|>|λ2|≥≥|λn|，则对任何非零向量v0(a10)，幂法的算式成立.又设A有n个线性无关的特征向量，λ1对应的r个线性无关的特征向量为x1,x2,,xr，则由(2.2)式有如果A的主特征值为实的重根,即λ1=λ2==λr,且|λr|>|λr+1|≥≥|λn|，第二十九页，共八十四页，2022年，8月28日为A的特征向量，这说明当A的主特征值是实的重根时，定理5的结论还是正确的.应用幂法计算A的主特征值λ1及其对应的特征向量时，如果|λ1|>1(或|λ1|<1)，迭代向量

vk的各个不等于零的分量将随k→∞

而趋向于无穷(或趋向于零)，这样在计算机实现时就可能“溢出”.为克服这个缺点，就需要将迭代向量加以规范化.第三十页，共八十四页，2022年，8月28日设有一向量v0，将其规范化得向量为其中max(v)表示v的绝对值最大的分量.即如果有则max(v)=vq，且q为所有绝对值最大的分量中的最小下标.在定理12的条件下幂法可这样进行：任取一初始向量v00(a10)，构造规范化向量序列为第三十一页，共八十四页，2022年，8月28日由(2.3)式第三十二页，共八十四页，2022年，8月28日第三十三页，共八十四页，2022年，8月28日收敛速度由比值r=|λ2/λ1|确定.总结上述结论，有同理，可得到第三十四页，共八十四页，2022年，8月28日定理13设A∈Rn×n有n个线性无关的特征向量，主特征值λ1满足|λ1|>|λ2|≥≥|λn|，则对任意非零初始向量v0=u0(a10)，有幂法计算公式为则有⑴⑵第三十五页，共八十四页，2022年，8月28日

例1用幂法计算矩阵的主特征值与其对应的特征向量.

解取v0=u0=(0,0,1)T

,则第三十六页，共八十四页，2022年，8月28日直到k=8时的计算结果见下表k12,4,1,40.5,1,0.2524.5,9,7.7590.5,1,0.861135.7222,11.4444,8.36111.44440.5,1,0.736045.4621,10.9223,8.230610.92230.5,1,0.753655.5075,11.0142,8.257611.01420.5,1,0.749465.4987,10.9974,8.249410.99740.5,1,0.750175.5002,11.0005,8.250111.00050.5,1,0.750085.5000,11.0000,8.250011.00000.5,1,0.7500从而

见书p303-例3.第三十七页，共八十四页，2022年，8月28日8.2.2幂法的加速方法1、原点平移法

由前面讨论知道，应用幂法计算A的主特征值的收敛速度主要由比值r=|λ2/λ1|来决定，但当r接近于1时，收敛可能很慢.这时，一个补救办法是采用加速收敛的方法.其中p为参数，设A的特征值为i，则对矩阵B的特征值为i-p，而且A,B的特征向量相同.引进矩阵B=A-pI.第三十八页，共八十四页，2022年，8月28日如果要计算A的主特征值1,只要选择合适的数p，使1-p为矩阵B=A-pI

的主特征值，且那么，对矩阵B=A-pI应用幂法求其主特征值1-p,收敛速度将会加快.这种通过求B=A-pI的主特征值和特征向量，而得到A的主特征值和特征向量的方法叫原点平移法.对于A的特征值的某种分布，它是十分有效的.第三十九页，共八十四页，2022年，8月28日

例4设A∈R4×4有特征值比值r=|λ2/λ1|≈0.9.做变换B=A-12I(p=12),则B的特征值为应用幂法计算B的主特征值μ1的收敛速度的比值为虽然常常能够选择有利的p值,使幂法得到加速,但设计一个自动选择适当参数p的过程是困难的.第四十页，共八十四页，2022年，8月28日下面考虑当A的特征值是实数时，怎样选择p使采用幂法计算λ1得到加速.且使收敛速度的比值设A的特征值都是实数，且满足则对实数p，使矩阵A-pI的主特征值为1-p或n-p时，当我们计算1及x1时，首先应选取p使第四十一页，共八十四页，2022年，8月28日显然，当2-p=-(n-p)时，即P=(2+n)/2＝P*

时ω为最小值，这时收敛速度的比值为当希望计算n时，应选取

p=(1+n-1)/2＝P*

使得应用幂法计算n得到加速.当A的特征值都是实数，满足且2,n能初步估计出来，我们就能确定P*的近似值.第四十二页，共八十四页，2022年，8月28日

例2用原点平移加速法求例1中矩阵A的主特征值与其对应的特征向量.对B应用幂法，仍取v0=(0,0,1)T

,则

解取p=-2.5,做平移变换B=A-pI，则第四十三页，共八十四页，2022年，8月28日迭代5步的计算结果见下表k12,4,3.540.5,1,0.87527,14,10.5625140.5,1,0.754536.76,13.5179,10.140613.51790.5,1,0.750746.7503,13.5007,10.125613.50070.5,1,0.750056.7500,13.5000,10.125013.50000.5,1,0.7500可得到B的主特征值为113.5000,主特征向量为v1

(0.5,1.0,0.7500)T

,因此，A的主特征值为1=1+p11.0000,主特征向量仍为x1=(0.5,1,0.7500)T.第四十四页，共八十四页，2022年，8月28日原点位移的加速方法，是一个矩阵变换方法.这种变换容易计算，又不破坏矩阵A的稀疏性，但p的选择依赖对A的特征值分布的大致了解.

见书p306-例5.第四十五页，共八十四页，2022年，8月28日设A∈Rn×n为对称矩阵，称为向量x的瑞利商，其中(x,x)=xTx为内积.由定理11知道，实对称矩阵A的特征值1及n可用瑞利商的极限值表示.下面我们将瑞利商应用到用幂法计算实对称矩阵A的主特征值的加速上来.2、瑞利商(Rayleigh)加速第四十六页，共八十四页，2022年，8月28日定理14设A∈Rn×n为对称矩阵，特征值满足对应的特征向量xi满足(xi,xj)=δij

(单位正交向量)

，应用幂法公式(2.9)计算A的主特征值1，则规范化向量uk的瑞利商给出1的较好的近似值为由此可见，R(uk)比μk更快的收敛于1.第四十七页，共八十四页，2022年，8月28日

证明由(2.8)式及得第四十八页，共八十四页，2022年，8月28日幂法的瑞利商加速迭代公式可以写为其中A为n阶实对称矩阵.对给定的误差限，当|μ

k–μk-1|<时，取近似值第四十九页，共八十四页，2022年，8月28日8.2.3反幂法

反幂法是用于求非奇异矩阵A的按模最小的特征值和对应特征向量的方法.而结合原点平移法的反幂法则可以求矩阵A的任何一个具有先了解的特征值和对应的特征向量。设矩阵A非奇异,其特征值i(i=1,2,,n),满足其相应的特征向量x1,x2,,xn线性无关，则A-1的特征值为1/i

，对应的特征向量仍为xi

(i=1,2,,n).第五十页，共八十四页，2022年，8月28日此时,A-1的特征值满足因此,对A-1应用幂法,可求出其主特征值1/n

μ

k

和特征向量

xn

uk.从而求得A的按模最小特征值

n

1/μk

和对应的特征向量

xn

uk，这种求A-1的方法就称为反幂法.第五十一页，共八十四页，2022年，8月28日为了避免求A-1,可通过解线性方程组Avk=uk-1得到vk，采用LU分解法，即先对A进行LU分解A=LU,此时反幂法的迭代公式为

反幂法的迭代公式为第五十二页，共八十四页，2022年，8月28日对给定的误差，当|μk–μk-1|<

时，得显然，反幂法的收敛速度取决于比值，比值越小，收敛越快.第五十三页，共八十四页，2022年，8月28日定理15设A∈Rn×n为非奇异矩阵，且有n个线性无关的特征向量，其对应的特征值满足

|1|≥|2|≥≥|n-2|>|n|>0，则对任意非零初始向量u0(an0)，由反幂法计算公式构造的向量序列{vk}，{uk}满足⑴⑵第五十四页，共八十四页，2022年，8月28日

在反幂法中也可以用原点平移法加速迭代过程，或求其它特征值与其对应的特征向量.

如果矩阵(A-pI)-1存在，显然其特征值为对应的特征向量仍然是x1,x2,,xn，现对矩阵(A-pI)-1应用幂法，得到反幂法的迭代公式第五十五页，共八十四页，2022年，8月28日

如果p是A的特征值j的一个近似值，且设j与其它特征值是分离的，即就是说1/(j-p)是矩阵(A-pI)-1的主特征值，可用反幂法(2.12)计算特征值及特征向量.

设A∈Rn×n有n个线性无关的特征向量x1,x2,,xn，则第五十六页，共八十四页，2022年，8月28日其中同理可得：第五十七页，共八十四页，2022年，8月28日定理16设A∈Rn×n有n个线性无关的特征向量，矩阵A的特征值及对应的特征向量分别记为i

及xi(i=1,2,,n)，而p为j的近似值，(A-pI)-1存在，且⑴⑵则对任意非零初始向量u0(aj0)，由反幂法计算公式(2.12)构造的向量序列{vk}，{uk}满足且收敛速度为第五十八页，共八十四页，2022年，8月28日由该定理知，对A-pI(其中p≈j)应用反幂法，可用来计算特征向量xj，只要选择p是j的一个较好的近似且特征值分离情况较好，一般r很小，常常只要迭代一二次就可完成特征向量的计算.反幂法迭代公式中的vk是通过解方程组求得的,为了节省工作量,可以先将A-pI进行三角分解于是求vk相对于解两个三角形方程组第五十九页，共八十四页，2022年，8月28日实验表明,按下述方法选择u0是较好的:选u0使用回代求解(2.13)即得v1,然后再按公式(2.12)进行迭代.

反幂法计算公式见书p311.前面已提到.

见书p311-例6.第六十页，共八十四页，2022年，8月28日8.3豪斯霍尔德方法(1)用初等反射矩阵作正交相似变换约化一般实矩阵A为上海森伯格矩阵.8.3.1引言

本节讨论两个问题(2)用初等反射矩阵作正交相似变换约化对称矩阵A为对称三对角矩阵.于是，求原矩阵特征值问题，就转化为求上海森伯格矩阵或对称三对角矩阵的特征值问题.第六十一页，共八十四页，2022年，8月28日8.3.2用正交相似变换

约化一般实矩阵为上海森伯格矩阵

设A∈Rn×n，下面来说明，可选择初等反射矩阵U1,U2,,Un-2使A经正交相似变换约化为一个上海森伯格阵.(1)设第六十二页，共八十四页，2022年，8月28日其中c1=(a21,,an1)T∈Rn-1

,不妨设c1≠0，否则这一步不需要约化.于是,可选择初等反射阵使，其中令第六十三页，共八十四页，2022年，8月28日则其中第六十四页，共八十四页，2022年，8月28日(2)第k步约化：重复上述过程，设对A已完成第1步,,第k-1步正交相似变换，即有或且第六十五页，共八十四页，2022年，8月28日其中为k阶上海森伯格阵，设ck≠0,于是可选择初等反射阵Rk使其中，Rk计算公式为第六十六页，共八十四页，2022年，8月28日令则第六十七页，共八十四页，2022年，8月28日其中为k+1阶上海森伯格阵，第k步约化只需计算及(当A为对称矩阵时，只需要计算).(3)重复上述过程，则有第六十八页，共八十四页，2022年，8月28日定理17(豪斯霍尔德约化矩阵为上海森伯格阵)设A∈Rn×n则存在初等反射矩阵U1,U2,,Un-2

使为上海森伯格矩阵.总结上述结论，有算法1(豪斯霍尔德约化矩阵为上海森伯格阵)设A∈Rn×n，本算法计算U0TAU0=H(上海森伯格型)，其中U0=U1U2Un-2为初等反射阵的乘积.1.U0←I第六十九页，共八十四页，2022年，8月28日2.对于k=1,2,,n-2(1)计算初等反射阵Rk使本算法约需要5n3/3次乘法运算，要明显形成U0还需要附加2n3/3次乘法.(2)约化计算(3)U0←U0Uk第七十页，共八十四页，2022年，8月28日

例7用豪斯霍尔德方法将矩阵约化为上海森伯格阵.

解选取初等反射阵R1使,其中c1=(2,4)T.

(1)计算R1:第七十一页，共八十四页，2022年，8月28日则有

(2)约化计算:则得到上海森伯格阵为第七十二页，共八十四页，2022年，8月28日8.3.3用正交相似变换

约化对称矩阵为对称三对角矩阵定理18(豪斯霍尔德约化对称矩阵为对称三对角矩阵)设A∈Rn×n为对称矩阵，则存在初等反射矩阵U1,U2,,Un-2使为对称三对角矩阵.第七十三页，共八十四页，2022年，8月28日证明由定理17,存在初等反射矩阵U1,U2,,Un-2

使为上海森伯格阵，且An-1亦是对称阵，因此，An-1为对称三对角阵.

由上面讨论可知，当A为对称阵时，由Ak←Ak+1=AkUkAk一步约化计算中只需计算Rk及RkA22(k)Rk

.又由于A的对称性，故只需计算RkA22(k)Rk的对角线以下元素.注意到第七十四页，共八十四页，2022年，8月28日引进记号则有对对称阵A用初等反射阵正交相似约化为对角三对角阵大约需要2n3/3次乘法.用正交矩阵进行相似约化有一些特点，如构造的，Uk容易求逆，且Uk的元素数量级不大，这个算法是十分稳定的.

算法2见书p318.第七十五页，共八十四页，2022年，8月28日8.4QR方法8.4.1QR算法

Rutishauser(1958)利用矩阵的三角分解提出了计算矩阵特征值的LR算法，Francis(1961,1962)利用矩阵的QR分解建立了计算矩阵特征值的QR算法.

QR方法是一种变换方法，是计算一般矩阵(中小型矩阵)全部特征值问题的最有效方法之一.第七十六页，共八十四页，2022年，8月28日(1)上海森伯格矩阵的全部特征值问题；(2)计算对称三对角矩阵的全部特征值问题，

目前QR方法主要用来计算：且QR方法具有收敛快，算法稳定等特点.

对于一般矩阵A∈Rn×n(或对称矩阵)，则首先用豪斯霍尔德方法将A化为上海森伯格阵B(或对称三对角阵)，然后再用QR方法计算B的全部特征值.

设A∈Rn×n，且A对进行QR分解，即A=QR，第七十七页，共八十四页，2022年，8月28日其中R为上三角阵,Q为正交阵,于是可得到一新矩阵B=RQ=QTAQ.显然，