第8章 矩阵特征值问题的数值解法

第8章矩阵特征值问题的数值解法8.1引言8.2幂法及反幂法8.3雅可比方法8.4QR算法8.1引言

工程技术中有多种振动问题，如桥梁或建筑物的振动，机械零件、飞机机翼的振动，及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题.

下面先复习一些矩阵的特征值和特征向量的基础知识.

定义1⑴已知n阶矩阵A=(aij)，则称为A的特征多项式.一般有n个根(实的或复的，复根按重数计算)称为A的特征值.用λ(A)表示A的所有特征值的集合.

A的特征方程⑵设λ为A的特征值，相应的齐次方程组

注：当A为实矩阵时，

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

定理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.

定理2

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

定理3设A∈Rn×n，则有

定义2

设n阶矩阵A=(aij)，令

下面讨论矩阵特征值界的估计.⑴；⑵集合称为复平面上以aii为圆心，以ri为半径的n阶矩阵A的n个Gerschgorin圆盘.

定理4(Gerschgorin圆盘定理)

特别地，如果A的一个圆盘Di是与其它圆盘分离(即孤立圆盘)，则Di中精确地包含A的一个特征值.⑴设n阶矩阵A＝(aij)，则A的每一个特征值必属于下面某个圆盘之中⑵如果A有m个圆盘组成一个连通的并集S，且S与余下n-m个圆盘是分离的，则S内恰包含A的m个特征值.或者说A的特征值都在n个圆盘的并集中.

证明只就⑴给出证明.设λ为A的特征值，即Ax=λx，其中x=(x1,x2,,xn)T0.或

记，考虑Ax=λx的第k个方程，即于是即这说明，A的每一个特征值必位于A的一个圆盘中，并且相应的特征值λ一定位于第k个圆盘中(其中k是对应特征向量x绝对值最大的分量的下标).

例2

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

解矩阵A的3个圆盘为

由定理8，可知A的3个特征值位于3个圆盘的并集中，由于D1是孤立圆盘，所以D1内恰好包含A的一个特征值λ1(为实特征值)，即A的其它两个特征值λ2,λ3包含在D2,D3的并集中.

定义3设A∈Rn×n为对称矩阵，对于任一非零向量x，称为对应于向量x的瑞利(Rayleigh)商.

定理5设A∈Rn×n为对称矩阵(其特征值次序记为λ1≥λ2≥≥λn)，则1.(对任何非零x∈Rn)；2.；3..

证明只证1，关于2,3自己作练习.

由于A为实对称矩阵，可将λ1,λ2,,λn

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

正交规范化，则有(xi,xj)=δij，设x0为Rn中任一向量，则有于是从而1成立.结论1说明瑞利商必位于λn和λ1之间.

关于计算矩阵A的特征值问题，当n＝2,3时，我们还可按行列式展开的办法求(λ)=0的根.但当n较大时，如果按展开行列式的办法，首先求出(λ)的系数，再求(λ)的根，工作量就非常大，用这种办法求矩阵的特征值是不切实际的，由此需要研究求A的特征值及特征向量的数值解法.

本章将介绍一些计算机上常用的两类方法，一类是幂法及反幂法（迭代法），另一类是正交相似变换的方法（变换法）.

幂法与反幂法都是求实矩阵的特征值和特征向量的向量迭代法，所不同的是幂法是计算矩阵的主特征值(矩阵按模最大的特征值称为主特征值，其模就是该矩阵的谱半径)和相应特征向量的一种向量迭代法，而反幂法则是计算非奇异(可逆)矩阵按模最小的特征值和相应特征向量的一种向量迭代法.下面分别介绍幂法与反幂法.8.2幂法及反幂法现讨论求λ1及x1的方法.

设实矩阵A=(aij)有一个完全的特征向量组，即A有n个线性无关的特征向量，设矩阵A的特征值为λ1,λ2,,λn,相应的特征向量为x1,x2,,xn.已知A的主特征值λ1是实根，且满足条件

8.2.1幂法（又称乘幂法）

幂法的基本思想是:任取非零的初始向量v0

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

表示vk的第i个分量，则当k充分大时，有即为A的主特征值1的近似值.由于

迭代公式实质上是由矩阵A的乘幂Ak与非零向量v0相乘来构造向量序列{vk}={Akv0}，从而计算主特征值λ1及其对应的特征向量，这就是幂法的思想.的收敛速度由比值来确定，r越小收敛越快，但当r≈1时收敛可能很慢.

定理6设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|，为A的特征向量，这说明当A的主特征值是实的重根时，定理6的结论还是正确的.

应用幂法计算A的主特征值λ1及其对应的特征向量时，如果|λ1|>1(或|λ1|<1)，迭代向量

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

而趋向于无穷(或趋向于零)，这样在计算机实现时就可能“溢出”.为克服这个缺点，就需要将迭代向量加以规范化.

