西安科技大学 研究生 数值分析课件7章矩阵特征值与特征向量的计算

1、7 矩阵特征值与特征向量的计算设为n阶方阵，所谓的特征值问题是求数和非零向量，使成立。数称作的一个特征值，非零向量称作与特征值对应的特征向量。求给定方阵的特征值与特征向量是先求解特征方程然后对应于每一个特征值，再求解退化的齐次线性方程组从而得到的特征值及对应的特征向量。但是这种方法计算机很大，计算过程复杂，因此有必要研究相对简单的数值解法。本章主要介绍三类计算特征值的方法：计算大型（稀疏）矩阵主特征的幂法与反幂法，计算中小型（实对称）矩阵全部特征值的Jacobi法，计算中小型矩阵全部特征值的QR法。7.1 特征值估计在矩阵特征值计算中，有时需要对特征值所在范围给出一个估计。这里介绍一种从矩阵的

2、元素出发，运用较简便的运算估计特征值的方法。定义7-1 设，称由不等式在复平面上确定的区域为矩阵的第个盖尔圆（Gerschgorin圆），并用表示。其中称为盖尔圆的半径。定理7-1 矩阵的一切特征值均落在它的个盖尔圆的并集中，即。证明 设是的任一特征值，是对应的特征向量。令，则。由，可得。即 于是有 这表明任一特征值，从而也在的第个盖尔圆的并集中。例7-1 估计矩阵的特征值范围。解 的4个盖尔圆为： 画在复平面上其区域如图7-1所示。图7-1 例7-1盖尔圆分布图于是的全部特征值就在这4个盖尔圆的并集中。为了更确切地知道某个特征值落在哪个或哪几个盖尔圆的并集中，给出如下第二盖尔圆盘定理。定理7

3、-2 若的个盖尔圆中，有个盖尔圆构成的一个连通域（所谓连通域，是指其中的任意两点都可以用位于该区域内的一条折线连接起来），且该连通域与其余个盖尔圆严格分离，则在该连通域中恰好有的个特征值（重特征值按重数重复计算）。特别地，每个孤立的盖尔圆恰有的一个特征值（证明从略）。由定理2可知，在例1中与中各有的一个特征值，而与构成的连通部分中有两个特征值，但不能确定这两个特征值具体落在哪个盖尔圆中。例7-2 估计矩阵的特征值范围。解 的两个盖尔圆为：，在复平面上的区域如图7-2所示。图7-2 例7-2盖尔圆分布图此时只能判断的两个特征值落在与的并集中，至于是每个盖尔圆中各有一个特征值还是两个特征值都落在其

4、中一个盖尔圆上则无法确定。实际上，由于，所以两个特征值都不会在盖尔圆中，而是落在盖尔圆中。对于某些矩阵，可利用相似变换矩阵具有相同特征值的性质得到更确切的特征值范围。设，取正数构成对角阵，对作相似变换，令，由于相似于，所以与的特征值完全相同，又由于与的主对角线元素对应相等，所以与的盖尔圆圆心相同。这表明，若适当选取正数，可以改变盖尔圆的半径，从而有可能将相交的盖尔圆分离得到仅含一个特征值的孤立盖尔圆。选取的一般方法是：欲使的第个盖尔圆的半径大而其余盖尔圆变小，就取，其余。例7-3 求矩阵的特征值范围。解 的3个盖尔圆为：，其中与相交，而孤立。记中所含的一个特征值为，如图7-3所示。为分离与，可

5、以让的第3行元素绝对值变大，第3列元素绝对值变小。现取，则图7-3 例3盖尔圆分布图 图7-4 例7-3分离后盖尔圆分布图其3个盖尔圆分别是：，显然，的盖尔圆是3个孤立的盖尔圆，如图7-4，注意，此情况下，的半径变大了。例7-4 设矩阵按行严格对角占优，则可逆。证明 由线性代数知，可逆的充分条件是，而（其中是的特征值），所以只要证明即可。设是的任一特征值，则必存在某个盖尔圆使。若，则有，而这与按行严格对角占优矛盾，故应有，由的任意性，得。7.2 幂法与反幂法在线性代数中，设是阶方阵，若存在个线性无关的特征向量，则称这个特征向量构成的一个完全的特征向量组。例如，对矩阵，通过求解特征方程，不难求出

6、的三个特征值为，的三个特征值为。方阵可以找到三个线性无关的特征向量，而方阵找不到三个线性无关的特征向量。我们称方阵可对角化，而不可对角化。7.2.1 幂法幂法的基本思想是构造一个向量序列使之逼近主特征值（矩阵的按模最大的特征值）对应的特征向量，然后求出主特征值。该方法简单易行，但收敛速度较慢。现设有一个完全的特征向量组，其对应的特征值是。已知的主特征值是单根，即特征值满足条件任取一个非零初始向量，由矩阵构造向量序列由于的完全特征向量组可以作为向量空间的一组基，因此可由线性表示，即有 （设）于是 其中。注意到，故当时，因此有由于是主特征值对应的特征向量，其乘上常数因子仍为的特征向量，故当充分大时

