数值分析矩阵特征值特征向量计算

合集下载

数值分析—矩阵特征值问题的数值计算

( vk +1 )i vk = α1 x1 , lim = λ1 . k k →∞ λ k →∞ ( v ) k i 1
lim
可见，当 k 充分大时， vk 近似于主特征向量（相差一个常数倍）， vk +1 与 vk 的对应非零分量的比值近似于主特征值。在实际计算中，需要对计算结果进行规范化。因为当 λ1 <1 时， vk 趋于零；当 λ1 >1 时， vk 的非零分量趋于无穷，从而计算时会出现下溢或上溢。为此，对向量
λ1 ≥ λ2 ≥ L ≥ λn ,
∑x y
i =1 i
n
i
为向量 x 和 y 的内积。
定理 8.3 设 A 为 n 阶实对称矩阵，其特征值都为实数，排列为对应的特征向量 x1 , x2 ,L , xn 组成正交向量组，则有 1) 对任何非零向量 x ∈ R n ，有 λn ≤ R( x) ≤ λ1 ， 2) λ1 = max R ( x) = R( x1 ) ，
Ζ = ( z1 , z2 ,L, zn )T ∈ R n ，记 max( Ζ) = zi ，其中 zi = Ζ ∞ ，这样，我们有如下乘幂法的实用
的计算公式：任取 v0 = u0 ≠ 0 ，对于 k = 1, 2,L 分别计算 vk = Auk −1 , uk = vk / max(vk ). 求出对应矩阵的主特征向量和特征值的近似值，有下面的定理。定理 8.4
m1 0 M = 0 0 0
称为质量矩阵，而
0 m2 0 0 0
0 0 m3 0 0
0 0 0 m4 0
0 0 0 0 m5 0 0 −k4 k 4 + k5 − k5 0 0 0 − k5 k5

数值分析第3章矩阵特征值与特征向量的计算

第三章矩阵特征值和特征向量计算
工程实践中有许多问题，如桥梁或建筑物的振动，机械
机件、飞机机翼的振动，及一些稳定性分析和相关分析可转化为求矩阵特征值与特征向量的问题。
设A (aij )nn是n阶方阵, 如果数和 n 维非零向量x满足 Ax x，则称为 A 的一个特征值， x称为矩阵A对应于的特征向量。
（2）在正交相似变换下，矩阵元素的平方和不变。设A (aij )nn ,
n
n
U为正交矩阵，记B U T AU (bij )nn ,则
ai2j
bi2j
i, j1
i, j1
Jacobi方法的基本思路：通过一次正交变换,将A中一对非
零的非对角元素化成零,并且使得非对角元素的平方和减少。
反复进行上述过程，使变换后的矩阵的非对角元素的平方和
2a pq
(
/
4)则有a
(1) pq
a(1) qp
1
0 1 0.93
1 0 0 0.93 1
0
0
1
0
0
0.93
1
0 1/ 0.93 1 0
0 0.93 1/ 0.93
按算法迭代3次， 3.0000954,与准确值3的误差小于10－4，u
(1, 0.9992431, 0.9991478)T 与准确值(1, -1,1)T比较，残差 r 0.001.
空间中的二维坐标旋转矩阵。
坐标旋转矩阵U pq ( )是正交矩阵.
设A为实对称矩阵，且apq aqp 0,若记
A(1)
U
T pq
AU
pq
(ai(j1) )
aaq((p1q1p))
a pp a pp
cos2 aqq sin2 sin2 aqq cos2

数值分析--第四章--特征值特征向量计算（乘幂法）

数值分析--第四章--特征值特征向量计算（乘幂法）
摘要：n阶⽅阵A满⾜AX=λx，λ为矩阵A的特征值，x为特征值对应的特征向量。

⼀.乘幂法（求模最⼤特征值及对应特征向量）
设矩阵A有n个相性⽆关的特征向量x1,x2,x3,.....xn，相应的特征值λ1,λ2,λ3,.....λn（由⼤到⼩排列）。

迭代法引⼊：上⼀章学了迭代法求解线性⽅程组Ax=b的解，给定任⼀的初始值v0，不断迭代可以得到Ax=b的解。

同理，给定任⼀⾮零的n维向量v0，不断迭代可以得到矩阵A的特征向量，
对于初始向量v0可以由A的n个线性⽆关的特征向量表⽰：
带⼊迭代⽅程中：
当迭代次数k趋近于⽆穷⼤时，可得到最⼤特征值λ1对应的特征向量a1x1(与x1线性相关)
同理，当迭代次数趋近于⽆穷⼤时，可得到绝对值最⼤的特征值，λ1
其中，m表⽰向量中的绝对值最⼤的那个元素值
如何利⽤迭代法求解按模最⼤特征值和特征向量
说明：
1.初始值：我们给定初始值x0=[1,1,1]^T,取特征值1
2.第⼀次迭代：
对应的近似特征值取：
3.第⼆次迭代：
⼆.改进乘幂法
这个规范化处理的⽬的：防⽌数据溢出或是数据消失
从上⾯可以看出，改进乘幂法即是每次迭代出的特征向量都进⾏⼀次规范化处理。

数值分析QR方法求矩阵特征值和特征向量

四.实验代码: function [H,B]=Hessenberg(A) n=length(A);B=eye(n);for k=1:n-2X=zeros(n-k,1);H=eye(n);for i=1:n-kX(i)=A(i+k,k);enda=max(abs(X));if a==0.0breakendX=X/a;c=X(1);b1=sqrt(sum(X.^2));if X(1)>=0b1=-b1;endX(1)=X(1)-b1;b=b1^2-b1*c;H0=eye(n-k)-X*X'/b;for i=1:n-kfor j=1:n-kH(i+k,j+k)=H0(i,j);end endA=H*A*H;B=B*H;endH=A;一.实验题目:QR方法求矩阵的特征和特征向量二.设计目的:学会利用镜面变换进行矩阵的QR分解及利用将幂法求特征值和特征向量,熟悉Matlab编程环境。

三.设计原理:利用镜像变换将A相似变换为Hessenberg B矩阵。

记录变换矩阵。

运用Householder矩阵进行QR分解,QR方法为:B1=BB1=Q1R1B2=R1Q1....Bm=QmRmBm+1=RmQmBm+1与Bm相似,从而特征值相等。

再利用原点位移的反幂法求B(或A)的特征向量。

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

设A∈R n×n为非奇异矩阵,A 的特征值依次记为|λ1|≥|λ2|≥|λ3|≥…≥|λn |,相应的特征向量为x1 ,x2,…,x n,则A-1的特征值为|1/λn|≥|1/λn-1|≥…≥|1/λ1 | ,相应的特征向量为x n ,x.所以计算A的按模最小的特征值λn的问题就是计算n-1,…,x1A-1的按模最大的特征值问题。

对于A-1应用幂法迭代(称为反幂法),可求得矩阵A-1的主特征值1/λn,从而求得A的按模最小的特征值λn。

西安科技大学研究生数值分析课件7章矩阵特征值与特征向量计算

