Transcript 第八章

第8章
矩阵特征值计算
1
8.1
特征值性质和估计
8.2
幂法及反幂法
8.3
正交变换与矩阵分解
8.4
QR方法
8.1
特征值性质和估计
8.1.1
特征值问题及其性质
设矩阵 A  R nn ，特征值问题是求   C和非零向量 x  R n ,
使
Ax  x,
（1.1）
其中 x 是矩阵 A 属于特征值  的特征向量.
求 A的特征值问题（1.1）等价于求 A 的特征方程
p( )  det( I  A)  0
（1.2）
的根.
2
定理1
设  为 A  R nn 的特征值, Ax  x，x  0,
则
（1） c 为 cA的特征值（ c 为常数 c  0 );
（2）   为 A  I 的特征值，即 ( A  I )  (   ) x;
（3）  k为 Ak 的特征值；
3
定理2 （1）设 A  R nn 可对角化，即存在非奇异矩阵
P, 使
 1


P 1 AP  



2






n 

的充要条件是 A 具有 n个线性无关的特征向量.
（2） 如果 A 有 m个 (m  n)不同的特征值 1 , 2 , , m ,
则对应的特征向量 x1 , x2 ,, xm线性无关.
4
定理3
设 A  R nn为对称矩阵，则：
（1） A的特征值均为实数；
（2） A有 n个线性无关的特征向量；
（3）存在一个正交矩阵 P使
 1


P 1 AP  



2




,

n 
且 i (i  1,, n)为 A特征值， 而 P  (u1 , u2 ,, un )的列
向量 u j 为 A 的对应于  j 的特征向量.
5
定理4
设 A  R nn 为对称矩阵（其特征值次序记为
1  2  n ), 则
( Ax, x)
1. n 
 1
( x, x )
( Ax, x)
2. 1  maxn
;
xR
( x, x )
x0
记 R( x) 
证明
(对任何非零向量x  R n );
( Ax, x)
n  minn
.
xR
( x, x )
x0
（1.3）
( Ax, x)
, x  0, 称为矩阵 A 的瑞利（Rayleigh）商.
( x, x )
只证 1.
由于 A为实对称矩阵，可将 1 , 2 ,, n 对应的特征
向量 x1 , x2 ,, xn 正交规范化，则有 ( xi , x j )   ij .
6
设 x  0 为 R n中任一向量，则有展开式
1
2
 n 2
x 2    i   0,
 i1

n
x   i xi ,
i 1
于是
n
( Ax, x )

( x, x )

i 1
n
2
i

i 1
i
.
2
i
从而1成立.
结论1说明瑞利商必位于 n 和 1 之间.
7
特征值估计与扰动
8.1.2
定义1
设 A  (aij ) nn .令
n
(1) ri   aij
(i  1,2,  , n);
j 1
j i
(2) 集合 Di  {z z  aii  ri , z  C} .
称复平面上以 aii 为圆心，以 ri 为半径的所有圆盘为 A的格
什戈林（Gerschgorin）圆盘.
8
定理5 （格什戈林圆盘定理） （1） 设 A  (aij ) nn，
则 A的每一个特征值必属于下述某个圆盘之中
n
  aii  ri   aij
(i  1,2,  , n).
（1.4）
j 1
j i
或者说， A的特征值都在复平面上 n 个圆盘的并集中.
（2） 如果 A 有 m个圆盘组成一个连通的并集 S ，
且 S 与余下 n  m个圆盘是分离的，则 S 内恰包含 A的 m个
特征值.
特别地，如果 A的一个圆盘 Di 是与其他圆盘分离的
(即孤立圆盘)，则 Di 中精确地包含 A 的一个特征值.
9
证明
只就(1)给出证明.
Ax  x,
设  为 A 的特征值，即
其中 x  ( x1 , x2 ,, xn )T  0.
记 xk  max xi  x
1i  n

 0, 考虑 Ax  x 的第 k 个方程，
即
n
a
j 1
kj
x j  xk ,
或
n
(  akk ) xk   akj x j ,
j k
于是
n
  akk xk   akj x j  xk
j k
n
a
j k
kj
,
10
即
n
  akk   akj  rk .
j k
这说明， A的每一个特征值必位于 A 的一个圆盘中，
并且相应的特征值  一定位于第 k 个圆盘中.
其中 k 是对应特征向量 x 绝对值最大的分量的下标.
11
利用相似矩阵性质，有时可以获得 A的特征值进一步
的估计，即适当选取非奇异对角阵
D 1
 11






并做相似变换 D
1
 21

 aij j
AD  
 
i






1 
n 


 .
 nn
适当选取  i (i 1,2,,n), 可使某些圆盘半径及连通
性发生变化.
12
例1
估计矩阵
4

A  1
1

1
0
1
0 

1 
 4

特征值的范围.
解
A的3个圆盘为
D1 :   4 1,
D2 :   2,
D3 :   4  2,
由定理5，可知 A的3个特征值位于3个圆盘的并集中，
由于 D1 是孤立圆盘，所以 D1内恰好包含 A的一个特征值
（ 1为实特征值），即 3  1  5
13
A的其他两个特征值 2 , 3包含在 D2 , D3 的并集中.
现选取对角阵
D 1
1




1



0.9 

做相似变换
 4

1
A  A1  D AD   1

 0.9

1
0
0.9
0 
10 
.

9 
4 

A1 的3个圆盘为
E1 :   4 1,
E2 :  
10
,
9
E3 :   4 1.8. 14
这样，3个圆盘都成为了孤立圆盘，每一个圆盘都包
含 A的一个特征值（为实特征值）且有估计
3  1  5,
19
19

 2  ,
9
9
 5.8  3  2.2.
15
下面讨论当 A有扰动时产生的特征值扰动，即 A有微
小变化时特征值的敏感性.
定理6 （Bauer-Fike定理） 设  是 A  E  R nn 的一
个特征值，且 P 1 AP  D  diag (1 , 2 ,, n )则有
min     P 1
 ( A)
p
P
p
E p,
（1.5）
其中  为矩阵 p的范数，p  1,2, .
p
证明 只要考虑    ( A).这时 D  I 非奇异，设 x
是 A  E 对应于  的特征向量，由 ( A  E  I ) x  0 左乘 P 1
可得
16
( D  I )( P 1 x)  ( P 1 EP)( P 1 x),
P 1 x  ( D  I ) 1 ( P 1 EP)( P 1 x),
P 1 x是非零向量.上式两边取范数有
( D  I ) 1 ( P 1 EP)
p
 1.
而对角矩阵 ( D  I ) 1 的范数为
( D  I )
1
 ,
m
1
p
m  min    ,
 ( A)
所以有
P 1
p
E
p
P
p
 m.
这就得到（1.5）式.这时总有  ( A) 中的一个  取到 m值. 17
由定理6可知 P 1 P  cond ( P) 是特征值扰动的放大系
数，但将 A对角化的相似变换矩阵 P不是唯一的，所以取
cond ( P) 的下确界
 ( A)  inf{ cond ( P) P 1 AP  diag (1 , 2 ,, n )},
（1.6）
称为特征值问题的条件数.
只要  ( A)不很大，矩阵微小扰动只带来特征值的微小
扰动.但是 ( A) 难以计算，有时只对一个 P，用 cond ( P)
代替 ( A).
18
特征值问题的条件数和解线性方程组的条件数是两个
不同的概念，对于一个矩阵 A，两者可能一大一小.
关于计算矩阵 A的特征值问题，当 n  2,3时，还可以
按行列式展开的方法求特征方程的根.但当 n 较大时，如果
按展开行列式的方法，首先求出 p ( )的系数，再求 p ( )
的根，工作量就很大，用这种方法求特征值是不切实际的,
需要研究求 A的特征值及特征向量的数值方法.
19
8.2
8.2.1
幂法及反幂法
幂法
幂法是一种计算矩阵主特征值（矩阵按模最大的特征
值）及对应特征向量的迭代方法，特别适用于大型稀疏矩
阵.
反幂法是计算海森伯格阵或三对角阵的对应一个给定
近似特征值的特征向量的有效方法之一.
设实矩阵 A  (aij ) nn 有一个完全的特征向量组，其特
征值为 1 , 2 ,, n ，相应的特征向量为 x1 , x2 ,, xn .
20
已知 A的主特征值是实根，且满足条件
1  2  3    n ,
（2.1）
现讨论求 1及 x1的方法.
幂法的基本思想是任取一个非零的初始向量 v0，由矩
阵 A构造一向量序列
 v1  Av0 ,