7、，迭代向量是的特征向量的近似向量。为了利用迭代向量求出主特征值的近似值，设表示的第个分量，则 于是 这表明两相邻迭代向量对应分量的比值收敛于主特征值，亦即当充分大时，可用两相邻迭代向量的分量比作为主特征值的近似值，即若主特征值是的重实特征值，即，对应的个线性无关特征向量为。则有当充分大时，即仍为主特征值对应的特征向量的近似向量，相邻两迭代向量的分量比仍为主特征值的近似值。综上所述，有定理7-3 设是阶实矩阵，具有完全的特征向量组，主特征值是重根，即则对任意非零初始向量，迭代向量满足 ，或 ，这样用非零初始向量及矩阵构造向量序列以计算的主特征值及相应的特征向量的方法称为幂法。不过从上面的讨论中可

8、以看到，如果或，迭代向量当时，其不为零的分量就会趋于无穷大或趋于零。为克服这个缺点，可以在每步迭代中加上对向量规范化的步骤，使迭代向量的数量级保持在一个稳定的量级上，归纳起来，幂法的计算步骤为：步骤1 给定非零初始向量，精度，令；令 ，； 步骤2 迭代，其中表示绝对值最大的分量；步骤3 规范化；步骤4 若且，则即为的近似特征向量，即为的近似值；否则，转步骤2继续迭代。例7-5 用幂法计算 的主特征值和相应的特征向量，结果见表7-1。表7-1K（规范化向量）0(1,1,1)11(0.9091,0.8182,1)2.75000005(0.7651,0.6674,1)2.55879110(0.749

9、4,0.6508,1)2.55879115(0.7483,0.6497,1)2.536625616(0.7483,0.6497,1)2.536584017(0.7482,0.6497,1)2.536559818(0.7482,0.6497,1)2.536559819(0.7482,0.6497,1)2.536537420(0.7482,0.6497,1)2.5365323而此题的准确值为 7.2.2 幂法的加速幂法的收敛速度由比值来确定，越小收敛越快，而当时收敛可能很慢。为了克服这一缺点，常采用原点平移法对幂法进行加速。 设，其中是待定参数。显然，若的特征值为，则的相应特征值为，且、的特征向量

10、相同。这是因为对有特征方程，而对有特征方程，所以另一方面，若是的对应的特征向量，即则 原点平移法的思想是引入矩阵，恰当的选择参数，使是的主特征值，且其速比，这样用幂法求的主特征值的收敛速度就快于用幂法求的主特征值，而一旦求出，由可得的主特征值，达到了加速的目的。但是为了选取恰当的选择参数，需要对的特征值的分布的大致了解。 例7-6 设4阶方阵有特征值其速比。作变换则的特征值为，其速比。设的实特征值满足 若的值可大致估计出，若要求，考察待定参数的选取。 在原点平移法通过变换后，不论如何选取，矩阵的主特征值也只能是在或 。若希望求，就应选择，使称为的主特征值，即 这时的收敛速比是比值和中的较大者，

11、即 显然依赖于的选取，记做。为了使应用幂法求的主特征值的收敛速度尽可能快，我们希望选择最佳参数，使由的表示式（求二者之间的较大值）和对的最小化要求，只有当时，达到最小。由于会有得到矛盾的结果（），所以只能是即 类似地，若用反幂法求最小特征值，若，可大致估计，取最佳平移参数例7-7 取，用原点平移法，计算例7-7中矩阵的主特征值。 解 作变换，则 对应用幂法，计算结果见表7-2。即，则的主特征值为与例7-5比较，上述结果比例7-5迭代15次还好。表7-2k（规范化向量）0(1,1,1)15(0.7516,0.6522,1)1.79140116(0.7491,0.6511,1)1.78884437

12、(0.7488,0.6501,1)1.78733008(0.7484,0.6499,1)1.78691529(0.7483,0.6497,1)1.786658710(0.7482,0.6497,1)1.78659147.2.3 反幂法设方阵按模最小的特征值是，且，则可逆。于是，由，可得，这表明是的主特征值。反幂法就是将幂法应用于，通过求出的主特征值得到的按模最小的特征值及其对应的特征向量。 定理7-4 设是阶实矩阵，具有完全的特征向量组，其特征值满足则对任意非零初始向量，按下述方式构造的迭代向量 ，满足 ， ，在实际计算中，可先对进行LU分解，通过求解 ，来求解方程组。反幂法的计算步骤为：步骤

