第8章 矩阵特征值问题计算

第8章矩阵特征值问题计算

8.1引言

物理、力学和工程技术中的很多问题在数学上都归结为求矩阵特征值问题。例如，振动问题（大型桥梁或建筑物的振动、机械的振动、电磁振荡等），物理学中某些临界值的确定，这些问题都归结为下述数学问题。

定义１　（1）已知，则称

，

为的特征多项式．

的特征方程

(8.1.1)

一般有个根（实的或复的，重根按重数计算．当时，为实系数次代数方程，其复根共轭成对出现），称为的特征值．用表示的所有特征值的集合，并称之为的谱集．

（２）设为的特征值，相应的齐次方程组

　　　　　　　　　　　　　(8.1.2)

的非零解称为的对应于的特征向量．

例１　求的特征值及特征向量，其中

，

解　矩阵的特征方程为

，

求得的特征值为：．对应的各特征向量是对应齐次方程

，．

的解．计算得到：

．

下面我们来看几个有关特征值的结论．

定理8.1.1　设为的特征值，且，其中０，则

（１）为特征值（为非零常数）；

（２）为的特征值，即；

（３）为的特征值；

（４）设为非奇异矩阵，那么且为的特征值，即．

定理8.1.2　设为阶矩阵的特征值，则

（１）（是矩阵的迹）；

（２）．

定理8.1.3　设，则

．

定理8.1.4　设为分块上三角阵，即

，

其中每个对角块均为方阵，则．

定理8.1.5　设与为相似矩阵（即存在非奇异矩阵使），则

（１）与有相同的特征值；

（２）如果是的特征向量，则是的特征向量．

定理8.1.5说明，一个矩阵经过相似变换，则的特征值不变．

定义8.1.2　设，如果有一个重数为的特征值且对应于矩阵的线性无关的特征向量个数少于，称为亏损矩阵．

一个亏损矩阵是一个没有足够特征向量的矩阵，亏损矩阵在理论上和计算上都存在困难．

定理8.1.6（１）可对角化，即存在非奇异矩阵使

的充要条件是具有个线性无关的特征向量．

（２）如果有个（不同的特征值，则对应的特征向量线性无关．

定理8.1.7（对称矩阵的正交约化）设为对称矩阵，则：

（１）的特征值均为实数；

（２）有个线性无关的特征向量；

（３）存在一个正交矩阵使

且为特征值，而的列向量为的对应于的特征向量．

下面讨论矩阵特征值界的估计．

定义8.1.3　设．令：

（１）；

（２）集合．称复平面上以为圆心，以为半径的所有圆盘为的Gerschgorin圆盘．

定理8.1.8（Gerschgorin圆盘定理）（１）设，则的每一个特征值必属于下述某个圆盘之中

．

或者说，的特征值都在复平面上个圆盘的并集中．

（２）如果有个圆盘组成一个连通的并集，且与余下－个圆盘是分离的，则内恰包含的一个特征值．

证明　只就（１）给出证明，设为的特征值，即有使

．

记　，考虑的第个方程，即

，

或

，

于是

，

即

．

这说明，的每一个特征值必位于的一个圆盘中，并且相应的特征值一定位于第个圆盘中（其中是对应特征向量绝对值最大的分量的下标）．

利用相似矩阵性质，有地可以获得的特征值进一步的估计，即适当选取非奇异对角阵

，

并做相似变换．适当选取可使某些圆盘半径及连通性发生变化．

例２　估计矩阵

特征值的范围．

解　的３个圆盘为

由定理8.1.8，可知的３个特征值位于３个圆盘的并集中，由于是弧立圆盘，所以内恰好包含的一个特征值（为实特征值），即

．

的其他两个特征值包含在，的并集中．

现选取对角阵

做相似变换

．

的３个圆盘为

显然，３个圆盘都是孤立圆盘，所以，每一个圆盘都包含的一个特征值（为实特征值）且有估计

下面给出理论上有关通过酉相似变换及正交相似变换可以约化一般矩阵到什么程序的问题．

定理8.1.9（Schur定理）　设，则存在酉阵使

（上三角阵）

其中为的特征值．

当时，如果限制用正交相似变换，由于有复的特征值，不能用正交相似变换约化为上三角阵．用正交相似变换能将约化到什么程度呢？