7 矩阵特征值与特征向量地计算设A 为n 阶方阵,所谓A 地特征值问题是求数λ和非零向量x ,使x Ax λ=成立.数λ称作A 地一个特征值,非零向量x 称作与特征值λ对应地特征向量.求给定方阵地特征值与特征向量是先求解特征方程()||0E A ϕλλ=-=然后对应于每一个特征值i λ,再求解退化地齐次线性方程组()0i E A x λ-=从而得到A 地特征值i λ及对应地特征向量x .但是这种方法计算机很大,计算过程复杂,因此有必要研究相对简单地数值解法.本章主要介绍三类计算特征值地方法：计算大型（稀疏）矩阵主特征地幂法与反幂法,计算中小型（实对称）矩阵全部特征值地Jacobi 法,计算中小型矩阵全部特征值地QR 法.7.1 特征值估计在矩阵特征值计算中,有时需要对特征值所在范围给出一个估计.这里介绍一种从矩阵地元素出发,运用较简便地运算估计特征值地方法.定义7-1 设()n m ij A a C ⨯=∈,称由不等式||ii i z a R -≤在复平面上确定地区域为矩阵A 地第i 个盖尔圆（Gerschgorin 圆）,并用i G 表示.其中1||ni ij j j i R a =≠=∑称为盖尔圆i G 地半径(1,2,,)i n =.定理7-1 矩阵()n m ij A a C ⨯=∈地一切特征值均落在它地n 个盖尔圆地并集中,即1(1,2,,)ni jj G i n λ=∈=.证明设λ是A 地任一特征值,12(,,,)T n x x x x =是λ对应地特征向量.令01||max ||i i i nx x ≤≤=,则00i x ≠.由Ax x λ=,可得001()ni j j i j a x x λ==∑.即∑≠==-ni j j j j i i i i x a x a 000001)(λ于是有 000000011i i jni j j ji ni j j i jji i i R x x ax x aa ≤≤=-∑∑≠=≠=λ这表明任一特征值0i G λ∈,从而也在A 地第n 个盖尔圆地并集中.例7-1 估计矩阵10.10.20.30.530.10.210.310.50.20.30.14A ⎡⎤⎢⎥⎢⎥=⎢⎥-⎢⎥---⎣⎦地特征值范围. 解 A 地4个盖尔圆为：1:|1|0.6G z -≤ 2:|3|0.8G z -≤ 3:|1| 1.8G z +≤ 4:|4|0.6G z +≤画在复平面上其区域如图7-1所示.图7-1 例7-1盖尔圆分布图于是A 地全部特征值就在这4个盖尔圆地并集中.为了更确切地知道某个特征值落在哪个或哪几个盖尔圆地并集中,给出如下第二盖尔圆盘定理.定理7-2 若A 地n 个盖尔圆中,有m 个盖尔圆构成地一个连通域（所谓连通域,是指其中地任意两点都可以用位于该区域内地一条折线连接起来）,且该连通域与其余n m -个盖尔圆严格分离,则在该连通域中恰好有A 地m 个特征值（重特征值按重数重复计算）.特别地,每个孤立地盖尔圆恰有A 地一个特征值（证明从略）.由定理2可知,在例1中2G 与4G 中各有A 地一个特征值,而1G 与3G 构成地连通部分中有两个特征值,但不能确定这两个特征值具体落在哪个盖尔圆中.例7-2 估计矩阵10.80.50A -⎡⎤=⎢⎥⎣⎦地特征值范围. 解 A 地两个盖尔圆为：1:|1|0.8G z -≤,2:|0|0.5G z -≤在复平面上地区域如图7-2所示.图7-2 例7-2盖尔圆分布图此时只能判断A 地两个特征值落在1G 与2G 地并集中,至于是每个盖尔圆中各有一个特征值还是两个特征值都落在其中一个盖尔圆上则无法确定.实际上,由于1,21(12λ=±,1,2||0.5λ=>,所以两个特征值都不会在盖尔圆2G 中,而是落在盖尔圆1G 中.对于某些矩阵,可利用相似变换矩阵具有相同特征值地性质得到更确切地特征值范围.设()ij n m A a ⨯=,取正数12,,,n d d d 构成对角阵12diag(,,,)n D d d d =,对A 作相似变换,令1()iij n n jd B DAD a d -⨯==,由于B 相似于A ,所以B 与A 地特征值完全相同,又由于B 与A 地主对角线元素对应相等,所以B 与A 地盖尔圆圆心相同.这表明,若适当选取正数12,,,n d d d ,可以改变盖尔圆地半径,从而有可能将相交地盖尔圆分离得到仅含一个特征值地孤立盖尔圆.选取12,,,n d d d 地一般方法是：欲使A 地第i 个盖尔圆i G 地半径大而其余盖尔圆变小,就取1i d >,其余1()j d j i =≠.例7-3 求矩阵2050.841011210A j ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地特征值范围. 解 A 地3个盖尔圆为：1:|20| 5.8G z -≤,2:|10|5G z -≤,3:|10|3G z j -≤其中1G 与2G 相交,而3G 孤立.记3G 中所含地一个特征值为3λ,如图7-3所示.为分离2G 与1G ,可以让A 地第3行元素绝对值变大,第3列元素绝对值变小.现取diag(1,1,2)D =,则12050.44100.52410B DAD j -⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦图7-3 例3盖尔圆分布图图7-4 例7-3分离后盖尔圆分布图其3个盖尔圆分别是：1:|20| 5.4G z '-≤,2:|10| 4.5G z '-≤,3:|10|6G z j '-≤ 显然,B 地盖尔圆是3个孤立地盖尔圆,如图7-4,注意,此情况下,3G '地半径变大了.例7-4 设矩阵()ij n n A a ⨯=按行严格对角占优,则A 可逆.证明由线性代数知,A 可逆地充分条件是||0A ≠,而1||nj j A λ==∏（其中j λ是A 地特征值）,所以只要证明0j λ≠即可(1,2,,)j n =. 设λ是A 地任一特征值,则必存在某个盖尔圆i G 使∑≠=≤-ij ij i ii a R a λ.若0j λ=,则有∑≠≤ij ij ii a a ,而这与A 按行严格对角占优矛盾,故应有0λ≠,由λ地任意性,得||0A ≠.7.2 幂法与反幂法在线性代数中,设A 是n 阶方阵,若A 存在n 个线性无关地特征向量,则称这n 个特征向量构成A 地一个完全地特征向量组.例如,对矩阵320230005A -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,110430102B -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦通过求解特征方程,不难求出A 地三个特征值为1231,5λλλ===,B地三个特征值为1232,1λλλ===.方阵A 可以找到三个线性无关地特征向量,而方阵B 找不到三个线性无关地特征向量.我们称方阵A 可对角化,而B 不可对角化. 7.2.1 幂法幂法地基本思想是构造一个向量序列使之逼近主特征值（矩阵地按模最大地特征值）对应地特征向量,然后求出主特征值.该方法简单易行,但收敛速度较慢.现设()ij n n A a ⨯=有一个完全地特征向量组12,,,n x x x ,其对应地特征值是12,,,n λλλ.已知A 地主特征值是单根1λ,即特征值满足条件12||||||n λλλ>≥≥任取一个非零初始向量0u ,由矩阵A 构造向量序列102210110k k k u Au u Au A u u Au A u++=⎧⎪==⎪⎪⎨⎪==⎪⎪⎩由于A 地完全特征向量组可以作为向量空间n R 地一组基,因此0u 可由12,,,n x x x 线性表示,即有01122n n u a x a x a x =+++ （设10a ≠）于是011122211111121()()k k k k k n n nn kk k i i i k i u A u a x a x a x a x a x a x λλλλλλελ===+++⎡⎤=+=+⎢⎥⎣⎦∑ 其中21()nk i k i i i a x λελ==∑.注意到),,2(11n i i=<λλ,故当k →∞时,0k ε→,因此有111k k u a x λ≈由于1x 是主特征值1λ对应地特征向量,其乘上常数因子11k a λ仍为1λ地特征向量,故当k 充分大时,迭代向量k u 是1λ地特征向量地近似向量.为了利用迭代向量求出主特征值1λ地近似值,设()k i u 表示k u 地第i 个分量,则1111111()()()[]()()()k i i k ik i i k iu a x u a x ελε+++=+ 于是 11()lim()k ik k iu u λ+→∞= 这表明两相邻迭代向量对应分量地比值收敛于主特征值,亦即当k 充分大时,可用两相邻迭代向量地分量比作为主特征值地近似值,即11()()k ik iu u λ+≈若主特征值是A 地r 重实特征值,即12(1)r r n λλλ===≤≤,对应地r 个线性无关特征向量为12,,,n x x x .则有01111()r nkk k i k i i i i i i r u A u a x a x λλλ==+⎡⎤==+⎢⎥⎣⎦∑∑当k 充分大时,11rkk i i i u a x λ=≈∑即k u 仍为主特征值对应地特征向量地近似向量,相邻两迭代向量地分量比仍为主特征值地近似值.综上所述,有定理7-3 设A 是n 阶实矩阵,具有完全地特征向量组,主特征值是r 重根,即112||||||||(1)r r n r n λλλλ++>≥≥≥≤≤则对任意非零初始向量0u ,迭代向量0k k u A u = 满足 111lim(0)rki ikk i u a x a λ→∞==≠∑ ,11()lim ()k ik k iu u λ+→∞= 或 11rk k i i i u a x λ=≈∑,11()()k ik iu u λ+≈ 这样用非零初始向量0u 及矩阵A 构造向量序列{}k u 以计算A 地主特征值1λ及相应地特征向量地方法称为幂法.不过从上面地讨论中可以看到,如果1||1λ>或11<λ,迭代向量k u 当k →∞时,其不为零地分量就会趋于无穷大或趋于零.为克服这个缺点,可以在每步迭代中加上对向量规范化地步骤,使迭代向量地数量级保持在一个稳定地量级上,归纳起来,幂法地计算步骤为：步骤 1 给定非零初始向量0u ,精度12,εε,令00v u =；令(0)10max()v λ=,1=k ；步骤 2 迭代1-=k k Av u ,()1max()k k u λ=,其中)max(k u 表示k u 绝对值最大地分量；步骤3 规范化max()kk k u v u =；步骤 4 若11k k v v ε--<且()(1)112||k k λλε--<,则k v 即为1λ地近似特征向量,()1k λ即为1λ地近似值；否则,1+=k k ,转步骤2继续迭代.例7-5 用幂法计算1.0 1.00.51.0 1.00.250.50.252.0A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地主特征值和相应地特征向量,结果见表7-1.表7-1而此题地准确值为1 2.5365258λ= 1(0.748221,0.649661,1.000000)T x =7.2.2 幂法地加速幂法地收敛速度由比值21r λλ=来确定,r 越小收敛越快,而当1r ≈时收敛可能很慢.为了克服这一缺点,常采用原点平移法对幂法进行加速.设B A pE =-,其中p 是待定参数.显然,若A 地特征值为12,,,n λλλ,则B 地相应特征值(1,2,,)i k i n =为12,,,n p p p λλλ---,且A .B 地特征向量相同.这是因为对A 有特征方程||0i A E λ-=,而对B 有特征方程|||()|0i i B k E A p k E -=-+=,所以,i i i i p k k p λλ=+=-另一方面,若i x 是A 地对应i λ地特征向量,即i i i Ax x λ=则 ()()i i i i i i Bx A pE x Ax px p x λ=-=-=-原点平移法地思想是引入矩阵B ,恰当地选择参数p ,使11k p λ=-是B 地主特征值,且其速比2211maxB A p r r p λλλλ-=<=-,这样用幂法求B 地主特征值1k 地收敛速度就快于用幂法求A 地主特征值1λ,而一旦1k 求出,由11k p λ+=可得A 地主特征值,达到了加速地目地.但是为了选取恰当地选择参数p ,需要对A 地特征值地分布地大致了解. 例7-6 设4阶方阵A 有特征值15(1,2,3,4)j jj λ=-=其速比210.9A r λλ=≈.作变换 (12)B A pEp =-=则B 地特征值为12k =,21k =,30k =,41k =-,其速比2112B A k r r k ==<. 设A 地实特征值满足121n n λλλλ->≥≥>若2,n λλ地值可大致估计出,若要求1λ,考察待定参数p 地选取. 在原点平移法通过变换pE A B -=后,不论p 如何选取,矩阵地B 主特征值也只能是在n p λ-或 1p λ-.若希望求1λ,就应选择p ,使1p λ-称为B 地主特征值,即1||||n p p λλ->-这时B 地收敛速比B r 是比值21||/||p p λλ--和1||/||n p p λλ--中地较大者,即211||||max ,||||n B p p r p p λλλλ⎧⎫--=⎨⎬--⎩⎭显然B r 依赖于p 地选取,记做()B r p .为了使应用幂法求B 地主特征值地收敛速度尽可能快,我们希望选择最佳参数*p ,使*()min ()B B r p r p =由B r 地表示式（求二者之间地较大值）和)(*p r B 对)(p r B 地最小化要求,只有当2||||n p p λλ-=-时,()B r p 达到最小.由于2n p p λλ-=-会有得到矛盾地结果（2n λλ=）,所以只能是2()n p p λλ-=--即 *22np λλ+=类似地,若用反幂法求最小特征值n λ,若1n λ-,1λ可大致估计,取最佳平移参数*112n p λλ-+=例7-7 取0.75p =,用原点平移法,计算例7-7中矩阵A 地主特征值.解作变换B A pE =-,则0.2510.510.250.250.50.25 1.25B ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦对B 应用幂法,计算结果见表7-2.即1 1.7865914k ≈,则A 地主特征值1λ为110.75 2.5365914k λ=+=与例7-5比较,上述结果比例7-5迭代15次还好.表7-27.2.3 反幂法设方阵A 按模最小地特征值是n λ,且0n λ≠,则A 可逆.于是,由n n n Ax x λ=,可得11n n nA x x λ-=,这表明1nλ是1A -地主特征值.反幂法就是将幂法应用于1A -,通过求出1A -地主特征值得到A 地按模最小地特征值及其对应地特征向量.定理7-4 设A 是n 阶实矩阵,具有完全地特征向量组,其特征值满足12||||||0n λλλ≥≥≥>则对任意非零初始向量00u v =,按下述方式构造地迭代向量11k k u A v --= ,max()kk k u v u =满足lim max()n k k n x v x →∞=, 1lim max()k k nu λ→∞= /max()k n n v x x ≈,1max()k nu λ≈在实际计算中,可先对A 进行LU 分解,通过求解1k k Ly v -= ,k k Uu y =来求解方程组1k k Au v -=.反幂法地计算步骤为：步骤1 预先取定非零向量00u v =；给定精度12,εε；取(0)0m a x ()nu μ=；步骤2 对矩阵A 作LU 分解,A LU =；令1=k ；步骤3 求解方程组1k k Ly v -= ,k k Uu y = 得到迭代向量k u ；步骤4 规范化max()kk k u v u =步骤5 若11k k v v ε--<且()(1)2||k k n n μμε--<,则k v 即为A 地对应于n λ地近似特征向量,()1k nμ即为n λ地近似值；否则,令1+=k k ,转步骤3继续迭代.7.3 矩阵地两种正交变换本节先介绍镜面(初等)反射变换和平面旋转变换,它们是QR 算法和Jacobi 算法地基础.7.3.1 豪斯荷尔德（House holder ）变换定义7-2 设有方阵B ,若当1i j >+时(,1,2,,)i j n =,0ij b =,则称B 是上Hessenberg 矩阵,即1112121222,1n n n n nn b b b b b b B b b -⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦定义7-3 设向量ω满足21ω=,矩阵2T H E ωω=- （ω是列向量）称为初等反射矩阵,又称House holder 矩阵,记为()H ω,即211212212221212222122()2212n n n n n H ωωωωωωωωωωωωωωωω⎡⎤---⎢⎥---⎢⎥=⎢⎥⎢⎥---⎢⎥⎣⎦其中(1,2,,)i i n ω=是ω地分量.可以证明初等反射阵是对称阵()T H H =.正交阵()T H H E =. 例7-8 设向量0α≠,试证矩阵222TH E ααα=- 是一个初等反射阵. 证明令2αωα=,则 222221||||||||1αωααα=== 由定义7-3,2222TTH E E ααωωα=-=-是初等反射阵.定理7-5 设,x y 是两个不相等地n 维列向量,且22||||||||x y =,则存在一个初等反射阵H,使得Hx y =证明令2||||x yx y ω-=-,由例7-8可知22()()22||||T T Tx y x y H E E x y ωω--=-=-- 是一个初等反射阵.由于22||||()()T T T T Tx y x y x y x x y x x y y y -=--=--+ 注意到22||||||||x y =,即T T x x y y =,又()T T T T x y x y y x == ,故22||||2()T Tx y x x y x -=-从而22()()2||||T T x y x x y x Hx x x y --=--y y x x =--=)(. 例7-9 设1(1,2,2),(1,0,0)T T x e ==,用Householder 变换将x 化为与1e 同方向地向量.解因为2||||3x =,可设13y e =,则22||||||||x y = 取21,1,1)||||T x y w x y -==--,构造Householder 矩阵[]11122212111,1,12123311221T H E ww -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=-=--=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦则13Hx e =推论设向量12(,,,)0T n x x x x =≠,12()||||r sign x x =,且1x r ε≠-,则存在初等反射阵1222||||T T uu H E E uu u ρ-=-=- 使1Hx r ε=- .其中,1(1,0,,0)T ε=,1u x r ε=+,22||||/2u ρ=.设12(,,,)T n u u u u =,则12(,,,)T n u x r x x =+22222122222112111||||[()]221(2)2()n n u x r x x r rx x x x r r x ρ==++++=+++++=+引入初等反射阵地目地,是设法用一系列初等反射阵将原始矩阵约化成上Hessenberg 阵.由于约化过程是逐列进行地,我们先给出计算Hx 地算法步骤,该算法算出H 及r ,使Hx r ε=-,u 地分量冲掉x 地分量.（1）1max ||i i nx η≤≤=；（2）(1,2,,)ii i x x u i n η←==,此步规范化是为避免计算r 时产生溢出；（3） 12211()()nii r sign x x ==∑；（4）11u u r ←+；（5） 1ru ρ=；（6） r r η←；于是初等反射阵1T H E uu ρ-=-,1Hx r ε=-.如果要将H 作用于矩阵A ,设i a 是A 地第i 列向量,则12(,,,)n A a a a =,12(,,,)n HA Ha Ha Ha = 其中,11()()(1,2,,)T T i i i i Ha E uu a a u a ui n ρρ--=-=-=.下面讨论用初等反射阵约化原始矩阵A 称为上Hessenberg 阵地步骤.11121(1)(1)2122211121(1)(1)212212n n n n nn a a a a a a a A A A a a a a a ⎡⎤⎢⎥⎡⎤⎢⎥===⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦步骤1 不妨设(1)210a ≠（否则这一步不需约化）,选择初等反射阵1R ,使(1)12111R a r ε=-,其中： 1(1)(1)2212112(1)1211112111121211111()(())(1)1()2ni i T r sign a a u a r n u r r a R E u u εερρ=-⎧=⎪⎪⎪=+-⎨⎪==+⎪⎪=-⎩∑是维单位坐标列向量令11100U R ⎡⎤=⎢⎥⎣⎦则(2)(2)(2)(1)111213111212111(2)(2)(1)(1)222312112210A a A a A R A U AU a A R a R A R ⎡⎤⎡⎤===⎢⎥⎢⎥⎣⎦⎣⎦其中,(2)11A 是21⨯阵,(2)22a 是2n -维列向量,(2)23A 是2n -阶方阵.步骤k 设对A 已进行了1k -步约化,即111(2)()()()()11121,111,11(2)()()()1222,12,2()()()1,1,()()()1,1,11,()()(),1(2,3,,1)k k k k k k k k k k k n k k k k kn k k k k kk k k k n k k k k k k k k nk k k nkn k nnA U A U k n a a a a a a r a a a a r a a a a a a a a a ----+--++++++==-⎡⎢-⎢=-⎣()()()111213()()22230k k k k k A a A a A ⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦⎡⎤=⎢⎥⎣⎦其中,()11k A 是(1)k k ⨯-阵,()22k a 是n k -维列向量,()23k A 是n k -阶方阵.设()220k a ≠,选初等反射阵()k R n k -阶,使()221k k k R a r ε=-,其中1ε是n k -维单位坐标向量,可得1()()221,1()221()1,1()(())()nk k k k k ik i k k kk k k k kk nT k k k k r sign a a u a r r r a R E u u ερρ+=++-⎧=⎪⎪⎪=+⎨⎪=+⎪⎪=-⎩∑ 令 00k k E U R ⎡⎤=⎢⎥⎣⎦则 ()()()1112131()()2223()()()111213()12300k k k k k k k k k k k k k k k k k k k k k A a A R A U A U R a R A R A a A R r R A R ε+⎡⎤==⎢⎥⎣⎦⎡⎤=⎢⎥-⎣⎦ 可见1k A +地左上角1k +阶子阵为上Hessenberg 阵,从而约化又进了一步.重复此过程,直到122112211(2)122(3)233(1)1n n n n n nn A U U U AU U U a r a r a r a -----=⨯⨯⨯⎡⎤⎢⎥-⨯⨯⎢⎥⎢⎥=-⨯⎢⎥⎢⎥⎢⎥-⎣⎦使原始矩阵A 在一系列初等反射阵地作用下,约化为上Hessenberg 阵.综上所述,有定理7-6.定理7-6 如果A 是n 阶实矩阵,则存在初等反射阵122,,,n U U U -,使221122n n U U U AU U U C --=（上Hessenberg 阵）例7-10 试证矩阵A 与其约化成为地上Hessenberg 阵C 有相同地特征值.证明记221n P U U U -=,由于初等反射阵是正交对称阵,故122T n P U U U -=,且P 是正交阵,故T PAP C =.于是||||||||||||T T C E PAP E P A E P A E λλλλ-=-=-=-其中T PP E =,||||1T P P =.这表明A 与C 具有相同地特征多项式,即两者有相同地特征值.进一步,设y 是C 地对应于特征值λ地特征向量,即Cy y λ=,则有T PAP y y λ= ()()T T A P y P y λ=这表明T P y 为A 地对应于λ地特征向量,于是求原始矩阵A 地特征值与特征向量可转化为求上Hessenberg 阵C 地特征值和特征向量.定理7-7 若A 是实对称矩阵,则存在初等反射阵122,,,n U U U -使2211221112211()n n n n n U U U AU U U c b b c b C b b c ----⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦对称三对角阵证明由定理7-6,存在初等反射阵可使A 约化为上Hessenberg 阵C ,当A 是对称矩阵时,C 亦为对称阵,即T C C =,且T C 亦为上Hessenberg 阵,故C 是对称三对角阵.例7-11 用豪斯荷尔德方法将下述矩阵化为上Hessenberg 阵.1437232427A A ---⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦解（1）对1k =,确定变换阵111000U R ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,(1)2124a ⎡⎤=⎢⎥⎣⎦ 其中1R 为初等反射阵,使(1)121110R a r ⎡⎤=-⎢⎥⎣⎦(1)12124.472136r a ==≈(1)12111 6.472136244u a r ε⎡⎡⎤+=+=≈⎢⎢⎥⎣⎦⎣⎦11121()2)28.94427r r a ρ=+≈[]1111110 6.4721361 6.472136401428.944270.4472070.8944230.8944230.447216TR E u u ρ-=-⎡⎤⎡⎤=-⎢⎥⎢⎥⎣⎦⎣⎦--⎡⎤=⎢⎥-⎣⎦（2）计算(1)122R A .记(1)221232(,)27A a a ⎡⎤==⎢⎥⎣⎦,于是 (1)1221112 3.1304967.155419(,) 1.788855 1.341640R A R a R a --⎡⎤==⎢⎥-⎣⎦其中,111111111()()(1,2)T T i i i i R a E u u a a u a u i ρρ--=-=-=（3）计算(1)121A R 及(1)1221()R A R ,即求 1(1)121211(1)1223373.1304967.1554191.788855 1.341640T T T b A R b R R R A b ⎡⎤--⎡⎤⎡⎤⎢⎥⎢⎥==--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦⎣⎦7.6026340.4472127.800030.3999990.399999 2.200000-⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦其中,11111()(1,2,3)T T T Ti i i b R b b u u i ρ-=-=（4）计算2111A U AU =.(1)12121(1)1221447.6026340.4472124.4721367.8000030.39999900.399999 2.2000000A R A r R A R ⎡⎤--⎡⎤⎢⎥⎢⎥⎢⎥==--⎢⎥⎢⎥-⎢⎥-⎢⎥⎣⎦⎢⎥⎣⎦为上Hessenberg 阵.7.3.2 平面旋转变换定义7-4 称矩阵111(,)111i j csi P i j scj ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第列第列第行第行为平面旋转矩阵,又称Givens 矩阵,其中cos c θ=,sin s θ=.平面旋转阵(,)P i j 是一个正交阵,与单位阵只有在(,),(,),(,i i i j j j j i四个位置上地元素不一样,用其左乘矩阵A 只改变A 地第i 行和第j 行元素.设12(,,,)T n x x x x =则平面旋转变换Px y =地结果为⎪⎩⎪⎨⎧≠=+-=+=ji k x y cx sx y sx cx y k kj i j j i i ,若令/i c x =,j s x =, 则平面旋转变换向量y 地第i个分两为22j i x x +,第j 个分量为0,其余分量即为x 对应地分量.和初等反射变换一样,用平面旋转变换也可以将一个方阵化为上Hessenberg 矩阵,也可以将将一个方阵化为上三角矩阵.7.4 QR 算法7.4.1 矩阵地QR 分解定理7-8 设A 是可逆矩阵,则存在正交矩阵121,,,n P P P -使121()n P P P A R -=上三角矩阵且R 地主对角元素0(1,2,,1)ii r i n >=-.证明若10(2,3,,)j a j n ==,则A 地第一列不需约化.若有某个 10(2)j a j n ≠≤≤,则可选择1(1,)j P j P =使A 地第一列中第j 个元素变为零.一般地,可设平面旋转矩阵12131,,,n P P P ,使(2)(2)11121(2)(2)222113122(2)(2)200nn nn nn r a a a a P P P A A a a ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦记111312nP P P P =,则12P A A =.同理,若(2)20(3,4,,)j a j n ≠=,可选取23242,,,n P P P 使(2)(2)(2)1112131(3)(3)22232(3)(3)2212323333(3)(3)3nn n n n n nn r a a a r a a P P P A A a a a a -⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦记2223nP P P =,则213P P A A =.重复上述过程,可得一系列正交阵121,,,n P P P -使11121222121n n n nn r r r r r P P P A R r -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦ 定理7-9 （矩阵地QR 分解）如果n 阶实矩阵A 可逆,则A 可分解为一正交阵Q 和上三角阵R 地乘积,即A QR =,且当R 地对角元素都为正数时分解唯一.证明由定理8知存在正交阵11,,n P P -使121n P P P A R -=为上三角阵,记121T n Q P P P -=,于是T Q A R =由于(1,2,,1)i P i n =-是正交阵,则T Q 亦为正交阵,故A QR =. 若A 有两种QR 分解,记为1122A Q R Q R ==其中12,R R 为上三角阵且主对角元素都为正数,12,Q Q 为正交阵,于是12121T Q Q R R -=注意121R R -是上三角阵地乘积,结果仍为上三角阵,而12,TQ Q 是正交阵,所以121R R -也应是正交阵.若记121D R R -=,由其上三角性T D 应是下三角阵,1D -应是上三角阵；由其正交性由1T D D -=,故D 只能是对角阵,且有2T D D D E ==.又因12,R R 地主对角元素都为正数,即有222212diag[,,,]diag[1,1,,1]n D d d d E ===故1(1,2,,)i d i n ==,则D E =,于是12R R =,12Q Q =.例7-12 求矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=212240130A 地QR 分解. 解方法1：利用初等反射阵进行QR 分解令(0)1(0,0,2)T a =,取(0)112||||2d a ==,则)2,0,2(81211)0(111)0(11-=--=e d ae d a u1110012010100TH E u u ⎡⎤⎢⎥=-=⎢⎥⎢⎥⎣⎦,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=1302402121A H 再令(0)2(4,3)T a =,取(0)222||||5d a ==,则(1)2212(2)22121,3)||||T a d e u a d e -==--2224312345TH E u u ⎡⎤=-=⎢⎥-⎣⎦令2210014305534055H H ⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦于是21212051002H H A R ⎡⎤⎢⎥=-=⎢⎥⎢⎥-⎣⎦故123405521243005155002100T TA H H R ⎡⎤-⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦方法2：利用平面旋转阵进行QR 分解. 取1202,0100221221=+==+=s c ,则130********T ⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦,132********T A ⎡⎤⎢⎥=-⎢⎥⎢⎥--⎣⎦再取53)3(43,54)3(44222222-=-+-==-+=s c ,则231004305534055T ⎡⎤⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎣⎦,2313212051002T T A R ⎡⎤⎢⎥=-=⎢⎥⎢⎥-⎣⎦ 故13233405521243005155002100T T A T T R ⎡⎤-⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦例7-13 求矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=110133044A 地QR 分解,使得R 地对角线元素为正数.解 A A =1地第一列T x ]0,3,4[1=,521=x .用1x 构造镜面反射阵1H ,使得T y x H ]0,0,5[111==,令T y x u ]0,3,1[111-=-=,有⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-=-=10005453053542221111u u u E H T ,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-==11054005355112A H A 2A 地第2列对角线以下为T x ]1,0[2=,122=x .用2x 构造镜面反射阵2~H ,使得T y x H ]0,1[~222==,令T y x u ]1,1[222-=-=,易得 ⎥⎦⎤⎢⎣⎡=-=01102~222222u u u E H T,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=010100001~122H H 于是有R A H A =⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-==54001105355333,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-==010540535305421H H Q容易验证,QR A =.请读者用平面旋转变换对本例地矩阵A 进行QR 分解.7.5.3 QR 算法QR 算法就是利用QR 分解构造一个矩阵序列{}k A ,当k 充分大时,k A 是近似地上三角矩阵,而该上三角阵地对角元素便是原始矩阵A 地全部特征值.设1()n n ij n n A A a R ⨯⨯==∈,对A 做QR 分解,即A QR =其中R 为上三角阵,Q 为正交阵.利用这个分解可得新矩阵（对QR 交换乘积）2T A RQ Q AQ == 由于2A 是1A 经过正交相似变换得到地,因此2A 与1A 有相同地特征值.再对2A 做QR 分解,按上述方式又可得新矩阵3A ,且3A 与2A 也具有相同地特征值.具体地说,其步骤为：设1A A =,做QR 分解111A Q R =求矩阵211111T A R Q Q A Q ==求得k A 后对k A 作QR 分解k k k A Q R =求矩阵1Tk k k k k k A R Q Q A Q +==只要A 可逆,由定理9可知,按上述方法可唯一确定矩阵序列{}k A ,且序列中任意k A 与原始矩阵有相同特征值.因此只要恰当选择正交相似变换阵12,,,k Q Q Q ,使1111111T T TT TT T k k k k k k k k k k k k k k A Q A Q Q Q A Q Q Q Q Q A Q Q Q +----====当k →∞时,逼近一个上三角阵,便可求出A 地全部特征值（为所逼近上三角阵地主对角元素）.可见,QR 算法地关键在于选择正交变换阵(1,2,)k Q k =.从定理7-8地证明看到,正交变换阵k Q 是一系列平面转换矩阵地乘积,这些平面旋转矩阵是用来将k A 地主对角线以下元素约化为零地.如果将QR 算法直接应用于原始矩阵,计算量很大,所以在实际计算中,总是先将原始矩阵用豪斯赫尔德方法约化为上Hessenberg 阵,而后再对上Hessenberg 阵应用QR 算法.可以证明,由上Hessenberg 阵用QR 算法生成地矩阵序列中地每个矩阵仍为上Hessenberg 阵.7.5 雅可比方法雅可比方法是用来计算实对称矩阵地全部特征值及特征向量地一种有效方法.它地基本思想是,通过一组正交相似变换对称矩阵A 化为对角矩阵,得其全部特征值.定理7-10 设A 为n 阶对称矩阵,T C PAP =,其中P 为正交矩阵,则22||||||||F F C A = 证明一方面2222111||||()()nnnFiji i j i A a tr A A λ======∑∑∑另一方面2221||||()()()nTFi i C tr C C tr C C λ====∑由假设()()i i A C λλ=,故22||||||||F F C A =.设n n A R ⨯∈为对称矩阵,(,)P i j 为一平面旋转矩阵,则T C PAP =（其中()ij n n C c ⨯=）地元素计算公式为：（1）22cos sin 2sin cos ii ii jj ij c a a a θθθθ=++22sin cos 2sin cos jj ii jj ij c a a a θθθθ=+-（2）1()sin 2cos 22ij ji jj ii ij c c a a a θθ==-+ （3）第i 行元素和第j 列元素cos sin (,)ik ki ik jk c c a a k i j θθ==+≠ （4）第j 行元素和第i 列元素 cos sin (,)jk kj jk ik c c a a k i j θθ==+≠（5）(,,)lk lkc a l k i j =≠这说明,当A 经过一初等正交相似变换化为C 时,只需按上述公式计算C 地第i 列.第j 列元素,由对称性可得第i 行和第j 行元素,C 地其余元素与A 地对应元素相同.设A 地非对角元素0ij a ≠,我们可选择平面旋转阵(,)P i j ,使T C PAP =地非对角元素0ij ji c c ==.由定理11可选择(,)P i j ,使sin 2cos 202jj iiij ji ij a a c c a θθ-==+=即选择θ,使22(||)4ij ii jja tg a a πθθ=≤-其中定理7-11 设n n A R ⨯∈为对称阵,0ij a ≠为A 地一个非对角元素,则可选择一平面旋转阵(,)P i j ,使T C PAP =地非对角元素0ij ji c c ==且T C PAP =与A 地元素满足下述关系（1）2222(,)ik jk ik jkc c a a k i j +=+≠（2）222222ii jj ii jj ij c c a a a +=++ （3）22(,,)iklk c a l k i j =≠证明由上面地计算ij c 公式直接计算可知（1）成立.由（1）及定理7-10可证（2）.如果用()S A 表示A 地非对角线元素地平方和,()D A 表示A 地对角线元素平方和,则2()()2ijD C D A a =+ ,2()()2ij S C S A a =- 这说明C 地对角线元素平方和比A 地对角线元素平方和增加了22ij a ,C 地非对角线元素平方和比A 地非对角线元素平方和减少了22ij a .下面介绍雅可比方法.首先在A 地非对角元素中选择绝对值最大地元素（称为主元素）,如11||max ||i j lk l ka a ≠=可设110i j a ≠,否则A 已经对角化了.由定理12,选择一平面旋转矩阵111(,)P i j ,使111TAP AP =地非对角元素11110i j j i c c ==. 再选(1)1()lkn n A a ⨯=地非对角元素中地主元素,如 22(1)(1)||max ||0i j lk l ka a ≠=≠由定理12,又可选择一平面旋转矩阵222(,)P i j ,使2212T A P A P =地非对角元素2222(2)(2)i j j i a a ==（注意上次消除了地主元素这次又可能变为不是零）. 继续这个过程,连续对A 实行一系列平面旋转变换,消去非对角线绝对值最大地元素,直到将A 地非对角元素全化为充分小为止,从而求得A 地全部（近似）特征值.定理7-12 （雅可比方法地收敛性）设()ij n n A a ⨯=为实对称矩阵,对A 施行上述一系列平面旋转变换1(1,2,)Tm m m mA P A P m -==则 lim ()m m A D→∞=对角矩阵证明记()()m m lk n n A a ⨯=,()2()m m lk l kS a ≠=∑由定理7-11地（2）可得()212()m m m ij S S a +=-其中 ()()||max ||m m ijlk l ka a ≠= 又由于()2()2()(1)()m m m lk ij l kS a n n a ≠=≤-∑即()2()(1)m m ij S a n n ≤- 由以上得12(1)(1)m m S S n n +≤-- 反复应用上式,即得1102(1)(2)(1)m m S S n n n ++≤->-故 lim 0m m S →∞= 可以证明()lim m ll m a →∞存在(1,2,,)l n =. 下面介绍特征向量地计算.由雅可比收敛定理知,当m 充分大时2112T TTmm P P P AP P P D ≈记12T T T T m m R P P P =,则T m R 地列向量就是A 地近似特征向量.计算Tm R 可采用累积地办法,用一数组R 保存Tm R ,开始时R E ←,以后对A 每进行一次平面旋转变换,就进行计算Tm R RP ←用初等正交阵T m P 右乘R 只需计算R 地两列元素,若记(,)m m P P i j =,则Tm RP 地计算公式为()()cos ()sin (1,,)()()sin ()cos li li lj li li lj l n θθθθ←+⎧⎪=⎨←+⎪⎩R R R R R R关于sin θ和cos θ地计算如下.由定理7-11知,当0ij a ≠时,可选θ满足2tg2ij ii jja a a θ=-方ii jj a a ≠时,由22tg 1tg21tg dθθθ=≡- 得到tg θ地二次方程2tg 2tg 10d θθ+-=解得tg θ=选取tg 0d d θ>=<由此得 |tg |1θ≤可由集合{},,ii jj ij a a a 来计算sin ,cos θθ,设0,||max ||ij ij lk l ka a a ≠≠=,则210tg ,()10cos sin cos ii jj ija a d a d t s d d c t ct sθθθθ-⎧=⎪⎪⎪≥⎧⎪=≡=⎨-<⎨⎩⎪⎪=≡⎪⎪=⋅=≡⎩如果jj ii ij a a a -<<,则12ij ii jja t d a a ≈=-,将c,s 代入定理7-9地（1）中可得ii ii ij jj jj ij ij ji c a ta c a ta c c ⎧=+⎪=-⎨⎪==⎩ 每迭代一次地主要工作是选m A 地非对角线元素中地主元素与计算T 111m m m +++=A P AP .首先计算sin ,cos ,θθ,只需计算1m +A 地第i 列,第j 列元素,再算对称元素,不用做3个矩阵地乘法.计算机计算时,需要两组工作单元,以便存储A （或m A ）和R .可用()2()m m lk l ka ε≠=<∑S 控制迭代终止,其中ε是要求地精度.例7-14 用雅可比方法计算对称阵210121012⎡⎤-⎢⎥--⎢⎥⎢⎥-⎣⎦A = 地特征值.解第1步0=A A ,选非对角线元素中地主元素121(1,2)a i j =-==0,1,1/0.7071068,1/0.7071068d t c s ======T 111100.7071068030.70710680.70710680.70710682⎡⎤-⎢⎥==-⎢⎥⎢⎥--⎣⎦A P AP第2步在1A 中选非对角元素地主元素(1)130.7071068(1,3)a i j =-==0.7071068,0.5176381,0.8880738,0.4597008d t c s ====T 22120.63397460.325057600.325057630.627963000.62798302.366025-⎡⎤⎢⎥=--⎢⎥⎢⎥-⎣⎦A P A P 第3步在2A 中选非对角元素地主元素(2)230.627930(2,3)a i j =-==0.5047869,0.6153960,0.8516540,0.5241045d t c s =-=-==-T 33230.63397460.27683660.17036420.27683663.38644600.170364201.979579⎡⎤--⎢⎥=-⎢⎥⎢⎥-⎣⎦A P A P 第4步在3A 中选非对角元素地主元素(3)120.2768366(1,2)a i j =-==4.971292,0.09958013,0.9950785,0.09909004d t c s ====T 44340.606407200.169525803.4140130.016881400.16952580.016881401.979579⎡⎤-⎢⎥=⎢⎥⎢⎥-⎣⎦A P A P 第5步在4A 中选非对角元素地主元素(4)130.1695258(1,3)a i j =-==4.050038,0.1216293,0.9926842,0.1207395d t c s ==== 2T 255450.58578790.20382521000.203825210 3.4140130.0167579000.016757902.000198--⎡⎤⨯⎢⎥=⨯⎢⎥⎢⎥⎣⎦A P A P 于是A 地特征值为1233.414013, 2.000198,0.5857879λλλ===A 地精确特征值为12(1 3.414214λ=≈,22λ=,32(10.585786λ=-≈ 且可逐步求出412345T T T T T T R P P P P P =地列向量,即得A 地近似特征向量.雅可比方法是一个求对称矩阵A 地全部特征值及特征向量地迭代方法,精确度较高,但计算量较大,对稀松带状矩阵经过平面旋转变换后其稀松带状将被破坏,所以很少使用.习题71．设911203111(2102113810A j j B ⎡⎤⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦试估计它们地特征值所在地范围.2．编写幂法程序,并求矩阵732341213A -⎡⎤⎢⎥=⎢⎥⎢⎥--⎣⎦地主特征值及对应地特征向量（准确到小数点后3位）.3．若p 是A 地特征值j λ地一个近似值,且||||()j i p p i j λλ-<-≠则1j pλ-是1()A pE --地主特征值.试用反幂法求矩阵134231111A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地最接近于6地特征值及对应地特征向量.4．设有向量(2,1,2)Tx=,试构造初等反射阵H,使(3,0,0)THx=.5．设(2,3,0,5)Tx=,(1,0,0,5)Te=,用Householder变换化x为与e同方向向量.6．设031042212A⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,求其QR分解.7．设221022212A⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,求其QR分解.8．利用初等反射阵将134312421A⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦正交相似约化为对称三对角阵.9．试用平面旋转变换阵对矩阵A作QR分解,其中111021245A⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦.10．按下列要求编写程序框图.（1）将一般矩阵用豪斯赫尔德方法约化称上Hessenberg阵.（2）对矩阵作QR分解.（3）对上Hessenberg阵应用QR算法求全部特征值及相应地特征向量.11．用QR算法求矩阵120211013A⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦地全部特征值.12．设A是对称矩阵,λ和(1)x x=是A地一个特征值及相应地特征向量.又设p是一个正交阵,使1(1,0,0,,0)Tpx e==证明T=是第一行和第一列除了λ外,其余元素均为零.B PAP。

数值分析第四章矩阵特征值与特征向量的计算

192.9996. 973
12
➢ 幂法的加速—原点移位法
应用幂法计算矩阵A的主特征值的收敛速度主要
由比值 r=|2/1|来决定, 但当r接近于1时, 收敛可能
很慢. 这时可以采用加速收敛的方法.
引进矩阵
B=A－0I
其中0为代选择参数. 设A的特征值为1, 2, …, n, 则B的特征值为1－0, 2－0, …, n－0, 而且A, B
10
2 1 0 例用幂法求矩阵 A 0 2 1
0 1 2
的按模最大的特征值和相应的特征向量.
取 x(0)=(0, 0, 1)T, 要求误差不超过103.
解 y 0 x 0 0 ,0 ,1 T ,
x 1 A 0 0 y , 1 , 2 T , 1 m x ( 1 ) ) a 2 , x
y(1)
x(1)
1
(0,0.5,1)T
x ( 2 ) A ( 1 ) 0 . 5 y , 2 , 2 . 5 T ,2 m x ( 2 ) ) 12 1a . 5 ,
y(2)
x(2) 2
(0.2,0.8,1)T
x ( 3 ) A ( 2 ) 1 . 2 y , 2 . 6 , 2 . 8 T ,3 m x ( 3 ) ) 2 a . 8 ,
x
(
k
1
)
Ax
(k )
A k1 x (0)
在一定条件下, 当k充分大时:
1
x ( k 1) i
x
( i
k
)
相应的特征向量为: x(k1) 4
➢ 幂法的理论依据
n
对任意向量x(0), 有 x(0) tiui ,
i1
x(k1) Ax(k) Ak1x(0)

特征值和特征向量计算的数值方法

特征值和特征向量计算的数值方法在数学和计算机科学领域中，特征值和特征向量是非常重要的概念。

特征值和特征向量的计算有许多不同的数值方法，本文将介绍其中一些常见的数值方法，并分析它们的优劣和适用范围。

一、特征值和特征向量的定义在矩阵理论中，给定一个n×n的矩阵A，如果存在一个非零向量v和一个标量λ，使得Av=λv，那么称v为矩阵A的特征向量，λ为矩阵A的特征值。

特征值和特征向量的计算可以帮助我们理解矩阵的性质以及解决一些实际问题。

二、幂法幂法是计算特征值和特征向量的常用数值方法之一。

幂法的基本思想是通过多次迭代，逐渐逼近矩阵的特征值和特征向量。

具体操作如下：1. 初始化一个非零向量b0；2. 进行迭代计算：bi+1 = A * bi / ||A * bi||；3. 取出近似特征向量的最后一列：v = bn；4. 进行迭代计算特征值：λ = (Av)T * v / (vT * v)。

幂法的主要优点是简单易懂，且只需要进行矩阵向量乘法和内积计算。

然而，幂法仅能求取具有最大特征值的特征向量，而且对于存在多个特征值相等的情况并不适用。

三、反幂法反幂法是幂法的一种改进方法，用于求取矩阵A的最小特征值和对应的特征向量。

反幂法的基本步骤如下：1. 初始化一个非零向量b0；2. 进行迭代计算：bi+1 = (A - μI)^-1 * bi / ||(A - μI)^-1 * bi||；3. 取出近似特征向量的最后一列：v = bn；4. 进行迭代计算特征值：λ = (Av)T * v / (vT * v)。

反幂法的改进之处在于引入了矩阵的逆运算，通过使用矩阵A减去一个合适的常数μ乘以单位矩阵来实现。

反幂法适用于矩阵A的特征值接近于μ的情况。

四、QR方法QR方法也是一种常用的特征值计算方法，它适用于求解所有特征值以及对应的特征向量。

QR方法的基本思想是将一个矩阵分解为正交矩阵Q和上三角矩阵R的乘积，然后迭代地将矩阵A转化为更接近上三角形的形式。

矩阵特征值的数值解法

矩阵特征值的数值解法矩阵的特征值是在矩阵与其特征向量之间的关系中的数值解。

特征值在各个领域中都有广泛应用，包括物理、工程、金融等。

在解决实际问题时，我们经常需要计算矩阵的特征值，因此研究如何求解矩阵特征值的数值方法是非常重要的。

1. 幂迭代法（Power Iteration）幂迭代法是求解矩阵特征值的一种简单而常用的数值方法。

它的基本思想是通过不断迭代矩阵与向量的乘积，使得向量趋近于该矩阵的一个特征向量。

具体步骤如下：（1）初始化一个非零的初始向量x。

（2）进行迭代计算，即$x^{(k+1)}=Ax^{(k)}/，Ax^{(k)}，$。

（3）当向量x的相对误差小于一些预设的精度要求时，停止迭代，此时的x即为矩阵A的一个特征向量。

（4）将x带入特征值的定义式$\frac{Ax}{x}$，计算出特征值。

幂迭代法的优点是简单易实现，计算速度较快，缺点是只能求解特征值模最大的特征向量，而且对于存在特征值模相近的情况，容易收敛到错误的特征值上。

2. QR迭代法（QR Iteration）QR迭代法是一种较为稳定的求解矩阵特征值的数值方法。

它的基本思想是通过不断进行QR分解，使得矩阵的特征值逐渐收敛。

具体步骤如下：（1）将矩阵A进行QR分解，得到正交矩阵Q和上三角矩阵R，令$A_1=RQ$。

（2）将$A_1$再次进行QR分解，得到新的矩阵$A_2=R_1Q_1$。

（3）重复步骤（2），直到得到收敛的矩阵$A_k$，此时$A_k$的对角线上的元素即为矩阵A的特征值。

QR迭代法的优点是对于特征值模相近的情况仍然能够收敛到正确的特征值上。

缺点是每次QR分解都需要消耗大量的计算量，迭代次数较多时计算速度较慢。

3. Jacobi迭代法（Jacobi's Method）Jacobi迭代法是一种通过对称矩阵的对角线元素进行迭代操作，逐步将非对角元素变为零的求解特征值的方法。

具体步骤如下：（1）初始化一个对称矩阵A。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

(1) 任取一个非零向量 v0，要求满足 (x1,v0) 0 (2) 对 k = 1, 2, ... ，直到收敛，计算
vk uk , vk
( uk , Auk ) pk ( uk , uk )
1
vk 1 A pk I uk
21
几点注记
带位移的反幂法中需要计算 vk 1 A pk I uk

v k 1
lim vk
k

1
12
改进的幂法
改进的幂法
(1) 任取一个非零向量 v0，要求满足 (x1,v0) 0 (2) 对 k = 1, 2, ... ，直到收敛，计算
vk uk , vk 1 Auk vk
定理：设 A 有 n 个线性无关的特征向量，其特征值满足 1 2 3 n
Ak u0 Ak 1u0 证： uk , vk 1 k max( A u0 ) max( Ak u0 )
k 1 k

( Auk , uk ) ( A u0 , A u0 ) k k ( uk , uk ) ( A u0 , A u0 )
2 2 k 1 i i i 1 n
10
幂法
幂法中存在的问题 , | 1 | 1 k vk 1 1 x1 0, | 1 | 1
改进方法：规范化
vk 1 Avk
vk uk , vk 1 Auk vk
x1 lim uk k x1
11
幂法
1 的计算
vk Aபைடு நூலகம் v0 uk k vk A v0
定理：设 A 有 n 个线性无关的特征向量，其特征值满足
1 2 n1 n 0
则由反幂法生成的向量满足 xn lim uk , lim vk k xn k

1
n
18
反幂法的加速
n 反幂法的收敛速度取决于 r 的大小 n 1
当 r 接近于 1 时，反乘幂法收敛会很慢！
A-1
的特征值为：
1
1

1
2

1
n1

1
n
对应的特征向量仍然为 x1, x2, ..., xn
反幂法：对矩阵 A-1 使用幂法
17
反幂法
反幂法
(1) 任取一个非零向量 v0，要求满足 (x1,v0) 0 (2) 对 k = 1, 2, ... ，直到收敛，计算
vk uk , vk 1 A1uk vk
Di C
| aii | aij j 1, j i
n
Gerschgorin 圆盘
i=1, 2, ... , n
定理：(Gerschgorin 圆盘定理) 设是 A 的特征值，则
Di
i 1 n
若有 m 个圆盘互相连通，且与其它圆盘都不相连，则这 m 个圆盘内恰好包含 m 个特征值。

1 x1
k 1
2 越小，收敛越快 1
8
幂法的收敛性
当 k 充分大时，有
vk 1k1 x1 vk1 1k11 x1
vk 1 1vk
又 vk 1 Avk
v k 1 j vk j
1 ( j =1, 2, ... , n )
Avk 1vk
可以使用原点平移法对反幂法进行加速
问题：如何选择参数 p ？
离 n 越近越好（但不能相等）
19
幂法的Rayleigh商加速
定理设 A 是 n 阶实对称矩阵，其特征值为
1 2 n
对应的特征向量 x1, x2, ..., xn 满足： ( xi , x j ) ij ，使用改进的乘幂法计算 A 的按模最大特征值 1 时， uk 的 Rayleigh商给出了 1 的较好的近似，即 2k 2 ( Auk , uk ) 1 O 1 ( uk , uk )
k 1 n i k 1 1 1 x1 i xi 1 i2 Ak 1v0 Auk k k n A v0 i k 1 1 x1 i xi i2 1
n

i 1
2 2k i i
2 k 1 O 2 1
20
Rayleigh 商加速
Rayleigh 商加速
xn lim uk k xn
( uk , Auk ) ( xn , Axn ) lim n k ( u , u ) ( xn , xn ) k k
用幂法计算矩阵 B 的主特征值：1 - p
15
举例
例：用带位移的幂法计算下面矩阵的主特征值和对应的特征向
量，取 p=0.75
1.0 1.0 0.5 A 1.0 1.0 0.25 0.5 0.25 2.0
16
反幂法
反幂法
计算矩阵的按模最小的特征值及其特征向量假设：(1) |1| |2| … |n-1| > |n| > 0 (2) 对应的 n 个线性无关特征向量为：x1, x2, ..., xn
2 幂法的收敛速度取决于 r 的大小 1
当 r 接近于 1 时，乘幂法收敛会很慢！
幂法的加速：原点平移法
选择适当的 p 满足：
带位移的幂法
令 B = A – pI，则 B 的特征值为：i - p (1) | 1 p || j p | ( j = 2, ... , n ) j p 2 (2) max 2 j n p 1 1 保持主特征值加快收敛速度
则由改进的幂法生成的向量满足 x lim uk 1 , lim vk k x1 k

1
13
举例
例：用改进的幂法计算下面矩阵的主特征值和对应的特征向量
1.0 1.0 0.5 A 1.0 1.0 0.25 0.5 0.25 2.0
14
幂法的加速
k k (2) Ax x A x x
1 1 B P AP , Ax x By y , y P x (3)
(4) 若 A 对称，则存在正交矩阵 Q，使得
QT AQ diag(1 , 2 ,, n )
4
圆盘定理
设 A=(aij)Rnn ，记
1
A pk I vk1 uk
带位移的反幂法可以用于计算任何一个特征值 k 将参数 p 取为 k 附近若已知特征值，计算特征向量时，可使用带位移的反幂法令 p 足够靠近 k
22
6
R( x )
幂法
幂法（乘幂法，幂迭代）
计算矩阵的主特征值（按模最大）及其特征向量假设：(1) |1| > |2| … |n| 0 (2) 对应的 n 个线性无关特征向量为：x1, x2, ..., xn 计算过程： (1) 任取一个非零向量 v0，要求满足 (x1,v0) 0
数值分析
第八章矩阵特征值计算
—— 幂法与反幂法
1
本章内容
特征值基本性质幂法与反幂法正交变换与矩阵分解 QR 方法
2
本讲内容
特征值基本性质幂法幂法的加速反幂法
3
特征值性质
特征值与特征向量
Ax = x
( C, x 0 )
性质 (1) Ax x ( A I ) x ( ) x
(2) 对 k = 1, 2, ... ，直到收敛，计算
vk Avk 1
7
幂法的收敛性
收敛性分析
设 v0 1 x1 2 x2 n xn (1 0)
v1 Av0 11 x1 22 x2 nn xn
k vk Avk1 11k x1 22k x2 nn xn k k 1k 1 x1 2 2 x2 n n xn 1 1
vk 为 1 的近似特征向量
9
幂法的收敛性
定理：设 A 有 n 个线性无关的特征向量，其特征值满足 1 2 3 n
则由幂法生成的向量满足 ( vk 1 ) j vk lim k 1 x1 , lim 1 k k ( v ) 1 k j
2 注：幂法的收敛速度取决于的大小 1
5
Rayleigh 商
定理：设 A 是 n 阶实对称矩阵，其特征值为 1 2 n
则对任意非零向量 x，有
( Ax , x ) n 1 ( x, x )
且
( Ax, x ) ( Ax , x ) 1 max , n min x 0 ( x, x ) x 0 ( x, x ) ( Ax, x ) 称为矩阵 A 关于 x 的 Rayleigh 商。 ( x, x )