Transcript 第七章对称特征值

第七章 对称特征值问题的计算方法
基本性质
对称QR方法
Jacobi方法
7.1 基本性质
对称矩阵 特征值均为实数，而且其特征向量可以
构成Rn的一组标准正交基，即有
定理1(谱分解定理）若A是n阶对称矩阵，则存
在正交矩阵Q使得
T
Q AQ    diag(1, , n )
定理2（极小极大定理）设A是n阶对称矩阵，并
假定A的特征值为λ1≥… ≥ λn,则有
uT Au
uT Au
i  maxn min T  min
max T .
n
 n i 1 0u u u
 i 0u u u
n

其中 k 表示 Rn 中所有k维子空间的全体。
定理3 设n阶对称矩阵A和B的特征值分别为
1   n 和 1   n
则有
i  i  A  B 2 , i  1, 2, , n.
该定理表明对称矩阵的特征值总是十分良态的。
定理4 设A和A+E是两个n阶实对称矩阵，并假定q1是A的一
个单位特征向量，Q=[q1, Q2]是n阶正交矩阵，QTAQ和QTEQ分
块如下
T

Q AQ  
0
T
若
0 T

, Q EQ  

D2 
e
e 

E22 
d
d  min     0, E 2  ,
 ( D2 )
4
则存在A+E的一个单位特征向量 q1 使得
T T 2
1 1
sin   1  q q
4
4
 e 2  E 2.
d
d
  arccos q1T q1 .
定理4表明，特征向量的敏感性依赖于对应特征值与
其他特征值之间的分离程序。
定理5 设A是m×n阶实矩阵，则存在正交矩阵U和V，
使得
 r 0 
U AV  

0
0


,  r ),  1   2    r  0.
T
其中  r  diag( 1 ,
注：定理5给出的是矩阵A的奇异值分解，数
1   2 
  r   r 1 
 n  0
称为A的奇异值。V的列向量称为A的右奇异向量；U
的列向量称为A的左奇异向量。
定理6 设A和B均为m ×n实矩阵，并假定它们有奇异值分别为
则有
1 
  n 和 1 
n
 i   i  A  B 2 , i  1, 2, , n.
该定理表明矩阵的奇异值也是十分良态的。
作业：习题7
2，3，4
7.2 对称QR方法
•实对称矩阵的三对角化
•隐式对称QR迭代
•隐式对称QR算法
7.2.1 三对角化
对称矩阵的三对角化就是把一个对称矩阵经一
系列正交约化使其成为对称的三对角矩阵。
显然，把第6章“矩阵的上Hessenberg化”算法
应用于实对称矩阵A，则由对称性可知，所得到的上
H矩阵就是对称的三对角矩阵。
注意到A的对称性，上Hessenberg化算法可以改
进，以减少运算量。这就是本节要介绍的实对称矩
阵的三对角化算法。即经n-2次Householder约化，把
A化成三对角矩阵T
H n2
H 2 H1 AH1 H 2
H n2  T .
第1步：把矩阵A分块写成
A  A(0)
 a11(0)

 v
vT 

A2 
依v构造n-1阶Householder变换 H1  I n1   vvT ,并令
1 0 
(0)
H1  
(
H
v

sgn(
a
1
21 ) v 2 )

0 H1 
用H1约化A(0)得到
(0)
1
0

a


11
A(1)  H1 A(0) H1  

0
H
1  v

vT  1 0   a11(0)
v T H1 




A2  0 H1   H1v H1 A2 H1 
显然，第1步约化的主要工作量是计算 H1 A2 H1 ，设
则有
1
u   A2 v, w  u   (v1T u )v,
2
H1 A2 H1  A2  vwT  wvT .
1
H1 A2 H1  A2  vwT  wvT , u   A2 v, w  u   (v1T u )v.
2
上式的推导过程：
H1 A2 H1  A2  vwT  wvT 的计算过程为：
由对称性，
T
(22 H
: n1 )( IA(2vv: nT ,)2A:2 (nI)v(2 :vvn,1)
)
Hu1 A
T
T
) T )(
u (2
) A
(
v
(2
:
n
)
u (2 : n) / 2)v(2 : n)
)
vv
A2 :n
( I(2 :nvv
w
2
n T   vvT A2   2 vvT A2vvT
 2A2:vv
A2 i 
 for
n  1  2 vT A vvT   A v  1  2vvT A v vT
 v j  vi T: A
 A2 for
  2 2
2 
2
2
2
A(i, j )  A(i, j )  v(i ) wT ( j )  w(i )v( j )
T
T
2
T
1
v
v
A
vv


v
A


v
A
 A2  vA(j ,Ai )2 vA12(i, 2jvv



2
2
2
) 2
 v  u  12  (vT u ) v    u  12  (vT u ) v  vT
 A2 endfor
endfor
 A2  vwT  wvT
T
算法7.2.1(实对称矩阵三对角化：Householder变换法)
Tridiagonal( A(1: n,1: n))
for k  1: n  2
[v,  ]  House( A( k  1: n, k ))
u (2 : n)   A(2 : n, 2 : n)v(2 : n,1)
w(2 : n)  u (2 : n)   (v(2 : n)T u (2 : n) / 2)v(2 : n)
for i  2 : n
for j  i : n
A(i, j )  A(i, j )  v(i ) w( j )  w(i )v( j )
A( j , i )  A(i, j )
endfor
endfor
endfor
4n 3
时间复杂度： .
3
7.2.2 隐式对称QR迭代
实现了对称矩阵的三对角化之后，下面的任务就
是选取适当的位移进行QR迭代。由于实对称矩阵的
特征值均为实数，只需进行带原点单步位移即可。
带原点位移的QR迭代格式：
Tk  k I  Qk Rk 
 k  0,1,
Tk 1  Rk Qk  k I 
; T0  T .
由于QR迭代保持上Hessenberg形和对称性的特点，
则迭代产生的Tk都是对称三对角矩阵。和非对称
QR方法一样，仍然假定迭代中所出现的Tk均是不
可约的，即次对角线元均不为零。
Tk  k I  Qk Rk 
 k  0,1,
Tk 1  Rk Qk  k I 
; T0  T .
（1）位移的选取：
第1种选法：选取 μk=T(n,n).
第2种：取μk为Tk右下角子矩阵
 n1
Tk (n  1: n, n  1: n)  
  n1
 n1 
 n 
的两个特征值中靠近αn的那一个，即取
k   n    sgn( )  2   n21 ,
其中   ( n 1   n ) / 2. 这就是著名的威尔金森(Wilkinson)
位移. Wilkinson证明了这两种位移最终都是三次收敛的，
并且说明了为什么后者比前者好[12].
（2）如何具体实现一次QR迭代
T   I  QR, T  RQ   I
当然，可以利用Givens变换直接实现 T   I 的QR分解，进而
完成一步QR迭代。但是，更漂亮的做法是以隐含的方式来实现
由 T 到T 的变换。
QR迭代实质是用正交变换将T约化为 T ，即T  QT TQ .因此,根
据定理6.4.3，T 本质上是由Q的第1列完全确定的。从利用二元正
交变换实现QR迭代的过程可知，Qe1=G1e1，其中G1=G(1,2,ϴ)，
且ϴ满足
 cos  sin   1    *
  sin  cos       0 .

 1   
因此，可以先对T用G1做一次正交约化，得到
   0
   0

B  G1TG1T  
    


0
0




然后调用算法7.2.1，把B正交约化为三对角阵，即得到 T .
算法7.2.2(带Wilkinson位移的隐式对称QR迭代)
WilkinsonQR( A(1: n,1: n))
  (T (n  1, n  1)  T (n, n)) / 2

  T (n, n)  T (n, n  1) 2   sgn( )  2  T (n, n  1) 2
x  T (1,1)   ; y  T (2,1)
for k  1: n  1
[c, s ]  Givens( x, y )
T  Gk TGkT
if k  n  1
x  T (k  1, k ); y  T (k  2, k )
endif
endfor
时间复杂度：10n.

7.2.3 隐式对称QR算法
算法7.2.3(计算实对称矩阵的谱分解：隐式QR方法)
（1）输入A;
（2）三对角化：T=QTAQ;
（3）收敛性判定：
i）将满足 ti 1,i  ti ,i 1   ti ,i  ti 1,i 1  u 的元素置零。
ii）确定最大的非负整数m和最小的非负整数l,使得
0
T11 0
T   0 T22 0  , 其中 T11 , T22 , T33均为三对角阵，且


( n l  m )( n l  m )
l l
mm
T

R
,
T

R
,
T

R
.
0 T33 
 0
22
11
33
iii）若m=n，输出有关信息，结束; 否则进入下一步.
（4）QR迭代：对T22调用算法7.2.2迭代一次，即
T22  GT22GT , G  G1G2 Gnl m1.
（5）Q  Q diag( I l , G, I m ), 然后返回（3）.
作业：习题7
9，11，13
7.3 Jacobi方法
Jacobi方法是求实对称矩阵全部特征值和特征
向量的最古老的方法之一，是由Jacobi于1846年
提出的。他依据实对称矩阵可通过正交相似变换
约化为对角阵这一特点，用一系列适当选取的平
面旋转变换将一个实对称矩阵逐步约化为对角阵。
尽管该算法收敛速度与QR算法无法相比，但该
算法具有编程简单、并行效率高等特点，近年来
又重新受到人们的重视。
7.3.1 经典Jacobi方法
设A=[αij]是n阶实对称矩阵。Jacobi方法的目标
是将A的非对角“范数”
E ( A) 
n
A F  a 
2
i 1
2
ii
n
n
2
a
 ij
i 1 j 1
j i
逐步约化为零. 所用的基本工具就是Givens正交变换：
G( p, q, )  I  (cos   1)(ep eqT  eq eTp )  sin  (e p eqT  eq eTp ).
其中，1 ≤ p<q ≤n,ep,eq均是单位向量.Jacobi方法每
迭代一步需要完成：
1）确定 p和q；确定旋转角ϴ，使得
c  s   pp  pq   c s    pp


 s c  



  qp  qq    s c   0
0 
.

 qq 
2）用G=G(p,q, ϴ)对A一次正交约化：A:=GTAG.
（1）Jacobi迭代法的正交约化过程：
显然，每迭代一步，倘若将A约化为B. 则A,B仅第p
行（列）与第q行（列）不一样，其余元素均相同。
 pk   kp  c kp  s kq , k  p, q;
且
 kq   qk  s kp  c kq , k  p, q;
 pp  c 2 pp  2sc pq  s 2 qq ,
 qq  s 2 pp  2 sc pq  c 2 qq ,
 pq   qp  (c 2  s 2 ) pq  sc( pp   qq ).
以上计算式子的推导，可借下两式：
 pn  c p1  s q1
c  s   p1


 s c  
 qn   s p1  c q1

  q1
1 p 1q 
 c1 p  s1q

  c s 

 s c   
 c  s
 np  nq  
nq
 np
c pn  s qn 
s pn  c qn 
s1 p  c1q 


s np  c nq 
（2）c和s的计算：从最后一个式子知c和s应满足
(c 2  s 2 ) pq  sc( pp   qq )  0
如果 αpq=0，则只需取c=1,s=0即可。如果αpq ≠ 0，则令
 qq   pp
s

, t  tan   ,
2 pq
c
于是有
解得
t 2  2 t  1  0
t    1 2
取t为绝对值较小者，即
t    sgn( ) 1   
2
sgn( )
  1 2
由于 1  tan 2   sec2   1/ cos 2  , 则有
c
1
1 t
2
, s  tc.
（3）p和q的选择：注意到Jacobi算法的初衷是要把
n
A F    ii2
E ( A) 
2
i 1
逐步约化为零。首先正交约化不改变矩阵的F范数，即
A
 B
F
F
B和A的对角线元素只有两个不同的，而且
2
2
2
2
2
2
 pp
  qq
 2 pq
  pp
 qq2  2 pq
  pp
  qq2
故有
E ( B)  B
2
2
F
n
n
2
2
    A F  ii2  2 pq
 E ( A)2  2 pq
.
i 1
2
ii
2
i 1
欲使 E ( A) 尽量小，只需取p和q，使
 
 pq  max  ij .
1i  j  n
算法7.3.1(对称矩阵特征值问题的Jacobi方法)
Jacobi ( A(1: n,1: n), G (1: n,1: n)  I )
while(true)
A( p, q )  max { A(i, j ) }
1i  j  n
if( A( p, q)   ) stop
  ( A(q, q )  A( p, p )) / (2 A( p, q))
t  sgn( ) / (   1   2 )
c  1/ 1   2 ; s  tc
c  s 
 c s
A(1: n,1 : n)  
A(1: n,1: n) 


s
c

s
c




 c s
G (1: n,1: n)  G (1: n,1: n) 


s
c


endwhile
作业：习题7
13，14