2
v

Av

A
v0 ,

1
 2



k 1
 vk 1  Avk  A v0 ,



（2.2）
21
称为迭代向量.
由假设，v0 可表示为
v0  1 x1   2 x2     n xn
(设1  0), （2.3）
于是
k
k
k
vk  Avk 1  Ak v0 11 x1  2 2 x2  n n xn
n
 1k [1 x1   i (i / 1 ) k xi ]  1k (1 x1   k ),
i2
其中
n
 k    i (i / 1 ) k xi .
i 2
22
由假设
i / 1 1 (i  2,3,,n),
故
lim  k  0,
k 
从而
lim
k 
vk

k
1
这说明序列
 1 x1.
vk
越来越接近 A的对应于 1的特征向量，

或者说当 k 充分大时
k
1
（2.4）
23
vk  11k x1 ,
（2.5）
即迭代向量 vk为 1的特征向量的近似向量（除一个因子外）.
再考虑主特征值 1 的计算，用 (vk )i表示 vk 的第 i个
分量，则
1 ( x1 ) i  ( k 1 ) i 
(vk 1 ) i
 1 
,
(vk ) i
 1 ( x1 ) i  ( k ) i 
（2.6）
故
(vk 1 ) i
 1 ,
k  (v )
k i
lim
（2.7）
也就是说两相邻迭代向量分量的比值收敛到主特征值.
24
这种由已知非零向量 v0 及矩阵 A 的乘幂 Ak构造向量
序列 {vk } 以计算 A的主特征值 1及相应特征向量的方法
称为幂法.
由（2.6）式知， (vk 1 ) i  1 的收敛速度由比值
( vk ) i
2
来确定，r 越小收敛越快，但当 r  2  1 时
1
1
收敛可能就很慢.
1 ( x1 ) i  ( k 1 ) i 
(vk 1 ) i
（2.6）
 1 
,
(vk ) i
 1 ( x1 ) i  ( k ) i 
r
25
定理7
设 A  R nn 有 n个线性无关的特征向量，主
特征值  满足
1
1  2  3    n ,
则对任何非零初始向量 v (1  0) ,(2.4),(2.7)式成立.
即
lim
k 
vk

k
1
 1 x1.
(vk 1 ) i
lim
 1 ,
k  (v )
k i
（2.4）
（2.7）
26
如果 A的主特征值为实的重根，即 1  2    r ，
且
r  r 1    n ,
又设 A有 n个线性无关的特征向量，1 对应的 r 个线性无
关特征向量为 x1 , x2 ,, xr ，则由（2.2）式
r
n
vk  A v0   { i xi   i (i / 1 ) k xi },
ir 1
 v1  Av0 , i1

r
r
2
v
vk 2 
 Av
A (v设

1
0,
lim
x
 i（2.2）
xi  0).



i
i
k  k
 1
i1 
i1

k 1

Av

A
v0 ,
这说明当
的主特征值是实的重根时，定理7的结论还
 vA
k 1
k

27


是正确的.
k
k
1
应用幂法计算 A的主特征值 1 及对应的特征向量时，
如果 1  1（或 1  1 ），迭代向量 vk 的各个不等于零
的分量将随 k  而趋向于无穷（或趋于零）.
这样在计算机实现时就可能“溢出”.
为了克服这个缺点，就需要将迭代向量加以规范化.
设有一向量 v  0，将其规范化得到向量
v
u
,
max( v)
其中 max( v)表示向量 v 的绝对值最大的分量， 即如果有
28
vi0  max vi ,
1i  n
则 max( v)  vi ，且 i0 为所有绝对值最大的分量中的最小
o
下标.
主特征值为单特征值的条件下幂法可这样进行：
任取一初始向量v0  0 (1  0) ，构造向量序列 max( v)
v1  Au0  Av0 ,
u1 
A2 v0
v2
u2 

,
2
max( v2 ) max( A v0 )
Av0
v1

,
max( v1 ) max( Av0 )
A2 v0
v2  Au1 
,
max( Av0 )
29
 
Ak v0
vk 
,
k 1
max( A v0 )
Ak v0
uk 
,
k
max( A v0 )
由（2.3）式
k
n




i
k
k
k
A v0   i i xi  1 1 x1   i   xi ,
i1


i2
 1 
n
（2.8）
k
（2.3）
n 0),

v0 1 x1  2 x2  n xkn ( 设1 
 i 

1 1 x1   i 
xi 

k

i 2

 1 
A v0


uk 
max ( Ak v0 )

 
n
 i
k
max 1 1 x1   i 

 
i 2
 1
 




 xi  



k
30
n

 i 

xi

 1 
1 x1   i 
i 2
k
n

 i
max 1 x1   i 


i 2
 1

k



 xi 




x1
( k  ).
max ( x1 )
这说明规范化向量序列收敛到主特征值对应的特征向量.
同理，可得到
k
n


 i 
k

1 1 x1   i 
xi 


i 2

 1 


vk 
,
k 1
n


 i 
k 1
max 1 1 x1   i 
 
 xi 

i 2

 1


31
k
n


 i 

1 max 1 x1   i 
 xi 


i 2

 1

  ,
max( vk ) 
1
k 1
n


 i 
max 1 x1   i 
 
 xi 

i 2

 1


( k  ),
收敛速度由比值 r 
2
确定.
1
32
定理8
设 A  R nn 有 n个线性无关的特征向量，主
特征值 1 满足 1  2  3    n ,则对任意非零
初始向量 v0  u0 (1  0)，按下述方法构造的向量序列
{uk },{vk } :







则有
v0  u0  0,
vk  Auk 1 ,
(k 1,2,),
（2.9）
k  max( vk ),
uk  vk /  k .
(1) lim u k 
k 
x1
,
max ( x1 )
(2) lim  k  1.
k 
33
例2
用幂法计算
 1.0

A   1.0
 0.5

1.0
1.0
0.25
0.5 

0.25 
2.0 
的主特征值和相应的特征向量. 计算过程为
v0  u0  (1,1,1)T ,
u1  v1 / 1.
v1  Au0 ,
v2  Au1 ,
1  max( v1 ),
2  max( v2 ),
u2  v2 / 2 .

结果如表8-1.
34
表8  1
k
T
u（
规范化向量）
k
1
max( vk )
0
(1
1)
1
5
(0.9091
(0.7651
0.8182
0.6674
1)
1)
2.7500000
2.5587918
10
(0.7494
0.6508
1)
2.5380029
15
(0.7483
0.6497
1)
2.5366256
16
17
(0.7483
(0.7482
0.6497
0.6497
1)
1)
2.5365840
2.5365598
18
(0.7482
0.6497
1)
2.5365456
19
(0.7482
0.6497
1)
2.5365374
20
(0.7482
0.6497
1)
2.5365323
35
表8-1的结果是用8位浮点数字进行运算得到的，u k的
分量值是舍入值.
于是得到
1  2.5365323
T
(
0
.
7482
,
0
.
6497
,
1
)
.
及相应的特征向量
1和相应的特征向量的真值（8位数字）为
1  2.5365258,
~
x1  (0.74822116 , 0.64966116, 1)T .
36
8.2.2
加速方法
原点平移法
由前面讨论，应用幂法计算 A的主特征值的收敛速度
主要由比值 r  2 来决定，但当 r 接近于1时，收敛可能
1
很慢.
一个补救的办法是采用加速收敛的方法.
引进矩阵
B  A  pI ,
其中 p为选择参数.
37
设 A 的特征值为 1 , 2 ,, n ，则 B的相应特征值为
1  p,2  p,,n  p,
而且 A, B 的特征向量相同.
如果要计算 A的主特征值 1，就要适当选择 p , 使 1  p
仍然是 B 的主特征值，且使
2  p

 2.
1  p
1
对 B应用幂法，使得在计算 B的主特征值 1  p的过程中
得到加速.
这种方法通常称为原点平移法.
38
例3
设 A  R 44 有特征值
 j  15  j
( j  1,2,3,4),
比值 r  2 / 1  0.9. 作变换
B  A  pI
( p  12),
则 B的特征值为
1  2, 2 1, 3  0, 4  1.
应用幂法计算 B的主特征值 1的收敛速度的比值为
2
2  p
2
1

 0.9.
 
1
1  p
1
2
39
选择有利的 p值，虽然能够使幂法得到加速，但问题
在于如何选择适当的参数 p .
设 A的特征值满足
1  2    n1  n ,
（2.10）
则不管 p如何， B  A  pI 的主特征值为 1  p或 n  p .
当希望计算 1及 x1时，首先应选择 p 使
1  p  n  p ,
且使收敛速度的比值

 2  p n  p
  max 
,

 1  p 1  p


  min .


40
2  n

p

p
n
2
显然，当
, 即p
 p *时

2
1  p
1  p
 为最小，这时收敛速度的比值为
2  n
  p*
2  p *

.
 n
1  p *
1  p * 21  2  n
当 A的特征值满足（2.10）且 2 , n 能初步估计时，
就能确定 p *的近似值.
当希望计算 n 时，应选择
1  2    n1  n ,
1  n 1
p
 p*,
2
使得应用幂法计算 n 得到加速.
（2.10）
41
例4
计算矩阵
 1.0

A   1.0
 0.5

1.0
1.0
0.25
0.5 

0.25 
2.0 
的主特征值.
作变换 B  A  pI , 取 p  0.75，则
 0.25

B   1.0
 0.5

1.0
0.25
0.25
0.5 

0.25 .
1.25 
对 B应用幂法，计算结果如表8-2.
42
表8  2
k
T
u（
规范化向量）
k
max( vk )
0
(1
1
1)
5
(0.7516
0.6522
1)
1.7914011
6
(0.7491
0.6511 1)
1.7888443
7
(0.7488
0.6501 1)
1.7873300
8
(0.7484
0.6499
1)
1.7869152
9
(0.7483
0.6497
1)
1.7866587
10
(0.7482
0.6497
1)
1.7865914
由此得 B的主特征值为 1 1.7865914
, A的主特征值 1为
1  1  0.75  2.5365914,
43
与例2结果比较，上述结果比例3迭代15次还好.
若迭代15次，1  1.7865258
（相应的 1  2.5365258
).
原点位移的加速方法，是一个矩阵变换方法.
这种变换容易计算，又不破坏矩阵 A的稀疏性，但 p
的选择依赖于对 A的特征值分布的大致了解.
44
瑞利商加速
定理9
设 A  R nn为对称矩阵，特征值满足
1  2  3    n ,
对应的特征向量满足 ( xi , x j )   ij ，应用幂法计算 A 的主
特征值 1，n则规范化向量u k的瑞利商给出
k1的较好的近似
n

 i 
