1、 矩阵计算矩阵计算-特征值特征值的数值解法的数值解法 2013113123 罗来星 特征值的数值解法 2013113123 罗来星 摘要摘要：对于,H RAR 矩阵特征值问题就是求C 及非零向量 x，使Axx,这样的 称为矩阵 A 的特征值，向量 x 称为对应特征值的特征向量。通过圆盘定理：在计算矩阵的特征值时，能够给出特征值大小的一个范围，用范数粗略估计特征值。（1）可以用乘幂法求矩阵的特征向量来求出特征值，它主要用来求矩阵按模最大的特征值及其对应的特征向量。其优点是算法简单，容易计算机实现，缺点是收敛速度慢，其有效性依赖于矩阵特征值的分布情况。（2）近代发展起来的 QR 方法是用来计算任意
2、矩阵的全部特征值和特征向量的最有效的方法之，其中 Householder 变换和 Givens 变换是 QR 方法的核心内容。（3）实际问题中，遇到较多的是实对称矩阵，它的特征值问题的计算比一般情况简单些。对于实对称矩阵我们一般用 Jacobi 方法来求解。用 Jacobi 方法求得的结果精度一般都比较高，特别是求得的特征向量正交性很好。所以 Jacobi 方法是求实对称矩阵全部特征值和特征向量的一个较好的方法。它的弱点是计算量大，对原矩阵是稀疏矩阵，旋转变换后不能保持其稀疏的性质，所以 Jacobi 方法只能应用于阶数不高的“满矩阵”。关键词关键词：特征值；圆盘定理；QR 算法；Jacobi
3、 方法。引言引言：特征值的概念是在 1855 年由凯莱提出的，它先于矩阵的提出，但直到 20 世纪中期特征值的名称才得到统一。这里将介绍几种常用的计算矩阵特征值和特征向量的数值方法，如相似变换法，包括 Householder 变换、Givens 变换、QR 方法、对称矩阵的 Jacobi 方法等。幂法的实际使用是在 20 世纪早期开始的，但在此之前它已被反复提到过，如 1931 年 Krylov利用幂法产生的序列求出了矩阵的特征多项式。反幂法是 1944 年由 Wielandt 提出的。而用Jacobi 方法计算对称矩阵特征值要追溯到 1845 年。QR 迭代方法是在 1961 年由 Fran
4、cis 和Kublanovskaya 分别独立提出的。在纯数学角度来看，特征值和特征向量是将二次型化为标准型的一个重要的工具，在科学研究和工程技术中的许多问题，如信息系统设计、非线性最优化、数量经济分析、微小振动研究等等，都涉及矩阵的特征值与特征向量的计算。对于,H RAR 矩阵特征值问题就是求C 及非零向量 x，使得Axx,这样的 称为矩阵 A 的特征值，向量 x 称为对应特征值的特征向量。在工程技术中许多实际问题最终都要归结为矩阵特征值和特征向量的计算问题。在线性代数中个，我们知道，矩阵 A 的特征值的求解可通过矩阵 A 的特征多项式 1112121221212detnnnnaaLaaaL
5、afIAMMMaaLa 的根来得到。上式中 I 是 n 阶单位阵，detIA 表示IA 的行列式，它是的n 次代数多项式。众所周知，5 次以上的代数多项式是没有求根公式的，它们的求解大部分是相当困难的。同时高次多项式的重根的计算往往精度较低等。因此，从数值计算的观点来看，用特征多项式来求矩阵特征值的方法并不可取，必须建立有效的数值方法。在实际应用中，求矩阵的特征值和特征向量通常分为变换法和迭代法两种。变换法是从原矩阵出发，用有限个正交相似变换将其化为便于求出特征值的形式，如对角矩阵、三角形矩阵等。这类方法有工作量小和应用范围广的优点，但由于舍入误差的影响，其精度往往不高。迭代法的基本思想是将特
6、征值及特征向量作为一个无限序列的极限来求得。这类方法对舍入误差的影响有较强的稳定性，但通常工作量较大。实际问题中提出的特征值问题的要求是不相同的，有些问题只要求计算绝对值最大或最小的特征值，但更多的则要求计算全部特征值和特征向量。必须针对问题特点进行具体分析，选择适当的方法。下面介绍几种目前在计算机上比较常用的矩阵特征值问题的数值方法。本论本论：1.矩阵特征值问题的有关理论 本节叙述一些与特征值有关的概念与结论。命题（1）设,(1,2,)H RijiAaRiL nA是 的特征值则有 1;niiDet A 11,nniuiitr Aa 式中 tr（A）为矩阵 A 的迹。在计算矩阵的特征值时，如果
7、能够给出特征值大小的一个范围，在很多情况下是非常有用的。在前面的范数讲到，可以用范数粗略估计特征值外，还可以利用著名的Gerschgorin 圆盘定理来估计特征值。定理 1（Gerschgorin 圆盘定理）设,H RijAaC 则设 A 的每一个特征值 1niiD 其中 1,1,2,niijijjj iDz zaaiL n 为第 i 个圆盘。定理 2 如果 A 的 n 个圆盘中的 m 个圆盘形成一个连通集 S，且 S 与其他 n-m 个圆盘是分离的。则在连通集 S 中恰好有 A 的 m 个特征值。特别的，当 m=1 时，即每个孤立圆盘中恰好有一个特征值。定义 1 设 A 为 n 阶实对称矩阵
8、，对任意非零向量 x，称 ,Ax xR xx x 为对应于向量 x 的 Rayleigh 商，其中1,niiix yx y 为向量 x 和 y 的内积。定理 3 设 A 为 n 阶实对称矩阵，其特征值都为实数，排列为 12,nL 对应的特征向量 x1,x2,L,xn 组成正交向量组，则有 （1）对任何非零向量 n1,nxRR x有 （2）11max()()R xR x （3）min()()nnR xR x 证明 设 x=0，则有表达式：211,0nniiiiixxx x 2221111(,)nmmniiiiiiiAx x 由此可见（1）成立。在 Rayleigh 商中分别取 x=x1和 x=x
9、n，可得到 Rayleigh 商的最大值和最小值，即（2）和（3）成立。2.乘幂法 （1）乘幂法和加速方法 乘幂法是通过求矩阵的特征向量来求出特征值的一种迭代法。它主要用来求矩阵按模最大的特征值及其对应的特征向量。其优点是算法简单，容易计算机实现，缺点是收敛速度慢，其有效性依赖于矩阵特征值的分布情况。设矩阵H RAR 的 n 个特征值满足 12nL，（1）对应的 n 个特征向量 x1，x2，K,xn 线性无关。称模最大的特征值1 为主特征值，对应的特征向量 x1 为主特征向量。乘幂法用于求主特征值和特征向量。它的基本思想是任取非零的初始向量 v0，由矩阵A 构造一个向量序列：11,2,kkko
10、vAvA vnL 由假设 v0 可表示为：01 122nnvxxLx （2）若记kiv 为向量 vk 的第 i 个分量，则有 01 111 1121knnkkkkikiiiiiikiivA vaxa xaxa x 111 111 1,kkiikkiivxvx 其中12knjkiiix。若10,0,ix 则由lim0kk 知 11 111lim,limkkkkkkvvxv 可见，当 k 充分大时，vk近似于主特征向量（相差一个常数倍），vk+1与 vk的对应非零分量的比值近似于主特征值。在实际计算中，需要对计算结果进行规范化。因为当11 时，vk趋于零；当11 时，vk的非零分量趋于无穷，从而计
11、算时会出现下溢或上溢。为此，对向量12,nnZz zL zR 记 max（Z）=zi，其中,izZ 这样，我们有如下乘幂法的实用的计算公式：任取000v，对于 k=1,2，L 分别计算 1,/max()kkkkkvAuuvv （3）求出对应矩阵的主特征向量和特征值的近似值，有下面的定理。定理 4.设H RAR 的特征值(1,2,)iiL n 满足上式并且对应的 n 个线性无关的特征向量1,2,ix iL n。给定初始向量01niiiva x，0i，则由上式生成的向量序列有 111lim,limmaxmaxkkkxvx 证明 若记 mk=max（vk），有 2120111kkkkkkkkkAuA
12、 uA uuLmm mm mLm 由于 uk的最大分量为 1，即 max（uk）=1,故 110maxkkkm mLmA u 又注意到假定00uv，从而 00100,maxmaxkkkkkkA vA vvuAvA v 可以推出 2011 111 12111 1011 101 111 11,maxmaxmaxmaxknkkkikikkkkkkkkkA vxxxxA vuxA vxxkxx 同理，可得到 1 111 111 1111 11maxmaxkkkkkkkxxvxx 1 1111 11maxmaxmaxkkkxvkx 定理得证。由定理的证明可见，乘幂法的收敛速度由21/的大小决定。若 A
13、的特征值不满足（1）式，将有不同的情况。如果12r，且1rr。可以做类似的分析，对初始向量（2）和计算公式（3）有 111lim,limmaxmaxriiikkrkkiiixuvx 可见，uk 仍收敛于一个主特征向量。对特征值的其他情况，讨论较为复杂，完整的乘幂法程序要加上各种情况的判断。-3.QR 算法 近代发展起来的 QR 方法是用来计算任意矩阵的全部特征值和特征向量的最有效的方法之一，对它的研究和应用日益扩大。这里仅对这一算法的基本思想和实现算法作一个简要的介绍。（1）Householder 变换和 Givens 变换 1.Householder 变换 定义 2 设向量2,1nR,则称
14、2THI (1)为 Householder 矩阵或 Householder 变换。有时也称为镜面反射矩阵。例如：1 2 2,3 3 3r，则2=1，于是 7-4-41=-2=-41-89-4-81TH I 为 Householder 矩阵。Householder 矩阵 HH 有如下性质：（1）H 是对称正交阵，即 H=HT=H-1.。事实上，显然有 HT=H,又由21T 得知 244TTTTH HHII (2)对任何22,y=xnxR记y=Hx,有 。（3）记 S 为与 垂直的超平面，则几何上 x 与 y=Hx 关于平面 S 对称。事实上，由2TyHxIx 得知 2Txyx 上式表明向量 x-
15、y 与 平行，注意到 y 与 x 的长度相等。它的几何解释是向量 x 经变换后的象 y=Hx 是 x 关于 S 对称的向量。对应于性质(2),有下面的定理。定理 7 设 ,nx yRxy 且22xy，则有 Householder 矩阵 H，使得 Hx=y。证明 令2,2TxyHIxy。则有21。由TTx xy y 可知 2222TTTTTxyxx xx yy yxyxyxy 从而可得 2222TTxyxyxHxxxxxxyyxy 定理得证。该定理的一个重要应用是对12,0Tnxx x L x 有 Householder 矩阵 H，使得 1Hxe (2)其中112,1,0,0Tsign xxeL
16、 。矩阵 H 的计算公式为 111TuxexHI （3）关于符号 的选取，是为了使分母的 尽量大，从而有利于数值计算的稳定性。(2)的意义是对向量做消元运算，即使向量的第二个分量到第 n 个分量全都化为零。当然，也可以构造相应的 Householder 矩阵，使向量的连续若干个分量为零。例子 1：对于向量 x=(3,5,1)T,构造 Householder 矩阵 H，使 121,0,0THxsign xx 解 1226,6xsign xx ，19,5,1,1Tuxe。而从2108,54u，按（3）得 12745994529551955315495153THI 2.Givess 变换 House
17、holder 变换可视一个向量的连续若干个分量化为零，有时希望把指定的一个分量化为零，此时可以使用 Givess 变换。定义 3 将 n 阶单位矩阵 In改变第 i，j 行和第 i，j 列的四个元素得到的矩阵 10cossin,0sincos01iJJ i jj 称为 Givens 矩阵或 Givens 变换。容易验证，J（i，j，）为正交矩阵，表示在 n 维空间中将互相正交的两个坐标轴在其所决定的平面上旋转一个角度，并保持正交坐标系的其它轴不动。所以 J（i，j，）也称为平面旋转矩阵，为旋转角。设nxR，对 x 进行 Givens 变换得到的向量记为 y，即 y=J（i，j，）x，则 y可以
18、看做为对向量 x 旋转角后所得的向量。若令 1212,TTnnxx x L xyy y L y 则向量 x 和 y 的分量关系为 cossinsincos,iijjijkkyxxyxxyx ki j 可以看出向量 x 和 y 仅有第 i，j 两个分量不同，且不难验证2222xy。从上式也可以看出，适当的选择角度，便可以使向量 x 的第 j 个分量 xj=0。为此，只需令：2222cos,sin,jiijijxxxxxx 若220ijxx 设H RAR,将矩阵 A 左乘以矩阵 J（i，j，），则矩阵 A 仅有 i，j 两行元素发生改变，其余元素不变。事实上，设,pqpqAaBJ i j k Ab
19、，便有：,cossin,1sincospqpqiqiqjqjqiqjqbapi jbaaqnbaa 类似地，,TAJ i j 只改变 A 的第 i 列和第 j 列的元素，,TJ i jAJ i j 只改变矩阵 A 的第 i 行，第 j 行，第 i 列和第 j 列的元素。设H RijAaR 为对称矩阵，J（i，j，）为一平面旋转矩阵，则 ,TijBJ i jAJ i jb 的元素的计算公式为 2222cossin2sincossincos2sincos1sin2cos22cossin,cossin,iiiijjijjjiijjijijjijjiiijikkiikjkjkkjjkikimmiimb
20、aaabaaabbaaabbaaki jbbaaki jbbali j mi j 而且，不难验证 22222222.iijjijiijjijbbbaaa 4.Jacobi 方法 实际问题中，遇到较多的是实对称矩阵，它的特征值问题的计算比一般情况简单些。我们知道，若矩阵H RAR 为对称矩阵，则存在正交阵 p，使 12,TnPAPdiagD D 的对角元1,2,iiL n 就是 A 的特征值，Tp的列向量就是对应于特征值的特征向量。于是求实对称矩阵 A 的特征值问题转为寻找正交矩阵 P，使TPAPD 为对角阵，而这个问题的关键是如何构造 P。Jacobi 方法是用来计算实对称矩阵的全部特征值及其
21、对应的特征向量的一种变换方法。Jacobi 方法基本思想是对矩阵作一系列正交相似变换，使其非对角线元素收敛到零。所用的正交相似变换即为 Givens 变换。定理 3 设H RAR 为对称矩阵，若 B=PAPT,P 为正交阵，则有FFBA。证明 设(1,2,)TiL n 为 A 的特征值，则 2222111nnnTijiFijiAatr A Atr A 。另一方面，矩阵 B 的特征值也为1,2,iiL n，2221niFiBtr B。因此，22FFBA，定理得证。设 A 的非对角线元素0a ,我们可选择矩阵 J=J(I,J,),使 B=JAJT 的非对角线元素 bij=bji=0。为此，由矩阵
22、B 元素的计算公式可知，可选择，使 2tan2,4ijiijjaaa （2）如果用 D(A)表示 A 的对角线元素的平方和，用 S(A)表示 A 的非对角线元素的平方和，这对 B=JAJT ,由（1）式和（2）式 可知 222,2ijijD BD AaS BS Aa 这说明 B 的对角线元素的平方和比 A 的对角线元素的平方和增加了22ija，而 B 的非对角线元素的平方和减少了22ija。这就是 Jacobi 方法求矩阵特征值和特征向量的出发点。下面说明 Jacobi 方法的计算过程。先在 0OijAAa 中选择非对角元中绝对值最大的 0ija。可设 0ija 0，否则已经对角化了。由（2）
23、式选择平面选择矩阵 J1，使111TijJ A JA 的元素 10ija。计算出A1，再类似地选择 J2，计算2212TAJ A J,继续这个过程，连续对 A 进行一系列平面旋转变换，消除非对角线绝对值最大的元素，直到将所有的非对角线元素全化为充分小为止。定理 4.设 H RAR 为对称矩阵，对 A 施行上述一系列平面旋转变换11,2,TmmmmAJ AJmL,则有lim0mmS A。证：设 max0,1,2,mmijiki kaamL,由于 21222,1mmmijmmmikiji kS AS AaS Aan na 则有 1211mmS AS An n 反 复 利 用 上 式，即 得1102
24、121nmS AS Ann n。故lim0mmS A，定理得证。设 m 充分大时，有 2112TTTmmmAJ LJ J AJ J LJD D 为对角阵，则mA 的对角线元素就是 A 的近似特征值，12TTTmMQJ J LJ 的列向量就是对应的近似特征向量。可用mS A 控制迭代终止，其中是所要求的精度。例 用 Jacobi 方法计算矩阵 210121012A 的 特征值和特征向量。解：先取（i，j）=（1,2），1cot20,sincos2，所以 11111110102221110,0322200111222TJAJ AJ 再取（i，j）=（1,3），有1cot2,cos0.88808,s
25、in0.459702 20.8880800.459700100.4597000.88808J 22120.633980.3250500.3250530.6279700.627972.36603TAJ AJ 这里我们看到经变换后所得的非对角线元素的最大绝对值逐次变小。继续做 下去可以得 90.585780.000000.000000.000002.000000.000000.000000.000003.41421A 91290.500000.707100.500000.707100.000000.707100.500000.707100.50000TTTQJ J LJ 矩阵 A 的近似特征值和特
26、征向量均已求出。用 Jacobi 方法求得的结果精度一般都比较高，特别是求得的特征向量正交性很好。所以 Jacobi 方法是求实对称矩阵全部特征值和特征向量的一个较好的方法。它的弱点是计算量大，对原矩阵是稀疏矩阵，旋转变换后不能保持其稀疏的性质，所以 Jacobi 方法只能应用于阶数不高的“满矩阵”。结论结论：矩阵特征值、特征向量理论有着悠久的发展历史和极其丰富的内容。作为一种基本的工具，矩阵特征值、特征向量在数学学科与其它科学技术领域都有广泛的应用。现代科学技术的发展为矩阵特征值、特征向量的应用开辟了更广阔的前景。从发展现状来看，矩阵计算的应用非常广泛，所以我们很有必要熟练的掌握关于矩阵计算的知识，要会用，更要深入理解，在社会发展过程中了解矩阵的实际应用是必不可少的，它能使我们有更好的发展前景。参考文献参考文献：1陈公宁.矩阵理论及其应用M.北京：科学出版社.2007,4.77-80.2同济大学数学系.线性代数M.上海：上海科学技术出版社,1999.74-76.3郭华，刘小明.特征值与特征向量在矩阵运算中的作用J.渝州大学学报.2000.6,17(2):72-75.4李迪.中国数学史简编 M.沈阳：辽宁人民出版社.1984,9.124-125.