定理8.1.10（实Schur分解）　设，则存在正交矩阵使

，

其中对角块为一阶或二阶方阵，且每一个一阶是的实特征值，每个二阶对角块的两个特征值是的两个共轭复特征值．

我们转向实Schur型的实际计算．

定义8.1.4　设为阶实对称矩阵，对于任一非零向量，称

为对应于向量的瑞利（Rayleigh）商．

定理8.1.11　设为对称矩阵（其特征值次序记为），则

（1） (对任何非零)；

（2）；

（3）．

证明　只证１，关于２，３留作习题．

由于为实对称矩阵，可将对应的特征向量正交规范化，则有．设为中任一向量，则有展开式

，

于是

．

从而１成立，结论１说明瑞利商必位于和之间．

关于计算矩阵的特征值问题，当时，我们还可按行列式展开的办法求的根．但当较大时，如果按展开行列式的办法，首先求出的系数，再求的根，工作量就非常大，用这种办法求矩阵特征值是不切实际的，由此需要研究求的特征值及特征向量的数值解法．

本章将介绍一些计算机上常用的两类方法，一类是幂法及反幂法（迭代法），另一类是正交相似变换的方法（变换法）．

8.2幂法及反幂法

8.2.1 幂法

幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，特别适用于大型稀疏矩阵。反幂法是计算海森伯格或三对角阵的对应一个给定近似特征值的特征向量的有效方法之一。

设实矩阵有一个完全的特征向量组，其特征值为，相应的特征向量为．已知的主特征值是实的，且满足条件

，　 (8.2.1)

现讨论求及的方法．

幂法的基本思想是任取一个非零的初始向量，由矩阵构造一向量序列

(8.2.2)

称为迭代向量．由假设，可表示为

，　　　　(8.2.3)

于是

其中，显然．

由假设知，故

． (8.2.4)

这说明序列越来越接近的对应于的特征向量，或者说当充分大时

， (8.2.5)

即迭代向量为的特征向量的近似向量．

下面再考虑主特征值的计算，用表示的第个分量，则

，　　　　　　　　(8.2.6)

因此

，　　　　　　　　　　　(8.2.7)

也就是说两相邻迭代向量分量的比值收敛到主特征值．

这种由已知非零向量及矩阵的幂构造向量序列{}以计算的主特征值（利用(8.2.7)式）及相应特征向量（利用(8.2.5)式）的方法称为幂法．

由(8.2.6)知，的收敛速度由比值来确定，越小收敛越快，当接近于１时收敛会很慢．

定理8.2.1　设有个线性无关的特征向量，主特征值满足

，

则对任何非零初始向量，(8.2.4), (8.2.7)式成立．

如果的主特征值为实的重根，即，且

，

又设有个线性无关的特征向量，对应的个线性无关特征向量为，则由(8.2.2)式

这说明当的主特征值是实的重根时，定理５的结论还是正确的．

应用幂法计算的主特征值及对应的特征向量时，如果＞１（或＜１）,迭代向量的各个不等于零的分量将随而趋向于无穷（或趋向于零），这样在计算机实现时就可能“溢出”．为了克服这个缺点，就需要将迭代向量加以规范化．

设有一向量，将其规范化得到向量

其中表示向量的绝对值最大的分量，即如果有

，

则＝，且为所有绝对值最大的分量中的最小下标．

在定理8.2.1的条件下幂法可这样进行：任取一初始向量，构造向量序列

由(8.2.3)式

　　(8.2.8)

这说明规范化向量序列收敛到主特征值对应的特征向量．

同理，可得到

收敛速度由比值确定．总结上述讨论，有

定理8.2.2　设有个线性无关的特征向量，主特征值满足，则对任何非零初始向量，按下述方法构造的向量序列：

　　　　　　　(8.2.9)

则有

（１），

（２）．

例３　用幂法计算

的主特征值和相应的特征向量．计算过程如表８－１．

[解]下述结果是用Visual C++6.0在微机上进行运算得到的，的分量值是舍入值．于是得到2.5365326，及相应的特征向量（0.74822631，0.64966657，1)．和相应的特征向量的真值（８位数字）为

表８-１　幂法计算结果

