第9章 矩阵特征值的数值解法

第9章 矩阵特征值的数值解法9.1 引言矩阵特征值问题有广泛的应用背景. 例如动力系统和结构系统中的振动问题、电力系统的静态稳定分析上、工程设计中的某些临界值的确定等，都归结为矩阵特征值问题. 数学中诸如方阵的对角化及解微分方程组等问题，都要用到特征值的理论. 本章介绍n 阶实矩阵n n ⨯∈R A 的特征值与特征向量的数值解法.定义9.1.1 已知n 阶实矩阵()n n ij a ⨯=∈R A ，如果存在常数λ和非零向量x ，使λ=Ax x 或 ()λ-=0A I x (9.1.1)那么称λ为A 的特征值(eigenvalue)，x 为A 的相应于λ的特征向量(eigenvector). 多项式111212122212()det()n n n n n nn a a a a a a p a a a λλλλλ-⎡⎤⎢⎥-⎢⎥=-=⎢⎥⎢⎥-⎣⎦A I (9.1.2)称为特征多项式(characteristic polynomial)，det()0λ-=A I (9.1.3)称为特征方程(characteristic equation).注 式(9.1.3)是以λ为未知量的一元n 次代数方程，()det()n p λλ=-A I 是λ的n 次多项式. 显然，A 的特征值就是特征方程(9.1.3)的根. 特征方程(9.1.3)在复数范围内恒有解，其个数为方程的次数(重根按重数计算)，因此n 阶矩阵A 在复数范围内有n 个特征值. 除特殊情况 (如2,3n =或A 为上(下)三角矩阵)外，一般不通过直接求解特征方程(9.1.3)来求A 的特征值, 原因是这样的算法往往不稳定. 在计算上常用的方法是幂法与反幂法和相似变换方法. 本章只介绍求矩阵特征值与特征向量的这两种基本方法. 为此将一些特征值和特征向量的性质列在此处.定理9.1.2 设n 阶方阵()ij n n a ⨯=A 的特征值为12,,,n λλλ，那么(1) 121122n nn a a a λλλ+++=+++;(2) 12det n λλλ=A .定理9.1.3 如果λ是方阵A 的特征值，那么(1) kλ是kA 的特征值，其中k 是正整数;(2) 当A 是非奇异阵时，1λ是1-A 的特征值. (3) ()n p λ是()n p A 的特征值，其中()n p x 是多项式2012()n n n p x a a x a x a x =++++.定义9.1.4 设,A B 都是n 阶方阵. 若有n 阶非奇异阵P ，使得1-=P AP B ，则称矩阵A 与B 相似(similar)，1-P AP 称为对A 进行相似变换(similarity transformation)，P 称为相似变换矩阵(similarity transformation matrix).定理9.1.5 若矩阵A 与B 相似，则A 与B 的特征值相同. 定理9.1.6 如果A 是n 阶正交矩阵，那么 (1) 1T -=A A ，且det 1=A 或1-; (2) 若=y Ax ，则22=yx , 即T T ⋅=⋅x x y y .定理9.1.7 设A 是任意n 阶实对称矩阵，则 (1) A 的特征值都是实数； (2) A 有n 个线性无关的特征向量.定理9.1.8 设A 是任意n 阶实对称矩阵，则必存在n 阶正交矩阵P ，使得1T -==P AP P AP Λ，其中12diag(,,,)n λλλ=Λ是以A 的n 个特征值12,,,n λλλ为对角元素的对角矩阵.定理9.1.9 (圆盘定理) 矩阵()ij n n a ⨯=A 的任意一个特征值至少位于复平面上的几个圆盘1,2,,n ，中的一个圆盘上。