k
k
k
A v0   i i xi  1 1 x1   i  2 k  xi , （2.8）
i2   1
i1 ( Au , u )


k
k
2



 1  O 
.
  1  
(uk , uk )


证明 由（2.8）式及
Ak u0
uk 
,
k
max( A u0 )
Ak 1u0
vk 1  Auk 
,
k
max( A u0 )
45
得
n
( Auk ,uk ) ( Ak 1u0 , Ak u0 )


k
k
(u k ,uk )
( A u0 , A u0 )

j 1
n

j 1
 
2
 1  O 

  1





2k

.


2
j
2j k 1
2
j
2j k
（2.11）
46
8.2.3
反幂法
反幂法用来计算矩阵按模最小的特征值及其特征向量，
也可用来计算对应于一个给定近似特征值的特征向量.
设 A  R nn 为非奇异矩阵， A 的特征值次序记为
1  2  3    n ,
相应的特征向量为 x1 , x2 ,, xn，则 A1 的特征值为
1
n

1
n 1

1
1
对应的特征向量为 xn , xn 1 ,, x1 .
,
47
因此计算 A 的按模最小的特征值n的问题就是计算 A1
的按模最大的特征值的问题.
对于 A1应用幂法迭代（称为反幂法），可求得矩阵 A1
的主特征值 1 / n ，从而求得 A的按模最小的特征值 n .
反幂法迭代公式为：
任取初始向量 v0  u0  0 , 构造向量序列
vk  A1u k 1

vk
u 
k

max( vk )

( k  1,2, ).
迭代向量 vk 可以通过解方程组 Avk  uk 1 求得.
48
定理10
设 A为非奇异矩阵且有 n 个线性无关的特征
向量，其对应的特征值满足
1  2    n1  n  0,
则对任何初始非零向量 u0 ( n  0) ，由反幂法构造的向量
序列 {vk },{uk }满足
(1) lim u k 
k 
xk
,
max ( xk )
( 2) lim max( vk ) 
k 
1
n
.
收敛速度的比值为 r  n .
n 1
49
反幂法中也可以用原点平移法来加速迭代过程或求其
他特征值及特征向量.
如果矩阵 ( A  pI ) 1 存在，其特征值为
1
1
1
,
,,
,
1  p 2  p
n  p
对应的特征向量仍然是 x1 , x2 ,, xn .
对矩阵 ( A  pI ) 1应用幂法，得到反幂法的迭代公式
 u0  v0  0, 初始向量


1
 vk  ( A  pI ) uk 1


 uk  vk / max( vk ) (k 1,2,).
（2.12）
50
如果 p是 A的特征值  j 的一个近似值，且设  j 与其
他特征值是分离的，即
 j  p  i  p
就是说
(i  j ),
1
是 ( A  pI ) 1的主特征值，
j  p
这时也可用反幂法计算特征值及特征向量.
51
设 A  R nn 有 n个线性无关的特征向量 x1 , x2 ,, xn，
则
n
u0   i xi ( j  0),
i1
( A  pI )  k u0
vk 
,
( k 1)
max (( A  pI )
u0 )
( A  pI )  k u0
uk 
,
k
max (( A  pI ) u0 )
其中
n
( A  pI )  k u0    i (i  p )  k xi ,
i 1
52
同理可得：
定理11
设 A  R nn 有 n个线性无关的特征向量，A
的特征值及对应的特征向量分别记为 i及 xi (i  1,2,, n),
而 p为  j 的近似值，( A  pI ) 1存在，且
 j  p  i  p
(i  j ).
则对任意的非零初始向量 u0 ( j  0) ，由反幂法迭代公式
（2.12）构造的向量序列 {vk },{uk } 满足
x, j 初始向量
u

v

0

0
0
(1) lim u k 
,

k 
max ( xj1)

vk  ( A  pI ) uk 1
（2.12）


1

(2
, k 1,2,).
) lim
k ) v )
ukmax(
 vk /vmax(
(
k
k
 p
j
53
即
1
p
 j
max (vk )
(当k  ),
且收敛速度由比值 r   j  p / min i  p 确定.
i j
由该定理知，对 A  pI (其中 p   j ) 应用反幂法，
可用来计算特征向量 x j .
只要选择的 p是  j 的一个较好的近似且特征值分离情
况较好，一般 r 很小，常常只要迭代一二次就可完成特征
向量的计算.
54
反幂法迭代公式中的 vk是通过解方程组
( A  pI )vk  uk 1
求得的.
为了节省工作量，可以先将 A  pI 进行三角分解
P( A  pI )  LU ,
其中 P 为某个排列阵.
于是求 vk 相当于解两个三角形方程组
Lyk  Puk 1 ,
Uvk  yk .
55
可以按下述方法选择 u0 ：选 u0 使
Uv1  L1 Pu0  (1,1, ,1)T
（2.13）
用回代求解（2.13）即得 v1，然后再按公式（2.12）进行迭
代.
反幂法计算公式
初始向量
 u0  v0  0,
1. 分解计算


vk  ( A  pI ) 1 uk 1
（2.12）

P( A 
 pI )  LU , 且保存 L, U及P信息.


uk  vk / max( vk ) (k 1,2,).
2. 反幂法迭代
(1)
解Uv1  (1,1,,1)T 求v1
56
1  max( v1 ), u1  v1 / 1
(2) k  2,3,
1)
解 Lyk  Puk 1 求 yk
解 Uvk  yk 求 vk
2) k  max( vk )
3)
计算 uk  vk / k
57
例5
用反幂法求
2

A  1
0

1
3
1
0