K

（规范化向量）

0

（1，

1，

1）

 

1

（0.90909091，

0.81818182,

1）

2.7500000

2

（0.83760684，

0.74358974,

1）

2.6590909

3

（0.79901559，

0.70303527,

1）

2.6047009

4

（0.77741499，

0.68033766,

1）

2.5752666

5

（0.76510819，

0.66740583,

1）

2.5587919

6

（0.75802535，

0.65996327,

1）

2.5494056

7

（0.75392531，

0.65565500,

1)

2.5440035

8

（0.75154396，

0.65315271,

1)

2.5408764

9

（0.75015815，

0.65169652,

1)

2.5390602

10

（0.74935078，

0.65084814,

1)

2.5380032

11

（0.74888009，

0.65035355,

1)

2.5373874

12

（0.74860558，

0.65006510,

1)

2.5370284

13

（0.74844545，

0.64989683,

1)

2.5368191

14

（0.74835202，

0.64979867,

1)

2.5366969

15

（0.74829751，

0.64974139,

1)

2.5366257

16

（0.74826571，

0.64970797,

1)

2.5365841

17

（0.74824715，

0.64968847,

1)

2.5365598

18

（0.74823632，

0.64967709,

1)

2.5365457

19

（0.74823000，

0.64967045,

1)

2.5365374

20

（0.74822631，

0.64966657,

1)

2.5365326

8.2.2 加速方法

1．原点平移法

由前面讨论知道，应用幂法计算的主特征值的收敛速度主要由比值来决定，但当接近于１时，收敛可能很慢．这时一个补救的办法就是采用加速收敛的方法．

引进矩阵

其中为选择参数．设的特征值为，则的相应特征值为，而且，的特征向量相同．

如果需要计算的主特征值，就要适当选择，使－仍然是的主特征值，且使

．

对应用幂法，使得在计算的主特征值－的过程中得到加速．这种方法通常称为原点平移法．对于的特征值的某种分布，它是十分有效的．

例４　设有特征值

，

比值，作变换

（＝12）,

则的特征值为

．

应用幂法计算的主特征值的收敛速度的比值为

．

虽然常常能够选择有利的值，使幂法得到加速，但设计一个自动选择适当参数的过程是困难的．

下面考虑当的特征值是实数时，怎样选择使采用幂法计算得到加速．

设的特征值满足

，　　 (8.2.10)

则不管如何，的主特征值为或．当我们希望计算及时，首先应选择使

，

且使收敛速度的比值

．

显然，当，即时为最小，这时收敛速度的比值为

．

当的特征值满足(8.2.10)且能初步估计时，我们就能确定的近似值．

当希望计算时，应选择

，

使得应用幂法计算得到加速．

例５　计算例３中矩阵的主特征值．

作变换，设＝0.75，则

．

对应用幂法，计算结果如表８.２．

表8-2　幂法计算结果

K

（规范化向量）

1

2

3

4

5

6

7

8

9

10

0.875

0.78333333

0.76834862

0.75382166

0.7516

0.74913041

0.7487941

0.74836879

0.74831858

0.74824505

0.75

0.7

0.66513761

0.65796178

0.65217778

0.65105965

0.65007315

0.64989782

0.64972854

0.64970127

1

1

1

1

1

1

1

1

1

1

2.0000000

1.8750000

1.8166667

1.8004587

1.7914013

1.7888444

1.7873301

1.7869153

1.7866589

1.7865914

由此得的主特征值为1.7865914，的主特征值为

＝2.5365914.

与例３结果比较，上述结果比例３迭代１５次还好．

原点位移的加速方法，是一个矩阵变换方法．这种变换容易计算，又不破坏矩阵的稀疏性，但的选择依赖于对的特征值分布的大致了解．

2．瑞利商加速

由定理8.1.11知，对称矩阵的及可用瑞利商的极值来表示．下面我们将把瑞利商应用到用幂法计算实对称矩阵的主特征值的加速收敛上来．

定理8.2.3　设为对称矩阵，特征值满足

，

对应的特征向量满足，应用幂法（公式(8.2.9)）计算的主特征值，则规范化向量的瑞利商给出的

．

证明　由(8.2.8)式及

，

得

．　　(8.2.11)

