8矩阵特征值的计算
矩阵特征值的计算
§1幂法与逆幂法
1.1幂法
幂法是一种计算矩阵的按模最大的特征值(称此特征值为矩阵的主特征值)及对应的特征向量的迭代方法,特别适用于大型稀疏矩阵。反幂法是计算海森伯格矩阵或三对角矩阵的对应一个给定近似特征值的特征向量的有效方法之一。
设n阶实矩阵A(aij)有n个线性无关的特征向量x1,x2,应的特征值为1,2,
,xn,对
,n,且主特征值1是实数,满足条件 |1||2||3|
|n|,
现讨论求1及x1的方法。
幂法的基本思想是任取一个初始向量z0,由矩阵A构造一向量序列
z1Az02zAzAz021
zAzAk1z
k0k1
称为迭代向量,又假设,z0可表示为
z0a1x1a2x2
于是
anxn (设a10)
zkAzk1a11kx1a22kx2
k
annkxn
k
2nk
xaxa 1[a1xn] 12n
11
则当k时,就有
lim
k
zk
k1
a1x1
这说明序列k充分大时
zk
k1
越来越接近A的对应于特征值的特征向量,或者说,当
zka11x1
即迭代向量zk为1的特征向量的近似向量(除一个因子外)。
若用(zk)i表示zk的第i个分量,则
k1
2
a1x1a2x2
(zk1)i11k(zk)i2
a1x1a2x2
1
k1
n
anxn
1
k
annxn
1i
故
lim
(zk1)i
1
k(z)ki
也就是说,两相邻迭代向量的分量的比值收敛到主特征值。 应用幂法计算矩阵A的主特征值1及对应的特征向量时,若迭代向量zk的分量会趋于无穷大;而
11,
11时,迭代向量zk的分量会趋于
零。无论哪种情形,在计算机实现时,都可能 “溢出”,为了克服这个缺陷就需要将迭代向量zk规范化,以免发生上溢或下溢的现象。 幂法的迭代格式为:
yAz;kk1
mk||yk||; k1,2,yzkk.
mk
在这种情况下,当k时,就有
x1
zk
||x1||
m
1k
例1 用幂法计算
332
A1044
362
的主特征值和相应的特征向量。
例1 用幂法计算
10.51.5
A11.50.25
0.50.252.5
的主特征值和相应的特征向量。 解 取z0(1,1,1),计算程序为
A=[3 3 2;10 4 4 ;3 6 2],z=[1 1 1]‘ >> for i=1:20
y=A*z
m=norm(y,inf) z=y/m
end
计算过程如下表所示
T
下述结果是用8位浮点数字进行运算得到的,zk得分量是舍入值,于是得到
1=3.03658162
相应的特征向量(0.7483 0.6497 1.0000),而1与相应的特征向量的真值(16位数字)为
T
1=3.[1**********]718
x1=(0.[1**********]438 0.[1**********]995 1.[1**********]000)T.
例1 用幂法计算
332
A1044
362
的主特征值和相应的特征向量。
解 取z0(0,0,1),计算过程如下表所示 ……
T
1.2幂法的加速方法
有前面的讨论可知,应用幂法计算矩阵A的主特征值的收敛速度主要由比值|
1
|来界定,若比值接近于1,收敛速度就非常慢。这时的补救办法2
是加速收敛的方法。 一、原点平移法
设1,2,征值为
,n是n阶矩阵A的特征值,则n阶矩阵BApI的特
1p,2p,,np
而且A,B有相同的特征向量。可以通过适当选择p,使得
2p2
1p1
从而用幂法计算矩阵BApI的按模最大特征值时,收敛速度就会显著提高。
例如,用幂法计算矩阵
310
A133
034
若用幂法迭代,发现有规律的摆动,我们取p=-4,再用幂法求BApI的按模最大特征值,迭代六次就得到精确到小数点后三位的结果1=5.125。 这种方法称为原点平移法,其关键是p值的选择。若我们能用某种方法估算按次大特征值2时,这时p的选择最好接近2。若对次大特征值的信息无法得到,则应注意在迭代过程中发现次大特征值。例如在计算的前后两次迭代有规律的摆动就可判明1和2接近符号相反,这时p
选为接近
p值,达到加速为止。
二、Rayleigh商加速技术
设A是n阶实对称矩阵,则可选他的特征向量x1,x2,
,xn满足
1 ij
(xi,xj)δij i,j1,2,
0 ij
应用Rayleigh商,则有
,n
(Azk,zk)
(zk,zk)
a
i1ni1
n
22k1ii
2k11
[aai2(i/1)2k1]
21
i2n
n
a
22kii
12k[a12ai2(i/1)2k]
i22k
1O2
1
易见,当A是n阶实对称矩阵,用zk的Rayleigh商改进近似特征值1,比原来精度提高一倍。
上述两种加速技术,其思想比较简单,但都难以自动完成,因此在编制程序进行实际计算时和幂法一样复杂。
1.3逆幂法
逆幂法用来计算矩阵按模最小的特征值与特征向量,也可用来计算对应于一个给定近似特征值的特征向量。
设n阶实矩阵A(aij)为非奇异矩阵,A特征值次序计为
|1||2|
和相应的特征向量为x1,x2,
|n1||n|0
,xn,则A1的特征值为 1|
|1||1|
|
1
n
||
n121
对应的特征向量为x1,x2,,xn。
1
因此,计算A的按模最小的特征值问题就是计算A的按模最大的特征值问题
逆幂法的迭代公式为
Ayz;
kk1
mk||yk
||; k1,2,yzkk.
mk
(I)
迭代公式(I)的第一个式子,每迭代一步都要求解一个线性方程组,而方程组得系数矩阵均为A,在实际计算时,(I)可改写为
Lykzk1;
Uykyk
mk||yk||; k1,2,
zyk.k
mk
(II)
其中ALU,L是单位下三角矩阵,U为上三角矩阵,A1的按模最大的特征值为
1
n
。
若|n1||n|0,则当k时,就有
x1
zk
||x1||
(II)
1mkn
其收敛速度为|
n1
|。 n
为了提高收敛速度,若同幂法,可以采取原点平移技术,如果已知A的某个特征值i的近似值,并且满足
0|i||j|, ij (IV)
则AI的特征值是i,特征向量仍是xi,而AI的按模最小的特征值是i,若用矩阵AI进行逆迭代
(AI)yz;
kk1
k1,2,mk||yk||;
yzkk.
mk
(V)
则可求得特征向量的近似值。若是A的特征值i的一个粗糙的近似值,还对进行改进。事实上,由于,当k时,有mk
1
,因此,i
当k时,
1
i,只要近似值比较接近i,其收敛速度是惊mk
人的,以至于迭代很少的几次就能得到满意的结果。
§2 双重步位移的QR算法
2.1 QR算法的基本思想
1958年Rutishauser利用矩阵的三角分解提出了计算矩阵特征值的LR算法,1962~1962年Francis等人利用矩阵的QR分解建立了计算矩阵特征值的QR方法。QR方法是电子计算机问世以来计算数学的重大进展之一,是目前计算中小型稠密矩阵全部特征值和特征向量的最有效的方法之一.
设A是n阶实矩阵,记A1A,对A1进行QR分解得,
A1Q1R1 (2.1)
作矩阵 A2R1Q1Q1A1Q 1 (2.2) 再对A2进行QR分解得, A2Q2R2 (2.3) 作矩阵 A3R2Q2Q2A2Q 2 (2.4) 若此继续下去,……,
当求得Ak后,对Ak进行QR分解得, AkQkRk (2.5)
TT
T
作矩阵 Ak1RkQkQkAkQk (2.6)
这样就产生一个矩阵序列Ak。只要A是可逆矩阵,由QR分解可知,
Qk是正交矩阵,Rk为上三角矩阵,且若要求Rk的主对角元均为正数,则
上述每一步分解是唯一的,则由上述方法就可以完全确定Ak。
定理1 (基本QR算法)设AA1是n阶实矩阵,构造QR分解算法:
T
AkQkRk, 其中QkQI, Rk为上三角矩阵;
(2.7)
AmRmQm, k1,2,
记QkQ1Q2Qk, RkRkRk1R1,则有
T
(1)Ak1相似于Ak,即Ak1QkAkQk; (2.8)
T
(2)Ak1(Q1Q2Qk)TA1(Q1Q2
k
T
Qk)QkA1Qk; (2.9)
(3)Ak的QR分解式为AQk Rk。 (2.10) 证明 (1),(2)显然,现证(3)
Qk RkQ1Q2Qk1(QkRk)Rk1R1Qk1 AkRk1
由(2)可得,QkAk1A1QkAQk,从而Qk1AkAQk1,所以
Qk RkQk1 AkRk1AQk1Rk1 (2.11)
由此进一步递推可得结果。
定理2 (基本QR算法的收敛性)设A是n阶实矩阵, (1)如果A的特征值满足:|1||2|
|n1||n|0;
1
,n),且P
1
(2)A有标准形APΛP,其中Λdiag(1,2,
1
有三角分解PLU(L为单位下三角矩阵,U为上三角矩阵)则由QR
算法产生的矩阵序列Ak本质上收敛于上三角矩阵,即 当k时,
1本质上2Ak–––R
若记Ak(aij),则 (i)limaiii;
k
(k)
(2.12) n
(k)
(ii)当ij时,limaij0;当ij时,limaij不一定存在。
k
k
(k)(k)
证明略
若A是对称矩阵,由定理2立得
定理3 如果对称矩阵A满足定理2的条件,则由QR算法产生的矩阵序列Ak收敛于对角矩阵
Λdiag(1,2,
,n)。 (2.13)
2.2 实Schur标准形
由于实际应用中所遇到的大量特征值问题都是关于实矩阵的,因此我们自然希望设计只涉及实运算的QR迭代,即给定ARnn,令AA1,构造迭代;
AkQkRk
k1,2,
Ak1R
kQk
其中Qk为正交矩阵,Rk为上三角矩阵.
(2.14)
然而,此时由于复共轭特征值的存在,我们自然不能期望(2.14)产生的Ak仍然逼近于一个上三角矩阵(即A的Schur标准形).那么,Ak将趋向于什么呢?这就涉及实矩阵在正交相似变换下的标准形问题.完全类似于
Schur分解定理的证明可以证明:
定理2.1 (实Schur分解) 设ARnn,则存在正交矩阵QRnn,使得
R11OQTAQOR12R22OR1mR2m, (2.15) Rmm
其中Rii或者是一个单元素(即是一个实数),或者是一个具有一对共轭复特征值的2阶实矩阵.
这个定理的证明留作作业.
(2.15)通常被称作实矩阵A的实Schur分解,而右边的拟上三角矩阵被称作实矩阵A的实Schur标准形.显然,只要求得一个实矩阵的实Schur标准形,我们就容易求得它的全部特征值.
从前面关于正交迭代法收敛性的粗略分析,不难想象迭代法(2.14)将逼近于A的实Schur分解.事实上也确实如此,在一定条件下,可证(2.14)产生Ak的将逼近于A的实schur标准形.
其次(2.14)作为一种实用的迭代法是没有竞争力的,其原因有二:一是每次迭代的运算量太大(大约是O(n3));二是收敛速度太慢(依赖于特征值的分离程度).因此,要使其成为一种高效的方法,就必须尽可能的减少每次迭代的运算量,提高其收敛性,这就是本节下面所讨论的主要内容.
在下面的讨论中,如无特别说明,我们总假定给定的n阶矩阵A是实矩阵.
2.3 上Hessenberg化
从矩阵的QR分解可知,对一般的方阵完成一次QR迭代所需的运算量是O(n),但是,如果在进行QR迭代之前,先对A进行适当的相似变换,使其具有较多的零元素,即适当选取非奇异矩阵Q0,使
1A0Q0AQ0, (2.16) 3
具有较多的零元素,然后再对A0进行QR迭代,则可望使每次迭代的运算量大为减少.
当然,我们只能利用矩阵计算的三个基本工具来实现Q0和A0的选取,首先来看利用Householder变换,可以得到什么样的A0.
对于给定的A(aij)Rnn,第一步,我们自然应该选取Householder变换H1,使得H1A的第一列尽可能多的零元素(至多只能是n-1个零元素).然而,为了保证对A进行相似变换,在对A进行了行变换之后,必须亦对A进行同样的列变换,即应将H1亦右乘于H1A上变为
H1AH1
这样,为了保证已在H1A的第一列所出现的零元素不致于在右乘H1时被破坏掉,我们应该选取H1具有如下形状
101, (2.17) H10H n11
利用形如(2.17)的Householder变换对A进行相似变换,先将A分块为
aA11
1
即有 T2。 A22
a11H1AH1H11
其中1(a21,a21,T, (2.18) H1A22H1,a1n),A22是A的右下角的T2H1T,an1),2(a12,a13,
n-1阶主子阵.由(2.18)易知,Householder变换H1的最佳选择应该使得
H11pe1, (2.19)
其中pR,这样一来,即可选取形如(2.17)e1是n-1阶单位矩阵的第一列,
的Householder变换H1,使得(2.18)的第一列有n-2个零元素.
然后,在对A22H1A22H1进行同样的考虑,又可找到Householder变换
10H2R(n1)(n1), 0H2
使得
(H2A22H2)e10. 0
从而令
10, H20H2
即有
. H2H1AH1H20O
如此进行n-2步,即可找到n-2个Householder变换H1,H2,,Hn2,使得
Hn2
其中H(hij)满足 H2H1AH1H2Hn2H.
hij0, ij1,
即
h11h21
0H0
0
现令 h12h22h3200h13h23h33h430h1,n1h2,n1h3,n1h4,n1hn,n1h1nh2nh3n, (2.20) h4nhnn通常称形如(2.20)的矩阵为上Hessenberg矩阵.
Q0H1H2
则有 Hn2
TQ0AQ0H (2.21)
即上述过程是将A正交相似变换成一个上Hessenberg矩阵H(具有(n-1)(n-2)/2个零元素).通常称分解式(2.21)为A的上Hessenberg分解.
现在,我们在考虑对上Hessenberg矩阵H进行一次QR迭代
HQR, HRQ
所需的运算量是多少呢?.基于H的特殊形状,H的QR分解可用Givens变换实现,即计算n-1个Givens变换GiG(i,i1,), i1,2,使 ,n1,
Gn1Gn2G1HR
为上三角矩阵.然后, H可按如下方式计算
T HRG1TG2
2TGn1. 容易算出完成这一过程所需的运算量是4n(减少了一个数量级).
此外,容易证明这样得到 H的仍是上Hessenberg矩阵,这样我们就可一直迭代下去,每次迭代的运算量都是4n.
前面利用Householder变换约化矩阵A为上Hessenberg形的方法,可总结为如下的实用算法. 2
算法2.1
(1)输入A(aij),k:=1;
(2)计算n-k阶Householder变换Hk,使
ak1,kak2,k0,A:HAH, Hkkkank0
其中Hkdiag(Ik,Hk).
(3)如果kn2,则k:k1转入(2)否则输出有关信息,结束 这一算法计算出A的上Hessenberg形就存放在A所对应的存储单元内,所需运算量是5n3/3;如果需要积累Q0H1H2
算量2n3/3.
此外,用算法2.1计算得到的上Hessenberg矩阵H满足 Hn2,则还需再增加运
HQT(AE)Q
其中Q是正交矩阵,||E||Fch||A||F,这里c是一常数(详见文献[71]) 当然,我们亦可用Givens变换将A约化为上Hessenberg形,一般需要运算量为10n/3.但是,如果A有较多的零元素,则适当安排Givens变换的次序,可使运算量大为减少.另外,为了节省运算量,也可采用列主元的Gauss消去法,将A利用非正交的相似变换约化为上Hessenberg形.不过这样做,虽然运算量少,但数值稳定性差.
尽管,一般来讲,上Hessenberg分解不是唯一的,然而,可以证明 定理2.2 设ARnn32有如下两个上Hessenberg分解
UTAUH, VTAVG (2.22) 其中U(u1,u2,,un)和V(v1,v2,,vn)是n阶正交矩阵,H(hij)和
G(gij)是上Hessenberg矩阵.若u1v1,而且H的次对角元素hi1,i均不为零,则存在对角元素均为1或-1的对角矩阵D,使得
UVD, HDGD, (2.23)
(即uivi, |hij||gij|, i,j1,2,,n).
证明 假定对某个m(1mn)已证
ujjvj , j1,2,m, (2.24)
其中11,j1或1,下面我们来证,存在m1为1或1,使得
um1m1vm1,
从(2.22)可得
AUUH, AVVG.
分别比较上面两个矩阵等式的第m列,可得
Aumh1mu1
Avmg1mv1hmmumhm1,mum1, (2.25) gmmvmgm1,mvm1, (2.26)
TT分别在(2.25)和(2.26)两边左乘ui和vi (i1,2,,m)可得
hi muiTAum
gi mvAvmTi i1,2,,m. (2.27)(2.28)
由(2.24),(2.27)和(2.28)可得
hi mimgi m, i1,2,,m. (2.29)
将(2.29)代入(2.25),并利用(2.24)和(2.26),可得
hm1,mum1m(Avm12g1mv1
m(Avmg1mv12mgmmvm) gmmvm)
mgm1,mvm1. (2.30)
由(2.30)即知
|hm1,m||gm1,m|,
而hm1,m0,故(2.30)蕴含着um1m1vm1 ,其中m1为1或1.
因此,利用归纳法原理即知(2.23)成立.
一个上Hessenberg矩阵H(hij),如果其次对角元素均不为零,即hi1,i0,i1,2,,n1,则称它是不可约的.定理2.2表明:如果QTAQH为不可约的上Hessenberg矩阵,则Q和H完全由Q的第一列确定(这里是在相差一个正负号的意义下的唯一)均不为零
2.4 双重步位移的QR迭代
下面我们来考虑如何加快QR迭代的收敛速度.假定A的特征值是
|1||n1||n|.
前面关于正交迭代法的收敛性分析表明,随着QR迭代的进行,Ak的右下角的对角元素ann将收敛到特征值(k)n,收敛速度是线性的.速率是|n/n1|.由此可见,如果对某一常数,将A的特征值重新编号,使其满足
|1||n1||n|,
而且比值|n|/|n1|很小,则用AI代替A进行QR迭代可望其收敛速度更快.因此,这就促使我们考虑如下的迭代:
TH1Q0AQ0 (上Hessenberg分解) (2.31) HkIQkRk (QR分解)
HRQI, k1,2,.kkk1
通常称迭代(2.31)为带位移的QR迭代,其中被称作位移.易见,这一迭代产生的Hk与A仍是正交相似的.
不失一般性,我们假定在迭代(2.31)中出现的上Hessenberg矩阵都是
不可约的.如若不然,在迭代的某一步,已有
(k)H11HkOO (k)H22
(k)(k)则我们可以分别对H11和H22进行QR迭代即可.
现在我们来考虑如何有效地选取.当然,大家从前面的分析已经明白,选取的与A的某个特征值越靠近越好,可是,问题是,在我们对A的特征值的信息知道甚少的前提下,如何选取才能使其与A的某个特征值较为靠近呢?
理论分析与实际计算的经验表明,QR迭代产生的矩阵序列右下角最先显露A的特征值.因此,我们可以利用QR迭代的这一特点来选取位移.如果显露的是A的实特征值,即是A的较好的近似特征值时,当然,此时
(k)就可简单地选取hnn作为第k+1次迭代的位移.然而,如果显露的是A
的复共轭特征值时,即
(k)hmmGk(k)hnm(k)hmn, mn1 (2.32) (k)hnn
的特征值是一对互相共轭的复数1和2,且与A的特征值比较接近时,就应该选择Gk的某一特征值i作为位移,但这样一来,就引进了复运算,而这是我们所不希望的.为了避免复运算的出现,人们想到用1和2连续作两次位移,即进行
H1IU1R1,
H1RU111I,
H12IU2R2,
H2R2U22I
这里我们记.对上面迭代产生的矩阵进行一些简单推算,可得
MQR, (2.33)
H2QHHQ, (2.34)
其中
M(H1I)(H2I), (2.35)
QU1U2, RR2R1, (2.36)
由(2.35)可得
MH2sHtI, (2.37)
其中
(k)(k)s22hmmhnnR,
t12det(Gk)R.
因此M是一个实矩阵,而且如果在迭代过程中选取R1和R2的对角元素均为实数,则由(2.33)可推知,Q亦是实的;从而由(2.34)知H2亦是实的,这就是说,在没有误差的情况下,用1和2连续作两次位移进行QR迭代产生H2的仍是实的上Hessenberg矩阵.但是,在实际计算时,由于舍入误差的影响,如此计算得到的H2一般并不是实的.
因此,为了确保计算得到的H2仍是实矩阵,根据(2.33)和(2.34),我们自然想到按如下的步骤来计算H2:
(1) 计算MHsHtI;
(2) 计算M的QR分解:MQR;
(3) 计算H2QHQ
然而,如此计算得第一步形成M就需运算量为O(n),这使得我们前面为减少每次迭代所所需运算量而所做的努力全部付之东流. 3T2
幸运的是,定理2.2告诉我们:不论采用什么样的办法求得正交矩阵Q,使QHQH2是上Hessenberg矩阵,只要保证Q的第一列于Q的第一列一样,则就H2与H2本质上是一样的(所有元素的绝对值相等).因此,我们可有很大的自由度去寻找更有效的方法来实现由H到H2的变换.
首先,我们从(2.33)知,Q的第一列于M的第一列共线,即Qe1由 Me1单位化得到.而由(2.37)容易算出 T
Me1(1,2,3,0,
其中
(k)2(k)(k)(k)1(h11)h12h21sh11t
)k()k()k 2h(
2k1)(h(kh)s, h112232h1,0)T, 32
其次,如果Householder变换P0将Me1变为ae1(即P0(Me1)ae1).其中aR,则Me1aP0e1.故,P0的第一列就与Me1共线,从而P0e1Qe1.而由第二章关于Householder变换的理论知,P0可以按如下方式确定:
P0diag(P0,In3)
其中 P0I3vv, 2/vv, TT
v
(1asign(1),2,3)T,a现在令BP0HP0,则我们只要能够找到第一列为e1的正交矩阵Q使HQBQ为上Hessenberg矩阵,那么H就是我们希望得到的H2.由前面所介绍的约化一个矩阵为上Hessenberg形的Householder方法可知,这是容易办到的.这只需确定n-1个Householder变换P1,P2,T,Pn1,使
(Pn1P2P1)B(P1P2Pn1)H
为上Hessenberg矩阵,即有QP1P2Pn1的第一列为e1,而且由于B所具有的特殊性,实现这一约化过程所需的运算量仅为O(n2).
事实上,由于用P0将H相似变换为B只改变了H的前三行前三列,故B具有如下形状
,
仅比上Hessenberg形矩阵多三个可能的非零元素“”.由B的这种特殊性,易知用来约化B为上Hessenberg矩阵的第一个Householder变换P1具有如下形状样 BP0HP0O
P1diag(1,P1,In4)
其中P1为3阶Householder变换,而且P1BP1具有如下形状
0 P1BP10
0 O,
如此递推地进行,不难推出,第k次约化所用的Householder变换Pk具有如下形状
Pkdiag(Ik,Pk,Ink3),
其中Pk为3阶Householder变换,k2,3,n,,而且Pn3PBP11Pn具有如下形状
Pn3P1BP1Pn3O,
因此,最后一次约化所用的Householder变换具有如下形状
Pn2diag(In2,Pn2)
其中Pn2为2阶Householder变换.
这样我们就找到了一种实现由H到H2的变换方法,它既避免了复运算的出现,又加少了运算量.当然,这一变换过程对Gk的两个特征值均为实数的情形也是可行的.因此,我们就不必在选取位移时区别显露的是实特征值还是复共轭特征值的情形,而只须取作Gk的两个特征值即可.
综述上面的讨论,就得到了著名的Francis双重步位移的QR迭代算法: 算法2.2
(1)输入不可约上Hessenberg矩阵H(hij)Rnn;
(2)m:n1, k:0, s:hmmhnn, t:hmmhnnhmnhnm,
2t:h11h12h21sh11t, y:h21(h11h22s), z:h21h32;
(3)如果kn2,则转入(5);否则,确定Householder矩阵PkR33,使
xPky0,
z0
H:PkHPk (Pkdiag(Ik,Pk,Ink3));
(4)x:hk2,k1, y:hk3,k1,z:
k:k1,转入(3);
(5)确定Householder变换Pn2R22hk4,k1, 当kn-3时, 0, 当kn-3时,,使
xPk, y0
H:Pn2HPn2 (Pn2diag(In2,Pn2))
(6)迭代结束.
这一算法的运算量是6n2;如果需要把正交变换积累起来,还需再增加运算量6n2.
2.5 双重步位移的QR算法
前面的讨论已解决了用QR方法求一个给定的实矩阵的实Schur分解的几个关键的问题.然而,作为一种实用的算法,还需给出一种有效的判定准则,来判定迭代过程中所产生的上Hessenberg矩阵的次对角元素何时可以忽略不计,一种简单而且适用的准则是:当
|hi1,i|(|hii||hi1,i1|), (2.38)
时就将hi1,i看作0.这样做的理由是,在前面约化A为上Hessenberg矩阵H时就已引进了量级为||A||F的误差.
将算法2.1和2.2与收敛准则(2.38)结合起来,就得到了现在流行的双重步位移的QR算法,这一算法是计算一给定的n阶实矩阵A的实Schur分解:QAQT,其中Q为正交矩阵,T为拟上三角矩阵(即对角块为1×1或2×2方阵的块上三角矩阵) T
算法2.3(QR算法)
(1)输入A.
(2)上Hessenberg化:用算法2.1计算A的上Hessenberg分解
THU0AU0; Q:U0.
(3)收敛性判定:
(i)把所有满足条件:
|hi1,i|(|hii||hi1,i1|)
的hi,i1置零;
(ii) 确定最大的非负整数m和最小的非负整数l,使
lH11H12H13HOH22H23nlm
Om, OH33
lnlmm
其中H33为拟上三角形,而H22为不可约的上Hessenberg矩阵;
(iii)如果mn,则输出有关信息,结束;否则进行下一步.
(4)双重步位移QR迭代,对H22用算法2.2迭代一次得H22:PTH22P,其中PP0P1Pnml2.
Q:Qdiag(Il,P,Im)
H12:H12P,H23:PTH23
(5)转入(3).
实际计算的统计表明,这一算法每隔离出一个1阶或2阶子矩阵,平均约需2次Francis迭代.因此,平均来讲,如果只计算特征值,则这一算法的运算量大约是8n;如果Q和T都需要,则为15n.
误差分析的结果表明,这一算法计算所得到的实Schur标准形T正交相似于一个非常靠近A的矩阵,即 33
QT(AE)QT,
其中QTQI,||E||2||A||2
T;计算所得到的Q是几乎正交的,即 QQIF,||F||2,
其中为机器精度,详见文献[37].
§3 特征向量和不变子空间
这一节,我们来讨论在用QR方法求得给定的矩阵的特征值之后,如何求其对应的特征向量和对应的不变子空间的一组标准正交基.
3.1 特征向量的计算
设ARnn,并假定已经利用QR方法求得A的特征值的一个近似值.现在,我们来讨论如何求A之对应于的近似特征向量.
目前,解决这一问题最好的方法就是反幂法,它是乘幂法的自然推广(将乘幂法用于A1上面得到的),常用的是带位移的反幂法,其基本迭代格式为如下:
(AI)vkzk, (3.1a)
zkvk/||vk||2,k1,2,, (3.1b)
其中是事先选定的常数,称作位移;z0是事先给定的向量,称作初始向量.
从(3.1)可以看出,每迭代一次就需解一个线性方程组,这要比乘幂法运算量大得多.但是,由于方程组的系数矩阵不随k的变化而变化,所以能够事先对它进行列选主元素的LU分解,然后每次迭代就只需解两个三角形方程组即可.
另外需要顺便指出的是,这里只是为了下面的分析方便,而在(3.1b)中||||2进行了规范化,在实际使用时是以||||进行规范化的.
假定A是非亏损的,既存在X(x1,x2,,xn)Cnn非奇异,使得
,n), X1AXΛdiag(1,2,
而且还可以假定||xi||1, i1,2,
展开 ,n.现将初始向量z0按x1,x2,,xn
z01x12x2nxnixi, (3.3)
i1n
再假定与A的特征值j最靠近,并有
0|j||i|, ij, (3.4a)
j0. (3.4b)
则从(3.1a),(3.2)和(3.3)可得
vkk(AI)kz0
ki(i)kxi
i1n
kj(j)k(xjuk),
(3.5)
其中k是一个正数,而
iiukxi0, k ijjj
其收敛速度依赖于|j|/min|i|的大小.将(3.5)代入(3.1b)即ijk
得
zkvkk(xjuk), ||vk||2||xjuk||2
其中k为满足|k|1的复数,从而有
Hdist(R(zk),R(xj))||zkzkxjxH
j||2
(xjuj)(xjuj)H
(xjuj)(xjuj)
ijHxjxHj20, k 收敛速度依赖于|j|/min|i|的大小.换句话说,就是zk将按方向
收敛于A的特征向量,与j越靠近就越快.
由此可见,从收敛速度的角度来考虑,用(3.1)进行迭代时,取得越靠近A的某个特征值越好.但是,当与A的特征值很靠近时,AI就与一个奇异矩阵很靠近,每迭代一步就需解一个非常病态的线性方程组.然而,实际计算的经验和理论分析的结果表明:AI的病态性,并不影响其收敛速度,而且当与A的某个特征值很靠近时,常常只需迭代一次就可得到相当好的近似特征向量.为了弄清这点,我们作下面的简要分析.
首先,我们来确立一个判定一个向量v是否作为A的近似特征向量的判定准则.设是A的一个特征值,v是Cn中的一个单位向量(||v||21),定义
r(AI)v, (3.6)
为向量v的剩余向量.从(3.6)易得
(ArvH)vAvrvHvAvrv,
即v是对应于Arv的特征值向量.如果||r||2很小,则亦||rv||2很小;从而如果A之对应于的特征向量是良态的话,当||r||2很小时,v就是A之对应于的很好的近似特征向量.因此,我们可用v的剩余向量的大小来衡量v是否可作为对应于的近似特征向量.
现在,假定 HH
||min||1, (3.7) (A)
其中1是一个很小的正数(1O()).再假定对给定的初始向量z0是用列主元消去法解方程组(3.1a)得到向量v1的计算值v1的,则由Gauss消去法的误差分析结果知
(AIE)v1z0, (3.8)
其中||E||22(2O()).。这样,由(3.1b)计算得
z1v1/||v1||2, (3.9)
这里为了避免符号上的麻烦,忽略了计算z1所引起的误差。
由(3.8)可得向量z1的剩余向量为
r(AI)z1()z1Ez1z0/||v1||2
于是
||r||212||v1||21, (3.10)
这里假定||z0||21。
由此可见,如果计算得到的v1有很大的范数,则由(3.1)迭代一次所得到的向量就有范数很小的剩余向量,从而在特征值问题不是十分病态的条件下,就得到到了很好的近似特征在向量。
现在来说明在条件(3.7)成立的前提下确有v1的范数很大。
设AIE的奇异值分解为
AIEUΣWH, (3.11)
其中U(u1,u2,,un),W(w1,w2,
n0。 ,wn)Un,Σdiag(1,2,,n),1
由是AI的特征值易得
n(AI)||1
其中n(AI)表示AI的最小奇异值;再由第一章的推论6.6得
nn(AI)||E||212, (3.12) 将z0按u1,u2,,un展开,有
z0ajuj,|aj|2||z0||21
j1j1nn
从而有
v(AIE)z0WΣU(ajuj)11H
j1j1nnajjwj, (3.13)
于是
na||v||2j
j1j21/2|an|n|an|。 12
这样一来,只要|an|不是很小(即z0在un方向上不是十分亏损)的话,则||v||2就很大。因此通常反幂法只需迭代一次就足够了。
这里还需指出一点的是,在比较病态时,利用反幂法在进行第二次迭代,一般不会得到更好的近似特征向量。这是由于v1的单位化向量z1可表示为
z1nwnjwj
j1n1
其中n非常接近1,j (j1,2,
按u1,u2,,n1)都是小量,第二次迭代时,将z0,un展开,其系数将主要由nwn决定,即
z1(nwun)un(jwH
juj)uj H
n
j1n1
而此n时很小,因而un和wn将与A对应于的左右特征向量很接近,从而|unwn|s()将很小,即此时z1在un方向上的投影几乎等于零,因此,H
第二次迭代得到的z2就不会有范数很小的剩余向量。现在再看一个具体的例子。
设
11A10 110
他有特征值10.99999和21.00001,以及对应的特征向量x1(1,105)T和x2(1,105)T,两个特征值的条件数均为105数量级。取1,z0(0,1)T应用反幂法迭代一次(在十位十进制的浮点系下进行)则z1(1,0)T,并且||Az1z1||21010。可是再迭代一次将产生
T,则有||Az2z2||21。 z2(0,1)
上述分析表明,利用反幂法球特征向量时,位移量取作较精确的近似特征值最好。此时,一般只需迭代一次就可得到很好的近似特征向量。因此,通常总是在用某种方法求得A的近似特征值之后,再用反幂法求对应的特征向量。连同QR方法一起使用反幂法求特征向量的基本步骤如下:
(1)计算A的Hessenberg分解:U0AU0H;
(2)使用双重步位移的QR方法求出H的特征值,而不积累变换矩阵;
(3)对每个计算的特征值,在(3.1)中取AH,进行迭代,求出z,使Hzz;
(4)计算xU0z(则x就是对应于特征值的近似特征向量)。 T注3.1 在上述方法的第三步用反幂法迭代时,其收敛性是用向量r(AI)z的大小来衡量的;迭代的初始向量,是采用随机选取的方法来选取的,在实际计算时,有专门实现随机选取的程序。
3.2不变子空间的计算
设
ARnn
,2
l
的
1
特征值
1,2,,n
满足
{1
,
,l},并假定集,L{,1,2,},l}关于复共n{
轭是封闭的(即若L,则必有L),则在Rn中存在唯一的l维子空间X使得AXX,而且A在X上的限制A|X的特征值正好是
1,2,,l,即由1,2,,l可唯一确定A的一个l维子空间,下面我们
将主要讨论如何求由1,2,,l确定的不变子空间X的一组标准正交基。
大家知道,如果已知A的实Schur分解为
TTl
QTAQ1112
OT22nl ,
lnl
而且T11的特征值正好是1,2,
则Q的前l列就是对应于1,2,,l,
,l
的A的l维不变之空间X的一组标准正交基。遗憾的是,如果我们使用双
重步位移的QR方法来实现A的实Schur分解的话,则所得到的拟上三角形的对角块是非常任意的,并不能按我们事先需要的次序排列,即若计算得到的Schur分解是
S12lS
QTAQS11
OS22nl
lnl
则我们并不能期望有(S11)L。但是如果我们能够找到一个正交矩阵Q使得
TTl
QTAQT1112
OT22nl
lnl
仍是拟上三角形矩阵,且(T11)L,则问题就解决了。这样一来,我们要求指定的1,2,
,l所对应的A的不变子空间的一组标准正交基的问
题,就归结为按指定次序来重排一个拟三角的对角块的问题。
先考虑A是22矩阵的情形。假定我们已计算出正交矩阵Q使
1t12
QAQT, 12
02
T
希望调换1和2的位置,即确定一个正交矩阵Q,使得
T
QTQ2
0
。 1
容易算出T对应于的特征向量是
x(t12,21)T,
即有Tx2x。再取正交矩阵Q为使Qx的第二个分量为零的Givens变换,则
QAQe12e1。
因此,令QQQ,则
T
2
QAQ
0
Tt12
。 1
这样,我们就解决了n=2时交换1和2的问题。
对于一般情形,A的特征值全为实数的条件之下,使用上述技巧有计划的交换相邻的特征值,就可把(A)的任意子集移到T之对角线的前面,从而求得所需的不变子空间的一组标准正交基。具体算法如下:
算法3.1
(1)
输入上三角矩阵Ttij,正交矩阵Q和指定的特征值子集L{1,2,
(2)
如果{t11,t22,
,l}。
,tll}L,则输出有关信息结束,否则,k:
=1。进行下一步。
(3)
若tkkL而tk1,k1L,则确定ccos,ssin,使
cstk,k1
t
sck1,k1tkk0
T
T:GkTGk,Q:QGk, (GkG(k,k1,));
否则,进行下一步。
(4)
如果kn1,则k:k1,转入(3);否则转入(2)。
注3.2 其中输入的T和Q是由已知矩阵A应用双重步位移
的QR方法计算得到的,即
QTAQT。
当A的特征值并非都是实数时,利用QR方法求得的实Schur标准形之对角块就必含有22的块,此时就会涉及到22块的交换问题。具体地讲,假
定
T11T12T
OT22
其中T11是1阶或2阶矩阵,T22是具有一对复共轭特征值的2阶矩阵,且
(T11)(T22)我们需要确定一个正交矩阵Q,使得
T11T12
, QTQ
OT22
T
满足T11是2阶矩阵,且(T11)(T22)。这可用前面类似的技巧完成,不过要梢微复杂一些,请同学们作为练习自己补出。
§4 对称QR方法
对称QR方法就是求解实对称矩阵的特征值问题的QR方法,是将QR方法应用于对称矩阵,并充分利用其对称性而得到的;是目前求解中小型稠密对称矩阵的特征值和对应的特征向量的最有
效方法之一。
设ASRnn,我们知道,为了减少每次QR迭代所需的运算量,QR方法的第一步就是将给定矩阵A约化为上Hessenberg矩阵,即计算正交矩阵U0,使得
U0AU0T, (4.1)
为上Hessenberg矩阵,而当A对称时,T亦是对称的,从而T是对称三对角矩阵,因此,对实对称矩阵而言,我们首先应该将给定的矩阵A三对角化。而且,在约化过程中再充分利用对称性,还可使约化的运算量大为减少。
将A作如下分块
T
a1v01A
n1 vA00
1n1
从约化一个矩阵为上Hessenberg矩阵的Householder方法不难推出,利用Householder变换将A约化为对称三对角矩阵的第k步为:
(I)计算Householder变换HkR(nk)(nk),使得
Hkvk1ke1, kR;
(II)计算
1ak1
nk1vk
Tvk
HkAk1Hk
Ak
1nk1
这里k1,2,,n1。
如果用上述约化过程产生的ak,k和Hk,定义
a11T
1
a2
n1
Ik, Hkn1O
an
O, Hk
QH1H2Hn1,
则
QTAQT。
从上述约化过程容易看出,第k步约化的主要工作量是计算
HkAk1Hk,设
HkIvvT, vRnk
则利用Ak1的对称性,易得
HkAk1HkAk1vwTwvT, (4.2)
其中
1
wu(vTu)v, uAk1v, (4.3)
2
利用(4.2)和(4.3),并注意到对称性,容易设计出运算量为2(nk)2的计算HkAk1Hk的算法。因此,完成整个约化所需的运算量为而非对称时把一个n阶矩阵约化为上Hessenberg矩阵所需2n3/3,
的运算量为5n3/3。
上述三对角化过程的详细细节,请读者作为练习自己补出,并在此基础上总结出一个实用的算法。
在完成了把A约化为对称三对角矩阵T的任务之后,我们要做的第二件事就是选取适当的位移进行QR迭代。由于此时A的特征值全是实数,因而再使用双重步位移是完全没有必要的,只需进行单步位移即可。
设T是对称三对角矩阵。我们来考虑单步位移的QR迭代:
TkIQkRk
k0,1,2,
TRQ
Ikkk1
(4.4)
其中T0T, Qk是正交矩阵,Rk是上三角矩阵,根据QR迭代保持上Hessenberg形和对称性不变的特点,我们立即知道迭代(4.4)
产生的Tk都是对称三对角矩阵。此外,容易算出每次迭代只需O(n)次运算即可完成,这比非对称的情形降低了一个数量级。 与非对称的QR迭代一样,这里我们亦可假定迭代所用Tk的均是不可约的,即
a1(k)
(k)Tk1
满足i(k)0, 1,2,
1(k)
a2(k)
n1(k)
(k)n1
(k)an
,n1。
现在再来看位移应该如何选取,从非对称QR方法的讨论
(k)
中,我们知道一种简单易行的取法是取an(通常称作Rayleigh
商位移)。但是,另一种更好的选法就是取为矩阵
(k)
an
1(k)n1
(k)
靠近an的特征值,即取
k)
n(1
(k)an
(k)ansign (4.5)
1(k)(k)
其中(an1an)。这就是著名的Wilkinson位移,Wilkinson
2
(参见[71])曾证明了这两种位移都是最终三次收敛的,并说明了为什么后者较前者更好一些。
下面再来考虑如何具体实现一步带位移的QR迭代:
TIQR
(4.6)
TRQI
其中Q为正交矩阵,R为上三角矩阵,
a1T1
1
a2
n1
, 0。
k
n1
an
当然,我们可以利用Givens变换来实现TI的QR分解,进而完成一步迭代。但是,更漂亮的方法是利用定理2.2以隐含的方式来实现由T到T的变换。
大家知道,迭代(4.6)的实质是将T用正交相似变换变为T,即TQTTQ。因而从定理2.2知道T本质上由Q的第一列完全确定。从利用Givens变换实现TI的QR分解过程易知Qe1G1e1,其中G1是将TI的第一列之第二个元素变为零的Givens变换,即
G1G(1,2,),
其中满足
cos
sin
sina1
。
cos10
T
令 BG1TTG1。
则B(例如,n=5)有如下形状
0000B0,
00000
仅比三对角矩阵多两个非零元“”。由此即知,只需用Givens变换将B和约化为三对角矩阵即可得到所需的三对角矩阵T。下面就n=5的情形,来说明这一约化过程。
首先用(2,3)坐标平面上的Givens变换G2将B之(1,3)和(3,1)位置上的元素消为零,这样在(2,4)和(4,2)位置上又会出现两个我们不希望有的非零元;接着再用(3,4)平面上的Givens变换G3将(2,4)和(4,2)位置上的元素消为零,又在(3,5)和(5,3)位置上出现两个非零元;最后利用(4,5)平面上的Givens变换G4将(3,5)和(5,3)位置上的元素消为零,而得到所需的三对角矩阵。这一过程用图示方式形象地表述如下:
× ×
×
×
G2
×
出境
从图示可以看出,整个约化过程就是把三对角线之外不受欢迎的非零元“”逐步“驱赶”出矩阵之外。因此,我们可以称这一约化方法为“驱逐出境法”。 对于一般的n,我们可以通过n-2次Givens变换G2G3
Gn1使
TGn1T
G2(G1TTG1)G2Gn1
是三对角矩阵,且(G1G2Gn1)e1G1e1Qe1。
综上所述,可得带Wilkinson位移的对称QR迭代算法如下:
算法4.1
(1) 输入T的对角元素a1,a2,
,an和次对角元素
1,2,,n,1
k:=1。
(2) :(an1an)/2,
:ann21/(sign, x:a1, y:1。
(3) 计算ccos, ssin和使
csx, scy0
akkcsak:ak1sck
k
s
, ak1sc
kc
k1:ck1。
(4) 如果k1,则k1:;否则进行下一步。 (5) 如果kn1,则
x:k, y:sk1, k:k1,
转入(3),否则迭代结束。
这一算法运算量为:乘除运算20n,开方n-1次;如果对给定
的正交矩阵Q,还需计算QG1G2 Gn1,则还需再增加运算量4n2。
至此,我们已经解决了对称QR方法的几个关键问题,可以归
结为如下算法:
算法4.2(对称QR算法)
(1) 输入ASRnn和误差限, (2) 三对角化A:
计算householder变换H1H2
Hn2,使
2
T
(H1H2Hn2)AH(H1Hn )2
a1
1
(3) 受敛性检验, (i)将所有满足
1
a2
n1
, n1
an
|j|(|ai||ai1|)
的j置零; (ii)如果i0,i1,2,
否则
0:0,确定正整数pq,使得
i p10,
i0,pp,
1q,, ,
,n1,则输出有关信息,结束;
q1n10。
(4) QR迭代
对对称三角矩阵
appapp1, Tqaqq1
应用算法4.1,然后转步(3)。
这一算法是数值代数中最漂亮的算法。误差分析的结果表明,由算法4.2计算所得到的特征值1,2,,n满足
,n), QT(AE)Qdiag(1,2,
其中QTQI,||E||2||A||2,从而由第一章的Weyl定理知
|ii|2||A||2, i1,2,,n
其中i是A的精确特征值,且与i排列次序一致,这就是说,对称QR算法计算得到的特征值是相当精确的,相对误差不超过机器精度,但值得注意的是,特征向量的计算值并不一定亦有这样的精度,它与i和其他特征值的分离程度有关。