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

Ⅴ.2 化实矩阵为准三角形的方法
矩阵特征值的QR算法，通常在准三角阵的基础上进行。任何矩阵都可经有限次相似变换化成准三角形，但这里只介绍对实矩阵应用豪斯雷尔德(Householder)变换的方法。
Ⅴ.2.1 准三角阵
下述形式的矩阵称为上准三角形矩阵，简称准三角阵：

例如，下列矩阵A和B都是准三角阵：

或把元素所在的对角线称为对角线，所在的对角线称为次对角线，则准三角阵就是次对角线以下的元素(j<i—1)全等于零的矩阵。
定义 如果准三角阵的次对角线元素都不等于零，则称之为是不可约的；反之，如果次对角线有零元素，则称之为是可约的。
例如，上述矩阵A是不可约的，而B是可约的。
对于矩阵的列可给出相应的定义，例如矩阵H的第i列

称为准三角列。≠0时，称不可约，当=0时，称可约。因此，不可约准三角阵的前n－1列都是不可约的，可约准三角阵在前n－1列中至少有一列是可约的。
Ⅴ.2.2 化准三角形的豪斯霍尔德方法
对于任意给定的实矩阵

可通过n－2次豪斯霍尔德相似变换，自左至右逐列化成准三角形。
先看矩阵的第一列。如果该列已是准三角列，则无需变换，否则可应用定理４的推论，构造一个初等对称正交矩阵，进行豪斯霍尔德变换=化的第一列为准三角列。
事实上，如果取第一个分量等于零的列向量

则按式(Ⅴ-5)构成的初等对称正交矩阵有对角分块形式：

其中，，而。将矩阵作相应的分块，并令

其中，，于是用对进行豪斯霍尔德变换得
==
由此可见，在矩阵的第一列中，对角线以下的部分为

按定理4的推论，若取

则必有，即的第一列为准三角列。
依此类推，设经i－1次豪斯霍尔德变换后得出前i－1列为准三角列的矩阵
（Ⅴ-10)
如果的第i列不是准三角列，则可构成矩阵，进行豪斯霍尔德变换=，化第i列为准三角列。这时，取构成的列向量

因为的前i个元素等于零，故有对角分块形式

其中，，而。将矩阵按式(V-10)作相应的分块，并令

于是豪斯霍尔德变换
==
从而矩阵的相应块。注意到仅最后一列非零，即知的前i－1列仍是准三角列，并且第i列对角线以下的部分

按定理4的推论，若取
（Ⅴ-11）
则有，即的第i列也是准三角列。
式(Ⅴ-11)中的正负号可任取，但为了减小计算误差，应使向量的第一个分量有较大的模，故通常取与异号。
综上所述，用豪斯霍尔德方法化实矩阵A=为准三角形的全过程是
=，=，…， （Ⅴ-12)
即……。其中，若令，并用sign（a）表示a的符号，则初等对称正交矩阵的计算公式是
（Ⅴ-13)
例Ⅴ-4 用豪斯霍尔德方法化下列矩阵为准三角形：