8.2.3 反幂法

反幂法用来计算矩阵按模最小的特征值及其特征向量，也可用来计算对应与一个给定近似特征值的特征向量。

设为非奇异矩阵，的特征值依次记为

，

相应的特征向量为则的特征值为

，

相应的特征向量为．

所以计算的按模最小的特征值的问题就是计算的按模最大的特征值问题。

对于应用幂法迭代(称为反幂法)，可求得矩阵的主特征值１/，从而求得的按模最小的特征值。

反幂法迭代公式为:

任取初始向量，构造向量序列

．

　迭代向量可以通过解方程组求得．

定理8.2.4　设为非奇异矩阵且有个线性无关的特征向量，其对应的特征值满足

，

则对任何初始非零向量，由反幂法构造的向量序列满足

（１），

（２）．

收敛速度的比值为．

在反幂法中也可以用原点平移法来加速迭代过程或求其他特征值及特征向量。

如果矩阵存在,对其应用幂法，得反幂法的迭代公式

　　　　　　　　　　　　　 （2.12）

如果p是的特征值的一个近似值，且设与其他特征值分离,即

(),

就是说是的主特征值，可用反幂法(2.12)计算特征值及特征向量。这样我们有下面的定理

　　　 定理8.2.5设有个线性无关的特征向量,的特征值及对应的特征向量分别记为及(i=1,2,…,n),而为的近似值,存在，且

()，

则对任意的初始非零向量，有反幂法迭代公式(2.12)构造的向量序列满足

(1) ;

(2) ,即

且收敛速度由比值

确定。

由该定理知，对(其中)应用反幂法，可用来计算特征向量.只要选择的是的一个较好的近似且特征值分离情况较好,一般很小，常常只要迭代一二次就可完成特征向量的计算。

反幂法迭代公式中的是通过解方程组

求得的。为节省工作量，可先将进行三角分解　

，

其中P为某个排列阵，于是求相当于解两个三角形方程组

实验表明，按下述方法选择较好的：选使

(2.13)

用回代求解(2.13)即得，然后再按公式(2.12)进行迭代。

反幂法迭代公式:

　　　　　 1. 分解计算

,且保存信息。

　　　　　 2. 反幂法迭代　

　　　　 (1) 解求

　

(2) k=2,3,…

1) 解，求得;
解，求得

2)

3) 计算

例6用反幂法求

的对应于特征值λ=1.2679(精确特征值为)的特征向量.

[解]用部分选主元的三角分解将 (其中)分解为

，

其中

由，得

由，得

对应的特征向量是

，

由此看出是的相当好的近似。

特征值，的真值为。

8.3豪斯霍8.4尔德方法

8.4.1 引言

本节讨论下面问题:

(1) 用初等反射阵作正交相似变换约化一般实矩阵为上海森伯格阵(次对角元素下方的元素均为零)。

(2) 用初等反射阵作正交相似变换约化对称矩阵为对称三对角阵。

于是，求原矩阵特征值问题,就转化为求上海森伯格阵或对称三对角阵的特征值问题。

8.4.2 用正交相似变换约化一般矩阵为上海森伯格阵。

　　 设。下面说明可选择初等反射阵使经正交相似变换约化一个上海森伯格阵。

　 （1）设

其中,

当时进行下一步。当时，对维列向量待定，令

这是阶初等反射阵()，使其满足,为维向量)　即　

易见

现在令

,

则

其中

　　　 (2) 第步约化：重复上述过程，设对已完成第1步，…，第－1步正交相似变换，即有

其中为阶上海森伯格阵，．

设，于是可选择初等反射阵使,其中的计算公式为

(8.3.2)

令

,

则

(8.3.3)

其中为+1阶上海森伯格阵。第k步约化只需计算及 (当为对称阵时可以不算)。

(3) 重复上述过程，有。

定理8.3.1 (矩阵的上Hessenberg化)设,则存在初等反射阵使　

（上海森伯格阵）

算法1(Householder约化上H-矩阵)设,本算　法计算 (上海森伯格型)，其中为初等反射阵的乘积。

　　　 1.

　　　 2. 对于

计算初等反射阵使，即取

　（2）约化计算

（3）

算法分析：　　

