矩阵特征值的计算
矩阵特征值的计算
物理、力学和工程技术中的许多问题在数学上都归结为求矩力学和工程技术中的许多问题在数学上都归结为求矩阵的特征值和特征向量问题。
�计算方阵A的特征值,就是求特征多项式方程:
|A−λI|=0
即
λ+p1λ
nn−1
+⋅⋅⋅+pn−1λ+pn=0
的根。求出特征值λ后,再求相应的齐次线性方程组:
(A−λI)x=0
的非零解,即是对应于λ的特征向量。这对于阶数较小的矩阵是可以的,但对于阶数较大的矩阵来说,求解是十分困难,所以用这种方法求矩阵的特征值是不切实际的。
�如果矩阵A与B相似,则A与B有相同的特征值。
把A化为最简单的形式。因此希望在相似变换下,一般矩阵(Jordan)标准形。由于在一般情况下,的最简单的形式是约当约当(用相似变换把矩阵A化为约当标准形是很困难的,所以设法对矩阵A依次进行相似变换,使其逐步趋向于一个约当标准形,从而求出A的特征值。
下面主要介绍求部分特征值和特征向量的幂法、反幂法,求实对称矩阵全部特征值和特征向量的雅可比(Jacobi)方法,求任意矩阵全部特征值的QR方法。
§1幂法与反幂法
一、幂法
特�幂法:是一种求任意矩阵A的按模最大特征值及其对应的按模最大特征值及其对应特征向量的迭代算法。
该方法最大的优点是计算简单,容易在计算机上实现,对稀疏矩阵较为合适,但有时收敛速度很慢。
为了讨论简单,我们假设:
(1)n阶方阵A的特征值λ1,λ2,⋅⋅⋅,λn按模的大小排列
|λ1|>|λ2|≥⋅⋅⋅≥|λn|
(3)v1,v2,⋅⋅⋅,vn线性无关。
(1)
(2)vi是对应于特征值λi的特征向量(i=1,2,⋅⋅⋅,n);
任取一个非零的初始向量x0,由矩阵A构造一个向量序列
xk=Axk−1 , k= 1, 2, ⋅⋅⋅
称为迭代向量。
(2)
由于v1,v2,⋅⋅⋅,vn线性无关,构成n维向量空间的一组基,所以,初始向量x0可唯一表示成:
x0=a1v1+a2v2+⋅⋅⋅+anvn
于是
(3)
xk=Axk−1=A2xk−2=⋅⋅⋅=Akx0
kk
=a1λ1v1+a2λkv+⋅⋅⋅+aλ22nnvn
λ2kλnkk
=λ1[a1v1+a2()v2+⋅⋅⋅+an()vn]
λ1λ1
(4)
|λi|
因为比值|λ|
1
xk
limk=a1v1k→∞λ1
当k充分大时,有
kxk≈λ1a1v1
(5)
(6)
从而
k+1xk+1≈λ1a1v1
(7)
这说明当k充分大时,两个相邻迭代向量xk+1与xk近似地相差一个倍数,这个倍数便是矩阵A的按模最大的特征值λ1。若用
(xk)i表示向量xk的第i个分量,则
(xk+1)iλ1≈
(xk)i
最大的特征值。
(8)
即两个相邻迭代向量对应分量的比值近似地作为矩阵A的按模
因为xk+1≈λ1xk,又xk+1=Axk,所以有Axk≈λ1xk,因此向量xk可近似地作为对应于λ1的特征向量。
这种由已知的非零向量x0和矩阵A的乘幂构造向量序列
{xk}以计算矩阵A的按模最大特征值及其相应特征向量的方法
称为幂法。
|λ2|
由(4)式知,幂法的收敛速度取决于比值|λ|的大小。比
1
|λ2|
值越小,收敛越快,但当比值|λ|接近于1时,收敛十分缓慢。
1
用幂法进行计算时,如果|λ1|>1,则迭代向量xk的各个不为零的分量将随着k无限增大而趋于无穷。反之,如果
|λ1|
机上计算时就可能溢出停机。为了避免这一点,在计算过程中,常采用把每步迭代的向量xk进行规范化,即用xk乘以一个常1。数,使得其分量的模最大为数,使得其分量的模最大为1
�幂法迭代公式:
⎧yk=Axk−1
⎪
⎨mk=max(yk) ( k= 1, 2, ⋅⋅⋅)⎪x=y/m
kk⎩k
其中mk是yk模最大的第一个分量。相应地取
(9)
⎧λ1≈mk
⎨
⎩v1≈xk ( 或 yk)
例1设
(10)
⎡ 2 −1 0⎤
⎥A=⎢−1 2 −1⎢⎥
⎢⎣ 0 −1 2⎥⎦
用幂法求其模为最大的特征值及其相应的特征向量(精确到小数点后三位)。
Tx=(1, 1, 1)解取0,计算结果如表1所示。
表1
k
1234567
123-2.5-2.428-2.416-2.414
ykT
0-2-43.53.4283.4163.414
123-2.5-2.428-2.416-2.414
mk
12-43.53.4283.4163.414
11-0.75-0.714-0.708-0.707-0.707
xk
0-111111
11-0.75-0.714-0.708-0.707-0.707
当k=7时,xk已经稳定,于是得到:
λ1≈m7=3.414
及其相应的特征向量v1为:
v1≈x7=(−0.707, 1, −0.707)T
�应用幂法时,应注意以下两点:
(1)应用幂法时,困难在于事先不知道特征值是否满足(1)式,以及方阵A是否有n个线性无关的特征向量。克服上述困难预的方法是:先用幂法进行计算,在计算过程中检查是否出现了在计算过程中检查是否出现了预期的结果。如果出现了预期的结果,就得到特征值及其相应特征向量的近似值;否则,只能用其它方法来求特征值及其相应的特征向量。
(2)如果初始向量x0选择不当,将导致公式(3)中v1的系数
a1等于零。但是,由于舍入误差的影响,经若干步迭代后,xk=Akx0。按照基向量v1,v2,⋅⋅⋅,vn展开时,v1的系数可能不
等于零。把这一向量xk看作初始向量,用幂法继续求向量序列
xk+1,xk+2,⋅⋅⋅,仍然会得出预期的结果,不过收敛速度较慢。
如果收敛很慢,可改换初始向量。
二、原点平移法
|λ2|
由前面讨论知道,幂法的收敛速度取决于比值|λ|的大
1
小。当比值接近于1时,收敛可能很慢,一个补救的方法是采用
原点平移法。
设矩阵
B=A−pI
其中p为要选择的常数。
(11)
A与B除了对角线元素外,其它元素都相同,而A的特征
并且相应的特值λi与B的特征值µi之间有关系µi=λi−p,并且相应的特征向量相同。
这样,要计算A的按模最大的特征值,就是适当选择参数
p,使得λ1−p仍然是B的按模最大的特征值,且使
λ2−p|λ2|
λ1−p|λ1|
对B应用幂法,使得在计算B的按模最大的特征值λ1−p
的过程中得到加速,这种方法称为原点平移法。
例2设4阶方阵A有特征值
λi=15−i , i=1, 2, 3, 4
|λ2|r==0.9,令p=12作变换比值|λ1|
B=A−pI
则B的特征值为
µ1=2 , µ2=1 , µ3=0 , µ4=−1
应用幂法计算B的按模最大的特征值µ1时,确定收敛速度的比值为:
µ2λ2−p|λ2|
==0.5
所以对B应用幂法时,可使幂法得到加速。
虽然选择适当的p值,可以使得幂法得到加速,但由于矩阵的特征值的分布情况事先并不知道,所以在计算时,用原点平移法有一定的困难。
下面考虑当A的特征值为实数时,如何选择参数p,以使得用幂法计算λ1时得到加速的方法。
设A的特征值满足:
λ1>λ2≥λ3≥⋅⋅⋅≥λn−1>λn
则对于任意实数p,B=A−pI的按模最大的特征值λ1−p或λn−p。
如果需要计算λ1及v1时,应选择p使
|λ1−p|>|λn−p|
且确定的收敛速度的比值:
⎧λ2−pλn−p⎫
r=max⎨,⎬=min
⎩λ1−pλ1−p⎭
λ2+λn
当λ2−p=−(λn−p),即p=时,r为最小。这时用
2
幂法计算λ1及v1时得到加速。
如果需要计算λn及vn时,应选择p使
|λn−p|>|λ1−p|
且确定收敛速度的比值:
⎧λ1−pλn−1−p⎫
r=max⎨,⎬=min
⎩λn−pλn−p⎭
λ1+λn−1
当λ1−p=−(λn−1−p),即p=时,r为最小。这时
2
用幂法计算λ1及v1时得到加速。
原点平移的加速方法,是一种矩阵变换方法。这种变换容易计算,又不破坏A的稀疏性,但参数p的选择依赖于对A的特征值的分布有大致了解。
三、反幂法
征�反幂法:用于求矩阵A的按模最小的特征值和对应的特的按模最小的特征值和对应的特征向量,及其求对应于一个给定的近似特征值的特征向量。
设n阶方阵A的特征值按模的大小排列为:
|λ1|≥|λ2|≥⋅⋅⋅≥|λn−1|>|λn|>0
相应的特征向量为v1,v2,⋅⋅⋅,vn。则A的特征值为:
−1
1111≤≤⋅⋅⋅≤
对应的特征向量仍然为v1,v2,⋅⋅⋅,vn。因此,计算矩阵A的按模
−1
A最小的特征值,就是计算的按模最大的特征值。−1
�反幂法的基本思想:把幂法用到A上。
任取一个非零的初始向量x0,由矩阵A构造向量序列:
−1
xk=A−1xk−1 , k=1, 2, ⋅⋅⋅
−1
−1
(12)
用(12)式计算向量序列{xk}时,首先要计算逆矩阵A。由于计算A时,一方面计算麻烦,另一方面当A为稀疏阵时,A不算一定是稀疏阵,所以利用A进行计算会造成困难。在实际计在实际计算12)式等价于:时,常采用解线性方程组的方法求xk。(。(1212)式等价于:
−1
−1
Axk=xk−1 , k=1, 2, ⋅⋅⋅
为了防止溢出,计算公式为
(13)
⎧Ayk=xk−1
⎪
⎨mk=max(yk) ( k= 1, 2, ⋅⋅⋅)⎪x=y/m
kk⎩k
相应地取
(14)
1⎧
⎪λn≈
mk⎨
⎪v≈y ( 或 x)
kk⎩n
(15)
(13)式中方程组有相同的系数矩阵A,为了节省工作量,可先对矩阵A进行三角分解:
A=LU
再解三角形方程组:
(16)
⎧Lzk=xk−1
k=1, 2, ⋅⋅⋅⎨
⎩Uyk=zk
(17)
当A是三对角方阵,或是非零元素较少且分布规律的方阵时,无论存储或计算都比较便。根据幂法的讨论,我们知道,在一定条件下,可求得A的按模最大的特征值和相应的特征向量,从而得到A的按模最小的特征值和对应的特征向量,称这种方法为反幂法。
反幂法也是一种迭代算法,每一步都要解一个系数矩阵相同每一步都要解一个系数矩阵相同的线性方程组。
−1
�原点平移的反幂法
−1
A−pI(A−pI)p设为任一实数,如果矩阵可逆,则的特
征值为
111
,,⋅⋅⋅,λ1−pλ2−pλn−p
对应的特征向量仍为v1,v2,⋅⋅⋅,vn。
如果p是矩阵A的特征值λi的一个近似值,且
|λi−p|
1−1
(A−pI)则λ−p是矩阵的按模最大的特征值。因此,当给
i
出特征值λi的一个近似值p时,可对矩阵A−pI应用反幂法,求出对应于λi的特征向量。反幂法迭代公式中的yi通过方程组
(A−pI)yk=xk−1 , k=1, 2, ⋅⋅⋅
求得。
例3用反幂法求矩阵
⎡⎢ 2 −1 0 0⎤
A= ⎢−1 2 −1 0⎥
⎢⎥
⎢ 0 −1 2 −1⎥
⎣ 0 0 −1 2⎥⎦
的对应于特征值λ=0.4的特征向量。
解取xT
0=(1, 1, 1, 1)解方程组
(A−0.4I)y1=x0
得
y1=(−40, −65, −65, −40)T
m1=max(y1)=−65
x188T
1=my1=(13, 1, 1, 13)
1
再解方程组
(A−0.4I)y2=x1
得
y[1**********]5T
2=(−13, −13, −13, −13
720m2=max(y2)=−13
18989Tx2=y2=(, 1, 1, )m2144144
x1与x2的对应分量大体上成比例,向所以对应于λ=0.4的特征的特征向量为:
v=( 89
144, 1,, 89T
144) 1
§2雅可比(Jacobi)方法
�雅可比方法:是用来计算实对称矩阵A的全部特征值及其相应特征向量的一种变换方法。
在介绍雅可比方法之前,先介绍方法中需要用到的线性代数知识与平面上的旋转变换。
一、预备知识
(1)如果n阶方阵A满足:
ATA=I ( 即 A−1=A )
则称A为正交阵。
(2)设A是n阶实对称矩阵,则A的特征值都是实数,并且有互相正交的n个特征向量。
(3)相似矩阵具有相同的特征值。
(4)设A是n阶实对称矩阵,P为n阶正交阵,则B=PTAP也是对称矩阵。
(5)n阶正交矩阵的乘积是正交矩阵。
(6)设A是n阶实对称矩阵,则必有正交矩阵P,使
⎡λ1 ⎤⎢ λ⎥
2⎥=ΛPTAP=⎢⎢ ⋱⎥⎢⎥⎣ λn⎦(18)
其中Λ的对角线元素的是A的n个特征值,正交阵P的第i列是A的对应于特征值λi的特征向量。
由(6)可知,对于任意的n阶实对称矩阵A,只要能求得一
T),则可得到A的全部个正交阵P,使PAP=Λ(Λ为对角阵为对角阵)的全部
特征值及其相应的特征向量,这就是雅可比方法的理论基础。
二、旋转变换
设
⎛a11 a12⎞A=⎜⎜a a⎟⎟⎝2122⎠
为二阶实对称矩阵,即a12=a21。因为实对称矩阵与二次型是一一对应的,设A对应的二次型为
2f(x1,x2)=a11x12+2a12x1x2+a22x2(19)
由解析几何知识知道,方程f(x1,x2)=C表示在x1,x2平面上的一条二次曲线。如果将坐标轴Ox1,Ox2旋转一个角度θ,使得旋转后的坐标轴Oy1,Oy2与该二次曲线的主轴重合,如图所示。
则在新的坐标系中,二次曲线的方程就化成:
2λ1y12+λ2y2=C(20)
这个变换就是:
⎛x1⎞⎛cosθ −sinθ⎞⎛y1⎞⎜⎜sinθ cosθ⎟⎟⎜⎜x⎟⎟=⎜⎜y⎟⎟⎠⎝2⎠⎝2⎠⎝(21)
变换(21)把坐标轴进行旋转,所以称为旋转变换。其中
⎛cosθ −sinθ⎞P=⎜⎜sinθ cosθ⎟⎟⎝⎠
的变换过程即(22)T称为平面旋转矩阵。显然有PP=I,所以P是正交矩阵。上面
⎛λ1 ⎞PAP=⎜⎜ λ⎟⎟2⎠⎝T
由于
22⎛acosθ+asinθ+a12sin2θ 0.5(a22−a11)sin2θ+a12cos2θ⎞1122T⎜⎟PAP=⎜22⎟⎝0.5(a22−a11)sin2θ+a12cos2θ a11sinθ+a22cosθ−a12sin2θ⎠
所以只要选择θ满足:
0.5(a22−a11)sin2θ+a12cos2θ=0
即:2a12tan2θ=a11−a22(23)
πθ=(当a11=a22时,可选取)4
PTAP就成对角阵,这时A的特征值为:
22⎧λ=acosθ+asinθ+a12sin2θ⎪11122⎨22⎪λ=asinθ+acosθ−a12sin2θ1122⎩2
相应的特征向量为:
⎛cosθ⎞⎛−sinθ⎞v1=⎜⎜sinθ⎟⎟ , v2=⎜⎜ cosθ⎟⎟⎝⎠⎝⎠
三、雅可比方法
阵�雅可比方法的基本思想:是通过一系列的由平面旋转矩是通过一系列的由平面旋转矩阵
到构成的正交变换将实对称矩阵逐步化为对角阵,从而得构成的正交变换将实对称矩阵逐步化为对角阵,从而得到
A的全部特征值及其相应的特征向量。
首先引进R中的平面旋转变换。变换n
⎧xi=yicosθ−yisinθ⎪⎨xj=yisinθ+yicosθ
⎪⎩xk=yk k≠i,j
记为x(24)=Pijy,其中:
⎛1⎞⎜⎟⎜ ⋱⎟⎜ cosθ ⋯ −sinθ⎟i⎜⎟⎜ ⋮ ⋮⎟Pij=⎜⎟j sinθ ⋯ cosθ⎜⎟⎜ 1⎟⎜⎟⎜ ⋱⎟⎜ 1⎟⎝⎠
i j(25)
x=(x1,x2,⋅⋅⋅,xn)T
y=(y1,y2,⋅⋅⋅,yn)T
则称x=Pijy为Rn中xi,xj平面内的一个平面旋转变换,Pij称为xi,xj平面内的平面旋转矩阵。容易证明Pij具有如下简单性质:
①Pij为正交矩阵。
②Pij的主对角线元素中除第i个与第j个元素为cosθ外,其它元素均为1;非对角线元素中除第i行第j列元素为−sinθ,第j行第i列元素为sinθ外,其它元素均为零。
T③PA只改变的第i行与第j行元素,AP只改变的第i列
T与第j列元素,所以PAP只改变的第i行、第j行、第i列、
第j列元素。
设A=(aij)n×n (n≥3)为n阶实对称矩阵,aij=aji≠0为一对非对角线元素。令
(1)A1=PTAP=(aij)n×n
则A1为实对称矩阵,且A1与A有相同的特征值。通过直接计算知:
(1)⎧aii⎪(1)⎪ajj⎪
(1)⎪a⎪ij⎨⎪a(1)
⎪ik
1)⎪a(
jk⎪(1)⎪⎩akl=aiicos2θ+ajjsin2θ+aijsin2θ=aiisin2θ+ajjcos2θ−aijsin2θ=a(1)ji(1)=aki1=(ajj−aii)sin2θ+aijcos2θ2=aikcosθ+ajksin2θ k≠i,j(26)(1)=akj=−ajksinθ+ajkcos2θ k≠i,j=akl k,l≠i,j
2aij当取θ满足关系式:tan2θ=aii−ajj(27)
π(当aii=ajj时,可选取θ=)4
(1)(1)a=a时,ijji=0,且
(1)21)22⎧(aik)+(a(
jk)=aik+a2
jk , k≠i,j⎪⎪(1)2(1)2222(a)+(a)=a+a+2a⎨iijjiijjij⎪(1)22(a)=a⎪kl , k,l≠i,j⎩kl(28)
由于在正交相似变换下,矩阵元素的平方和不变,所以若用D(A)表示矩阵A的对角线元素平方和,用S(A)表示A的非对角线元素平方和,则由(28)式得:
2⎧D(A)=D(A)+2a1ij⎪⎨2(29)S(A)=S(A)−2a⎪1ij⎩
这说明:用Pij对A作正交相似变换化为A1后,A1的对角线元素平方和比A的对角线元素平方和增加了2aij;A1的非对角线2
22aA元素平方和比1的非对角线元素平方和减少了ij,且将事先
选定的非对角线元素消去了(即aij(1)=0)。因此,只要逐次地用这种变换,就可以使得矩阵A的非对角线元素平方和趋于的非对角线元素平方和趋于零,也即使得矩阵A逐步化为对角阵。
这里需要说明一点:并不是对矩阵A的每一对非对角线非零元素进行一次这样的变换就能得到对角阵。因为在用变换消去aij的时候,只有第i行、第j行、第i列、第j列元素在变化,如果aik或Pkj为零,经变换后又往往不是零了。
雅可比方法:就是逐步对矩阵A进行正交相似变换,消去非对角线上的非零元素,直到将A的非对角线元素化为接近于零为止,从而求得A的全部特征值,把逐次的正交相似变换矩阵乘起来,便是所要求的特征向量。
雅可比方法的计算步骤归纳如下:
第一步在矩阵A的非对角线元素中选取一个非零元素
aij。一般说来,取绝对值最大的非对角线元素;
第二步由公式tan2θ=a−a求出θ,从而得平面旋转矩
iijj
阵P1=Pij;
TA=P第三步11AP1,A1的元素由公式(26)计算。2aij
第四步以A1代替A,重复第一、二、三步求出A2及P2,继续重复这一过程,直到Am的非对角线元素全化为充分小(即小于允许误差)时为止。
第五步Am的对角线元素为A的全部特征值的近似值,P=P1P2⋅⋅⋅Pm的第j列为对应于特征值λj(λj为Am的对角线上第j个元素)的特征向量。
例4用雅可比方法求矩阵
⎡⎢ 2 −1 0⎤
A=⎢−1 2 −1⎥
⎢ 0 −1 2⎥
⎣⎥⎦
的特征值与特征向量。
解首先取i=1,j=2,由于a11=a22=2,故取θ=π
4,所以
⎡ 1/2 −1/2
PP⎢0⎤1=12=⎢⎢1/2 1/2 0⎥⎥
⎢⎣ 0 0 1⎥⎥⎦
⎡ 1 0 −1/2
A=PT⎢⎤⎥
11AP1=⎢⎢0 3 −1/2⎥
⎢ −1/2 2⎥⎣ −1/2⎥⎦
再取i=1,j=3,由
2×(−1
tan2θ=2)
1−2=2
得
sinθ≈0.45969 , cosθ≈0.88808
所以
⎡ 0.88808 0 −0.45969⎤⎢0 1 0⎥P=1⎢⎥⎢⎣ 0.45969 0 0.88808⎥⎦
⎡ 0.63398 −0.32505 0⎤
⎥A2=P2TAP2=⎢−0.32505 3 −0.62797⎢⎥
⎢⎣ 0 −0.62797 2.36603⎥⎦
继续做下去,直到非对角线元素趋于零,进行九次变换后,得
⎡0.58758 0.00000 0.00000⎤
⎥A9=⎢0.00000 2.00000 0.00000⎢⎥
⎢⎣0.00000 0.00000 3.41421⎥⎦
A9的对角线元素就是A的特征值,即
λ1≈0.58758 , λ2≈2.00000 , λ3≈3.41421
相应的特征向量为
⎡0.50000
v1=⎢⎢0.70710
⎢⎣0.50000⎤⎡ 0.70710
⎥ , v=⎢ 0.00000
2⎥⎢⎥⎢⎦⎣−0.70710⎤⎡ 0.50000
⎥ , v=⎢−0.70710
3⎥⎢⎥⎢⎦⎣ 0.50000
⎤
⎥⎥⎥⎦
相应的特征值的精确值
λ1=2−
相应的特征向量为
2 , λ2=2 , λ3=2+
2
⎡ 1/2⎤⎡ 1/2⎤⎡ 1/2⎤
⎢⎥⎢⎥⎢⎥v1=⎢1/2⎥ , v2=⎢ 0⎥ , v3=⎢−1/2⎥
⎢ 1/2⎥⎢−1/2⎥⎢ 1/2⎥⎣⎦⎣⎦⎣⎦
由此可见,雅可比方法变换九次的结果已经相当精确了。 用雅可比方法求得的结果精度都比较高,特别是求得的特征值向量正交性很好,所以雅可比方法是求实对称矩阵的全部特征雅可比方法是求实对称矩阵的全部特征值及其对应特征向量的一个较好的方法。但由于上面介绍的雅可比方法,每次迭代都选取绝对值最大的非对角线元素作为消去对象,花费很多机器时间。另外当矩阵是稀疏矩阵时,进行正交相似变换后并不能保证其稀疏的性质,所以对阶数较高的矩阵不宜采用这种方法。
调目前常采用一种过关雅可比方法。这种方法是选取一个单这种方法是选取一个单调
an=0)减小而趋于零的数列{an}(即a1>a2>⋅⋅⋅,且limn→∞
作为限值,这些限值称为“关”,对矩阵的非对角线元素规定一个顺序(例如先行后列、自左至右的顺序)。首先对限值a1按规定的顺序逐个检查矩阵的非对角线元素,碰到绝对值小于a1的元素就跳过去,否则就作变换将其化为零。重复上述过程,直到所有的非对角元素的绝对值都小于a1为止。再取a2,a3,⋅⋅⋅类似处理,直到所有的非对角线元素的绝对值都小于am时,迭代停止。这时的am应小于给定的误差限ε。
实际运算中常用如下的办法取限值:对于矩阵A,计算
A=A0的非对角线元素平方和S(A0),任取N≥n,取
ak=
S(Ak−1)
, k= 1 , 2 , ⋅⋅⋅N
§3
QR算法
QR算法:也是一种迭代算法,是目前计算任意实的非奇异
矩阵全部特征值问题的最有效的方法之一。该方法的基础:是构造矩阵序列{Ak},并对它进行QR分解。
由线性代数知识知道,若A为非奇异方阵,则A可以分解为正交矩阵Q与上三角形矩阵R的乘积,即A=QR,而且当R的对角线元素符号取定时,分解式是唯一的。
若A为奇异方阵,则零为A的特征值。任取一数p不是A的特征值,则A−PI为非奇异方阵。只要求出A−PI的特征值,就很容易求出A的特征值,所以假设A为非奇异方阵,并不妨碍讨论的一般性。
设A为非奇异方阵,令A1=A,对A1进行QR分解,即把
A1分解为正交矩阵Q1与上三角形矩阵R1的乘积
A1=Q1R1
令
A2=R1Q1=Q1TA1Q1
继续对A2进行QR分解并定义
A2=Q2R2
T
A3=R2Q2=Q2A2Q2
一般地,递推公式为
⎧A1=A=Q1R1⎨TA=RQ=QkkkAkQk , k=2, 3, ⋅⋅⋅⎩k+1
QR算法就是利用矩阵的QR分解,按上述递推公式构造矩阵按上述递推公式构造矩阵
序列{Ak}。只要A为非奇异方阵,则由QR算法就完全确定
{Ak}。这个矩阵序列{Ak}具有下列性质:
性质1所有Ak都相似,它们具有相同的特征值。性质1证明因为
Ak+1=RkQk=QkTAkQk=QkTQkT−1Ak−1Qk−1Qk
=⋅⋅⋅=QkTQkT−1⋅⋅⋅Q1TAQ1⋅⋅⋅Qk−1Qk
~~
若令Qk=Q1Q2⋅⋅⋅Qk,则Qk为正交阵,且有
~T~Ak=QkAQk
因此Ak与A相似,它们具有相同的特征值。
性质2Ak的QR分解式为性质2
~~
A=QkRk
k
~~Q=QQ⋅⋅⋅Q , R其中k12kk=RkRk−1⋅⋅⋅R1
证明用归纳法。显然当k=1时,有
~~
A=A1=Q1R1=Q1R1
~~k−1k−1
假设A有分解式A=Qk−1Rk−1
于是
~~~~QkRk=Q1Q2⋅⋅⋅Qk−1(QkRk)Rk−1⋅⋅⋅R1=Qk−1AkRk−1~T~
因为Ak=Qk−1AQk−1,所以
~~~~
QkRk=AkQk−1Rk−1=Ak
~~
因为Q1,Q2,⋅⋅⋅,Qk都是正交阵,所以Qk也是正交阵,同样Rk
也是上三角形阵,从而A的QR分解式为
k
~~Ak=QkRk
~T
由由前面的讨论知Ak+1=QkAQk。这说明QR算法的收敛性算法的收敛性由
~{Q正交矩阵序列k}的性质决定。
~~{Q}RQ1如果k收敛于非奇异矩阵∞,k为上三角形矩定理定理1
Ak存在并且是上三角形矩阵。阵,则limk→∞
~
证明因为{Qk}收敛,故下面极限存在
~T~T
limQk=limQk−1Qk=Q∞Q∞=Ik→∞k→∞
~T~TT
R∞=limRk=limAk+1Qk=limQkAQk−1=Q∞AQ∞
k→∞
k→∞
k→∞
由于Rk (k=1,2,⋅⋅⋅)为上三角形矩阵,所以R∞为上三角形矩阵。又因为
A∞=limAk=limQkRk=R∞
k→∞
k→∞
Ak存在,并且是上三角形矩阵。所以limk→∞
2(QR算法的收敛性)设A为n阶实矩阵,且定理定理2算法的收敛性)(1)A的特征值满足:|λ1|>|λ2|>⋅⋅⋅>|λn|>0(2)A=A1=XDX
−1
,其中D=diag(λ1,λ2,⋅⋅⋅,λn)且设
X−1有三角分解式X−1=LU(L为单位下三角阵,U为上三角
阵)。
则由QR算法得到的矩阵序列{Ak}本质上收敛于上三角形矩阵。
(k)
A=(a即kij)n×n满足
(k)
aii→λi,当k→∞
(k)aij→0,当i>j, k→∞(k)aij (j>i)的极限不一定存在。
T
证明因为Ak+1=QkAQk,矩阵Qk决定{Ak}的收敛性。又
~~~~
A1k=QkRk,我们利用A1k求Qk,然后讨论{Qk}的收敛性。
−1
A=A=XDX由定理条件,得1
A1k=XDkX−1=XDkLU=X(DkLD−k)(DkU)
令
DkLD−k=I+Bk
(k)
bB(i,j)其中k的元素ij为
(k)
bij
⎧λik
⎪lij(λ) , i>j , lij∈L=⎨j
⎪0 , i≤j⎩
于是
A1k=X(I+Bk)(DkU)
λi
由假设,当i>j时,λ
j
limBk=0
k→∞
设方阵X的QR分解式为
X=QxRx
由
A1k=(QxRx)(I+Bk)(DkU)=Qx(I+RxBkRx−1)(RxDkU)
Bk=0知,对充分大的k, I+RxBkRx−1非奇异,它应有唯一由limk→∞
的QR分解式kk,并且
limk=I, limk=I
k→∞
k→∞
于是
A1k=(Qxk)(kRxDkU)
k
RD但上三角阵kxU的对角线元素不一定大于零。为此,引入
对角矩阵
Dk=diag(±1,±1,⋅⋅⋅,±1)
kkRD以便保证(kxU)的对角线元素都是正数,从而得到A1
的QR分解式
Ak1=(QxkDk)(DkkRxDkU)
由Ak
1的QR分解式的唯一性得到
⎧⎪⎨Q~
k=QxkDk⎪⎩R~k
=DkkkRxDU从而
A~TQ~TTk+1=QkAk=Dkk(QxAQx)kDk
由于A=XDX−1=QxRxDR−1TxQx,所以
QTR−1
xAQx=xDRx
从而AT−1k+1=Dkk(RxDRx)kDk
其中
⎡⎢λ1 ∗ ∗ ⋅⋅⋅ ∗⎤
RR−1 λ ∗ ⋅⋅⋅ ∗⎥
20=xDRx=⎢
⎢⎥⎢ ⋱ ⋮⎥⎣ λ⎥n⎦
于是
A~T
k+1=DkQkR0kDk
因为R0为上三角阵,limk→∞
k=I , Dk为对角阵,且元素为1或
-1,所以
(k)aii→λi,当k→∞
(k)aij→0,当i>j, k→∞(k)aij (j>i)的极限不一定存在。
例5用QR算法求矩阵
⎡ 5 −2 −5 −1⎤
⎢ 1 0 −3 2⎥
⎥A=⎢
⎢ 0 2 2 −3⎥⎢⎥⎣ 0 0 1 −2⎦
的特征值。A的特征值为−1 , 4 , 1+2i , 1−2i。
解令A1=A,用施密特正交化过程将A1分解为
A1=Q1R1
⎡0.9806 −0.0377 0.1923 −0.1038⎤⎡5.0992 −1.9612 −5.4912 −0.3922⎤⎢0.1961 0.1887 −0.8804 −0.4192⎥⎢ 0 2.0381 1.5852 −2.5288⎥
⎥×⎢⎥=⎢
⎢ 0 0.9813 0.1761 0.0740⎥⎢ 0 0 2.5242 −3.2736⎥⎢⎥⎢⎥ 0 0 0.3962 −0.8989 0 0 0 0.7822⎣⎦⎣⎦
将Q1与R1逆序相乘,求出A2
⎡4.6517 5.9508 1.5299 0.2390
⎢0.3997 1.9401 −2.5171 1.5361
A2=R1Q1=⎢
⎢ 0 2.4770 −0.8525 3.1294⎢
⎣ 0 0 0.3099 −0.7031
用代替A1重复上面过程,计算11次得
⎤⎥⎥⎥⎥⎦
A12
⎡4.0000 ∗ ∗ ∗⎢ 0 1.8789 −3.5910 ∗=⎢
⎢ 0 1.3290 0.1211 ∗⎢
⎣ 0 0 0 −1.0000
⎤⎥⎥⎥⎥⎦
由A12不难看出,矩阵A的一个特征值是4,另一个特征值是-1,其他两个特征值是方程
.8789−λ −3.5910
.3290 0.1211−λ
=0
的根。求得为1±2i。