解 四阶矩阵应作两次豪斯霍尔德变换。第一次令i=1，由式(Ⅴ-13)和=，依次算出
，，
。
第二次令i=2，由式(Ⅴ-13)和=，即可依次算出最终结果
，，
。
Ⅴ.3 QR算法的基本原理
本节简要地介绍QR算法的一般原理及改进收敛性和提高计算速度的两项措施——先化原始矩阵为准三角形，然后再进行原点位移。具体计算方法在下节讨论。
Ⅴ.3.1 一般原理
定理1指出，分块三角形矩阵的特征值就是各对角块的特征值。如果对角块都是1×1或2×2阶矩阵，则可直接求出它的全部特征值。由此设想：对于任意给定的矩阵A=，如果对它连续作这样的相似变换
(V-14)
使所得矩阵序列{A}收敛为分块三角形矩阵，并且对角块都是1×1或2×2阶矩阵，那么就可以求出原始矩阵A的全部特征值。
实现这一设想的早期做法是LR算法。对于过程(V-14)中的每一次迭代=，首先对进行三角分解：=。然后取变换矩阵=，从而得==。其中，是单位下三角矩阵，是上三角矩阵。LR算法的缺点是数值计算不稳定，并且除某些特殊类型的矩阵外收敛性很差。
为了克服这些缺点，人们提出了QR算法：每次迭代首先把分解成U矩阵和上三角矩阵的乘积，即QR分解
= (Ⅴ-15)
然后取变换矩阵=，从而由和上式得
== (Ⅴ-16)
矩阵到的这种U相似变换，称为QR变换。
由于U矩阵的行和列都是单位向量，所以QR算法的显著优点是数据计算稳定，但是就收敛性和每次迭代的计算量来说，并不比LR算法好。解决这一矛盾的方法是迭代前先把原始矩阵化成准三角形，迭代的每一步进行原点位移。对于这两者，后面将分别加以讨论。
从上可知，QR算法的基础是矩阵的QR分解。与三角分解不同，QR分解是五条件的，对此有：
定理5 对于任何矩阵A (一般地可是复数矩阵)，存在U矩阵Q和对角线元素为非负实数的上三角矩阵R，使分解式A＝QR成立。如果A非奇异，则分解是惟一的。
本定理的证明基于下述引理：
引理 对于任何复元素向量，存在U矩阵M，使线性变换。
证明 如果b是零向量，则M可以是任何U矩阵，例如M＝I。如果b非零，则可取 M=TD，其中D为与向量b相应的对角U矩阵D＝diag。因为，故矩阵T有两种可能的取法：当时取T=I；当至少有一个为非零元素时，根据定理4的推论，取T为初等对称正交矩阵。注意到U矩阵的乘积仍是U矩阵，引理得证。
对于n阶矩阵A，可根据引理取n－1个U矩阵，使得用它们连续左乘A时，自左至右逐列消去对角线以下的元素，从而把A化成对角线元素为非负实数的上三角矩阵
(Ⅴ-17)
事实上，如果经i－1次变换已化A为，并且的前i－1列为三角列

则下次变换令，其中的n－i+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,…，n－1）
如果A不可约，则它的前n－1列都不可约。由=不可约，推知≠0，并且是不可约列；由不可约，进一步推知≠0，并且是不可约列；如此推下去，最后知≠0，并且是不可约列。这就是说，当A不可约时Q也是不可约准三角阵，并且R的前n－1个对角线元素都不等于零。显然，当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变换得到的矩阵序列{}常常收敛于分块三角形，并且每个对角块对应于模相等的特征值。特别地，对于具有多重复共轭特征值的实矩阵，{}的极限形式是分块三角形，并且对角块为1×1或2×2阶矩阵，其中2×2阶对角块对应于一对复共轭特征值。
由此可见，对于一般的矩阵，特别是对于具有复数特征值的实矩阵，采用加速最后一个对角线元素收敛的实数位移=是不适宜的。事实上，这时的不可能接近于任何复数特征值。考虑到{}的极限形式中最后一个对角块的阶数常常是2×2，作为改善办法，可取右下角子矩阵

的两个特性值和中模值靠近的一个当作位移值。实践表明，这种原点位移方法对加速QR算法收敛是很有效的。
V.3.4 二重QR算法
对于实矩阵A=，当取位移值为右下角2×2阶子块的特征值时，可能是复数，从而按式(V-20)和式(V-23)进行的QR变换需作复数运算。为了减轻计算量，避免实矩阵的复数运算，可采用把连续两次迭代合并为一次的二重QR算法，它的一般原理如下：
设经k－1次二重迭代后，原始矩阵已相似变换成。对于第k次二重迭代，首先从右下角2×2阶子矩阵求出位移值和，然后用它们按式(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分解的惟一性，二重QR变换可按式(V-24)~式(V-26)进行，即先算出并对它作QR分解，然后用所得Q对作U相似变换。实际上当和中有的特征值时，二重QR变换也可以这样做。可以证明，在这种情况下，不可约准三角阵二重QR变换的结果是可约的。
现在设是实矩阵，这时的位移值和若是复数则必为共轭对，从而它们的和与积都是实数。注意到式(V-24)的展开式：
=－ (+)+I
可见是实矩阵，于是由它分解出来的Q可取正交矩阵，按式(V-26)得到的仍是实矩阵。因此，对于实矩阵采用二重QR算法，可避免因原点位移导致的复数运算。

１２