本算法约需要次乘法运算，要明显形成还需加次乘法。

在计算的过程，的计算可能会引起溢出现象．为避免溢出现象，可以考虑向量的规范化，即，由计算得到．但应注意不能将保存到中．

在实际计算时，不必要得出矩阵形式，只需要得出算法中交待的相关参数及向量即可．

　　　 例7用豪斯霍尔德方法将

矩阵约化为上Hessenberg阵。

　　 解选取初等反射阵使　R₁ , 其中．

计算:（规范化）

则有　．

　(2) 约化计算：令

则

8.4.3 用正交相似变换约化对称阵为对称三对角阵

定理8.3.2(豪斯霍尔德约化对称阵为对称三对角阵)设,则存在初等反射阵使　

证明　由定理8.3.1，存在初等反射阵使，注意到，显然（对称矩阵），因此若对称，则必有对称，从而是对称三对角阵．

定理证毕

根据定理１７的结论，我们容易对算法１加以改进得之如下算法２．

算法２（豪斯霍尔得德约化对称阵为对称三对角阵）设为对称矩阵，本算法确定初等反射阵使为对称三对角阵．Ｃ的对角元存放在数组内，Ｃ的次对角元存放在数组内，数组的初始值可用来存放及，确定中向量的分量存放在Ａ的相应位置．冲掉，约化Ａ的结果冲掉Ａ，数组Ａ的上部分元素不变，如果第不需要变换则置为零．

1. 对于

(1)

(2) 确定变换

1) 计算

2) 如果，则，转4）

3) 计算

4)

5)

6)

7)

(3) 应用变换

1)

2) 计算及

对于

(a)

(b)

3) 计算

4) 计算对角线以下部分

对于

5) 继续循环

2.

对对称矩阵用初等反射阵正交相似约化为对称三对角阵大约需要次乘法。

用正交矩阵进行相似约化有一些特点，如构造的容易求逆，且的元素数量级不大，这个算法是十分稳定的。

8.5QR 方法

QR方法是一种变换方法，是计算一般矩阵(中小型矩阵)全部特征值问题的最有效方法之一。目前QR方法主要用来计算：(1)上海森伯格阵的全部特征值问题;(2)计算对称三对角矩阵的全部特征值问题。QR方法具有收效快，算法稳定等特点。

对于一般矩阵 (或对称矩阵),则首先用豪斯霍尔德方法将化为上海森伯格阵(或对称三对角阵),然后再用QR方法计算的全部特征值。

8.5.1 　基本QR算法

设，基本QR迭代算法如下：

　　　　　　（8.4.1）

其中为上三角阵，为正交阵。

定理8.4.1（基本QR算法的性质定理）记，则有

（1）；

（2）；

（3）。

[证明](1);

(2) ;

(3) 事实上，

由（2）知，代入上式得

，

从而便有

．

［证毕］

定理?????　设阶矩阵的个特征值满足，并设阶非奇异矩阵的第行是对应于的左特征向量．如果有ＬＵ分解，则由(8.4.1)产生的矩阵的对角线以下的元素趋向于零，同时对角元素趋向于

［证明］令：，则有．假定的LU分解为，其中是单位下三角矩阵，是上三角矩阵．这样，我们有

???????

其中．由于是单位下三角矩阵，而，故必有

． ???????

现令的ＱＲ分解为：．由的非奇异，可要求的对角元素均为正负．将这一分解代入???????可得

． ???????

当充分大时，是非奇异的，故它有如下的ＱＲ分解

???????

其中的对角元素均为正数．从(8.4.6)和(8.4.8)知

于是有

(8.4.9)

将???????代入?????７?，得

即已经找到了的一个ＱＲ分解．为保证分解中的上三角矩阵的对角元素均为正数，令

其中是的第个对角元素．于是我们有

将上式与(8.4.4)比较，并注意到ＱＲ分解的唯一性，就有

将上式代入(8.4.3)即有

再注意到

代入上式便得

．

由于，故

而是上三角矩阵，加之和是对角矩阵，因此，主对角线以下的元素趋于０．

定理证毕

8.5.2 上Hessenberg化

实际计算时，为了减少每次迭代所需的运算量，总是先将原矩阵Ａ经相似变换约化为一个上准三角矩阵——上Hessenberg矩阵，然后对约化后的矩阵进行ＱＲ迭代．

