CH矩阵特征值问题计算实用

会计学1CH矩阵特征值问题计算实用第1页/共41页第2页/共41页第3页/共41页计算矩阵的主特征根及对应的特征向量Waitasecond,whatdoesthatdominanteigenvalue

mean?Thatistheeigenvaluewiththelargestmagnitude.WhyintheearthdoIwanttoknowthat?Don’tyouhavetocomputethespectralradiusfromtimetotime?原始幂法条件：A有特征根|1|>|2|…|n|0，对应n个线性无关的特征向量思路：从任意出发………|i/1|<1当k

充分大时，有这是A关于1的近似特征向量§1乘幂法和反幂法第4页/共41页

§1.1

乘幂法

乘幂法是用来求矩阵A按模最大的特征值和相应的特征向量的方法.

设A是单构矩阵,即A有n个线性无关的特征向量.

A的n个特征值为

|12n对应的特征向量为

x1,x2,…xn

线性无关.我们要求1

和x1.第5页/共41页

乘幂法的基本思想是取初始向量v(0)Rn,作迭代

v(k+1)=Av(k)=Ak+1v(0),k=0,1,2,…产生迭代序列v(k).

由于x1,x2,…xn

线性无关,从而有

v(0)=a1x1+a2x2+…+anxn

故有

v(1)=Av(0)

v(2)=Av(1)

v(k)=Akv(0)=a11kx1+a22kx2+…+annkxn(8.1)……………=a11x1+a22x2+…+annxn

=a112x1+a222x2+…+ann2xn

第6页/共41页1.设|1>2n,这时,(8.1)式可写成若a10,则对充分大的k有因而有

或取

而特征向量x1v(k).乘幂法的收敛速度取决于|2/1|的大小.第7页/共41页

求矩阵A的按模最大的特征值

解取v(0)=(1,0)T,计算v(k)=Av(k-1),结果如下例2kv1(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.第8页/共41页

对非零向量v,用max(v)表示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.

乘幂法的规范化计算公式为:

任取初始向量u(0)=v(0)0,计算由于第9页/共41页所以又由第10页/共41页其收敛速度由比值|2/1|来确定,其值越小收敛越快.所以因此,当|k-k-1|<时,可取:1k,x1u(k).第11页/共41页

如用规范化乘幂法解例2,仍取u(0)=v(0)=(1,0)T,则有

故可取10.412627,x1(1,0.813138)T.k01234k0.250.410.4126020.412627u1(k)11111u2(k)00.80.8130080.8131360.813138用乘幂法求A的按模最大的特征值和相应特征向量.例3设

解取初值u(0)=v(0)=(1,1,1)T,计算得第12页/共41页kku(k)0123…101112107.26.5…6.0033526.0016756.000837(1,1,1)T(1,0.8,0.1)T(1,0.75,-0.111)T(1,0.730769,-0.188034)T…..(1,0.714405,-0.249579)T(1,0.714346,-0.249790)T(1,0.714316,-0.249895)T

可取16.000837,x1(1,0.714316,-0.249895)T.

实际上,A的3个特征值分别为1=6,2=3,3=2.第13页/共41页2.设1=2==r,且|1>r+1n,这时,(8.1)式可写成若a1,a2,…,ar不全为零,则对充分大的k有由于a1x1+a2x2+…+arxr是对应1的特征向量,若仍记为x1,则有:v(k)1kx1,故前面的结论仍然成立.3.设1=-2,且|1=|2|>3n,这时,(8.1)式可写成则对充分大的k有第14页/共41页v(2i)12i(a1x1+a2x2),v(2i+1)12i+1(a1x1-a2x2)于是有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于是有第15页/共41页x1k+1u(k+1)+1u(k),x2k+1u(k+1)-1u(k)的按模最大特征值和相应的特征向量。例4

用乘幂法求矩阵

解取初始向量u(0)=v(0)=(1,1,2)T,计算可得第16页/共41页Kku(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.465201,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)T第17页/共41页§1.2加速技术由于所以,乘幂法收敛速度取决于比值|2/1|,当|2/1|1时,收敛是很慢的.

1.Aitken加速方法由(8.2)式可知x2=13u(13)-1u(12)=(0,0.631924,0.631924)T.x1=13u(13)+1u(12)=(4.315961,8.631924,8.631924)T,

实际上,A的特征值为1=4,2=-4,3=1.第18页/共41页可见,序列k线性收敛于1.会达到加速收敛的目的.

构造Aitken序列

如把Aitken加速方法用于例3,则有第19页/共41页ku(k)107.26.5…6.0033526.0016756.000837(1,1,1)T(1,0.8,0.1)T(1,0.75,-0.111)T(1,0.730769,-0.188034)T…..(1,0.714405,-0.249579)T(1,0.714346,-0.249790)T(1,0.714316,-0.249895)Tk0123…1011126.266667………6.0000176.0000036.0000002.原点位移法

作矩阵B=A-pI,则B的特征值为mi=i-p(i=1,2,…,n),而且对应的特征向量相同.第20页/共41页则对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)=v(0)=(1,1,1)T,由规范化计算公式:例5

