化实矩阵为准三角形的方法.doc

.2 化实矩阵为准三角形的方法矩阵特征值的QR算法，通常在准三角阵的基础上进行。任何矩阵都可经有限次相似变换化成准三角形，但这里只介绍对实矩阵应用豪斯雷尔德(Householder)变换的方法。.2.1 准三角阵下述形式的矩阵称为上准三角形矩阵，简称准三角阵：例如，下列矩阵A和B都是准三角阵：或把元素所在的对角线称为对角线，所在的对角线称为次对角线，则准三角阵就是次对角线以下的元素(ji1)全等于零的矩阵。定义 如果准三角阵的次对角线元素都不等于零，则称之为是不可约的；反之，如果次对角线有零元素，则称之为是可约的。例如，上述矩阵A是不可约的，而B是可约的。对于矩阵的列可给出相应的定义，例如矩阵H的第i列称为准三角列。0时，称不可约，当=0时，称可约。因此，不可约准三角阵的前n1列都是不可约的，可约准三角阵在前n1列中至少有一列是可约的。.2.2 化准三角形的豪斯霍尔德方法对于任意给定的实矩阵可通过n2次豪斯霍尔德相似变换，自左至右逐列化成准三角形。先看矩阵的第一列。如果该列已是准三角列，则无需变换，否则可应用定理的推论，构造一个初等对称正交矩阵，进行豪斯霍尔德变换=化的第一列为准三角列。事实上，如果取第一个分量等于零的列向量则按式(-5)构成的初等对称正交矩阵有对角分块形式：其中，而。将矩阵作相应的分块，并令其中，于是用对进行豪斯霍尔德变换得=由此可见，在矩阵的第一列中，对角线以下的部分为按定理4的推论，若取则必有，即的第一列为准三角列。依此类推，设经i1次豪斯霍尔德变换后得出前i1列为准三角列的矩阵 （-10)如果的第i列不是准三角列，则可构成矩阵，进行豪斯霍尔德变换=，化第i列为准三角列。这时，取构成的列向量因为的前i个元素等于零，故有对角分块形式其中，而。将矩阵按式(V-10)作相应的分块，并令于是豪斯霍尔德变换=从而矩阵的相应块。注意到仅最后一列非零，即知的前i1列仍是准三角列，并且第i列对角线以下的部分按定理4的推论，若取 （-11）则有，即的第i列也是准三角列。式(-11)中的正负号可任取，但为了减小计算误差，应使向量的第一个分量有较大的模，故通常取与异号。综上所述，用豪斯霍尔德方法化实矩阵A=为准三角形的全过程是=，=， （-12)即。其中，若令，并用sign（a）表示a的符号，则初等对称正交矩阵的计算公式是 （-13)例-4 用豪斯霍尔德方法化下列矩阵为准三角形：解 四阶矩阵应作两次豪斯霍尔德变换。第一次令i=1，由式(-13)和=，依次算出，。第二次令i=2，由式(-13)和=，即可依次算出最终结果，。.3 QR算法的基本原理本节简要地介绍QR算法的一般原理及改进收敛性和提高计算速度的两项措施先化原始矩阵为准三角形，然后再进行原点位移。具体计算方法在下节讨论。.3.1 一般原理定理1指出，分块三角形矩阵的特征值就是各对角块的特征值。如果对角块都是11或22阶矩阵，则可直接求出它的全部特征值。由此设想：对于任意给定的矩阵A=，如果对它连续作这样的相似变换 (V-14)使所得矩阵序列A收敛为分块三角形矩阵，并且对角块都是11或22阶矩阵，那么就可以求出原始矩阵A的全部特征值。实现这一设想的早期做法是LR算法。对于过程(V-14)中的每一次迭代=，首先对进行三角分解：=。然后取变换矩阵=，从而得=。其中，是单位下三角矩阵，是上三角矩阵。LR算法的缺点是数值计算不稳定，并且除某些特殊类型的矩阵外收敛性很差。为了克服这些缺点，人们提出了QR算法：每次迭代首先把分解成U矩阵和上三角矩阵的乘积，即QR分解 = (-15)然后取变换矩阵=，从而由和上式得= (-16)矩阵到的这种U相似变换，称为QR变换。由于U矩阵的行和列都是单位向量，所以QR算法的显著优点是数据计算稳定，但是就收敛性和每次迭代的计算量来说，并不比LR算法好。解决这一矛盾的方法是迭代前先把原始矩阵化成准三角形，迭代的每一步进行原点位移。对于这两者，后面将分别加以讨论。从上可知，QR算法的基础是矩阵的QR分解。与三角分解不同，QR分解是五条件的，对此有：定理5 对于任何矩阵A (一般地可是复数矩阵)，存在U矩阵Q和对角线元素为非负实数的上三角矩阵R，使分解式AQR成立。如果A非奇异，则分解是惟一的。本定理的证明基于下述引理：引理 对于任何复元素向量，存在U矩阵M，使线性变换。证明 如果b是零向量，则M可以是任何U矩阵，例如MI。如果b非零，则可取 M=TD，其中D为与向量b相应的对角U矩阵Ddiag。因为，故矩阵T有两种可能的取法：当时取T=I；当至少有一个为非零元素时，根据定理4的推论，取T为初等对称正交矩阵。注意到U矩阵的乘积仍是U矩阵，引理得证。对于n阶矩阵A，可根据引理取n1个U矩阵，使得用它们连续左乘A时，自左至右逐列消去对角线以下的元素，从而把A化成对角线元素为非负实数的上三角矩阵 (-17)事实上，如果经i1次变换已化A为，并且的前i1列为三角列则下次变换令，其中的ni+1阶U矩阵，按引理对列向量选取，使得线性变换。于是，=的前i列都是上三角列。由式(-17)可得，若令，则得A=QR，这就证明了存在性。如果A非奇异，则由和易知R的对角线元素恒为正实数。下面证明惟一性。设A有两个QR分解式和，则由=得= （-18)其中是U矩阵，从而也是U矩阵。于是 （-19)因为上三角矩阵的逆和积都是上三角矩阵，而它的共轭转置则是下三角矩阵，故由式 (V-19)可知是对角矩阵。注意到式(V-18)，应是对角U矩阵，因为和的对角线元素为正实数，所以只能是单位矩阵=I，从而=。可见，非奇异矩阵A的QR分解是惟一的。应该指出，从上述证明可知，如果A是实矩阵，则Q可取正交矩阵。此外，关于R的对角线元素的规定，在QR算法中并无必要，实际上对于非奇异矩阵，无论怎样进行QR分解，可以证明Q和R各元素的模都是惟一确定的。 V.3.2 准三角阵的QR算法 准三角阵有两个有利于QR算法的性质： 第一，可约准三角阵可分割成分块准三角阵。例如，矩阵就是这样，其中的对角块都是不可约准三角阵。因此，对于准三角阵的QR算法，可归结为对不可约准三角阵的算法。 第二，QR变换保持准三角形不变。具体地说，有 定理6 如果A是不可约准三角阵，则经QR变换A=QR，B=Q*AQ=RQ，所得矩阵B仍是准三角形，并且当A非奇异时B不可约，当A奇异时B可约。其中的U矩阵Q，总是不可约准三角阵。 按此定理，如果初始不可约准三角阵奇异，则经一次QR变换即可约简成较小的不可约准三角阵；如果非奇异，则用QR算法将得到不可约准三角阵序列。对于后者，只要有某个次对角线元素收敛为零，即可约简成较小的不可约准三角阵。因此，基于准三角阵的QR算法，一般来说收敛性比满矩阵为好。准三角阵更突出的优点是计算量小。事实上，由定理6可知，基于准三角阵的 QR算法，每次迭代都是对不可约准三角阵按式(V-15)和式(V-16)作QR变换。因为和都是准三角阵，所以其中的QR分解，以及U相似变换的计算量都比较小。QR分解就是化为上三角形(=)，U相似变换就是矩阵乘积=。可以推算出：满矩阵QR变换的计算量与n3成正比，而准三角阵QR变换的计算量仅与n2成正比。由此可见，当矩阵的阶数n较大时，基于准三角阵的QR算法，可使计算量大幅度下降。定理6的证明如下：令矩阵A和Q的列分别为和，由分解式A=QR可得 （i=1,2,，n1）如果A不可约，则它的前n1列都不可约。由=不可约，推知0，并且是不可约列；由不可约，进一步推知0，并且是不可约列；如此推下去，最后知0，并且是不可约列。这就是说，当A不可约时Q也是不可约准三角阵，并且R的前n1个对角线元素都不等于零。显然，当A非奇异时0，当A奇异时=0。 由于上三角矩阵和准三角阵的乘积仍是准三角阵，故B=RQ是准三角阵。注意到B的次对角线元素，由Q不可约，即0可知：当A奇异时不可约；当A奇异时B可约，并且它的次对角线元素仅=0。证毕。 V.3.3 收敛性和原点位移可以证明：如果非奇异矩阵A的各个特征值其模都不相同，则在QR变换下对角线以下的元素将趋于零，对角线以上元素的模趋于确定的值，而对角线元素则趋于各特征值。若是准三角阵，则在通常情况下=其中，和是与迭代次数k无关的常数；。这就是说次对角线元素以因子线性地收敛为零，对角线元素以因子线性地收敛为特征值，其中最后一个对角线元素的收敛因子是，并且特征值按模从大到小在对角线上顺序排列。在实际计算中，线性收敛是不满意的。但是，从上面的分析可以看出，如果每次迭代前将的特征值都减去一个尽量接近于的数值 (这相当于将特征值的坐标原点位移)，并在迭代后加回去，那么收敛因子将减小为。特别是最后一个收敛因子将变得更小，即，从而使迅速趋向零，迅速趋向。这就是用特征值原点位移加速算法收敛的基本思想。在上述特定的情况下，作为的估计数值，显然应取位移值=。由于矩阵A+I的特征值与A的特征值满足关系=+(i=1,2,，n)，所以带原点位移的QR算法每次迭代分两步：第一步，对位移后的矩阵I按式(V-15)和式(V-16)作QR变换，即先QR分解I= (V-20) 然后用所得对I作U相似变换 （I）= (V-21)其中，矩阵的特征值都比的小。第二步，恢复特征值的原点，即取矩阵=+I (V-22)显然，与的特征值相同。从式(V-21)可知= I ，代入式(V-22)得= (V-23)这就是说，带原点位移的QR变换，可先按式(V-20)QR分解矩阵I，然后用所得对作U相似变换，即按式(V-23)求得。下面进一步讨论原点位移的方法。可以证明：对于具有等模特征值的矩阵A=，由QR变换得到的矩阵序列常常收敛于分块三角形，并且每个对角块对应于模相等的特征值。特别地，对于具有多重复共轭特征值的实矩阵，的极限形式是分块三角形，并且对角块为11或22阶矩阵，其中22阶对角块对应于一对复共轭特征值。由此可见，对于一般的矩阵，特别是对于具有复数特征值的实矩阵，采用加速最后一个对角线元素收敛的实数位移=是不适宜的。事实上，这时的不可能接近于任何复数特征值。考虑到的极限形式中最后一个对角块的阶数常常是22，作为改善办法，可取右下角子矩阵 的两个特性值和中模值靠近的一个当作位移值。实践表明，这种原点位移方法对加速QR算法收敛是很有效的。 V.3.4 二重QR算法 对于实矩阵A=，当取位移值为右下角22阶子块的特征值时，可能是复数，从而按式(V-20)和式(V-23)进行的QR变换需作复数运算。为了减轻计算量，避免实矩阵的复数运算，可采用把连续两次迭代合并为一次的二重QR算法，它的一般原理如下： 设经k1次二重迭代后，原始矩阵已相似变换成。对于第k次二重迭代，首先从右下角22阶子矩阵求出位移值和，然后用它们按式(V-20)和式 (V-23)连续作两次QR变换。第一次，对矩阵用位移值作QR变换I=+ I =第二次，对矩阵用位移值作QR变换 I=+ I=如果把这两次变换合并为一次，并且令Q=，R=(I)( I) (V-24)则由上列各式可推知=QR (V-25)=Q (V-26)其中，式(V-25)可用下法推证：QR= ( I) =( I) =( I) =( I) (I)= (I) ( I)如果和都不是的特征值，则I、 I和都是非奇异矩阵。因此，根据QR分解的惟一性，二