定义8.4.1若，当时．则称为上Hessenberg矩阵或上准三角矩阵．

下面我们希望讨论这样一个问题：对，存在上Hessenberg矩阵与之相似．也就是说，总存在一个非奇异矩阵使得

对于给定的，我们经次正交变换，就可将其约化为上Hessenberg矩阵．

首先，我们将矩阵分块写成

，

其中

第一步：构造Householder变换，使其满足

，

其中．是阶单位矩阵的第一列．再令

， (8.4.10)

于是

显然

．

也就是说经第一步正交变换后，我们得到与相似的矩阵形如

．

第二步：先对，进行类似于第一步的正交变换，将写成

构造阶Householder变换，使得是阶单位矩阵的第一列．令

使得

显然

从而令

于是

如此进行步，就可找到个Householder变换，使得

．

其中

．

再令

则有

(8.4.11)

通常称该式为A的上Hessenberg分解。

注意：根据以上Hessenberg约化过程可以发现，正交矩阵的第一列为，即．

算法8.4.1（计算上Hessenberg分解：Householder变换）

这一算法计算出的上Hessenberg矩阵就存放在Ａ所对应的存储单元内．运算量为；如果需要累积，则还需要再增加运算量．

此外，上述算法计算得到的上Hessenberg矩阵满足

其中是正交矩阵，，这里是一常数，是机器精度．

当然，我们亦可用Givens变换将上Hessenberg化，一般所需要的运算量大约是算法8.4.1的二倍．但是，如果有较多的零元素，则适当安排Givens变换的次序，可使运算量大为减少．另外，为了节省运算量，也可采用列主元的Gauss消去法将约化为上Hessenberg矩阵．不过这样做，虽然运算量少，但数值稳定性差．

尽管一般来讲上Hessenberg分解是不唯一的，然而我们可以证明

定理8.4.3　设有如下两个上Hessenberg分解：

（8.4.12）

其中和是阶正交矩阵，和是上Hessenberg矩阵．若，而且的次对角元素均不为零，则存在对角元素均为１或－１的对角矩阵使得

(8.4.13)

[证明]　假定对某个已证

(8.4.14)

其中．下面来证明存在使得

从（8.4.12）可得

分别比较上面两个矩阵等式的第列，可得

（8.4.15）

（8.4.16）

分别在(8.4.15)和(8.4.16)两边左乘，可得

再利用(8.4.16)就有

(8.4.17)

将(8.4.17)代入(8.4.15)，并利用 (6.4.16) 和(8.4.18)，可得

(8.4.18)

注意到，则

而，故(8.4.18)蕴含着

其中．

定理证毕

定义8.4.2一个上Hessenberg矩阵，如果其次对角元素均不为零，即，则称它是不可约的．

上述定理表明：如果这不可约的上Hessenberg矩阵，其中为正交矩阵，则和完全由的第一列确定（为里是在相差一个正负号的意义下的唯一）．

现在假定是上Hessenberg矩阵，我们来考虑对进行一次ＱＲ迭代的具体实现问题．第一步是计算的ＱＲ分解．由于的特殊性，这一步可用个平面旋转变换来完成，其具体计算过程如下：

第１次约化：作Givens变换

则

其中仅有前两行的元素与的元素不同，其余各行元素与均相同，并且

，　

直至第k－１次约化毕，我们得到

第k次约化：令

作Givens变换

，

那么

，

其中仅第k行和第k+1的元素与的元素不同，其余各行元素与均相同，并且

，　

这样，我们便知

是上三角矩阵．令，则，即这样就已完成了的ＱＲ分解．要完成一次ＱＲ迭代，我们还须计算

．

由于是（１，２）坐标平面内的旋转变换，因此仅有前两列与不同，而的前两列由的前两列的线性组合构成．又是上三角矩阵，故必有如下形状：

同样，是（２，３）坐标平面内的旋转变换，仅有第二列和第三列与不同，它们是的和第二列和第三列的线性组合，故有如下形状：

如此进行下去，最后我们得到的仍是上Hessenberg矩阵．而且不难算出，这样进行的一次ＱＲ迭代的运算量是．注意，对一般方阵进行一次ＱＲ迭代的运算量是．

