第九章矩阵特征值问题的数值方法

第9章矩阵特征值问题的数值方法9.1 特征值与特征向量9.2 Hermite矩阵特征值问题9.3 Jacobi方法9.4 对分法9.5 乘幂法9.6 反幂法9.7 QR方法9.1 特征值与特征向量设A是n阶矩阵，x是非零列向量. 如果有数λ存在，满足， (1)那么，称x是矩阵A关于特征值λ的特征向量.如果把(1)式右端写为 ，那么(1)式又可写为: x λ()0I A x λ-=||0I A λ-=即1110()||...n n n f I A a a a λλλλλ--=-=++++记 它是关于参数λ的n 次多项式，称为矩阵A 的特征多项式， 其中a 0=(-1)n ｜A ｜.(2)显然，当λ是A的一个特征值时，它必然是的根. 反之，如果λ是的根，那么齐次方程组(2)有非零解向量x，使(1)式成立. 从而，λ是A的一个特征值.A的特征值也称为A的特征根.()0 fλ=()0 fλ=矩阵特征值和特征向量有如下主要性质：定理9.1.1 n阶矩阵A是降秩矩阵的充分必要条件是A有零特征值.定理9.1.2 设矩阵A与矩阵B相似，那么它们有相同的特征值.定理9.1.3 n阶矩阵A与A T有相同的特征值.定理9.1.4 设λ≠λj是n阶矩阵A的两个互异特i征值，x、y分别是其相应的右特征向量和左特征向量，那么，x T y=0 .9.2 Hermite矩阵特征值问题•设A为n阶矩阵，其共轭转置矩阵记为A H. 如果A=A H，那么，A称为Hermite矩阵.9.2.1 Hermite 矩阵的有关性质设是Hermite 矩阵A 的n 个特征值. 有以下性质:•全是实数. 12,,...,n λλλ12,,...,n λλλ•有相应的n 个线性无关的特征向量，它们可以化为一组标准酉交的特征向量组 ,即 12,,...,n λλλ12,,...,n u u u H i j u u ijδ=•是酉空间中的一组标准酉交基. 12,,...,n u u u•记U=( ),它是一个酉阵，即U H U=UU H =I ，那么即A 与以 为对角元的对角阵相似. 12,,...,n u u u 1H n U AU D λλ⎛⎫ ⎪== ⎪ ⎪⎝⎭12,,...,n λλλ•A 为正定矩阵的充分必要条件是 全为正数.12,,...,n λλλ定理9.2.1 设 是Hermite 矩阵A 的n个特征值，那么证: 12,,...,n λλλ21max i i n A λ≤≤=21n i F i A λ==∑2222()()(())H A A A A A ρρρ===由21max ii n A λ≤≤=因此2221()()n H i F i A tr A A tr A λ====∑又由21n iF i A λ==∑得设x 是一个非零向量，A 是Hermite 矩阵，称 为矩阵A 关于向量x 的Rayleigh 商，记为R(x).H H x Ax x x 定理9.2.2 如果A 的n 个特征值为其相应的标准酉交的特征向量为 那么有12...n λλλ≥≥≥12,,...,n u u u 1()n R x λλ≥≥定理9.2.3 设A 是Hermite 矩阵 ，那么 100min ()min ()k n k k k x C x x C x R x R x λλ-+∈≠∈≠==且且或9.2.2 极值定理定理9.2.4(极值定理) 设Hermite 矩阵的n 个特征值为 ，其相应的标准酉交特征向量为 .用C k 表示酉空间C n 中任意的k 维子空间，那么12...n λλλ≥≥≥12,,...,n u u u 1100max min ()max min ()k k n k n k k C x C x k C x C x R x R x λλ-+-+∈≠∈≠==且且或9.2.3 Hermite矩阵特征值问题的性态矩阵特征值问题与求解线性方程组问题一样，都存在当矩阵A的原始数据有小变化(小扰动)时，引起特征值问题的变化有大有小的问题，如果引起的变化小，称该特征值问题是良态的.反之，称为病态的.矩阵特征值问题的性态是很复杂的，通常分别就单个特征值或整体特征值给出状态数进行分析. 对于Hermite矩阵，由于其特征值问题的特殊性质，其特征值都是良态的.下面先证明Hermite矩阵特征值的扰动定理.定理9.2.5 设矩阵A ，E ，A+E 都是n 阶Hermite 矩阵，其特征值分别为那么, 证 设矩阵A 关于特征值λ1，λ2，…，λn 的标准酉交特征向量为u 1，u 2，…，u n ， 是由u i ，u i+1，…，u n 生成的n-i+1维子空间. 对 中任意非零向量x,由极值定理,有 12n λλλ≥≥≥12n εεε≥≥≥12n μμμ≥≥≥1i n i i λεμλε+≤≤+1n i C -+1n i C -+111000()max max max n i n i n i Hi H x C x H H H H x C x x C x x A E x x xx Ax x Ex x x x x μ-+-+-+∈≠∈≠∈≠+≤=+且且且由定理9.2.3， 又由定理9.2.2，对任意x≠0，有从而有 另一方面, A=(A+E)-E. 记 为矩阵-E 的特征值，那么， 重复上面的过程,可得 从而有 10max n i Hi H x C x x Ax x xλ-+∈≠=且110max n i H n H x C x x Ex x x εε-+∈≠≥≥且1i i μλε≤+12nδδδ≥≥≥1i n i δε-+=-1i i λμδ≤+i i n μλε≥+定理9.2.5通常又称为Hermite 矩阵特征值的扰动定理定理9.2.6 设矩阵A 和A′=A+E 都是n 阶Hermite 矩阵，其特征值分别为和 ，那么这个定理表明，扰动矩阵E 使A 的特征值的变化不会超过 ‖E‖2. 一般‖E‖2Hermite 矩阵特征值是良态的.12n λλλ≥≥≥12nλλλ'''≥≥≥222i i E E λμλ-≤≤+9.3 Jacobi方法理论上，实对称矩阵A正交相似于以A的特征值为对角元的对角阵. 问题是如何构造这样的正交矩阵呢? Jacobi方法就是通过构造特殊的正交矩阵序列，通过相似变换使A 的非对角线元素逐次零化来实现对角化的.9.3.1 平面旋转矩阵与相似约化先看一个简单的例子.设 是二阶实对称矩阵，即a 21=a 12，其特征值为λ1，λ2. 令使得 记 容易验证B T =B, 且11122122a a A a a ⎛⎫= ⎪⎝⎭cos sin sin cos R θθθθ⎛⎫-= ⎪⎝⎭11T R AR λλ⎛⎫= ⎪⎝⎭11122122T b b B R AR b b ⎛⎫== ⎪⎝⎭22111112222212212211122222111222cos 2sin cos sin ()sin cos (cos sin )sin 2sin cos cos b a a a b b a a a b a a a θθθθθθθθθθθθ=++==-+-=-+解之得:当 时 当 时并规定 1122a a ≠12112221arctan 2a a a θ=-1122a a =21arctan 2pq pp qq a a a θ=-4πθ≤9.3.2 经典的Jacobi 方法 设A 是实对称矩阵，记A 1=A.Jacobi 方法的基本思想是用迭代格式A k+1=Q T k A k Q k , k=1,2,…构造一个相似矩阵序列，使｛A k ｝收敛于一个对角阵. 其中 Q k 为平面旋转矩阵，其旋转角θk 由使A k 的绝对值 最大元a (k)pq =a (k)qp =0 或按列依次使A 的非对角元 零化来确定.定理9.3.1 设A 是n 阶实对称矩阵，那么由Jacobi 方法产生的相似矩阵序列｛A k ｝的非对角元收敛于0. 也就是说，｛A k ｝收敛于以A 的特征值为对角元的对角阵.记其中E k 是A k 除主对角元外的矩阵.由平面旋转矩阵的性质中,对于 ,有 因此, ()()k k ii k A diag a E =+1T k pq k pq A R A R +=,i p q ≠(1)2(1)2()2()2()()()()k k k k ip iq ip iq a a a a +++=+22()212()k k k pq F F E E a +=-又由假设,因此, 这样,便有从而,当 22()()()()1,,1max (1)n k k k k pq ij pq pq i j n i j i j i j a a a n n a ≤≤≠=≠==-∑且且且22()(1)k k pqF E n n a ≤-222112222(1)(1)k k k F F F E E E n n n n +≤-≤---1,0k F k E +→∞→9.3.3 实用的Jacobi 方法•循环Jacobi 方法必须一次又一次扫描，才能使｛A k ｝收敛于对角阵 ，计算量很大. 在实际计算中，往往用一些特殊方法来控制扫描次数，减少计算量. 下面介 绍一种应用最为广泛的特殊循环Jacobi 方法——阈Jacobi 方法. 阈Jacobi 方法首先确定一个阈值δ，在对非对角元零化的一次扫描中，只对其中绝对值 超过阈值的非对角元进行零化. 当所有非对角元素的绝对值都不超过阈值后，将阈值减少， 再重复下一轮扫描，直至阈值充分小为止. 减少阈值的方法通常是先固定一个正整数M≥n ，扫描一次后，让 . 而阈值的下界是根据实际问题的精度要求选定的.M δδ→9.3.4 用Jacobi 方法计算特征向量 •假定经过k 次迭代得到A k+1=R T k …R T 1AR 1…R k ，(15) 这时A k+1是满足精度要求的一个近似的对角阵. 如果记Q k =R 1R 2…R k =Q k-1R k ，(16)那么，Q k 是一个正交矩阵，且(15)式又可表示为A k+1=Q T k AQ k .当A k+1的非对角元素充分小，Q k 的第 j 列q j 可以看成是近似特征值a (k+1)jj 相应的特征向量了.在实际计算中，可以按(16)式在迭代过程中形成Q k ，把Q k 看成是Q k-1右乘一个平面旋转矩阵得到. 不妨记 Q 0=I ，Q k 的元素按下式计算：()(1)(1)()(1)(1)()(1)cos sin ,sin cos ,,,,1,2,k k k ip ip k ip k k k k iq ip k iq k k k ij ijqq q qq q q q j p q i nθθθθ-----=+=-+=≠=9.4 对分法理论上，一个实对称矩阵正交相似于一个以其特征值为对角元的对角阵. 但是，经典的结果告诉我们，一个大于4次的多项式方程不可能用有限次四则运算求根. 因此，我们不可能期望只用有限次相似变换将一个实对称矩阵约化为一个对角阵.下面先介绍将一个实对称矩阵相似约化为实对称三对角矩阵的方法，再讨论求其特征值的对分法.9.4.1 相似约化为实对称三对角矩阵 将一个实对称矩阵正交相似约化为一个实对称三对角矩阵的算法，可归纳如下：记A (1)=A,对k=1,2,…,n -2①按(4)式、(5)式和(8)式计算； ②按(9)～(12)式，计算A (k+1).,k k k u σβ'和()1,()2,()(4)k k k k k k k k k nk a a u a σ++⎛⎫- ⎪ ⎪'= ⎪ ⎪ ⎪⎝⎭()()21,1()()(5)nk k k k k ik i k sign a a σ+=+=-∑()1,()(8)k k k k k k a βσσ+=-1()1(1)(),(9),(10)1,(11)2(),(12)k k k k T k k kk k k k k k k T T k k k k g A u u g y g u A A u y y u βεβε--+===-=-+9.4.2 Sturm 序列的性质 •设实对称三对角矩阵为其中βi ≠0 (i=1,2,…，n-1) 其特征矩阵为T-λI. 记T-λI 的第i 阶主子式为111221n n a T a a ββββ-⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭111211()i i i i a a p a λββλλββλ----=-•这是关于λ的i 次多项式，当i=n 时， p n (λ)=｜T-λI ｜是矩阵T 的特征多项式. 令p 0(λ)≡1，则有p 1(λ)=α1-λ，p i (λ)=(αi -λ)p i-1(λ)-β2i-1p i-2(λ)，i=2,3，…，n.(15)•多项式序列｛p i (λ)｝ (i=0,1,…，n)称为Sturm 序列111211()i i i i a a p a λββλλββλ----=-定理9.4.1｛p i (λ)｝ (i=1,2,…，n)的根都是 实根. 证 由(14)式，p i (λ)是i 阶实对称矩阵的特征多项式，因此，｛p i (λ)｝ (i=1,2,…，n)的根全是实根. 定理9.4.2 定理9.4.2 设α是p i (λ)的一个根，那么 ①p i-1(α)p i+1(α)≠0，即相邻的两个多项式无公共根； ②p i-1(α)p i+1(α)<0，即p i-1(α)与p i+1( α)反号. lim ()0,i p λλ→-∞>lim ()0,lim ()0,i i p p λλλλ→+∞→+∞<>当i 为奇数当i 为偶数定理9.4.4 p i (λ)的根都是单根，并且将p i+1(λ)的根严格隔离.9.4.3 同号数和它的应用•定义1 设p 0(λ)≡1，｛p i (λ)｝(i=1,2,…，n) 是一个Sturm 序列，称相邻的两个数中符号一致的数目为同号数，记为a i (λ). 若某个p i (λ)=0，规定与p i-1(λ)反号.•定理9.4.5 设两个实数x<y ，那么，形如(13)式的实对称 三对角矩阵T 的特征多项式在区间(x,y ］上根的数目为a(x)-a(y).9.4.4 求Hermite矩阵特征值的对分法•对分法的计算可归纳为以下4个部分①确定(13)式的矩阵T的全部特征值的分布区间.②在区间［a,b］中，用区间对分的方法找出只含T的一个特征值的子区间.③在只含一个特征值的子区间上的对分法.④同号数的计算.9.5 乘幂法乘幂法是适用于求一般矩阵按模最大特征值及相应特征向量的算法. •设A 是n 阶矩阵，其n 个特征值按模从大到小排序为 又假设关于λ1，λ2，…，λn 的特征向量v 1，v 2，…，v n 线性无关.123n λλλλ>≥≥≥201112211[()()]k k k k n n A x a v a v v λλλλλ=+++x k →λk 1a 1v 1 (k→∞).因此，x k 可看成是关于特征值λ1的近似特征向量.迭代格式为 1,0,1,k k k k k x z x x Az k ∞+===按模最大特征值λ1及其相应的特征向量v 1的乘幂法的计算公式：111,0,1,kk k k kT k k k T kk x z x x Az z Az k z z λ∞++====9.5.2 收缩方法•设矩阵A 的n 个特征值按模从大到小排序为，其相应的n 个线性无关特征向量为v 1，v 2，…，v n . 在计算A 的最大特征值λ1及相应特征向量v 1后，可以通过收缩方法，继续用乘幂法计算λ2及其相应的特征向量v 2.123n λλλλ>>≥≥•定义n 阶矩阵•把去掉A 1的第1行和第1列的n-1阶矩阵记为2121112221122211111111211211000n n T n n n n nn n n a v a a v a a v a AA v a v a a v a a v a α⎛⎫ ⎪--- ⎪=-= ⎪ ⎪---⎝⎭222112232113221132311233311333112112311311n n n n n n n n nn n n a v a a v a a v a a v a a v a a v a B a v a a v a a v a ---⎛⎫ ⎪--- ⎪= ⎪ ⎪---⎝⎭•那么，B 有与A 1除λ1外的相同的n-1个特征值 ｜λ2｜>｜λ3｜≥…≥｜λn ｜，可以用乘幂法计算λ2及其相应的 特征向量. 在计算λ1和v １后，按(15)式形成n-1阶矩阵B 的计算过程称为收缩方法.9.6 反幂法•反幂法可以求一个非奇异矩阵A 的逆矩阵A -1的按模最小的特征值及相应的特征向量，又可以求A 的一个近似特征值相应的特征向量. •9.6.1 求按模最小特征值及相应特征向量的反幂法,又称为反迭代法. 1(1)1,0,1,k k k k k T k k k Tk kx z x LUx z z x k z z ξ∞+++====9.6.2 求近似特征值的特征向量的反幂法•先对矩阵进行LU 分解，记 那么， (7)下面介绍一种选取特殊的初始向量x 0的反幂法——半迭代法.l I A λ-l I A LU λ-=1,0,1,k k k k k k k x z x Ly z Ux y k ∞+====•假设 ,选取初始向量x 0满足‖x 0‖∞=1，这时z 0=x 0.对照(7)式中的第二个式子.可把z 0看成满足Le=z 0.(8)这里，e=(1,1，…，1)T ，而z 0的各个分量的取值多少是无关重要的.这样，在第一个迭代步的计算中，只需求解(7)式中的上三角方程组Ux 1=e. “半迭代法”的命名也由此而得. l I A λ-9.7 QR 方法定理9.7.1设A 是n 阶矩阵，其n 个 特征值为 .那么存在一个酉矩阵U ，使 U ＨAU 是以为 对角元的上三角矩阵. 12,,,n λλλ12,,,n λλλ9.7.1 两个基本定理 定理9.7.2设A 是n 阶实矩 阵，那么，存在一个正交矩阵Q ，使Q ＴAQ 为一个准上三角矩阵，它的对角元是A 的一个特征值，对角元上的二阶块矩阵的两个特征值是A 的一对共轭复特征值.9.7.2 相似约化为上Hessenberg矩阵对一般n阶矩阵，QR算法的每一个迭代步需要O(n３)次乘法运算.如果矩阵阶数稍大，这个算法几乎没有实际的应用价值.通常采用的方法是先将矩阵相似约化为上Hessenberg形式的矩阵，在此基础上应用QR 迭代.这时，一个QR迭代步的乘法运算次数只需O(n２）次.所谓上Hessenberg 矩阵是指一个n 阶矩阵A ，如果当i>j+1时，a ij =0，称A 为上Hessenberg 矩阵.例如：一个5阶的上Hessenberg 矩阵具有如下的形式：**********0****00***000**A ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭下面介绍QR 方法时，都假设矩阵A 是一个上Hessenberg 矩阵.9.7.3 QR算法•设A是n阶矩阵且有QR分解A＝QR，(2)这里，Q是酉矩阵，R是上三角矩阵.如果A 是满秩并规定R有正对角元，这个分解是惟一的.一、QR 算法的基本思想 •记A ＝A １且有A １＝Q 1R 1.将等号右边两个矩阵因子的次序交换，得A ２＝R １Q １，且，(3)即A ２～A １. 不难证明: 即A k+1～A k ～…～A １，矩阵序列｛A k ｝有相同的特征值.记 12111A Q AQ -=111111k k k k kk k A Q A Q Q Q A Q Q ---+==12k k Q Q Q Q =12k kR R R R =容易得到 是A k 的一个QR 分解 k k k A Q R 如果A 是一个满秩的上Hessenberg 矩阵，可以证明，经过一个QR 迭代步得到的A 2＝Q －１１A 1Q １仍然是上Hessenberg 矩阵. 因为上Hessenberg 矩阵次对角线以下的元素全为0，因此，只要证明，当k→∞时，由迭 代格式(4)产生的矩阵A k 的次对角元趋向于零就可以了.二、 QR算法的收敛性•定理9.7.3设n阶矩阵A的n个特征值满足|λ１|>|λ２|>…>|λn|>0，其相应的n个线性无关特征向量为x1，x2，…，xn.记X＝(x1,x2,…,x n), Y= X－１.如果Y存在LU分解，那么，由(4) 式产生的矩阵A k基本收敛于上三角矩阵R.这里，基本收敛的含义指｛Ak ｝的元素中除对角线以下的元素趋于零外，可以不收敛于R的元素.三、 QR 算法的迭代过程•1. 一个QR 迭代步的计算①对l=1，2，…，n-1，构造n-1个平面旋转矩阵P l,l+1,使A 1的次对角元全部零化，实现A 1的QR 分解的计算，这里, ,1,1,1,11,1,,1,2,l l l l ll lj l l l l l l l l c s a a j l l n s c a a ++++++⎛⎫⎛⎫⎛⎫→=++ ⎪⎪ ⎪-⎝⎭⎝⎭⎝⎭1,,1,122221,1,,l l ll l l l l ll l l ll l l a a c s a a a a +++++==++221,ll l lr a a +=+②用P l,l+1右乘(24)，所得结果也放回矩阵A 相应的元素中.,1,1,,1,1,1,1(,)(,),1,2,,1l l l l i l i l il i l l l l l c s a a a a s c i l ++++++-⎛⎫→ ⎪⎝⎭=+2. QR 算法的迭代控制 当迭代步数k 充分大时，由迭代格式(4)产生的A k 的次对角元趋于0.在 实 际计算中，控制迭代次数常用的一种办法是，预先给定一个小的正数ε，在一个迭代步的计 算结束后，对l=n-1, n-2,…，1，依次判别次对角元的绝对值是否满足 或更严格的准则是 或不太严格的准则是 如果上面三个不等式中有一个成立，把 看做实际上为零. 1,l l a Aε+≤1,1,1min{,}l l ll l l a a a ε+++≤1,1,1{}l l ll l l a a a ε+++≤+1,l l a +。