1
4 
的对应于计算特征值   1.2679
（精确特征值为 3  3  3 )
的特征向量（用5位浮点数进行运算）.
解
用部分选主元的三角分解将 A  pI 分解为
P( A  pI )  LU ,
其中 p  1.2679
58
 1

L 
0
 0.7321

1

U  0
0

0

0 ,
1 
0
1
 0.26807
0.7321
1
0
0

P  0
1

1
0
0
0

1 .
0 


2.7321
,
0.2940510 3 
1
由 Uv1  (1,1,1)T，得
v1  (12692,  9290.3, 3400.8)T ,
u1  (1,  0.73198, 0.26795)T ,
由 LUv2  Pu1，得
59
v2  (20404,  14937, 5467.4)T ,
u2  (1,  0.73206, 0.26796)T ,
3 对应的特征向量是
x3  (1, 1  3, 2  3 )T  (1,  0.73205, 0.26795)T ,
由此看出 u2 是 x3 的相当好的近似.
特征值 3  1.2679  1 / 2  1.26794901
，3 的真值为
3  3  3  1.26794912.
60
8.3
正交变换与矩阵分解
61
8.3.1 豪斯霍尔德变换
定义2
设向量 w  R n 且 wT w  1 ，
H ( w)  I  2wwT
为初等反射阵(或称为豪斯霍尔德变换).
如果记 w  (1 , 2 , , n )T , 则
 1 212

  221
H ( w)  


  2 
n 1

 212
1 222



 2n2

 21n 

 22n 
. （3.1）



2 
1 2n 
62
定理12
设有初等反射阵 H  I  2wwT ,其中
wT w 1, 则
(1) H 是对称矩阵，即 H T  H .
(2) H 是正交矩阵，即 H 1  H .
(3) 设 A为对称矩阵，那么 A1  H 1 AH  HAH 亦是对
称矩阵.
证明
只证 H的正交性，其他都可通过验证得到.
H T H  H 2  ( I  2wwT )( I  2wwT )
 I  4wwT  4w( wT w) wT  I .
63
设向量 u  0 ,
H  I 2
uu T
u
2
2
是一个初等反射阵.
初等反射阵的几何意义.
考虑以 w为法向量且过原点 O的超平面 S : wT x  0 .
设任意向量 v  R n ， 则 v  x  y , 其中 x  S , y  S .
于是
Hx  ( I  2wwT ) x
 x  2wwT x  x.
( wT x  0)
64
对于 y  S  ，Hy  ( I  2wwT ) y  y  2wwT y   y.
从而对任意向量 v  R n ，总有
Hv  x  y  v,
其中 v为 v 关于平面S的镜面反射(见图8-1).
图8-1
65
定理13
设 x, y为两个不相等的 n维向量， x 2  y 2 ,
则存在一个初等反射阵 H，使 Hx  y.
证明
令 w  x  y ，则得到一个初等反射阵
x y
2
H  I  2ww  I  2
T
x y
x y
T
T
(
x

y
),
2
2
而且
Hx  x  2
 x2
x y
x y
T
T
(
x

y
)x
2
2
( x  y )( xT x  y T x)
x y
2
2
.
66
因为
x y
T
T
T

(
x

y
)
(
x

y
)

2
(
x
x

y
x),
2
2
所以
Hx  x  ( x  y )  y.
w是使 Hx  y成立的唯一长度等于1的向量(不计符号).
67
定理14 (约化定理)
设 x  ( x1 , x2 , , xn )T  0 ，
则存在初等反射阵 H ,使 Hx  e，其中
1
H  I   1uu T ,
u  x e1 ,
证明
x
2
  sgn( x1 ) x 2 ,
1
 u
2
2
2
（3.2）
  (  x1 ).
记 y  e1 ,设 x  y ,取    x ，则有
2
 y 2 , 于是由定理13存在 H变换
H  I  2wwT ,
其中 w  x  e1 , 使 Hx  y  e1.
x  e1 2
68
记 u  x  e1  (u1 , u2 ,, un )T . 于是
H  I 2
uu T
u
2
 I   1uu T ,
2
其中 u  ( x1   , x2 , , xn )T ,   1 u 22 .
2
显然

1
u
2
2
2

1
(( x1  ) 2  x22  x22 )  (  x1 ).
2
如果  和 x1 异号，那么计算 x1  时有可能出现两相
近数相减的情况，有效数字可能损失.
69
取  和 x1 有相同的符号，即取
1/ 2
 n 2
  sgn( x1 ) x 2  sgn( x1 )  xi 
 i1

.
在计算  时，为了避免上溢或下溢，将 x
d x
,

x 
x
d
(设d  0)，
则有 H 使 H x   e1 , 其中
 H   I  (  ) 1 u u T ,

2
    / d , u   u / d ,     / d ,
H   H .

70
例6
设 x  (3,5,1,1)T，则 x
u  x  ke1  (9,5,1,1, )T , u
H  I   1uu T
2
2
2
 6 .取 k  6,
 108,  
  27

1   45

54   9

 9

 45
29
5
5
1
u
2
9
5
53
1
2
2
 54,
 9

 5
,

1

53 

可以验证 Hx  (6,0,0,0)T .
71
8.3.2
吉文斯变换
2
设 x, y  R ,
则变换
 y1   cos 

  sin 
y 

 2 
sin   x1 



,



cos   x2 
或 y  Px
是平面上向量的一个旋转变换，其中
 cos 
P( )  
  sin 

sin  

cos  

为正交矩阵.
72
R n中变换： y  Px
其中 x  ( x1 , x2 ,, xn )T , y  ( y1 , y2 , , yn )T , 而
i
1

 

1

cos 



P  P(i, j , )  



 sin 





j
sin 
1

1
cos 













1

 

1
i
（3.3）
j
73
称为 R n中平面 {xi , x j }的旋转变换，也称吉文斯变换.
P  P(i, j ,  )  P(i, j ) 称为平面旋转矩阵.
显然， P(i, j ,  )具有性质：
(1) P与单位阵 I 只是在 (i, i ), (i, j ), ( j , i ), ( j , j ) 位置
元素不一样，其他相同.
(2) P为正交矩阵 ( P 1  PT ).
(3) P(i, j ) A (左乘)只需计算第 i 行与第 j 行元素，
即对 A  (aij ) mn
74
 ail   c

 
 a   
 jl    s
s  ail 





c  a jl 

(l  1,2, , n).
其中 c  cos  , s  sin  .
(4) AP(i, j )(右乘)只需计算第 i 列与第 j 列元素
a ,
li
alj   ali ,
alj
 cs

s

c

(l  1,2, , m).
利用平面旋转变换，可使向量 x 中的指定元素变为零.
75
定理15 (约化定理)
设 x  ( x1 , xi , x j , xn )T ,
其中 xi , x j不全为零，则可选择平面旋转阵 P(i, j ,  ),使
i
j
Px  ( x1 ,  xi,  0,  xn )T
其中 xi  xi2  x 2j ,   arctan( x j / xi ).
证明
取 c  cos   xi / xi, s  sin   x j / xi .
由 P(i, j, ) x  x  ( x1, x2 ,, xi,, xj ,, xn )T ,
利用矩阵乘法，显然有
 xi  cxi  sx j , xj   sxi  cx j ,

 xk  xk (k  i, j ).
76
由 c, s 的取法得
xi 
xi2  x 2j , xj  0.
77
8.3.3
矩阵的QR分解与舒尔分解
定理16 设 A  R nn非奇异，则存在正交矩阵 P 使
PA  R, 其中 R为上三角阵.
证明 先用吉文斯变换给出构造 P的方法.
（1）第1步约化，由设有 j (1  j  n) 使 a j1  0，则可
选择吉文斯变换 P(1, j ) ，将 a j1处的元素化为零.若
a j1  0( j  2, , n) ，则存在 P(1, j ) 使得
 r11


P(1, n)P(1,2) A  



r12

( 2)
a22



an( 22)

r1n 

( 2)
a2 n 
( 2)

A
,



( 2) 
ann 
78
可简记为 P1 A  A( 2 ) ,
其中 P1  P(1, n)  P(1,2).
(2) 第 k 步约化：设上述过程已完成第1步至第 k  1步，
于是有
 r11



Pk 1P2 P1 A  