设有一向量v0，将其规范化得向量为其中max(v)表示v的绝对值最大的分量.即如果有则max(v)=vq，且q为所有绝对值最大的分量中的最小下标.在定理6的条件下幂法可这样进行：任取一初始向量v00(a10)，构造规范化向量序列为由(2.3)式收敛速度由比值r=|λ2/λ1|确定.总结上述结论，有同理，可得到

定理7设A∈Rn×n有n个线性无关的特征向量，主特征值λ1满足|λ1|>|λ2|≥≥|λn|，则对任意非零初始向量v0=u0(a10)，有幂法计算公式为则有⑴⑵

例1

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

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

,则直到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从而8.2.2幂法的加速方法1、原点平移法

由前面讨论知道，应用幂法计算A的主特征值的收敛速度主要由比值r=|λ2/λ1|来决定，但当r接近于1时，收敛可能很慢.这时，一个补救办法是采用加速收敛的方法.其中p为参数，设A的特征值为i，则对矩阵B的特征值为i-p

，而且A,B的特征向量相同.

引进矩阵B=A-pI.

如果要计算A的主特征值1,只要选择合适的数p，使1-p为矩阵B=A-pI

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

例4

设A∈R4×4有特征值比值r=|λ2/λ1|≈0.9.做变换B=A-12I(p=12),则B的特征值为应用幂法计算B的主特征值μ1的收敛速度的比值为

虽然常常能够选择有利的p值,使幂法得到加速,但设计一个自动选择适当参数p的过程是困难的.

定理

设A∈Rn×n为对称矩阵，特征值满足对应的特征向量xi满足(xi,xj)=δij

(单位正交向量)

，应用幂法公式(2.9)计算A的主特征值1，则规范化向量uk的瑞利商给出1的较好的近似值为

由此可见，R(uk)比μk更快的收敛于1.2、瑞利商加速

证明由(2.8)式及得

幂法的瑞利商加速迭代公式可以写为其中A为n阶实对称矩阵.

对给定的误差限，当|μ

k–μk-1|<时，取近似值反幂法

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

的特征值为1/i

，对应的特征向量仍为xi

(i=1,2,,n).此时,A-1的特征值满足因此,对A-1应用幂法,可求出其主特征值1/n

μ

k

和特征向量

xn

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

n

1/μk

和对应的特征向量

xn

uk

，这种求A-1的方法就称为反幂法.为了避免求A-1,可通过解线性方程组Avk=uk-1得到vk，采用LU分解法，即先对A进行LU分解A=LU,此时反幂法的迭代公式为

反幂法的迭代公式为对给定的误差，当|μk–μk-1|<

时，得显然，反幂法的收敛速度取决于比值，比值越小，收敛越快.

定理8设A∈Rn×n为非奇异矩阵，且有n个线性无关的特征向量，其对应的特征值满足

|1|≥|2|≥≥|n-2|>|n|>0，则对任意非零初始向量u0(an0)

，由反幂法计算公式构造的向量序列{vk}，{uk}满足⑴⑵

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

如果矩阵(A-pI)-1存在，显然其特征值为对应的特征向量仍然是x1,x2,,xn.

如果p是A的特征值j的一个近似值，且设j与其它特征值是分离的，即就是说1/(j-p)是矩阵(A-pI)-1的主特征值，可用反幂法应用于矩阵(A-pI)-1计算j及对应的特征向量.现对矩阵(A-pI)-1应用反幂法，得到如下迭代公式

定理9设A∈Rn×n有n个线性无关的特征向量，矩阵A的特征值及对应的特征向量分别记为i

及xi(i=1,2,,n)，而p为j的近似值，(A-pI)-1存在，且⑴⑵则对任意非零初始向量u0(aj0)

，由反幂法计算公式(2.12)构造的向量序列{vk}，{uk}满足且收敛速度为

由该定理知，对A-pI(其中p≈j)应用反幂法，可用来计算特征向量xj，只要选择p是j的一个较好的近似且特征值分离情况较好，一般r很小，常常只要迭代一二次就可完成特征向量的计算.反幂法迭代公式中的vk是通过解方程组求得的,为了节省工作量,可以先将A-pI进行三角分解于是求vk相对于解两个三角形方程组实验表明,按下述方法选择u0是较好的:选u0使用回代求解(2.13)即得v1,然后再按公式(2.12)进行迭代.例5

求矩阵A最接近于1.9的特征值和相应的特征向量取作迭代，结果如表：8.3Jacobi方法Jacobi

方法的基本思路矩阵的旋转变换Jacobi

迭代法收敛定理

说明解：

例6：

8.4QR算法8.4.1.

化矩阵为Hessenberg形8.4.2QR算法及其收敛性（1）镜面反射矩阵（2）化一般实矩阵为上Hessenberg矩阵（3）化对称矩阵为三对角矩阵