13、1 预先取定非零向量；给定精度；取；步骤2 对矩阵作分解，；令；步骤3 求解方程组 ，得到迭代向量；步骤4 规范化步骤5 若且，则即为的对应于的近似特征向量，即为的近似值；否则，令，转步骤3继续迭代。7.3 矩阵的两种正交变换 本节先介绍镜面(初等)反射变换和平面旋转变换，它们是QR算法和Jacobi算法的基础。7.3.1 豪斯荷尔德（House holder）变换定义7-2 设有方阵，若当时，则称是上Hessenberg矩阵，即定义7-3 设向量满足，矩阵 （是列向量）称为初等反射矩阵，又称House holder矩阵，记为，即其中是的分量。可以证明初等反射阵是对称阵、正交阵。例7-8 设向

14、量，试证矩阵是一个初等反射阵。证明 令，则由定义7-3，是初等反射阵。定理7-5 设是两个不相等的维列向量，且，则存在一个初等反射阵H，使得证明 令，由例7-8可知是一个初等反射阵。由于注意到，即，又 ，故从而。例7-9 设，用Householder变换将化为与同方向的向量。解 因为，可设，则 取 ，构造Householder矩阵则 推论 设向量，且，则存在初等反射阵 使 。其中，。设，则引入初等反射阵的目的，是设法用一系列初等反射阵将原始矩阵约化成上Hessenberg阵。由于约化过程是逐列进行的，我们先给出计算的算法步骤，该算法算出及，使，的分量冲掉的分量。（1）；（2），此步规范化是为避

15、免计算时产生溢出；（3） ；（4）；（5） ；（6） ； 于是初等反射阵，。如果要将作用于矩阵，设是的第列向量，则，其中，。下面讨论用初等反射阵约化原始矩阵称为上Hessenberg阵的步骤。步骤1 不妨设（否则这一步不需约化），选择初等反射阵，使，其中： 令 则其中，是阵，是维列向量，是阶方阵。步骤 设对已进行了步约化，即其中，是阵，是维列向量，是阶方阵。设，选初等反射阵，使，其中是维单位坐标向量，可得令 则 可见的左上角阶子阵为上Hessenberg阵，从而约化又进了一步。重复此过程，直到使原始矩阵在一系列初等反射阵的作用下，约化为上Hessenberg阵。综上所述，有定理7-6。定理7-

16、6 如果是阶实矩阵，则存在初等反射阵，使（上Hessenberg阵）例7-10 试证矩阵与其约化成为的上Hessenberg阵有相同的特征值。证明 记，由于初等反射阵是正交对称阵，故，且是正交阵，故。于是 其中，。这表明与具有相同的特征多项式，即两者有相同的特征值。 进一步，设是的对应于特征值的特征向量，即，则有 这表明为的对应于的特征向量，于是求原始矩阵的特征值与特征向量可转化为求上Hessenberg阵的特征值和特征向量。 定理7-7 若是实对称矩阵，则存在初等反射阵使 证明 由定理7-6，存在初等反射阵可使约化为上Hessenberg阵，当是对称矩阵时，亦为对称阵，即，且亦为上Hesse

17、nberg阵，故是对称三对角阵。例7-11 用豪斯荷尔德方法将下述矩阵化为上Hessenberg阵。解 （1）对，确定变换阵，其中为初等反射阵，使 （2）计算。记，于是其中，（3）计算及，即求 其中，（4）计算。为上Hessenberg阵。7.3.2 平面旋转变换定义7-4 称矩阵为平面旋转矩阵，又称Givens矩阵，其中，。平面旋转阵是一个正交阵，与单位阵只有在四个位置上的元素不一样，用其左乘矩阵只改变的第行和第行元素。设则平面旋转变换的结果为 若令， 则平面旋转变换向量的第i个分两为，第j个分量为0，其余分量即为对应的分量。和初等反射变换一样，用平面旋转变换也可以将一个方阵化为上Hesse

18、nberg矩阵，也可以将将一个方阵化为上三角矩阵。7.4 QR算法7.4.1 矩阵的QR分解定理7-8 设是可逆矩阵，则存在正交矩阵使且的主对角元素。证明 若，则的第一列不需约化。若有某个，则可选择使的第一列中第个元素变为零。一般地，可设平面旋转矩阵，使记，则。同理，若，可选取使记，则。重复上述过程，可得一系列正交阵使定理7-9 （矩阵的分解）如果阶实矩阵可逆，则可分解为一正交阵和上三角阵的乘积，即，且当的对角元素都为正数时分解唯一。 证明 由定理8知存在正交阵使为上三角阵，记，于是由于是正交阵，则亦为正交阵，故。若有两种分解，记为其中为上三角阵且主对角元素都为正数，为正交阵，于是注意是上三角