r12

r1k

r22


r2 k

(k )
akk

(k )
ank



r1n 

r2 n 
 
  A( k ) .
(k )
akn 
 

(k ) 
ann 
由设有 j (k  j  n) 使 a (jkk )  0 ， 若 a (jkk )  0( j  k  1,, n) ，
则可选择吉文斯变换 P(k , j )( j  k  1,, n)，使
79
Pk A( k )  P(k , n)P(k , k 1) A( k )
 Pk Pk 1P1 A  A( k 1) ,
其中 Pk  P(k , n)  P(k , k  1).
(3) 继续上述约化过程，最后则有
Pn 1  P2 P1 A  R
(上三角阵),
令 P  Pn 1  P2 P1 ，它是一个正交阵，有 PA  R.
也可以用豪斯霍尔德变换构造正交阵 P，记 A( 0)  A，
)
它的第一列记为 a1( 0.不妨设
a1( 0 )  0，可按公式（3.2）找
到矩阵 H1  R nn , H1  I  11u1u1T ，使
80
H1a1( 0)   1e1 ,
e1  (1,0,,0)T  R n .
于是
A(1)  H1 A( 0 )  ( H1a1( 0 ) H1a2( 0 )
  1
( 0)
 H1an )  
 0

b (1) 
,
(1) 
A 
其中 A (1)  (a1(1) , a2(1) , an(1)1 )T  R ( n 1)( n 1) .
一般地，设
A( j 1)
 D ( j 1)

 0

B ( j 1) 
,
( j 1) 
A

81
其中 D ( j 1)为 j  1 阶方阵，其对角线以下元素均为0，A ( j 1)
为 n  j  1阶方阵，设其第一列为 a1( j 1) ，可选择 n  j  1
的豪斯霍尔德变换 H j  R ( n j )( n j )，使
H j a1( j 1)   j e1 ,
e1  (1,0,,0)T  R n  j 1.
根据 H j构造 n n阶的变换矩阵 H j为
 I j 1
Hj 
 0

0 
,
Hj

于是有
A( j )  H j A( j 1)
 D( j )

 0

B( j ) 
,
( j) 
A 
82
A( j ) 和 A( j 1)有类似的形式，只是 D ( j )为 j 阶方阵，其对角
线以下元素是0，这样经过 n  1步运算得到
H n 1  H 2 H1 A  A( n 1)  R,
其中 R  A( n 1)为上三角阵，P  H n 1  H 2 H1为正交矩阵.
从而有 PA  R.
83
定理17
(QR分解定理) 设 A  R nn 为非奇异矩阵，
则存在正交矩阵 Q与上三角阵 R ，使 A有分解
A  QR ,
且当 R的对角元为正时，分解是唯一的.
证明
从定理16知，只要令 Q  PT 就有 A  QR，
下面证分解的唯一性，设有两种分解
A  Q1R1  Q2 R2 ,
其中 Q1 , Q2为正交阵， R1 , R2为对角元均为正的上三角阵，
则
84
AT A  R1T Q1T Q1 R1  R1T R1 ,
AT A  R2T Q2T Q2 R2  R2T R2 .
由假设及对称正定矩阵 AT A 的楚列斯基分解的唯一性，
则得 R1  R2 .从而可得 Q1  Q2 .
定理16保证了 A可分解为 A  QR.若 A 非奇异，则 R也
非奇异.如果不规定 R 的对角元为正，则分解不是唯一的.
一般按吉文斯变换或豪斯霍尔德变换方法作出的分解
A  QR， R 的对角元不一定是正的，设上三角矩阵 R  ( rij )，
只要令
rnn
r11
D  diag (
, ,
),
r11
rnn
85
则 Q  QD为正交矩阵，R  D 1 R 为对角元是 rii 的上三
角阵，这样 A  Q R 便是符合定理17的唯一QR分解.
86
87
88
89
例7
设用豪斯霍尔德变换作矩阵 A 的QR分解
2

A  1
1

解
2
1
3
3 

1 .
 1

按（3.2）式找豪斯霍尔德矩阵 H1  R 33.使
 2  *
   
H1  1    0 ,
 1  0
   
则有
  0.816497

H1    0.408248
  0.408248

 0.408248
0.908248
 0.0917517
 0.408248 

 0.0917517 ,
0.908248 
 90
  2.44949

H1 A  
0

0

0
1.44949
3.44949
 2.44949 

 0.224745 .
 2.22474 

再找 H 2  R 22 ，使 H 2 (1.44949,3.44949)T  (*,0)T，
得
1
H2  
0

1
0  

 0

H2  
0
  2.44949

H 2 ( H1 A)  
0

0

0
 0.387392
 0.921915
0
 3.74166
0


 0.921915 ,
0.387392 

0
 2.44949 

2.13809 .
 0.654654 

91
这时一个上三角阵，但对角元皆为负数，只要令 D   I ，
则有 R   H 2 H1 A 是对角元为正的上三角阵.取
 0.816497

T
Q  ( H 2 H1 )   0.408248
 0.408248

 0.534522
0.267261
0.801783
 0.218218 

0.872872 ,
 0.436436 

则得 A  QR.
92
除了QR分解，矩阵的舒尔（Schur）分解也是重要的工
具，它解决矩阵 A  R nn可约化到什么程度的问题，对复
矩阵 A  C nn,则存在酉矩阵 U，使 U T AU为一个上三角
矩阵 R，其对角元素就是 A 的特征值，A  URU H 称 A 的
舒尔分解，对于实矩阵 A，其特征值可能有复数，不能用
正交相似变换约化为上三角阵，但它可以约化为以下形式.
93
定理18 (实舒尔分解）设 A  R nn ，则存在正交矩阵
Q，使
 R11


Q T AQ  



R12

R22


R1m 

R2 m 
,
 

Rmm 

（3.4）
其中对角块 Rii (i  1,2,, m)为一阶或二阶方阵.
且每个一阶 Rii是 A的实特征值，每个二阶对角块 R jj
的两个特征值是 A 的两个共轭复特征值.
94
8.3.4
用正交相似变换约化一般矩阵为
上海森柏格阵
设 A  (aij )  R nn .
我们的目标是选择初等反射阵 U1 ,U 2 ,,U n 2 ，
使 A经正交相似变换约化为一个上海森伯格阵.
95
（1） 设
 a11

 a21
A


a
 n1
a12

a22


an 2

a1n 

a2 n    a11
c


 1

ann 
其中 c1  (a21, , an1 )T  R n 1，不妨设 c1  0,
(1)

A12
,
(1) 
A22 
否则这一
步不需要约化.
选择初等反射阵
R1  I  11u1u1T
使
R1c1   1e1
96
其中
1/ 2
n



 1  sgn ( a21 )  ai21  ,
 i 2




 u1  c1  1e1 ,


 1  1 ( 1  a21 ).
（3.5）
令
1
U1  


,
R1 
则
97
 a11
A2 U1 A1U1  
 R1c1
 a11

  1
 0

 

 0
( 2)
 A11
 
 0 c2
(1)
A12
R1 

(1)
R1 A22 R1 
( 2)
a12
( 2)
a13

( 2)
a22
( 2)
a32
( 2)
a23
( 2)
a33




an( 22)
an( 23)

a1(n2 ) 

( 2)
a2 n 
a3( 2n) 
 
( 2) 
ann

( 2)

A12
,
( 2) 
A22 
( 2)
( 2)
其中 c2  (a32
, , an( 22) )T  R n  2 , A22
 R ( n  2 )( n  2) .
98
(2) 第 k 步约化：
重复上述过程，设对 A已完成第1步，…，第 k  1步
正交相似变换， 即有 Ak  U k 1 Ak 1U k 1 ,
或
且
Ak  U k 1 U1 A1U1 U k 1 ,
(1)
 a11

  1


Ak  





( 2)
a12

a1(,kk11)
a1(kk )
a1(,kk)1

( 2)
a22


a2( k,k11)

 k 1
a2( kk )

(k )
akk
ak( k1),k
a2( k,k)1

ak( k,k)1
ak( k1),k 1



(k )
ank
an( k,k)1



a1(nk ) 

(k )
a2 n 

 
(k )

akn

ak( k1),n 
 

( k )  99
ann 
k
(k )
 A11
 
 0 ck
nk
(k )
 k
A12

(k ) 
A22  n  k ,
(k ) T
其中 ck  (ak( k1),k ,,ank
) R nk ,
(k )
(k )
为 k 阶上海森伯格阵， A22
A11
 R ( n k )( n  k ) .
设 ck  0，于是可选择初等反射阵 Rk 使 Rk ck   k e1 ，
其中, Rk 计算公式为
 k  sgn (a

(k )
)  aik
 ik 1
n
(k )
k 1,k


2
1/ 2



,
100
uk  ck  k e1 ,
 k  k ( k  ak( k1),k ),
（3.6）
Rk  I   k1u k u kT .
令
1
U k  


,
Rk 
则
( k 1)
 A11
Ak 1 U k AkU k  
 0 Rk ck
(k )
A12
Rk 

(k )
Rk A22 Rk 
101
( k 1)
 A11
 
 0 ck 1
( k 1)

A12
,
( k 1) 
A22 
（3.7）
( k 1) 为
其中 A11
k  1阶上海森伯格阵.
(k )
(k )
第 k步约化只需计算 A12
Rk .
Rk及 Rk A22
(k )
当 A为对称阵时，只需计算 Rk A22
Rk .
102
（3） 重复上述过程，则有
U n2 U 2U1 AU1U 2 U n2
 a11

  1











( 2)
a22
 2

( 3)
a33







 n2
an( n1,2n)1
 n1



 
 

 
 
( n1) 
ann

 An1.
总结上述讨论，有
103
定理19
（豪斯霍尔德约化矩阵为上海森伯格阵）设
AR nn , 则存在初等反射阵 U1 ,U 2 ,,U n 2 使
U n2 U 2U1 AU1U 2 U n2  U 0T AU0  H (上海森伯格阵).
本算法约需要
5 3
n 次乘法运算，要明显形成 U 0还需
3
2 3
n 次乘法.
要附加
3
104
例8
用豪斯霍尔德方法将
 4

A  A1   2
 4

3
3
2
 7

2 
7 
矩阵约化为上海森伯格阵.
解 选取初等反射阵 R1 , 使 R1c1   1e1 ， 其中
c1  (2,4)T
(1) 计算 R1 :
  max( 2,4)  4, c1  c1  (0.5, 1)T(规范化 )
105
  1.25 1.118034,
u1  c1 e1  (1.618034, 1)T ,
1  (  0.5) 1.809017,
 1   4.472136,
R1  I  11u1u1T .
则有
R1c1   1e1
(2) 约化计算：
令
1
U1  
0
0
,
R1 
106
则
4


A2 U1 A1U1    4.472136

0

7.602631
7.799999
 0.399999
 0.447214 

 0.400000   H .
2.200000 
如果 A是对称的，则 H  U 0T AU 0也对称，这时 H是
一个对称三对角阵.
107
定理20（豪斯霍尔德约化对称阵为对称三对角阵）设
A  R nn为对称矩阵，则存在初等反射阵 U1 ,U 2 ,,U n 2
使
U n  2 U 2U1 AU1U 2 U n  2
 c1 b1

 b1 c2 b2

 

bn  2




cn 1
bn 1



  C.

bn 1 
cn 
108
证明
由定理17，存在初等反射阵 U1 ,U 2 ,,U n 2 使
U n2 U 2U1 AU1U 2 U n2  H  An1为上海森伯格阵,且
An 1亦是对称阵，因此， An 1 为对称三对角阵.
由上面讨论可知，当 A为对称阵时,由 Ak  Ak 1  U k AkU k
(k )
的一步约化计算中只需计算 R 及 Rk A22
Rk .
k
(k )
又由于 A的对称性，故只需计算 Rk A22
Rk 的对角线以下
元素.
注意到
(k )
(k )
(k )
Rk A22
Rk  ( I   k1uk ukT )( A22
  k1 A22
uk ukT ).
109
引进记号
(k )
rk   k1 A22
uk  R n  k ,
t k  rk 
 k1
2
(ukT rk )uk  R n  k ,
则
(k )
(k )
Rk A22
Rk  A22
 uk t kT  t k ukT
(i  k 1,,n, j  k 1,,i ).
对对称阵 A用初等反射阵正交相似约化为对称三对角
2 3
阵大约需要 n 次乘法.
3
110
8.4
QR 方 法
8.4.1
QR算法
QR方法是一种变换方法，是计算一般矩阵（中小型矩
阵）全部特征值问题的最有效方法之一.
QR方法主要用来计算：
（1）上海森伯格阵的全部特征值问题，
（2）计算对称三对角矩阵的全部特征值问题，且QR
方法具有收敛快，算法稳定等特点.
111
对于一般矩阵 A  R nn（或对称矩阵），则首先用豪斯
霍尔德方法将 A 化为上海森伯格阵 B（或对称三对角阵），
然后再用QR方法计算 B 的全部特征值.
设 A  R nn ，且对 A进行QR分解，即
A  QR ,
其中 R为上三角阵， Q 为正交阵.
于是可得到一新矩阵
B  RQ  QT AQ.
显然， B 是由 A 经过正交相似变换得到，因此 B与 A特征
值相同.
112
再对 B进行QR分解，又可得一新的矩阵，重复这一过
程可得到矩阵序列：
设 A  A1
将 A1 进行QR分解 A1  Q1R1
作矩阵 A2  R1Q1  Q1T A1Q1
求得 Ak后将 Ak 进行分解

Ak  Qk Rk
形成矩阵 Ak 1  Rk Qk  QkT Ak Qk

QR算法，就是利用矩阵的QR分解，按上述递推法则
构造矩阵序列 { Ak } 的过程.
113
只要 A为非奇异矩阵，则由QR算法就完全确定 { Ak } .
定理21 （基本QR方法）设 A  A1  R nn . 构造QR算法：



Ak  Qk Rk (其中QkT Qk  I ,
Rk为上三角阵);
Ak 1  Rk Qk (k 1,2,),
~
~
记Q
，则有

Q
Q

Q
,
R
k
1 2
k
k  Rk  R2 R1
（4.1）
（1） Ak 1 相似于 Ak ，即 Ak 1  QkT Ak Qk ;
（2） Ak 1  (Q1Q2 Qk )T A1 (Q1Q2 Qk )
~
~
 QkT A1Qk ;
~ ~
k
（3） Ak 的QR分解式为 A  Qk Rk .
114
证明 （1）,（2）显然，现证（3）.
~ ~
用归纳法，显然，当 k  1 时有 A1  Q
，
1 R1  Q1 R1
设 Ak 1 有分解式
于是
~ ~
Ak 1  Qk 1 Rk 1 ,
~ ~
Qk Rk  Q1Q2 (Qk Rk )R1
 Q1Q2 Qk 1 Ak Rk 1R1
~
~
 Qk 1 Ak Rk 1
~ ~
 AQk 1Rk 1  Ak .
~T
~
其中利用了 Ak  Q
k 1 AQk 1.
115
由定理17知，将 Ak 进行QR分解，即将 Ak 用正交变换
（左变换）化为上三角矩阵
QkT Ak  Rk ,
T
其中 Qk  Pn 1  P2 P1 , 故
Ak 1  QkT Ak Qk  Pn1P2 P1 Ak P1T P2T PnT1.
这就是说 Ak 1 可由 Ak 按下述方法求得：
（1） 左变换Pn 1  P2 P1 Ak  Rk (上三角阵）；
（2） 右变换 Rk P1T P2T  PnT1  Ak 1.
116
定理22 （QR方法的收敛性）设 A  (aij )  R nn ,
(1) 如果 A的特征值满足： 1  2    n  0；
(2) A 有标准型 A  XDX 1其中 D  diag (1 , 2 ,, n ) ，
且设 X 1有三角分解 X 1  LU ( L 为单位下三角阵，U
为上三角阵），则由QR算法产生的 { Ak } 本质上收敛于上三
角矩阵，即
 1


本质上
Ak  R  



*
2
*

*
（当k  时）

 

n 


117
若记 Ak  (aij(k ) ) ，则
(1) lim aii( k )  i ;
（4.2）
(2) 当i  j时, lim aij( k )  0;
（4.3）
k 
k 
(k )
a
i

j
当
时 ij 极限不一定存在.
118
定理23
如果对称矩阵 A满足定理20的条件，则由QR
算法产生的 { Ak }收敛于对角阵 D  diag (1 , 2 ,, n ) .
关于QR算法收敛性的进一步结果为：
设 A  R nn ，且 A有完备的特征向量集合，如果 A 的
等模特征值中只有实重特征值或多重复的共轭特征值，则
由QR算法产生的 { Ak } 本质收敛于分块上三角矩阵（对角
块为一阶和二阶子块）且对角块中每一个2×2子块给出 A的
一对共轭复特征值，每一个一阶对角子块给出 A的实特征值，
即
119
 1



Ak  






*


m
* *  * *


 
* *  * *
,
B1  * * 
  
Bl 
其中 m  2l  n, Bi 为2×2子块，它给出 A一对共轭特征值.
120
8.4.2
带原点位移的QR方法
(k )
定理22中 lim ann
 n 的速度依赖于比值 rn  n / n1，
k 
当 rn 很小时，收敛较快，
如果 s 为 n 的一个估计，且对 A  sI 运用QR算法，
则 (n, n  1) 元素将以收敛因子 (n  s) /( n 1  s) 线性收
敛于零，(n, n) 元素将比在基本算法中收敛更快.
为了加速收敛，选择数列 {sk }，按下述方法构造矩阵
序列 { Ak } ，称为带原点位移的QR算法.
121
设 A  A1  R nn
对 A1  s1I 进行QR分解 A1  s1I  Q1R1
形成矩阵 A2  R1Q1  s1I  Q1T ( A1  s1I )Q1  s1I
 Q1T A1Q1
求得 Ak 后，将 Ak  sk I 进行QR分解
Ak  sk I  Qk Rk , k  3,4,
（4.4）
形成矩阵
Ak 1  Rk Qk  sk I  QkT Ak Qk
（4.5）
122
~
~
~T ~ ，
如果令 Q
，则有
Ak 1  Qk AQk
k  Q1Q2 Qk , Rk  Rk R2 R1
并且矩阵 ( A  s1I )( A  s2 I )( A  sn I )   ( A)有QR分解式
~ ~
 ( A)  Qk Rk .
在带位移QR方法中，每步并不需要形成 Q 和 R ，可按
下面的方法计算：
首先用正交变换（左变换）将 Ak  sk I 化为上三角阵，
即
Pn 1  P2 P1 ( Ak  sk I )  Rk
当 A 为上海森伯格阵或对称三对角阵时，Pi 可为平面旋转阵，
123
则
Ak 1  Pn 1  P2 P1 ( Ak  sk I ) P1T P2T  PnT1  sk I .
下面考虑用QR方法计算上海森伯格阵的特征值.
设 B 为上海森伯格阵，即
 b11

 b21
B



b12

b22



bn ,n 1
b1n 

b2 n 
.



bnn 
如果 bi 1,i  0(i  1,2, , n  1)，则称 B 为不可约上海森伯
格阵.
124
设 A  R nn ，由定理17可选正交阵 U 0使 H  U 0T AU 0
为上海森伯格阵，对 H 应用QR算法.
QR算法：H  H1
对于 k  1,2,
H k  Qk Rk (QR分解)

H k 1  Rk Qk

（4.6）
假设由（4.6）迭代产生的每一个上海森伯格阵H k 都是不可
约的，否则，若在某步有
p
H k 1
 H11
 
 0
n p
H12  p

.
H 22  n  p
125
于是，这个问题就分离为 H11与 H 22两个较小的问题.
当 p  n  1 或 n  2时，有
n 1
H k 1
 H11
 
 0
1
H12  n  1

( k 1) 
hnn  1
或
n2
2
 H11


 0

H12 
n 2
*
*
H k 1
 2 ,
* * 
( k 1)
由此可得到 H 的特征值 n  hnn
，或由 H k 1右下角二阶
阵的特征值求出 n 1 , n.
126
对降阶的 H11，用类似的方法可求出 H 的其余特征值.
实际上，每当 H k 1 的次对角元适当小时，就可进行分离.
例如，如果
h p 1, p   ( h pp  h p 1, p 1 ),
就把 h p 1, p 视为零.
一般取   10 t ，其中 t 是计算中有效数字的位数.
127
8.4.3 用单步QR方法计算上海森伯格阵特征值
上海森伯格阵的单步QR方法：选取 sk 并设
 h11

 h21
H 



h12

h22



hn ,n 1
h1n 

h2 n 
 H1 (设H为不可约阵).



hnn 
对于 k  1,2,  （用位移来加速收敛）
 H k  sk I  Qk Rk

 H k 1  Rk Qk  sk I
由 H k  H k 1 实际计算为
128
（1）左变换：Pn1,n  P23 P12 ( H1  s1I )  R1 (上三角阵).
T
（2）右变换：H 2  R1P12T P23
 PnT1,n  s1I .
其中 Pk ,k 1  P(k , k  1)为平面旋转阵.
（1） 左变换计算
hkk  hkk  s1
(k  1,2,, n),
确定平面旋转阵 P12  P(1,2) 使
 r11

0
P12 ( H 1  s1 I )   0



( 2)
h12
( 2)
h13

( 2)
h22
( 2)
h23

h32
h33



hn , n 1
h1(n2 ) 

( 2)
h2 n 

h3n .
 

hnn 
129
设已完成第1次，  第 k  1次左变换，即有
Pk 1,k P23 P12 ( H1  s1 I ) 
 r11











h1(,2k )1
h1(k2 )



rk 1,k 1
hk( k1),k
(k )
hkk
hk 1,k





hn ,n1
h1(n2 ) 

 

hk( k1),n 
(k )
.
hkn

hk 1,n 
 

hnn 
（4.7）
第 k 次变换的工作就是要确定平面旋转阵 Pk ,k 1  P(k , k  1)，
使 hk 1, k 变为0，且完成第 k 次左变换
130
Pk ,k 1 ( Pk 1,k  P12 ( H1  s1 I ))
这时只需计算（4.7）阵第 k 行及第 k  1行元素.
这是因为平面旋转阵 Pk ,k 1 只改变矩阵的 k 行和 k  1行.
继续这一过程，最后有
Pn 1,n  P23 P12 ( H1  s1 I )  R1 (上三角阵).
（2） 右变换计算
T
H 2  R1P12T P23
 PnT1,n  s1I ,
在第 k 次右变换 ( R1 P12T ) PkT,k 1中，只需计算 R1P12T  PkT1,k
第 k 列及第 k  1列元素.
131
hk ,k  hk ,k  s1
(k  1,2, , n).
最后
H 2  R1P12T PnT1,n  s1I
*

*




* 
* 
 
*
*

*
(为上海森伯格阵).



* 
由上述讨论指出，如果 H  R nn为上海森伯格阵，
则用QR算法产生的 H 2 , H 3 ,, H k , 亦是上海森伯格阵.
即上海森伯格阵在QR变换下形式不变.
132
讨论一个极端的情况
定理24
（2）
设：（1） H  R nn为不可约上海森伯格阵；
 为 H  H1 一个特征值. 则QR方法



H1  I  QR
(QR分解)
H 2  RQ  I
( 2)
 .
中 hn( 2,n)1  0, hnn
证明
记
 r11  r1n 


R
   (上三角阵).

rnn 

133
由设 H1为不可约阵，则上海森伯格阵 H1  I 亦为不可约.
由将上海森伯格阵 H1  I 约化为上三角阵 R 的平面旋转
变换的取法可知
rii  hi 1,i  0
(i  1,2, , n  1),
又因为 QT ( H1  I )  R 为奇异矩阵， 从而得到 rnn  0 .
因此，H 2 的最后一行为 (0,0,,0,  )， 即
( 2)
hn( 2,n)1  0, hnn
 .
(k )
这样在QR方法迭代中，参数 sk 可选为 hnn
，即 H k 的 ( n, n)
元素. 通常可以作为特征值的最好近似.
134
算法1 （上海森伯格矩阵的QR算法）
给定 H  R nn 为上海森伯格阵，本算法计算



H1  sI  Q1 R1
(QR分解) (取s  hnn )
H 2  R1Q1  sI
且 H 2覆盖 H ( H  H1 )
1.
h11  h11  s
2.
对于k  1,2,, n  1
(1) hk 1,k 1  hk 1,k 1  s
(2)
确定P(k , k  1)
使
135
 ck

  sk
(3)
sk  hkk   rkk
  

ck  hk 1,k   0



左变换
对于j  k ,, n
 hkj   ck

  
h

k

1
,
j

   sk
3.
sk  hkj 


ck  hk 1, j 
对于k  1,2,, n  1
(1)
右变换
对于i  1,2,, k  1
136
 ck
（hik , hi ,k 1 ) 
（hik , hi ,k 1 )
 sk
 sk 

ck 
(2) hkk  hkk  s
4. hnn  hnn  s
(k )，反复应用算法3就产生正
如果用不同的位移 sk  hnn
交相似的上海森伯格阵序列 H1 , H 2 ,, H k , .
当 hn( k, n)1充分小时，可将它置为零就得到 H 的近似特征
(k.
)
值 n  hnn
再将矩阵降阶，对较小矩阵连续应用算法.
137
例9
用QR方法计算对称三对角矩阵
2

A  A1   1
0

1
3
1
0

1 .
4 
的全部特征值.
(k )
解 选取 sk  ann
，则 s1  4.
P23 P12 ( A1  s1 I )  R
 2.2361




1.342
1.0954
0.4472 

 0.3651,
0.81650 
138
T
A2  RP12T P23
 s1 I
 1.4000

  0.4899

0

0.4899
3.2667
0.7454
 1.2915

A3   0.2017

0

0.2017
 1.2737

A4   0.0993

0

0.0993
3.0202
0.2724
2.9943
0.0072


0.7454 .
4.3333 
0


0.2724 ,
4.6884 
0


0.0072 ,
4.7320 
0
139
 1.2694

A5   0.0498

0

0.0498
2.9986
0


0 ,
4.7321
0
 1.2694
~
A5  
 0.0498
0.0498 
.
2.9986 
~
A
现在收缩，继续对 5 的子矩阵 A5  R 22 进行变换，得到
~
~
A5  P12 ( A5  s5 I ) P12T  s5 I
 1.2680
 
5

4

10

 410 5 
,
3.0000 
140
故求得 A 近似特征值为
3  4.7321, 2  3.0000, 1 1.2680.
而 A 的特征值是
3  3  3  4.7321, 2  3.0, 1  3  3  1.2679.
(k )
算法1是在实数中进行选择位移 sk  hnn
,
不能逼近一
个复特征值，所以算法3不能用来计算 H 的复特征值.
141
8.4.4*
双步QR方法（隐式QR方法）
第3节中将 A  R nn 经过正交相似变换化为上海森伯格
矩阵 H，即 U 0T AU 0  H，其中 H不是唯一的.
但是，如果规定了正交矩阵 U 0的第一列，则 U 0和 H除
差±1因子外唯一.
定理25 （隐式Q定理）设A  R nn ,
且：
（1） Q  (q1 , q2 ,, qn ) 及 V  (v1 , v2 ,, vn ) 都是
正交阵，且有 QT AQ  H , V T AV  G都是上海森伯格阵.
（2） H为不可约上海森伯格阵，且 q1  v（即
Q与 V
1
第1列相同）.
142
则：
（1） vi   qi ,且 hi ,i 1  g i ,i 1 (i  2, , n) ;
（2） G  D 1 HD,其中 D  diag (1,1,,1) ，
即 H和 G在 G  D 1 HD意义上“本质上相等”.
算法1不能用来求 H 的一个复特征值.
当 H的按模最小特征值是复数时，位移参数 sk , sk 1 可
取为某步 H k 右下角的二阶矩阵
 hn 1,n 1
G  
 hn ,n 1
的特征值.
hn 1,n 

hnn 
（4.8）
143
当 G 的特征值 s1与 s2为复数时，如果应用算法1就要
引进复数运算，这对于实矩阵 H是不必要的.
在某些条件下，可以用正交相似变换将 H约化为实舒
尔型.
隐式位移的QR方法，即用 s1与 s2作位移连续进行二次
单步的QR迭代，使用复位移，又避免复数运算.
（1） 设 H  H1  R nn 为上海森伯格阵，取共轭复
数 s1 , s2 作两步位移的QR方法，即
H1  s1I  Q1R1 ,
H 2  R1Q1  s1 I  Q1T H1Q1 ,
144
H 2  s2 I  Q2 R2 ,
H 3  R3Q3  s3 I  Q2T Q1T H1Q1Q2  QT H1Q,
其中Q  Q1Q2 ,
R  R2 R1.
（4.9）
显然 M  ( H1  s1I )( H1  s2 I ) 有QR分解
M  QR.
（4.10）
事实上，由(4.9)式并利用
( H 2  s2 I )  Q1T ( H1  s2 I )Q1  Q2 R2
有
M  ( H1  s2 I )Q1R1  (Q1Q2 R2Q1T )Q1 R1
145
 Q1Q2 R2 R1  QR.
且 M 阵为实矩阵，这是因为（即使 G 特征值为复数）
M  H12  ( s1  s2 ) H1  s1s2 I
（4.11）
其中
s1  s2  hn1,n1  hnn  s, s1s2  hn1,n1hnn  hn,n1hn1,n  t
为实数.
于是，（4.10）式为实矩阵 M 的QR分解，并且可以
选取 Q1 和 Q2使 Q  Q1Q2为实的正交阵.
由此得出
M  QR.
（4.10）
146
H 3  (Q1Q2 )T H1 (Q1Q2 )  QT H1Q
是实矩阵.
如果用下述算法就能保证 H 3是实矩阵
(a) 直接形成实矩阵 M  H12  sH1  tI
(b) 计算 M 阵的实QR分解 M  QR
(c) 令 H 3  QT H1Q
但是(a)需要 O(n 3 ) 次乘法运算，不实用.
（2） 根据隐式Q定理，如果按下述算法进行，就有可
能用 O(n 2 ) 次运算来实现从 H1到 H 3的转换.
(a′) 求与 Q 有相同第一列的正交阵 P0
147
(b′) 应用豪斯霍尔德方法将 P0T H1 P0 化为一个上海森
伯格阵，即
Pn  2  P2 P1 ( P0T H1 P0 ) P1 P2  Pn  2  H .
记 Q  P0 P1  Pn 2 , 上式为
QT H1Q  H .
显然， Q  的第一列与 P0 的第一列相同，即 Q 与 Q 第一
列相同（ Qe1  P0e1  Qe1 ）.
若 QT H1Q与 QT H1Q 两者都是不可约上海森伯格阵，
则由隐式Q定理 H 与 H 3 本质上相等.
（3） 如何寻求正交阵 P0 .
148
由于 M  QR (为 M 的QR分解），则
Me1  Q Re1  r11Qe1.
说明 Q 第一列即是 M 第一列的一个倍数.
于是，对 M 阵的第一列（非零）寻求初等反射阵 P0 ,
使
P0 ( Me1 )  r11e1 (其中r11   )
即
Me1  r11P0 e1.
这说明 P0 与 Q 具有相同的第一列.
由于 M  ( H  s1I )( H  s2 I ) ,
则
149
Me1  ( x, y, z,0,,0)T
其中
x  (h11  s1 )(h11  s2 )  h12h21  h112 h21  h12h21  sh11  t ,
y  (h11  s2 )h21  (h22  s1 )h21  h21 (h11  h22  s),
（4.12）
z  h21h32 .
双步QR方法：设 H  H1  R nn 为不可约上海森伯格阵.
(a) 计算 M 阵的第一列. 即按（4.12）式计算
Me1  ( x, y, z,0,,0)T ;
(b) 确定初等反射阵 P0使
P0 ( Me1 )   e1 ,
150
即确定初等反射阵 R0  R 33, 使
 x
 
R0  y    e1 ;
z
 
 R0
P0  

 3

;
I n3
(c) 计算初等反射阵 P1 , P2 ,, Pn 2, 使
Pn  2  P2 P1 ( P0 H1 P0 ) P1 P2  Pn  2  H 
为上海森伯格阵.
则 Q  Q1Q2与 Q  P0 P1  Pn 2 第一列相同,
且 H   H 3.
151
这样上面的算法就完成了从 H1到 H 3的变换，但没有
明显的应用到位移 s1 和 s2 .
152
153