8.5.3 带原点位移的QR迭代

从定理8.4.2已经知道，基本的QR算法是线性收敛的，其收敛速度取决于特征值之间的分离程度。为了加速其收敛速度，类似于反幂法，可引进原点位移。设第步迭代的位移为，则带原点位移的QR迭代如下：

这里是给定的上Hessenberg矩阵。

现在来讨论位移的选取。由于为上Hessenberg矩阵，故其最后一行仅有两个非零元素，若QR算法收敛，则当充分大时，就很小，因而就接近于的一个特征值。这样，根据从反幂法所得到的经验，我们可以选取位移为。事实上，对于这样选取的位移，可以证明，若很小的话，则经一次带原点位移QR迭代后，成立

（8.4.19）

显然，这只需考虑右下角的子矩阵的变化即可。由前面的讨论可知，将约化成上Hessenberg阵有步，现假定前面步已完成，此时的变为

，

这是因为前步不改变的最后一行。约化的第步就是要消去，即确定使得

．

从平面旋转变换的确定方法易知，此处

．

这样，通过简单的计算可知

，

即(6.4.24)成立．我们看到通过原点位移，特征值的渐近收敛速度从线性收敛加速面变成了二次收敛．

8.5.4 双重步位移的QR迭代

上面所讨论的带原点位移的ＱＲ迭代，存在严重的缺点：若具有复共轭特征值，则实位移一般并不能起到加速的作用．为了克服这一缺点，我们下面来介绍双重步位移的QR迭代，其基本思想是将带原点位移的ＱＲ迭代合并为一步，以避免复数运算．

设，考虑如下的迭代：

（8.4.20）

不失一般性，我们可以假定迭代（8.4.20）中出现的上Hessenberg矩阵都是不可约的．因若不然，迭代的某一步，已有

，

则我们可以分别对和进行QR迭代即可．

大家已经知道，在一定条件下取位移可起到加速收敛的作用；然而，大家也熟知实矩阵可以有复特征值．这样，假如的尾部子矩阵

有一对复共轭特征值和时，我们就不能期望最终收敛于Ａ的某个特征值，因而此种情形再取位移为就完全起不到加速收敛的作用．为了加速收敛，此时我们自然应该取或作位移．但这样一来就必须涉及复数运算，而这又是我们所不希望的．为了避免复运算的出现，人们想用和连续作两次位移，即进行

这里我们记．

但是，这里的和都是复数，要避免复数运算肯定不能按照上面给出的迭代直接进行，它只能是一个算法思想。

非常巧妙的是，上面的算法可以通过实运算得到．为了证实这一点，我们令

， (8.4.21)

, (8.4.22)

则

另外

即有

， (8.4.23)

, (8.4.24)

另一方面，由(8.4.23)可得

(8.4.25)

其中

因此是一个实矩阵；而且如果和均不是的特征值，并假定在迭代过程中选取和的对角元素均为正数，则由(8.4.23)可推知，Ｑ亦是实的；从而由(8.4.24)知也是实的．这也就是说，在没有误差的情况下，用和连续作两次位移进行QR迭代产生的仍是上Hessenberg矩阵．但是，实际计算时，由于舍入误差的影响，如此得到的一般并不一定是实的．这样为了确保计算得到的仍是实的，根据(8.4.23)和(8.4.25)，我们自然想到按如下的步骤来计算：

算；

计算的QR分解：；

计算．

然而，如此计算的第一步形成的运算量就是．当然，这是我们不希望的．幸运的是，定理8.4.3告诉我们：不论采用什么样的方法去求正交矩阵使是上Hessenberg矩阵，只要保证的第一列与的第一列一样，则就与本质上是一样的（所有元素的绝对值都相等）．当然，这需要是不可约来加以保证的．因此是不可约的，则我们就可有很大的自由度去寻求更有效的方法来实现由到的变换．下面的定理给出不可约的条件．

定理8.4.4若是不可约的上Hessenberg矩阵，且和均非的特征值，则也是不可约上Hessenberg矩阵．

［证明］用反证法．记，并假定存在使得，而．比较等式的两边矩阵的前列，得

由此可得

(8.4.26)

其中

由得