19、阵的乘积，结果仍为上三角阵，而是正交阵，所以也应是正交阵。若记，由其上三角性应是下三角阵，应是上三角阵；由其正交性由，故只能是对角阵，且有。又因的主对角元素都为正数，即有故，则，于是，。例7-12 求矩阵的分解。解 方法1：利用初等反射阵进行QR分解令，取，则，再令，取，则 令于是故方法2：利用平面旋转阵进行QR分解。取，则，再取，则，故例7-13 求矩阵的QR分解，使得R的对角线元素为正数。解 的第一列，。用构造镜面反射阵，使得，令，有 ，的第2列对角线以下为，。用构造镜面反射阵，使得，令，易得 ，于是有 ， 容易验证，。请读者用平面旋转变换对本例的矩阵A进行QR分解。7.5.3 算法算法就

20、是利用分解构造一个矩阵序列，当充分大时，是近似的上三角矩阵，而该上三角阵的对角元素便是原始矩阵的全部特征值。 设，对做分解，即其中为上三角阵，为正交阵。利用这个分解可得新矩阵（对交换乘积）由于是经过正交相似变换得到的，因此与有相同的特征值。再对做分解，按上述方式又可得新矩阵，且与也具有相同的特征值。具体地说，其步骤为：设，做分解求矩阵求得后对作分解求矩阵 只要可逆，由定理9可知，按上述方法可唯一确定矩阵序列，且序列中任意与原始矩阵有相同特征值。因此只要恰当选择正交相似变换阵，使当时，逼近一个上三角阵，便可求出的全部特征值（为所逼近上三角阵的主对角元素）。可见，算法的关键在于选择正交变换阵。从定

21、理7-8的证明看到，正交变换阵是一系列平面转换矩阵的乘积，这些平面旋转矩阵是用来将的主对角线以下元素约化为零的。如果将算法直接应用于原始矩阵，计算量很大，所以在实际计算中，总是先将原始矩阵用豪斯赫尔德方法约化为上Hessenberg阵，而后再对上Hessenberg阵应用算法。可以证明，由上Hessenberg阵用算法生成的矩阵序列中的每个矩阵仍为上Hessenberg阵。7.5 雅可比方法雅可比方法是用来计算实对称矩阵的全部特征值及特征向量的一种有效方法。它的基本思想是，通过一组正交相似变换对称矩阵化为对角矩阵，得其全部特征值。定理7-10 设为阶对称矩阵，其中为正交矩阵，则证明 一方面另一

22、方面由假设，故。设为对称矩阵，为一平面旋转矩阵，则（其中）的元素计算公式为：（1）（2）（3）第行元素和第列元素（4）第行元素和第列元素（5）这说明，当经过一初等正交相似变换化为时，只需按上述公式计算的第列、第列元素，由对称性可得第行和第行元素，的其余元素与的对应元素相同。设的非对角元素，我们可选择平面旋转阵，使的非对角元素。由定理11可选择，使即选择，使定理7-11 设为对称阵，为的一个非对角元素，则可选择一平面旋转阵，使的非对角元素且与的元素满足下述关系（1）（2）（3）证明 由上面的计算公式直接计算可知（1）成立。由（1）及定理7-10可证（2）。如果用表示的非对角线元素的平方和，表示的

23、对角线元素平方和，则 ，这说明的对角线元素平方和比的对角线元素平方和增加了，的非对角线元素平方和比的非对角线元素平方和减少了。下面介绍雅可比方法。首先在的非对角元素中选择绝对值最大的元素（称为主元素），如可设，否则已经对角化了。由定理12，选择一平面旋转矩阵，使的非对角元素。再选的非对角元素中的主元素，如由定理12，又可选择一平面旋转矩阵，使的非对角元素（注意上次消除了的主元素这次又可能变为不是零）。继续这个过程，连续对实行一系列平面旋转变换，消去非对角线绝对值最大的元素，直到将的非对角元素全化为充分小为止，从而求得的全部（近似）特征值。定理7-12 （雅可比方法的收敛性）设为实对称矩阵，对施行上述一系列平面旋转变换则 证明 记，由定理7-11的（2）可得其中 又由于即由以上得反复应用上式，即得故 可以证明存在。下面介绍特征向量的计算。由雅可比收敛定理知，当充分大时记，则的列向量就是的近似特征向量。计算可采用累积的办法，用一数组保存，开始时，以后对每进行一次平面旋转变换，就进行计算用初等正交阵右乘只需计算的两列元素，若记，则