9.2 幂法与反幂法9.2.1 幂法及其加速9.2.1.1 幂法幂法是计算矩阵按模最大特征值(largest eigenvalue in magnitude)及相应特征向量的迭代法. 该方法稍加修改，也可用来确定其他特征值. 幂法的一个很有用的特性是：它不仅可以求特征值，而且可以求相应的特征向量. 实际上，幂法经常用来求通过其他方法确定的特征值的特征向量. 下面探讨幂法的具体过程.设矩阵n n⨯∈RA 的n 个特征值满足1230n λλλλ>≥≥≥≥, (9.2.1)且有相应的n 个线性无关的特征向量12,,,n x x x ，则12,,,n x x x 构成n 维向量空间nR 的一组基, 因此1,,1,2,,nni i i i i n αα=⎧⎫==∈=⎨⎬⎩⎭∑R R z z x .在n R 中选取某个满足10α≠的非零向量01ni i i α==∑z x .用矩阵A 同时左乘上式两边，得011n ni i i i i i i ααλ====∑∑Az Ax x .再用矩阵A 左乘上式两边，得2201ni i i i αλ==∑A z x .这样继续下去，一般地有01,1,2,.nkk i i i i k αλ===∑A z x (9.2.2)记10,1,2,kk k k -===z Az A z ，则由式(9.2.2)得0111121, 1,2,.knn k k ki k i i i i i i i k λαλλααλ==⎡⎤⎛⎫⎢⎥===+= ⎪⎢⎥⎝⎭⎣⎦∑∑z A z x x x (9.2.3)由假设(9.2.1)，结合式(9.2.3)，得111lim,kkk αλ→∞=z x (9.2.4)于是对充分大的k 有111.k k λα≈z x (9.2.5)式(9.2.4)表明随着k 的增大，序列{}1kkλz 越来越接近A 的对应于特征值1λ的特征向量1x 的1α倍, 由此可确定对应于1λ的特征向量1x . 当k 充分大时，可得1x 的近似值.上述收敛速度取决于比值21λ. 事实上，由式(9.2.3)知，321122331111kkkkn nn kλλλααααλλλλ≤++++z x x x x . (9.2.6)再由式(9.2.1)得321111nλλλλλλ>≥≥≥. (9.2.7) 结合式(9.2.6)和式(9.2.7)知，序列{}1k kλz 收敛速度取决于比值21λλ.下面计算1λ. 由式(9.2.3)知1111011121,k n k k ik k i i i A A λλααλ++++=⎡⎤⎛⎫⎢⎥===+ ⎪⎢⎥⎝⎭⎣⎦∑z z z x x当k 充分大时, 11111k k λα++≈z x . 结合式(9.2.5)，得11k k λ+≈z z .这表明两个相邻向量大体上只差一个常数倍，这个倍数就是A 的按模最大特征值1λ. 记()T(1)(2)(),,,n k k k k z z z =z , 则有()11()lim ,1,2,,j k j k kz j n z λ+→∞==， (9.2.8)即两个相邻的迭代向量所有对应分量的比值收敛到1λ.定义9.2.1 上述由已知非零向量0z 及矩阵A 的乘幂k A 构造向量序列{}k z 来计算A 的按模最大特征值1λ及相应特征向量的方法称为幂法(power method)，其收敛速度由比值21γλλ=来确定，γ越小，收敛越快.注 由幂法的迭代过程(9.2.3)容易看出，如果11λ>（或11λ<），那么迭代向量k z 的各个非零的分量将随着k →∞而趋于无穷（或趋于零），这样在计算机上实现时就可能上溢(或下溢). 为了克服这个缺点，需将每步迭代向量进行规范化: 记()T(1)(2)()1,,,n k k k k k y y y -==y Az .若存在k y 的某个分量0()j ky ，满足0()()1max j j kk j ny y ≤≤=，则记0()max()j k k y =y . 将k y 规范化max()k k k =z y y ,这样就把k z 的分量全部控制在[1,1]-中.例如，设T(2,3,0,5,1)k =--y ，因为k y 的所有分量中，绝对值最大的的是5-，所以max()5k =-y ，故Tmax()(0.4,0.6,0,1,0.2)k k k ==--z y y .综上所述，得到下列算法：算法9.2.1 (幂法) 设A 是n 阶实矩阵，取初始向量0n∈R z ，通常取T 0(1,1,,1)=z ，其迭代过程是：对1,2,k =，有1,max(),.k k k k kk k m m -=⎧⎪=⎨⎪=⎩y Az y z y (9.2.9) 定理9.2.1 对式(9.2.9)中的序列{}k z 和{}k m 有()11lim max k k →∞=x z x , 1lim k k m λ→∞=， (9.2.10)其收敛速度由21γλλ=确定.证明 由迭代过程(9.2.9)知211111k k k k k k k k k k k m m m m m m -----====z y A z A y A z001max()kkk k ii m====∏A z A z A z , (9.2.11)其中01ni i i α==∑z x . 若10α≠，则由(9.2.3)知：011121kn k k ii i i λλααλ=⎡⎤⎛⎫⎢⎥=+ ⎪⎢⎥⎝⎭⎣⎦∑A z x x , 代入式(9.2.11)得11211121max kni i ii k kn i i i i λααλλααλ==⎛⎫+ ⎪⎝⎭=⎛⎫⎛⎫⎪+ ⎪ ⎪⎝⎭⎝⎭∑∑x x z x x ， (9.2.12) 故()11lim max k k →∞=x z x . (9.2.13)而()()()21211120max max max k k k k k k k k A A ------=====Ay A z A z y Az y z z 111211111121max kn k ii i i k n k i i i i a a a a λλλλλλ=--=⎡⎤⎛⎫+⎢⎥⎪⎢⎥⎝⎭⎣⎦=⎛⎫⎛⎫ ⎪+ ⎪ ⎪⎝⎭⎝⎭∑∑x x x x , (9.2.14) 于是1121111121max max()max kn i i i i k k k n i i i i a a m a a λλλλλ=-=⎛⎫⎛⎫⎪+ ⎪ ⎪⎝⎭⎝⎭==⎛⎫⎛⎫⎪+ ⎪ ⎪⎝⎭⎝⎭∑∑x x y x x , (9.2.15)故1lim k k m λ→∞=. (9.2.16)由式(9.2.6)和式(9.2.7)知：上述收敛速度由21γλλ=确定. 证毕！例9.2.1 用幂法求方阵123213335⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦A按模最大特征值及相应的特征向量，要求3110k k m m ---<.解 选取初始向量T0(1,1,1)=z ，按式(9.2.9)迭代，结果见表9.2.1.表9.2.1因此，所求按模最大特征值18.3589λ≈，相应特征向量T1(0.559817,0.559817,1)≈x .事实上，A 的按模最大特征值148.3588989λ=≈，相应特征向量T 1(0.55981649,0.55981649,1)=x , 故所得结果具有较高的精度.9.2.1.2 幂法的加速从上面的讨论可知，由幂法求按模最大特征值，可归结为求数列{}k m 的极限值，其收敛速度由21γλλ=确定. 当21γλλ=接近1时, 收敛速度相当缓慢. 为了提高收敛速度,可以利用外推法进行加速.因为序列{}k m 的收敛速度由21γλλ=确定，所以若{}k m 收敛，当k 充分大时，则有()121k k m O λλλ⎡⎤-=⎢⎥⎣⎦，或改写为()121kk m c λλλ-≈,其中c 是与k 无关的常数. 由此可得11211k k m m λλλλ+-≈-, (9.2.17)这表明幂法是线性收敛的. 由式(9.2.17)知1121111k k k k m m m m λλλλ+++--≈--.由上式解出1λ，并记为2k m +，即2221122121()22k k k k k k k k k k k k km m m m m m m m m m m m m ++++++++--==+-+-+， (9.2.18) 这就是计算按模最大特征值的加速公式.将上面的分析归结为如下算法：算法9.2.2 (幂法的加速) 设A 是n 阶实矩阵，给定非零初始向量0n∈R z ，通常取T 0(1,1,,1)=z . 对1,2k =，用迭代式1,max(),k k k k kk k m m-=⎧⎪=⎨⎪=⎩y Az y z y 求出12,m m 及12,z z . 再对3,4,k =，迭代过程为1212212,max(),(),2.k k k k k k k k k k k k k k m m m m m m m m m ------=⎧⎪=⎪⎪-⎨=-⎪-+⎪⎪=⎩y Az y z y (9.2.19) 当1k k m m ε--<(0ε>是预先给定的精度) 时，迭代结束，并计算k z ；否则继续迭代，直至满足迭代停止条件1k k m m ε--<.有关加速收敛技术，读者请参考文献[11].8.2.2 反幂法及其加速反幂法是计算矩阵按模最小特征值及相应特征向量的迭代法, 其基本思想是：设矩阵n n ⨯∈R A 非奇异，用其逆矩阵1-A 代替A , 矩阵A 的按模最小特征值就是矩阵1-A 的按模最大特征值. 这样用1-A 代替A 做幂法，即可求出1-A 的按模最大特征值，也就是矩阵A 的按模最小特征值；这种方法称为反幂法(inverse power method).因为矩阵A 非奇异，所以由i i i λ=Ax x 可知：11i i iλ-=A x x . 这说明：如果A 的特征值满足1210n-n λλλλ≥≥≥>>, (9.2.20)那么1-A 的特征值满足1211111nn-λλλλ>≥≥≥, (9.2.21)且对A 的应于特征值i λ的特征向量i x 也是1-A 的对应于特征值1i λ的特征向量.由上述分析知：对1-A 应用幂法求按模最大的特征值1n λ及相应的特征向量n x ，就是求A 的按模最小的特征值n λ及相应的特征向量n x .算法9.2.3 (反幂法) 任取初始非零向量0n∈R z ，通常取T 0(1,1,,1)=z . 为了避免求1-A ，对1,2,k =，将迭代过程(9.2.9)改写为:1,max(),.k k k k kk k m m -=⎧⎪=⎨⎪=⎩Ay z y z y (9.2.22) 仿定理9.2.1的证明，可得：定理9.2.2 对式(9.2.22)中的序列{}k z 和{}k m 有1lim ,lim max()nk k k k n nm λ→∞→∞==x z x , (9.2.23)其收敛速度由1n n γλλ-=确定.注 按(9.2.22)进行计算，每次迭代都需要解一个方程组1k k -=Ay z . 若利用三角分解法求解方程组，即=A LU ，其中L 是下三角矩阵，U 是上三角矩阵，这样每次迭代只需解两个三角方程组11,.k k --=⎧⎨=⎩L z Uy v ν 9.2.3 原点平移法为了提高收敛速度，下面介绍加速收敛的原点平移法.设矩阵p =-B A I ，其中p 是一个待定的常数，A 与B 除主对角线上的元素外，其他元素相同. 设A 的特征值为12,,,n λλλ，则B 的特征值为12,,,n p p p λλλ---,且A 与B 的特征向量相同.9.2.3.1 原点平移法的幂法设1λ是A 的按模最大的特征值，选择p ，使12,3,4,,i p p p i n λλλ->-≥-=，及2211p p λλλλ-<-. (9.2.24) 对B 应用幂法，得 算法9.2.4 对1,2,k =，有1(),max(),,k k k k kk k p m m -=-⎧⎪=⎨⎪=⎩y A I z y z y (9.2.25) 且()111lim ,lim ,max k k k k m p λ→∞→∞=-=x z x (9.2.26)其收敛速度由21()()p p λλ--确定.由式(9.2.24)知：在计算B 的按模最大特征值1p λ-的过程(9.2.25)中，收敛速度得到加速；算法9.2.4又称为原点平移下的幂法(power method with shift).9.2.3.2 原点平移下的反幂法设n λ是A 的按模最小的特征值，选择p ，使1,1,2,,2n n i p p p i n λλλ--<-≤-=-. (9.2.27)及11n n n n p p λλλλ---<-. (9.2.28)若矩阵p =-B A I 可逆，则1-B 的特征值为12111,,,,n p ppλλλ---且有1111,1,2,,2n n i i n p p pλλλ->≥=----. (9.2.29)对B 应用反幂法，得： 算法9.2.5 对1,2,k =，有1(),max(),,k k k k kk k p m m --=⎧⎪=⎨⎪=⎩A I y z y z y (9.2.30) 且()111lim ,lim ,max k k k k n m p λ→∞→∞==-x z x (9.2.31)其收敛速度由1())n n p p λλ---确定.由式(9.2.28)知：在计算1-B 的按模最大特征值1n pλ-的过程(9.2.30)中，收敛速度得到加速. 算法9.2.5又称为原点平移下的反幂法 (inverse power method with shift).定义9.2.8 原点平移下的幂法与原点平移下的反幂法统称为原点平移法.注 有的资料上的原点平移法专指原点平移下的反幂法；而有的资料上的反幂法指的就是原点平移下的反幂法.原点平移下的反幂法(算法9.2.7)的主要应用是：已知矩阵的近似特征值后，求矩阵的特征向量. 其主要思想是：若已知A 的特征值m λ的近似特征值为m λ，则m A I λ-的特征值就是(1,2,,)i m i m λλ-=，其中(1,2,,)i i m λ=是A 的特征值. 而按模最小的特征值是m m λλ-，相应的特征向量与A 的特征向量相同.利用公式(9.2.30)可求出{}{},k k m z ，并且有1lim ,lim max()mk k k k m mm m λλ→∞→∞==-x z x ，其收敛速度由11(min )m mi m i m i m i mλλλλλλλλ≤≤---=--确定. 于是当k 充分大时，可取m k ≈x z ， 1m m km λλ≈+, 只要近似值m λ适当好，收敛速度就很快，往往迭代几次就能得到满意的结果.例9.2.2 已知特征值λ的近似值0.3589λ=-，用原点平移下的反幂法求方阵123213335⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦A得对应特征值λ的特征向量.解 取0.3589p λ==-，对矩阵1.3589230.35892 1.3589333 5.3589A pI A I ⎡⎤⎢⎥-=+=⎢⎥⎢⎥⎣⎦.迭代公式(9.2.30)中的k y 是通过解方程组1()k k p --=A I y z求得的. 为了节省工作量，可先将p -A I 进行LU 分解.在LU 分解中尽量避免较小的rr u 当除数，通常可以先对矩阵p -A I 的行进行调换后再分解. 为此，我们可用001010100⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦P 乘p -A I 后再进行LU 分解，即 6001 1.35892333 5.3589()0102 1.358932 1.3589310033 5.3589 1.35892310033 5.35890.66671000.64110.5726,0.45301100 3.0710p -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥⎢⎥--⨯⎣⎦⎣⎦P A I LU 1( 1.2679)k k --=P A I y Pz ， 即 1k k -=LUy Pz .令11,.k k k k --==Lv Pz Uy v 选取0z ，使1T10(1,1,1)-==Uy L Pz ，得T 1(290929.45,290927.56,325732.90)=-y ，11max()325732.90m ==-y ， T 111(0.8932,0.8931,1)m ==--z y .由121-=Uy L Pz 得T 2(845418.49,845418.49,946558.42)=--y ，22max()946558.42m ==y ， T 222(0.8932,0.8932,1)m ==--z y .由于1z 与2z 的对应分量几乎相等，故A 的特征值为2110.35890.35889894354117946558.42m λλ≈+=-+=-, 相应的特征向量为T2(0.8932,0.8932,1)==--x z .而矩阵A 的一个特征值为40.358898943540674λ=-=，相应的特征向量为T (0.89315,0.89315,1)--，由此可见得到的结果具有较高的精度.9.3 QR 算法上一节我们介绍了求矩阵特征值的幂法和反幂法. 幂法主要用来求矩阵的按模最大特征值，而反幂法主要用于求特征向量. 本节将介绍幂法的推广和变形——QR 算法，它是求一般中小型矩阵全部特征值的最有效的方法之一，其基本思想就是利用矩阵的QR 分解. 矩阵n n ⨯∈R A 的QR 分解就是：用Householder 变换将矩阵A 分解成正交矩阵Q 与上三角矩阵R 的乘积，即=A QR . 下面首先介绍Householder 变换.9.3.1 豪斯霍尔德变换定义9.3.1 设()ij n n b ⨯=B 是n 阶方阵，若当1i j >+时，0ij b =，则称矩阵B 为上Hessenberg 矩阵(Hessenberg matrix )，又称为准上三角矩阵，它的一般形式为1112121222323,1n n n n n nn b b b bb b b b b b -⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦B . (9.3.1) 下面讨论如何将矩阵A 用正交相似变换化成(9.3.1)的形式. 为此先介绍一个对称正交矩阵——Householder 矩阵.定义9.3.2 设向量n∈R u 的欧氏长度21=u，I 为n 阶单位矩阵，则称n 阶方阵T ()2==-H H u I uu (9.3.2)为Householder 矩阵(Householder matrix). 对任何n∈R x ，称由()=H H u 确定的变换=y Hx (9.3.3)为镜面反射变换(specular reflection transformation)，或Householder (Householder transformation).注 Householder 变换，最初由A.C Aitken 在1932年提出. Alston Scott Householder 在1958年指出了这一变换在数值线性代数上的意义. 这一变换将一个向量变换为由一个超平面反射的镜像，是一种线性变换. 运用线性代数的知识，很容易证明：定理9.3.3 (9.3.2)式定义的矩阵H 是对称正交矩阵；对任何n∈R x ，由线性变换H =y x 得到y 的欧氏长度满足22=y x .反之，有下列结论：定理9.3.4 设,n∈R x y ，≠x y .若22=x y ，则一定存在由单位向量确定的镜面反射矩阵()H u ，使得H =x y .证 令2-=-x yu x y ，显然21=u . 构造单位向量u 确定的镜面反射矩阵T ()2==-H H u I uu ，()T T T T2222()()2()()22⎡⎤----⎢⎥=-=-=---⎢⎥⎣⎦x y x y x y x x y x Hx I uu x I x x x y x y . 又因为22=xy ，即T T =x x y y ，所以————————————————————海森伯格(Karl Adolf Hessenberg 1904年9月8日–1959年2月22日) 是德国数学家和工程师.豪斯霍尔德(Alston Scott Householder 1904年5月5日–1993年7月4日) 是美国数学家，在数学生物学和数值分析等领域卓有建树.2T T T 2T T T T T T T T T T ()()()()2(),-=--=--=--+=--+=-x y x y x y x y x y x x x y y x y y x x x y x y x x x x y x 于是T T 222()()()--=-=--=-x y x x y x Hx x x x y y x y.证毕！由定理9.3.4得： 算法9.3.1 若T 12(,,,)n x x x =x ，其中2,,n x x 不全为零，则由12T 11212T1T22sgn(),,(1,0,,0),(),()2n x x σσρσσρ-⎧=⎪=+=∈⎪⎪=+⎨⎪⎪==-=-⎪⎩R 12其中=x u x e e u uu H H u I I uu u (9.3.4) 确定的镜面反射矩阵H ，使得1σ=Hx e ，其中1,0,sgn()0,0,1,0.a a a a >⎧⎪==⎨⎪-<⎩例9.3.1 设T(1,2,2)=--x ，按式(9.3.4)的方法构造镜面反射矩阵H ，使得T (,0,0)*=Hx (*表示某非零元素).解12sgn()(3x σ==-=-x；T T T 1(1,2,2)(3,0,0)(4,2,2)σ=-=---=--u x e , 其中T 1(1,0,0)=e ；212()3[3(1)]12x ρσσ==+=--+-=12u ， 则所求镜面反射矩阵为1T 1004132210102(4,2,2)2211200122123ρ----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=-=---=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦H I uu ，且1323213232312021220---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦Hx .可以证明：定理9.3.6 对任意n 阶方阵()ij n n a ⨯=A ，存在正交矩阵Q ，使得T =B Q AQ为形如式(9.3.1)的上Hessenberg 阵.证 记(1)(1)(1)1112111121(1)(1)(1)21222212221(1)(1)(1)1212n nn n n n nn n n nn a a a a a a a a a a a a a a a a a a ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦A A ，(1)21(1)311(1)1n a a a ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦x . 由式(9.3.4)可构造1n -阶对称正交矩阵1H ：12121121122T 111111211112121T 1111sgn()sgn(),,(1,0,,0),(),,n i i n a a a a σσρσσρ=--⎧⎛⎫==-⎪ ⎪⎝⎭⎪⎪=+=∈⎨⎪=+⎪⎪=-⎩∑R 12其中=x u x e e u H I u u (9.3.5) 使得1111σ=H x e .记111⎡⎤=⎢⎥⎣⎦I Q H ，且1n n⨯∈R Q ，1I 是1阶单位矩阵. 显然1Q 是对称正交矩阵. 用1Q 对A 作相似变换，得(1)(2)(2)11121(2)(2)12221(2)(2)111112323(2)(2)20nn nn nn a a a a a a a a a σ-⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦记Q AQ Q A Q A . (9.3.6)记(2)(2)(2)T 2232422(,,,)n n a a a -=∈Rx ，同理可构造2n -阶对称正交矩阵2H ，使得2221σ=H x e (T 21(1,0,,0)n -=∈R e ).记222⎡⎤=⎢⎥⎣⎦I Q H ，其中2I 为2阶单位矩阵，则2Q 仍是对称正交矩阵，用2Q 对2A 作相似变换，得(1)(2)(3)(3)1112131(2)(3)(3)122232(3)(3)123332222223(3)(3)434(3)(3)30000nn n n n nn a a a a a a a a a aa a a σσ-⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦记Q A Q Q A Q A . (9.3.7)依此类推，经过k 步对称正交相似变换，得1111111k k k k k k -------=Q A Q Q A Q(1)(2)(3)(1)()()()1112131,111,11(2)(3)(1)()()()122232,122,12(3)(1)()()()2333,133,133(1)()()()1,11,1,11,000000k k k k k k k n k k k k k k k n k k k k k kk nk k k k k k k k k k k na a a a a a a a a a a a a a a a a a a a a a σσσ--+--+--+-----+-=()()()1,1()()()1,1,11,()()(),100000000k k k k k kk k k kn k k k k k k k k n k k k nk n k nna a a a a a a a a σ-++++++⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦记A . (9.3.8)重复上述过程，则有(1)(2)(3)(1)()1112131,11(2)(3)(1)()122232,12(3)(1)()12333,1322222213(1)()1,11,()1n n n n n n n n n n n nn n n n n n n n n n n n nn n nn a a a a a a a a a a a a a a a σσσσ-------------------⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦记Q A QQ A Q A . (9.3.9)由式(9.3.6)-(9.3.9)，可得122223332231132n n n n n n n n n n n n n -------------===A Q A Q Q Q A Q Q Q Q Q AQ Q Q .若记1n -=B A ，122n -=Q Q Q Q ，则Q 为正交矩阵，且有T =B Q AQ . 证毕！注1 由定理9.3.6知：因为任意n 阶方阵A 与n 阶上Hessenberg 矩阵B 相似，所以求A 的阵特征值的问题，就可转化为求上Hessenberg 矩阵B 的特征值的问题.注2 若A 是对称矩阵，则B 也是对称矩阵. 再由B 的形式(9.3.1)知，此时B 一定是对称三对角阵. 于是，求对称矩阵A 的阵特征值的问题，便可转化为求对称三对角阵B 的特征值问题.例9.3.2 设矩阵121222*********1⎡⎤⎢⎥-⎢⎥=⎢⎥-⎢⎥⎣⎦A 试用镜面反射变换求正交矩阵Q ，使TQ AQ 为上Hessenberg 矩阵.解 第1步 记1=A A ，(1)(1)(1)T T 1213141(,,)(2,1,2)a a a ==x ，利用式(9.3.4)构造三阶镜面反射阵1H ：112sgn(2)3σ===x ；T T T 1111(2,1,2)(3,0,0)(5,1,2)σ=+=+=u x e , 其中T 1(1,0,0)=e ；2(1)1111212()3(32)15a ρσσ==+=+=12u ； 则所求镜面反射矩阵为1T 111110050.66670.33330.666710101(5,1,2)0.33330.93330.1333,1500120.66670.13330.7333ρ----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=-=-=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦H I u u111100000.66670.33330.6667,00.33330.93330.133300.66670.13330.7333⎡⎤⎢⎥---⎡⎤⎢⎥==⎢⎥⎢⎥--⎣⎦⎢⎥--⎣⎦I Q H T 21111113003 2.33330.46670.066700.4667 1.5733 1.346700.0667 1.34670.0933-⎡⎤⎢⎥--⎢⎥===⎢⎥⎢⎥-⎣⎦A Q A Q Q AQ . 第2步 记(2)(2)T T23242(,)(0.4667,0.0667)a a ==-x , 利用式(9.3.4)构造2阶镜面反射阵2H ：222sgn(0.4667)0.4714σ===x ；T T T 2221(0.4667,0.0667)(0.4714,0)(0.9381,0.0667)σ=+=-+=-u x e , 其中T1(1,0)=e ；2(2)2222322()0.4714(0.47140.4667)0.4422a ρσσ==+=+=12u ； 则所求镜面反射矩阵为1T2222100.93810.99010.14151(0.9381,0.0667),010.06670.14150.98990.4422ρ--⎡⎤⎡⎤⎡⎤=-=--=⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦H I u u 222100100,000.99010.1415000.14150.9899⎡⎤⎢⎥⎡⎤⎢⎥==⎢⎥⎢⎥-⎣⎦⎢⎥⎣⎦I Q H T 322222213003 2.33330.4714000.4714 1.5733 1.500000 1.50000.5000-⎡⎤⎢⎥--⎢⎥===⎢⎥--⎢⎥-⎣⎦A Q A Q Q A Q . 由此得正交矩阵12100000.66670.23570.707100.33330.94290.000100.66670.23570.7070⎡⎤⎢⎥--⎢⎥==⎢⎥--⎢⎥-⎣⎦Q Q Q ， 使得T 313003 2.33330.4714000.4714 1.1667 1.500000 1.50000.5000-⎡⎤⎢⎥--⎢⎥==⎢⎥--⎢⎥-⎣⎦Q AQ A 为上Hessenberg 矩阵。