用原点位移法求例3中矩阵A的按模最大的特征值和特征向量.第21页/共41页计算可得kkU(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.249981)T(1,0.714287,-0.249995)T第22页/共41页这是因为|2/1|=1/2,而|m2/m1|=1/7,故对B应用乘幂法远比对A应用乘幂法收敛的快.

反幂法是求矩阵按模最小的特征值和相应特征向量的方法.取,16+2.5=6.000102,x1u(6)=(1,0.714287,0.249995)T§1.3反幂法

设A是n阶非奇异矩阵,其特征值为|1||2|…|n-1|>|n|>0对应的特征向量为x1,x2,…,xn,则有A-1的特征值为对应的特征向量为xn,xn-1,…,x1.

要想求n和xn只需对A-1应用乘幂法,任取初始向量u(0)=v(0)0,作第23页/共41页也可将上式改写成式(8.3)称为反幂法.显然有每一步求v(k)需要求解线性方程组,可采用LU分解法求解.第24页/共41页

反幂法还可结合原点位移法应用.设已求得矩阵A的特征值i的某个近似值对B应用反幂法可求出精度更高的i和xi.

设已求得例3中矩阵A的特征值的近似值16.003,和相应的特征向量x1(1,0.714405,-0.249579)T,试用带原点位移的反幂法求1和x1的更精确的值.

作原点位移,令B=A-E,

则B的特征值为例6

解取p=6.003,作矩阵B=A-6.003I,则第25页/共41页取初始向量u(0)=(1,0.714405,-0.249579)T,对B用反幂法计算可得:可见收敛速度非常快,这是因为B的3个特征值为1=-4.003,2=-3.003,3=-0.003,|3/2|0.000999很小.第26页/共41页反幂法

/*InversePowerMethod*/Ch.5PowerMethod–InversePowerMethod若

A有|1||2|…>|

n|，则A1有对应同样一组特征向量。11111lll…>-nnA1的主特征根A的绝对值最小的特征根Q:Howmustwecomputeineverystep?A:SolvealinearsystemwithAfactorized.若知道某一特征根i

的大致位置p

，即对任意j

i

有|

ip|<<|

jp|，并且如果(A

pI)1存在，则可以用反幂法求(A

pI)1的主特征根1/(ip

)，收敛将非常快。思路第27页/共41页Jacobi方法是求实对称矩阵全部特征值和特征向量的一种矩阵变换方法。§2Jacobi方法-正交变换法-QR分解

实对称矩阵A具有下列性质：

（1）A的特征值均为实数；

（2）存在正交矩阵R，使RTAR=diag(1,2,…,n),而第28页/共41页R的第i个列向量恰为i的特征向量;

直接求正交矩阵R是困难的.Jacobi提出用一系列所谓平面旋转矩阵逐次将A约化为对角矩阵.平面解析几何中的平面坐标旋转变换表示平面上坐标轴旋转角的变换.

（3）若记A1=RTAR,则A1仍为对称矩阵.§2.1平面旋转矩阵

在三维空间直角坐标系中,ox1y1平面绕着oz1轴旋转角的坐标变换为第29页/共41页Rpq()具有下列性质:

一般地,在n维向量空间Rn中,沿着xpyq平面旋转角的变换矩阵为称Rpq()为平面旋转矩阵.第30页/共41页

设实对称矩阵A=(aij)nn,记B=RpqT()ARpq()=(bij)nn则它们元素之间有如下关系:

（1）Rpq()为正交矩阵,即Rpq-1()=RpqT();

（2）如果A为对称矩阵,则RpqT()ARpq()也为对称矩阵,且与A有相同的特征值.

（3）RpqT()A仅改变A的第p行与第q行元素,ARpq()仅改变A的第p列与第q列元素.第31页/共41页所以有从而有(8.5)、(8.6)式可得

如果apq0,适当选取角,使第32页/共41页只需角满足从而

如果取|apq|=若记于是

则上式可记为第33页/共41页

由式(8.7),令t=tan,则t满足方程t2+2t-1=0

经典Jacobi算法是对A(0)=A施行一系列平面旋转变换:为保证||/4,取绝对值较小的根,有于是§2.2Jacobi方法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)(称为主元素)作为歼灭对象,构造平面旋第34页/共41页是给定的精度要求,则A的特征值可取为iaii(k),i=1,2,…,n.转矩阵Rk=Rpq(),经变换得到A(k)=(aij(k))nn,