将代入(8.4.26)，并注意到也是的多项式，就有

(8.4.27)

其中

记，并注意到是不可约的上Hessenberg矩阵，直接计算可知的第个分量为

这也就是说(8.4.27)有非零解，而这与和均非的特征值蕴含着非奇异矛盾．

定理证毕

基于定理8.4.2和6.4.4，我们可以从另外的途径来实现到的变换．首先，我们从(8.4.21)知，的第一列与的第一列共线（其实的第一列就相当于由的第一列单位化而得到的）．而由(8.4.25)容易算出

其中

其次，如果Householder变换使：，从而有，即的第一列就与共线，从而的第一列就可作为的第一列，即，而由关于Householder变换的理论知，可以按如下方式确定：

，

其中

现令

则我们只要能够找到第一列为的正交矩阵使为上Hessenberg矩阵，那么就是我们希望得到的．由本节所介绍的约化一个矩阵为上Hessenberg矩阵，的方法可知，这是容易办到的．这只需确定个Householder变换使

为上Hessenberg矩阵，则的第一列就为．而且由于所具有的特殊性，实现这一约化过程所需的运算量仅为．

事实上，由于用将相似变换为只改变了的前三行和前三列，故有如下形状

，

仅比上Hessenberg形多三个可能的非零元素＂＋＂．由的这种特殊性易知，用来约化为上Hessenberg形的第一个Householder变换具有如下形状

，

其中为３阶Householder变换，而且具有如下形状

．

一般地，第k次约化所用的Householder变换具有如下形状

其中为３阶Householder变换，，而且具有如下形状

．

因此，最后一次约化所用的Householder变换具有如下形状

其中为２阶Householder变换．

综述上面的讨论，就得到了著名的Francis双重步位移的ＱＲ迭代算法：

算法8.4.2（双重步位移ＱＲ迭代）

m=n-1

s=H(m,m)+H(n,n)

t=H(m,m)H(n,n)-H(m,n)H(n,m)

x=H(1,1)H(1,1)+H(1,2)H(2,1)-sH(1,1)+t

y=H(2,1)(H(1,1)+H(2,2)-s)

z=H(2,1)H(3,2)

for k=0:n-3

q=max{1,k}

H(k+1:k+3,q:n)= H(k+1:k+3,q:n)

R=min{k+4,n}

H(1:r,k+1:k+3)=H(1:r,k+1:k+3)

x=H(k+2,k+1)

y=H(k+3,k+1)

if k<n-3

z=H(k+4,k+1)

end

end

H(n-1:n,n-2:n)= H(n-1:n,n-2:n)

H(1:n,n-1:n)= H(1:n,n-1:n)

该算法的运算量为；如果需要累积正交变换投影还需要增加运算量．

8.5.5 隐式QR算法

前面的讨论已经解决了用QR方法求一个给定的实矩阵的Schur标准形的几个关键的问题．然而，作为一种实用的算法，还需给出一种有效的判定准则，来判定迭代过程中所产生的上Hessenberg矩阵的次对角元何时可以忽略不计．一种简单而实用的准则是，当

时，就将看作零．这样做的理由是，在前面约化为上Hessenberg矩阵时就已经引进了量级为的误差．

综合上面的讨论，就得到如下的隐式QR算法．该算法是计算一个给定的阶实矩阵的实Schur分解：，其中是正交矩阵，为拟上三角矩阵，即对角块为或方阵的块上三角矩阵，而且每个的对角块必为有一对复共轭特征值．

算法8.4.3（计算实矩阵的实Schur分解：隐式QR算法）

输入．

上Hessenberg化：用算法8.4.1计算的上Hessenberg分解，得．

收敛性判定：

(i)把所有满足条件

的置零；

(ii)确定最大的非负整数和最小的非负整数，使

其中为拟上三角形，而为不可约的上Hessenberg形；

（iii）如果，则输出有关信息，结束；否则进行下一步．

QR迭代：对用算法8.4.2迭代一次得

．

计算

，

然后转步（３）．

实际计算的统计表明，这一算法每分离出一个或子矩阵平均约需２次QR迭代．因此，如果只计算特征值，则运算量平均约为；如果都需要，则运算量平均约为．

关闭提示 关闭

确 认取 消