当A为实矩阵，如果限制用正交相似变换，由于A有复的特征值，A不能用正交相似变换约化为上三角阵.用正交相似变换能约化到什么程度呢？(Schur定理)

设A∈Rn×n，则存在酉矩阵U使其中rii(i=1,2,,n)为A的特征值.

下面给出理论上有关通过酉相似变换及正交变换可以约化一般矩阵A到什么程度的问题.其中Rii(i=1,2,,m)为一阶或二阶方阵，且每个一阶Rii是A的实特征值，每个二阶对角块Rii的两个特征值是A的两个共轭复特征值.

(实Schur分解)

设A∈Rn×n，则存在正交矩阵Q使63（1）镜面反射矩阵8.4.1.

化矩阵为Hessenberg形镜面反射矩阵的意义是“成批”消去向量的非零元素。例：解:#(1)用初等反射矩阵作正交相似变换约化一般实矩阵A为上Hessenberg矩阵.化矩阵为上Hessenberg矩阵讨论两个问题(2)用初等反射矩阵作正交相似变换约化对称矩阵A为对称三对角矩阵.

于是，求原矩阵特征值问题，就转化为求上Hessenberg矩阵或对称三对角矩阵的特征值问题.（2）化一般实矩阵为上Hessenberg矩阵

设A∈Rn×n，下面来说明，可选择初等反射矩阵U1,U2,,Un-2使A经正交相似变换约化为一个上Hessenberg阵.

(1)设其中c1=(a21,,an1)T∈Rn-1

,不妨设c1≠0，否则这一步不需要约化.于是,可选择初等反射阵令使其中则其中

(2)第k步约化：重复上述过程，设对A已完成第1步,,第k-1步正交相似变换，即有或且其中为k阶上Hessenberg阵，设ck≠0,于是可选择初等反射阵Rk使其中，Rk计算公式为令则其中为k+1阶上Hessenberg阵，第k步约化只需计算及(当A为对称矩阵时，只需要计算).

(3)重复上述过程，则有

定理11(Household约化矩阵为上Hessenberg阵)设A∈Rn×n则存在初等反射矩阵U1,U2,,Un-2

使为上Hessenberg矩阵.总结上述结论，有

例7用Household方法将矩阵约化为上Hessenberg阵.

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

(1)计算R1:则有

(2)约化计算:则得到上Hessenberg阵为（3）化对称矩阵为对称三对角矩阵推论(Household约化对称矩阵为对称三对角矩阵)设A∈Rn×n为对称矩阵，则存在初等反射矩阵U1,U2,,Un-2使为对称三对角矩阵.8.4.2QR算法及其收敛性

QR算法是一种变换方法，是计算一般矩阵(中小型矩阵)全部特征值问题的最有效方法之一.(1)上Hessenberg矩阵的全部特征值问题；(2)计算对称三对角矩阵的全部特征值问题，

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

对于一般矩阵A∈Rn×n(或对称矩阵)，则首先用Household方法将A化为上Hessenberg阵B(或对称三对角阵)，然后再用QR算法计算B的全部特征值.

设A∈Rn×n，且A对进行QR分解，即A=QR，其中R为上三角阵,Q为正交阵,于是可得到一新矩阵B=RQ=QTAQ.显然，B是由A经过正交相似变换得到，因此B与A的特征值相同.再对B进行QR分解，又可得一新矩阵,重复这一过程可得到矩阵序列：

设A=A1

将A1进行QR分解A1=Q1R1

作矩阵A2=R1Q1=Q1TR1Q1

QR算法，就是利用矩阵的QR分解，按上述递推法则构造矩阵序列{Ak}的过程.只要A为非奇异矩阵，则由QR算法就完全确定{Ak}.

定理13(基本QR算法)设A=A1∈Rn×n,构造QR算法:于是

证明(1),(2)显然，现证(3).用归纳法，显然，当k=1时有,

设Ak-1有分解式

定理14(QR算法的收敛性)设A=(aij)∈Rn×n,(1)如果A的特征值满足:(2)A有标准型A=XDX-1,其中D=diag(1,2,n),且设X-1有三角分解X-1=LU(L为单位下三角阵，U为上三角阵)，则由QR算法产生的{Ak}本质上收敛于上三角矩阵，即若记Ak=(aij(k))，则(1)(2)当i>j时，当i<j时aij(k)的极限不一定存在.

推论如果对称矩阵A满足定理14的条件,则由QR算法产生的{Ak}收敛于对角阵D=diag(1,2,n).

设A∈Rn×n，且A有完备的特征向量集合，如果A的等模特征值中只有实重特征值或多重共轭复特征值，则由QR算法产生的{Ak}本质收敛于分块上三角矩阵(对角块为一阶和二阶子块)，且对角块中每一个2×2子块给出A的一对共轭复特征值，每一个一阶对角子块给出A的实特征值，即