1、Ch8 矩阵特征值问题计算矩阵特征值问题计算引言引言1110102()()31140.定理设 为的特征值且，其中，则定理设 为的特征值且，其中，则（）为的特征值（为常数）；（）为的特征值（为常数）；（）为的特征值，即；（）为的特征值，即；（）为的特征值；（）为的特征值；（）设 为非奇异矩阵，那么且为的特征值，即（）设 为非奇异矩阵，那么且为的特征值，即n nkkARAxxxccAccpApIApI xp xAAAA xx 111(1,2,)()(1)();(2)det().定理设为 阶矩阵的特征值，则定理设为 阶矩阵的特征值，则iijn nnniiiiiniiinnAaatr AA 定义设，如
2、果 有一个重数为 的特征值 且对应定义设，如果 有一个重数为 的特征值 且对应矩阵 的线性无关的特征向量个数少于，则称 为亏损矩阵.矩阵 的线性无关的特征向量个数少于，则称 为亏损矩阵.n nARAkAkA 1()|(1,2,)2|,.定义设，（1）令；定义设，（1）令；（）称复平面上以为圆心，（）称复平面上以为圆心，以 为半径的圆为 的圆盘以 为半径的圆为 的圆盘nijn niijjj iiiiiiiiAarainDzzar zCarAGerschgorin 1(),|1,2,).(2).定理（G erschgori n定理）定理（G erschgori n定理）（）设则 的每个特征值必属于
3、下述某个圆盘中（）设则 的每个特征值必属于下述某个圆盘中（如果 有个圆盘组成一个连通的并集，且 与余下个如果 有个圆盘组成一个连通的并集，且 与余下个圆盘是分离的，则 内恰包含 的个特征值圆盘是分离的，则 内恰包含 的个特征值ijn niiiAaAarinAmSSnmSAm 11121222 (1,2,).定理（Schur定理）设，则存在酉矩阵使定理（Schur定理）设，则存在酉矩阵使其中为 的特征值其中为 的特征值n nnnHnniiARUrrrrrUAURrr inA 11121222,.定理（是Schur定理）设，则存在正交矩阵使得定理（是Schur定理）设，则存在正交矩阵使得其中对角块
4、为一阶或二阶方阵，每个一阶是 的是特征值，其中对角块为一阶或二阶方阵，每个一阶是 的是特征值，每个二阶对角块的两个特征值是 的两个共轭复特征值每个二阶对角块的两个特征值是 的两个共轭复特征值n nmmTmmiiiijjARQRRRRRQ AQRRRARA (,)()(,)定义设 为 阶实对称矩阵，对于任一非零向量，称定义设 为 阶实对称矩阵，对于任一非零向量，称为对应向量 的瑞利商（R ayl ei gh）.为对应向量 的瑞利商（R ayl ei gh）.AnxAx xR xx xx 12111,0,0(,)(1),0;(,)(,)(2)max;(,)(,)(3)min.(,)定理设为对称矩阵
5、，记为定理设为对称矩阵，记为的特征值，则的特征值，则nnn nnnx Rxnx RxARAAx xxx xAx xx xAx xx x 计算矩阵的主特征根及对应的特征向量计算矩阵的主特征根及对应的特征向量Wait a second,what does that dominant eigenvalue mean?That is the eigenvalue with the largest magnitude.Why in the earth do I want to know that?Dont you have to compute the spectral radius from time
6、 to time?原始幂法原始幂法条件：条件：A 有特征根有特征根|1|2|n|0，对应对应n个线性无个线性无关的特征向量关的特征向量nxx,.,1思路：思路：从任意从任意 出发出发0)0(0,11)0(niiix )0()1(A niiiix1 )1()2(A niiiix12 niikiikniikiikkxxA1111)1()(|i/1|2n,这时,(8.1)式可写成若a10,则对充分大的k有)()(11222111)(nknkkknaaaxxxv111)(xvakk)(11111)1(,kkkavxv因而有 nivvkiki,2,1/)()1(1或取)(1)()1()(2)1(2)(1
7、)1(11knknkkkkvvvvvvn而特征向量 x1 v(k).乘幂法的收敛速度取决于|2/1|的大小.求矩阵A的按模最大的特征值 解 取v(0)=(1,0)T,计算v(k)=Av(k-1),结果如下例例261515141Akv1(k)v2(k)v1(k)/v1(k-1)v2(k)/v2(k-1)01010.250.220.102500.0833330.410.4166530.0422920.0343890.412600.4126740.0174510.0141900.412630.41263 可取0.41263,x1(0.017451,0.014190)T.对非零向量v,用max(v)表
8、示v的按绝对值最大的分量,称向量u=v/max(v)为向量v的规范化向量.例如,设向量v=(2,1,-5,-1)T,则max(v)=-5,u=(-0.4,-0.2,1,0.2)T.可见规范化向量u总满足u=1.乘幂法的规范化计算公式为:)1()(kkAuv 任取初始向量u(0)=v(0)0,计算)max()(kkv,3,2,1,)()(kkkkvu由于1)0()1(Avu21)0(22)1()2(vAAuu，)(max()(21121111niikiniikiiiaaaaxxxxkkkkkkkk.21)0(1)2(2)1()(vAuAAuu)max()0()0(vAvAkk所以)max()ma
9、x(lim111111)(xxxxuaakk又由)(max()(2111211111niikiniikiiiaaaaxxxx)max()0(1)0(1vAvAAuvkkkk）（）（其收敛速度由比值|2/1|来确定,其值越小收敛越快.)(max)(max2111211111niikiniikiiiaaaaxxxx)max()(kkv所以1limkk因此,当|k-k-1|r+1n,这时,(8.1)式可写成)()(1111122111)(nknrkrrrkknraaaaaxxxxxv若a1,a2,ar不全为零,则对充分大的k有 22111)(rrkkaaaxxxv由于a1x1+a2x2+arxr 是
10、对应1的特征向量,若仍记为x1,则有:v(k)1kx1,故前面的结论仍然成立.3.设1=-2,且|1=|2|3 n,这时,(8.1)式可写成)()()1(1133122111)(nknkrkkknaaaaxxxxv则对充分大的k有 v(2i)12i(a1x1+a2x2),v(2i+1)12i+1(a1x1-a2x2)1(22111)(xxvaakkk于是有nivvkiki,2,1,/)()2(1 x1v(k+1)+1v(k),x2v(k+1)-1v(k)对于规范化的幂法，由于 u(k+2)=v(k+2)/k+2=Au(k+1)/k+2 =Av(k+1)/k+1k+2=A2u(k)/k+1k+2
11、于是有，211kk12 x1k+1u(k+1)+1u(k),x2k+1u(k+1)-1u(k)的按模最大特征值和相应的特征向量。例例4 用乘幂法求矩阵13162216114A 解 取初始向量u(0)=v(0)=(1,1,2)T,计算可得K ku(k)012345678910111213113.5536284.6792043.4611244.6354653.4526554.6321163.4543154.6319293.4542914.6319203.4542884.631924(1,1,2)T(0.454545,0.909091,1)T(0.537222,0.972116,1)T(0.4652
12、01,0.994041,1)T(0.539392,0.998269,1)T(0.465721,0.999627,1)T(0.539487,0.999892,1)T(0.465890,0.999975,1)T(0.539495,0.999993,1)T(0.465893,0.999999,1)T(0.539495,1,1)T(0.465893,1,1)T(0.539495,1,1)T(0.465893,1,1)T1.2 加速技术加速技术由于21()1max()()(8.2)kkkvo所以,乘幂法收敛速度取决于比值|2/1|,当|2/1|1时,收敛是很慢的.1.Aitken 加速方法由(8.2)式
13、可知 x2=13u(13)-1u(12)=(0,0.631924,0.631924)T.4,4631924.4454288.3213121 x1=13u(13)+1u(12)=(4.315961,8.631924,8.631924)T,实际上,A的特征值为1=4,2=-4,3=1.可见,序列k线性收敛于1.会达到加速收敛的目的.0lim12111kkk 构造Aitken序列kkkkkkk12212)(如把Aitken加速方法用于例3,则有 ku(k)107.26.56.0033526.0016756.000837(1,1,1)T(1,0.8,0.1)T(1,0.75,-0.111)T(1,0.
14、730769,-0.188034)T.(1,0.714405,-0.249579)T(1,0.714346,-0.249790)T(1,0.714316,-0.249895)Tk01231011126.2666676.0000176.0000036.000000k 2.原点位移法 作矩阵B=A-pI,则B的特征值为mi=i-p(i=1,2,n),而且对应的特征向量相同.则对B应用乘幂法可达到加速收敛的目的。解 由于A的特征值为1=6,2=3,3=2,故取p=2.5,则B的特征值为m1=3.5,m2=0.5,m3=-0.5,则如果选取p,使m1仍然是B的按模最大特征值，且满足取初始向量u(0)=
15、v(0)=(1,1,1)T,由规范化计算公式:121212ppmm例例5 用原点位移法求例3中矩阵A的按模最大的特征值和特征向量.5.00105.1050145.65.2 EAB计算可得,3,2,1,)()(kkkkvu)max()(kkv)1()(kkBuvkkU(k)01234567.53.7666623.5353963.5050023.5007183.500102(1,1,1,)T(1,0.733333,-0.2)T(1,0.716814,-0.238938)T(1,0.714643,-0.249061)T(1,0.714337,-0.249777)T(1,0.714293,-0.249
16、981)T(1,0.714287,-0.249995)T这是因为|2/1|=1/2,而|m2/m1|=1/7,故对B应用乘幂法远比对A应用乘幂法收敛的快.反幂法是求矩阵按模最小的特征值和相应特征向量的方法.取,16+2.5=6.000102,x1u(6)=(1,0.714287,0.249995)T1.3 反幂法反幂法 设A是n阶非奇异矩阵,其特征值为|1|2|n-1|n|0对应的特征向量为x1,x2,xn,则有A-1的特征值为1211111nnn对应的特征向量为xn,xn-1,x1.要想求n和xn只需对A-1应用乘幂法,任取初始向量u(0)=v(0)0,作也可将上式改写成,3,2,1,)()
17、(kkkkvu)max()(kkv)1(1)(kkuAv)1()(kkuAv)max()(kkv()()(8.3),1,2,3,kkkkuv式(8.3)称为反幂法.显然有)max(/lim,1lim)(nnkknkkxxu每一步求v(k)需要求解线性方程组,可采用LU分解法求解.反幂法还可结合原点位移法应用.设已求得矩阵A的特征值i的某个近似值对B应用反幂法可求出精度更高的i和xi.设已求得例3中矩阵A的特征值的近似值16.003,和相应的特征向量x1(1,0.714405,-0.249579)T,试用带原点位移的反幂法求1和x1的更精确的值.iip,取 作原点位移,令B=A-E,i 则B的特
18、征值为 )(,221ijijiiiniii且必有例例6 解 取p=6.003,作矩阵B=A-6.003I,则003.4010997.65014003.10B取初始向量u(0)=(1,0.714405,-0.249579)T,对B用反幂法计算可得:可见收敛速度非常快,这是因为B的3个特征值为1=-4.003,2=-3.003,3=-0.003,|3/2|0.000999很小.T)250000.0,714286.0,1(,00000167.6003.61)1(1uT)250000.0,714286.0,1(,000000007.6003.61)2(2u 反幂法反幂法 /*Inverse Power
19、 Method*/Ch.5 Power Method Inverse Power Method若若 A 有有|1|2|n|，则，则 A 1 有有对应同样一组特征向量。对应同样一组特征向量。11111 nnA 1 的主特征根的主特征根 A的绝对值最小的特征根的绝对值最小的特征根Q:How must we compute in every step?)(1)1(kkA A:Solve a linear system with A factorized.)()1(kkA 若知道某一特征根若知道某一特征根 i 的大致位置的大致位置 p，即对任意，即对任意 j i 有有|i p|j p|，并且如果，并且
20、如果(A pI)1存在，则存在，则可以用反幂法求可以用反幂法求(A pI)1的主特征根的主特征根 1/(i p)，收，收敛将非常快。敛将非常快。思思路路 Jacobi方法是求实对称矩阵全部特征值和特征向量的一种矩阵变换方法。2 Jacobi 方法方法-正交变换法正交变换法-QR分解分解 实对称矩阵A具有下列性质：（1）A的特征值均为实数；（2）存在正交矩阵R，使RTAR=diag(1,2,n),而 R的第i个列向量恰为i的特征向量;直接求正交矩阵R是困难的.Jacobi提出用一系列所谓平面旋转矩阵逐次将A约化为对角矩阵.平面解析几何中的平面坐标旋转变换表示平面上坐标轴旋转角的变换.（3）若记A
21、1=RTAR,则A1仍为对称矩阵.2.1 平面旋转矩阵平面旋转矩阵 1122cossinsincosyxyx 在三维空间直角坐标系中,ox1y1平面绕着oz1轴旋转角的坐标变换为1112221000cossin0sincoszyxzyx Rpq()具有下列性质:一般地,在n维向量空间Rn中,沿着xpyq平面旋转角的变换矩阵为行第行第qppq11cossin11sincos11)(R称Rpq()为平面旋转矩阵.设实对称矩阵A=(aij)nn,记B=RpqT()ARpq()=(bij)nn则它们元素之间有如下关系:（1）Rpq()为正交矩阵,即Rpq-1()=RpqT();（2）如果A为对称矩阵,
22、则RpqT()ARpq()也为对称矩阵,且与A有相同的特征值.（3）RpqT()A仅改变A的第p行与第q行元素,ARpq()仅改变A的第p列与第q列元素.222212cossinsin2sincossin2()sin2cos2(8.4)cossinsincos(,)ppppqqpqqqppqqpqpqqpqqpppqjppjjpjqjqqjjpjqijijbaaabaaabbaaabbaabbaabai jp q 222222(8.5)22ppqqpqppqqpqbbbaaa所以有从而2222qipiqipiaabb),(2222qpiaabbiqipiqip22FFAB2211,(8.6)n
23、nijijijijba，即有(8.5)、(8.6)式可得222222pqjiijpqjiijaabb21221222pqniiipqniiiaabb 如果apq0,适当选取角,使02cos2sin)(21pqppqqqppqaaabb只需角满足从而 如果取|apq|=若记cot2,|(8.7)24ppqqpqaaajiijpqjiijjiijaaab22222niiipqniiiniiiaaab12212122|maxijjia于是jiijpqanna2211）（，则jiijjiijannb22)1(21(,)(2jiijaA 则上式可记为2()(1)()(8.8)(1)n nBA 由式(8.
24、7),令t=tan,则t满足方程 t2+2t-1=0 经典Jacobi算法是对A(0)=A施行一系列平面旋转变换:为保证|/4,取绝对值较小的根,有于是2sgn()/(|1),0(8.9)1,0t,)1(cos212tsincos(8.10)t2.2 Jacobi 方法方法 A(1)=R1TA(0)R1,A(2)=R2TA(1)R2,A(k)=RkTA(k-1)Rk,每一步变换选择A(k-1)=(aij(k-1)nn 的非对角线元素中绝对值最大者apq(k-1)(称为主元素)作为歼灭对象,构造平面旋 是给定的精度要求,则A的特征值可取为iaii(k),i=1,2,n.转矩阵Rk=Rpq(),经
25、变换得到A(k)=(aij(k)nn,且apq(k)=0,这时由(8.8)式有从而由此递推得到 当k充分大时,或者(A(k),或者)()1(21()()1()(kknnAA)()1(21()()(AAkknn)(,0k),(lim21)(nkkdiagDA 另外,由于 A(k)=RkTA(k-1)Rk=RkTRk-1TR1TAR1R2Rk=RTAR,|max)(kijjia的全部特征值.解 记 A(0)=A,取p=1,q=2,apq(0)=a12(0)=2,于是有因此,R=R1R2Rk 的列向量xj(j=1,2,n)为A的近似特征向量.例例7 用Jacobi 方法计算对称矩阵612152224
26、A25.02)0(12)0(22)0(11aaa780776.0)1|/(|)sgn(,2t788206.0)1(cos212t615412.0cossin,t从而有所以 再取p=2,q=3,apq(1)=a23(1)=2.020190,类似地可得以下依次有1000788206.0615412.00615412.0788206.01000cossin0sincos)(1pqRR6020190.2961.0020190.2561552.60961.00438448.21)0(1)1(RARAT241166.40724794.00320386.8631026.0724794.0631026.043
27、8448.2)2(A496424.4209614.00209614.0320386.8595192.00595192.0183185.2)3(A496424.4208653.0020048.0208653.0377576.80020048.00125995.2)4(A485239.40020019.00388761.8001073.0020019.0001073.0125995.2)5(A485401.4000009.0001072.0000009.0388761.800001072.0125825.2)6(A485401.4000009.00000009.0388761.8000125825
28、.2)7(A从而A的特征值可取为 12.125825,28.388761,34.485401 为了减少搜索非对角线绝对值最大元素时间,对经典的Jacobi方法可作进一步改进.1.循环Jacobi方法:按(1,2),(1,3),(1,n),(2,3),(2,4),(2,n),(n-1,n)的顺序,对每个(p,q)的非零元素apq作Jacobi变换,使其零化,逐次重复扫描下去,直至(A)为止.2.过关Jacobi方法:取单调下降收敛于零的正数序列k,先以1为关卡值,依照1中顺序,将绝对值超过1的非对角元素零化,待所有非对角元素绝对值均不超过1时,再换下一个关卡值2,直到关卡值小于给定的精度.Jacobi方法具有方法简单紧凑,精度高,收敛较快等优点,是计算对称矩阵全部特征值和相应特征向量的有效方法,但计算量较大,一般适用于阶数不高的矩阵.
