第3章 矩阵特征值与特征向量的计算
一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。
3.1 特征值的估计
较粗估计ρ(A ) ≤ ||A ||
欲将复平面上的特征值一个个用圆盘围起来。
3.1.1 盖氏图
定义3.1-1 设A = [a ij ]n ⨯n ，称由不等式z -a ii ≤记为G i ，i = 1，2，…，n 。
∑a
j =1
j ≠i
n
ij
所确定的复区域为A 的第i 个盖氏图，
j =1j ≠i
n
n
定理3.1-1 若λ为A 的特征值，则λ∈
G
i =1
i
证明：设Ax = λx (x ≠ 0)，若k 使得x k =max x i =x
1≤i ≤n
∞
因为
∑a
j =1
n
kj
x j =λx k
⇒(λ-a kk ) x k =
∑a
j ≠k
n
kj
x j
⇒λ-a kk =
∑
j ≠k n
n
a kj x j x k
≤∑a kj
j =1j ≠k
n
x j x k
≤∑a kj
j =1j ≠k
n
⇒λ∈G k ⊂
G
i =1
i
0. 10. 2⎡1⎢0. 530. 1⎢例1 估计方阵A =⎢10. 3-1⎢
⎣0. 2-0. 3-0. 1
解：
0. 3⎤0. 2⎥⎥特征值的范围 0. 5⎥⎥-4⎦
G 1 = {z ：|z – 1|≤ 0.6}；G 2 = {z ：|z – 3|≤ 0.8}； G 3 = {z ：|z + 1|≤ 1.8}；G 4 = {z ：|z + 4|≤ 0.6}。
注：定理称A 的n 个特征值全落在n 个盖氏圆上，但未说明每个圆盘内都有一个特征值。
3.1.2 盖氏圆的连通部分
称相交盖氏圆之并构成的连通部分为连通部分。
孤立的盖氏圆本身也为一个连通部分。
定理3.1-2 若由A 的k 个盖氏圆组成的连通部分，含且仅含A 的k 个特征值。 证明: 令D = diag(a 11，a 12，…，a nn ) ，M = A – D ，记
⎛a 11⎫⎛0a 12 ⎪
a ⎪ a 21022
A (ε) =D +εM = ⎪+ε
⎪ ⎪a nn ⎭ ⎝⎝a n 1a n 2 a 1n ⎫
⎪
a 2n ⎪
(0≤ε≤1)
⎪
⎪
0⎪⎭
则显然有A (1) = A ，A (0) = D ，易知A (ε) 的特征多项式的系数是ε的多项式，从而A (ε) 的特征
值λ1(ε) ，λ2(ε) ，…，λn (ε) 为ε的连续函数。 A (ε) 的盖氏圆为：G i (ε) ={z :z -a ii ≤
∑|εa
j =1j ≠i
n
ij
|=ε∑|a ij |}⊂G i , (0≤ε≤1)
j =1j ≠i
n
因为A (0) = D 的n 个特征值a 11，a 12，…，a nn ，恰为A 的盖氏圆圆心，当ε由0增大到1时，λi (ε) 画出一条以λi (0) = a ii 为始点，λi (1) = λi 为终点的连续曲线，且始终不会越过G i ；
不失一般性，设A 开头的k 个圆盘是连通的，其并集为S ，它与后n – k 个圆盘严格分离，显然，A (ε) 的前k 个盖氏圆盘与后n – k 个圆盘严格分离。
当ε = 0时，A (0) = D 的前k 个特征值刚好落在前k 个圆盘G 1，…，G k 中，而另n – k 个特征值则在区域S 之外，ε从0变到1时，
。连续曲 G (ε) 与 G (ε) 始终分离（严格）
i
i
i =1
i =k +1
k
n
线始终在S 中，所以S 中有且仅有A 的k 个特征值。
注：1) 每个孤立圆中恰有一个特征值。
2) 例1中G 2，G 4为仅由一个盖氏圆构成的连通部分，故它们各有一个特征值，而G 1，G 3构成的连通部分应含有两个特征值。
3) 因为例1中A 为实方阵，所以若λ为A 的特征值，则λ也是A 的特征值，所以G 2，G 4
中各有一个实特征值。
3.1.3 盖氏圆与相似变换
由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换将盖氏圆的半径变小，以得到更好的估计。
原理：取对角阵作相似变换阵：P = diag(b 1，b 2，…，b n ) 其中b i > 0，i = 1，2，…，n 则B =P AP =diag (
-1
111
, , , ) A diag (b 1, b 2, , b n ) 与A 有相同特征值. b 1b 2b n
而B 的第i 个盖氏圆为：{z :z -a ii ≤
∑a ij
j =1
j ≠i
n
b j b i
，
⎛1 b 1 B =
⎝
1b 2
⎫⎪
⎪⎛a 11a 12⎪ a a 22⎪ 21
⎪ ⎪
1⎪⎝a n 1a n 2b n ⎪⎭
a 1n ⎫⎛b 1⎫
⎪ ⎪
a 2n ⎪ b 2⎪
⎪ ⎪
⎪ ⎪⎪ a nn ⎭⎝b n ⎪⎭
⎛
a 11
b 1 a = 21b
2 a b 1 n 1b
n ⎝
a 12
b 2
b 1
a 22 b a n 22
b n
b n b 1b
a 24n
b 2
a 1n
a nn
⎫
⎪⎪⎪⎪ ⎪⎪⎪⎪⎭
适当选取b 1，b 2，…，b n 就有可能使B 的某些盖氏圆的半径比A 的相应盖氏圆的半径小。 1) 欲缩小G i ，可取b i 最大。
2) 欲缩小除G i 外的圆，可取b i 最小。
⎛0. 90. 010. 12⎫ ⎪
例2，估计A = 0. 010. 80. 13⎪的特征值范围。
0. 010. 020. 4⎪⎝⎭
解：A 的三个盖氏圆分别为：
G 1 = {z ：| z – 0.9|≤ 0.13}；G 2 = {z ：| z – 0.8|≤ 0.14}； G 3 = {z ：| z – 0.4|≤ 0.03} λ3 ∈ G 3，较好。
为了更好地估计另外两个特征值可取b 3最小：
⎛1⎫⎛0. 90. 010. 012⎫ ⎪ ⎪-1
1取b 1 = b 2 = 1，b 3 = 0.1即P = ⎪，则B =P AP = 0. 010. 80. 013⎪ 0. 100. 20. 1⎪0. 4⎪⎝⎭⎝⎭
所以G 1' = {z ：| z – 0.9|≤ 0.022}；G 2' = {z ：| z – 0.8|≤ 0.023}；G 3' = {z ：| z – 0.4|≤ 0.3}
三个盖氏圆分离，故有λ1 ∈ G 1' ，λ2 ∈ G 2' ，λ3 ∈ G 3。
3.2 幂法与反幂法
幂法是求方阵的最大特征值及对应特征向量的一种迭代法。
3.2.1 幂法
设A n 有n 个线性相关的特征向量v 1，v 2，…，v n ，对应的特征值λ1，λ2，…，λn ，满足 |λ1| > |λ2| ≥ …≥ |λn | 1. 基本思想
因为{v 1，v 2，…，v n }为C 的一组基，所以任给x ≠ 0，x 所以有
n
(0)
(0)
(3.2.1)
=∑a i v i —— 线性表示
i =1
n
A x
k (0)
=A (∑a i v i ) =∑a i A k v i
k
i =1
i =1
n n
λk
=∑a i λi k v k =λ1[a 1v 1+∑(i ) k a i v i ]i =1i =2λ1
n
n
(3.2-2)
若a 1 ≠ 0，则因
λi
max(λ1a 1v 1) λ1max(a 1v 1) max(A k x (0) ) ≈==λ1
max(A k -1x (0) ) max(λ1k -1a 1v 1) λ1k -1max(a 1v 1)
k
k
另一方面，记max(x ) = x i ，其中|x i | = ||x ||∞，则当k 充分大时，
若a 1 = 0，则因舍入误差的影响，会有某次迭代向量在v 1方向上的分量不为0，迭代下去可求得λ1及对应特征向量的近似值。 2. 规范化
在实际计算中，若|λ1| > 1则|λ1k a 1| →∞，若|λ1|
⎧(k ) x (k ) ⎪y =
max(x (k ) ) ， k = 0，1，2，… (3.2-4) ⎨(k +1) ⎪=Ay (k ) ⎩x
(0)
定理3.2-1 任给初始向量x
v 1⎧(k )
lim y =⎪max(v 1) ≠0有，⎨k →∞
⎪lim max(x (k ) ) =λ1⎩k →∞
特征向量
(3.2-5)
特征值
证明：
y
(k )
x (k ) Ay (k -1)
==
max(x (k ) ) max(Ay (k -1) )
x (k -1)
A
A k x (0) max(x (k -1) )
== =(k -1) k (0)
x max(A x )
max(A )
max(x (k -1) )
λi k
) v i ](3. 2-2)
λi =21
=n
⎧k ⎫λ
max ⎨λ1[a 1v 1+∑a i (i ) k v i ]⎬
λ1i =2⎩⎭
λ1[a 1v 1+∑a i (
k
n
λi k
) v i ]
k →∞λa 1v 1v 1i =21
=→=
n
max(a v ) max(v ) ⎧⎫λ111
max ⎨[a 1v 1+∑a i (i ) k v i ]⎬
λ1i =2⎩⎭
[a 1v 1+∑a i (
n
而
max(x (k ) ) =max(Ay (k -1) )
k →∞
→max(A
v 1max(Av 1) ) =
max(v 1) max(v 1)
=
max(λv 1) =λ1
max(v 1)
注：若A 的特征值不满足条件(3.2.1)，幂法收敛性的分析较复杂，但若λ1 = λ2 = … = λ r且|λ1| > |λ r +1| ≥ …≥ |λn |则定理结论仍成立。
此时不同初始向量的迭代向量序列一般趋向于λ1的不同特征向量。 3. 算法
求maxa(x ) 的流程，设数组x (n ) 数向量x 的n 个分量
幂法流程：
⎛246⎫ ⎪
例1，用幂法求A = 3915⎪的最大模特征值及对应特征向量
41636⎪⎝⎭
见P312
function y = maxa(x) k=1;n=length(x); for i=2:n
if (abs(x(i))>abs(x(k)),k=i; end; end; y=x(k);
A=[2,4,6;3,9,15;4,16,36]; x0=[1;1;1]; y=x0/maxa(x0) x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) x1=A*y end; y
maxa(x1)
3.2.2 加速方法
⎧(k ) x (k ) ⎪y =
幂法的迭代公式：⎨max(x (k ) )
(k +1) ⎪=Ay (k ) ⎩x
当k →∞时y
(k )
→
v 1
，max(x (k ) ) → λ1，其中|λ1| > |λ2| ≥ …≥ |λn |
max(v 1)
注：幂法的收敛速度取决于比值|λ2| / |λ1|，考虑收敛加速
1. 特征值的Aitken 加速法
(1) 思想：由定理3.2.1的证明知
max(x (k ) ) =max(Ay (k -1) )
n
⎧⎫λ
max ⎨A [a 1v 1+∑a i (i ) k -1v i ]⎬
λ1i =2⎩⎭
=n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
i =2
λ1
λi k
) v i ]λi =21
=λ1
n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
max[a 1v 1+∑a i (
n
n
λi k λ
max[a 1v 1+∑a i () v i ]-max[a 1v 1+∑a i (i ) k -1v i ]
λ1λ1i =2i =2
max(x (k ) ) -λ1=λ1
n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
n
=λ1(
⇒
λ2k -1
) λ1
max[∑a i (
λi k -1λi
) (-1) v i ]λλ1i =22n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
n
(3.2.6)
=λ1(
λ2k -1
) λ1
n
λi λλ
max[a 2(-1) v 2+∑a i (i ) k -1(i -1) v i ]
λ1λ2λ1i =2
n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
≈λ1(
λ2k -1
) M →0(当k 充分大时) λ1
(≈
max(x (k +2) ) -λ1max(x (k +1) ) -λ1⇒≈(k +1)
max(x ) -λ1max(x (k ) ) -λ1
解之得
λ2
) λ1
max(x (k +2) ) -max(x (k ) ) -[max(x (k +1) )]2
λ1≈
max(x (k +2) ) -2max(x (k +1) ) +max(x (k ) ) =max(x (k +2) ) -
[max(x ) -max(x )](k +1)
=λ1
max(x (k +2) ) -2max(x (k +1) ) +max(x (k ) )
∆
(k +2)
(k +1)
2
(3.2.7)
使用λ1(k +2)作为λ1的近似值的算法称为Aitken 加速法。 (2) Aitken加速法
设{x k }线性收敛到x *，即存在c ，|c |
k →∞
2
x k x k +2-x k -x *+1
=0 令=则lim *k →∞x k +2-2x k +1+x k x k -x
算法：
x (0) →max(x (0) ) →y (0) →x (1) →max(x (1) ) →y (1)
→ →x
计算λ1
(k +2)
(k +2)
→max(x
(k +2)
(k +2)
) →y
(k +2)
=max(x
[max(x (k +2) ) -max(x (k +2) )]2
) -
max(x (k +2) ) -2max(x (k +1) ) +max(x (k ) )
流程图
例2 用幂法求方阵A 的最大模特征值，并用Aitkem 加速法
⎛-4140⎫ ⎪A = -5130⎪
-102⎪⎝⎭
解：见（P314）
x0=[1;1;1];y0=x0/maxa(x0) x1=A*y0;y1=x1/maxa(x1) x2=A*y1;y2=x2/maxa(x2)
l0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0)) while (abs(l1-l0))>0.01 x0=x1;x1=x2;l1=l0; x2=A*y2
maxk=maxa(x2) y2=x2/maxk
l0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0)) end;
2. 原点平移法
思想：由矩阵论知，若λ为A 的特征值则λ – a 为A – aI 的特征值，且特征向量相同。 若λ1 – a 为A – aI 的最大模特征值，且
λk -a λ2
λ1-a λ1
则对A – aI 计算λ1 – a 及对应的特征向量比对A 计算收敛得快，此即为原点平移法。
⎧(k ) x (k ) ⎪y =
计算λ1 – a 及特征向量的迭代公式⎨max(x (k ) )
(k +1) ⎪=(A -aI ) y (k ) ⎩x
特征向量：y
(k )
→
v 1
，max(x (k ) ) → λ1– a ，⇒ a + max(x (k ) ) → λ1。
max(v 1)
注：a 的选取较为困难。
⎛-310⎫⎛0⎫ ⎪ ⎪(0)
例3 设A = 1-3-3⎪，x = 0⎪，求最大模特征值及特征向量。
0-34⎪ 1⎪⎝⎭⎝⎭
解：（P315）
幂法：
A=[-3,1,0;1,-3,-3;0,-3,4]; x0=[0;0;1];k=1; y=x0/maxa(x0) x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) k=k+1 x1=A*y end; y
maxa(x1)
原点平移法：
A=[-3,1,0;1,-3,-3;0,-3,4]; x0=[0;0;1];k=1; y=x0/maxa(x0) x1=(A+4*eye(3))*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) k=k+1
x1=(A+4*eye(3))*y end; y
maxa(x1)-4
3. 对称矩阵的Rayleigh 商加速法
x T Ax
定义 设A 对称，x ≠ 0，则称R (x ) =T 为x 关于A 的Rayleigh 商
x x
思想：A 对称⇒特征值λ1，λ2，…，λn 均为实数，且存在特征向量v 1，v 2，…，v n 为标准正
交基。 设x
(0)
=∑a i v i ，a 1 ≠ 0，则x
i =1
n
(k )
=Ax
(k -1)
=... =A x
n
k (0)
=∑λi k a i v i
i =1
n
R (x (k ) ) =
(x ) A (x )
=(k ) T (k )
(x ) (x ) a
(k ) T (k )
(∑λa i v i ) A (∑λi k a i v i )
k i
T
n
(∑λi k a i v i ) T (∑λi k a i v i )
i =1
i =1
n
2
i =1n i =1n
=
∑λ
i =1n i =1
n
2k +12i i
λ
=
2k +12
11
a +∑λi 2k +1a i
i =2
n i =2
22k
λa ∑i i 2k +1211
n
222k
λ1a 1+∑λi 2k a i
λi 2k +1a i 2
) () ]λa 1i =21
=n
λa 22k
λ1a 1[1+∑(i ) 2k (i ) 2]
a 1i =2λ1
λ
a [1+∑(
λi 2k +1a i 2
) () λa 1i =21
=λ1
n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
1+∑(
n
λi 2k +1a i 2n λi 2k a i 2() () -∑() () ∑λa 1a 1i =2λ1
R (x (k ) ) -λ1=λ1i =21
n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
n
λi 2k a i 2λi
) () [-1]λa 1λ1
=λ1i =21n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
∑(
n
=λ1(
λ22k
) λ1
∑(
λi 2k a i 2λi
) () [-1]
a 1λ1i =2λ2
n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
n
≈λ1(
λ22k
) M ' λ1
当k 充分大时，M ' 与k 无关)
注；此比Aitken 加速中的(3.2-6)更快
⎧
y (k ) =⎪
⎪(k +1)
=公式⎨x
⎪(k )
⎪R (y ) =⎩
其中R (y
(k )
x (k )
max(x (k ) ) Ay (k ) 称为Rayleigh 商加速法。 (y (k ) ) T x (k +1) (y (k ) ) T (y (k ) )
v 1
max(v 1)
) →λ1, y (k ) →
注：有了R (x (k ) ) ，R (x (k +1)) ，R (x (k +2)) ，的值，可再用Aitken 加速法得到λ1的一个更好的近似值：
R (x (k +2) ) -λ1R (x (k +1) ) -λ1
≈因为
R (x (k +1) ) -λ1R (x (k ) ) -λ1
所以λ1≈R (x
(k +2)
[R (x (k +2) ) -R (x (k +1) )]2(k +2)
) -=λ 1(k +2) (k +1) (k )
R (x ) -2R (x ) +R (x )
⎛621⎫
⎪
例4 设A = 231⎪，用Rayleigh 商加速法求A 的最大模特征值及特征向量，并与幂法
111⎪⎝⎭
相比较。
解：（P317） 幂法：
A=[6,2,1;2,3,1;1,1,1]; x0=[1;1;1];k=1; y=x0/maxa(x0) x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) x1=A*y k=k+1 end; y
maxa(x1)
Rayleigh 商加速法： A=[6,2,1;2,3,1;1,1,1]; x0=[1;1;1];k=1;r=0; y=x0/maxa(x0) x1=A*y
while(abs(r1-r))>0.001 x0=x1;r1=r;
y=x0/maxa(x0) x1=A*y
r = y'*x1/(y'*y) k=k+1 end; y
maxa(x1) r
3.2.3 反幂法
——用A 代替A 作幂法，即反幂法 1. 求最小模特征值及相应的特征向量 若A 可逆，|λ1| > |λ2| ≥ …≥ |λn |为其特征值，则
-1
1
λn
为A -1的最大模特征值。
迭代公式：x (k +1) = A –1x (k ) ，k = 0，1，2，…，
但A –1不易求，通常可解方程组Ax (k +1) = x (k ) 来求x (k +1)
⎧x (k ) (k ) ⎪y =即有⎨max(x (k ) )
(k +1) ⎪=y (k ) ⎩Ax
k =0, 1, 2,... (3.2.12)
v n ⎧(k )
y →⎪⎪max(v n )
当k → ∞时有⎨
1
⎪max(x (k ) ) →⎪λn ⎩
1
或→λn
(k )
max(x )
注：为解(3.2-12)中的方程组。对A 作LR 分解（带行交模）P A = LR 则有
⎧x (k ) (k )
⎪y =max(x (k ) ) ⎪
⎪(k +1)
=Py (k ) ⎨Lz
⎪
⎪(k +1) (k +1) Rx =z ⎪⎩
k =0, 1, 2,...
2. 求任一特征值及相应特征向量 ——反幂法结合原点平移法
ˆ为λj 的近似值，则(A -λˆ) 的特征值是 思想：若已知λj j
-1
1λ1-λj
,...,
1λj -λj
,...,
1λn -λj
ˆλj -λj
而显然非常大（最大）比值很小
λi -λj λj -λj
1
迭代公式：
⎧(k ) x (k )
⎪y =
max(x (k ) ) ⎨
ˆI ) x (k +1) =y (k ) ⎪(A -λj ⎩
k =0, 1, 2,...
v j ⎧(k )
y →⎪max(v j ) ⎪
当k → ∞时有⎨
1
⎪max(x (k ) ) →
λj -λ⎪j ⎩
1ˆ+或λ→λn j
max(x (k ) )
ˆI ) =LR 注：(1) 若有分LR 解P (A -λj
则迭代公式
⎧x (k ) (k )
⎪y =max(x (k ) ) ⎪
⎪(k +1)
=Py (k ) ⎨Lz
⎪
⎪(k +1) (k +1) Rx =z ⎪⎩
k =0, 1, 2,... (3.2-16)
(2) 在(3.2-16)中直接取z (1) = (1，…，1) T 作初值开始迭代称为半次迭代法
1⎤⎡-12
⎢⎥ =-6.42，用带原点平移的反幂例5 设A =2-41的一个特征值的λ的近似值λ⎢⎥
1-6⎥⎢⎣1⎦
法求λ及相应的特征向量
见[P320] format long;
A=[-1,2,1;2,-4,1;1,1,-6]; x0=[1;1;1];
B=A+6.42*eye(3); C=lu(B); R = triu(C,0);
L =eye(3)+ tril(C,-1); y=x0/maxa(x0); z=[1,1,1]'; x1=inv(R)*z
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) z=inv(L)*y
x1=inv(R)*z end;
-6.42+1/maxa(x1)
预备知识：矩阵论
1. 矩阵QR 分解
定理 设A ∈R n ⨯n 可逆，则存在正交阵Q 与上三角阵R 使A=QR 注：方法 1) 使用史密斯正交变换
2) 使用Householder 变换（反） 3) 使Givens 变换
2. 矩阵Schur 分解
⎛R 11
n ⨯n T
定理 设A ∈R ，则存在正交阵Q 使Q AQ =
⎝R 12 R 1m ⎫
⎪
R 22 R 2m ⎪
——实Schur 型 ⎪ ⎪R mm ⎭
其中R ii 至多2阶。若1阶，其元素即A 的特征值若2阶其特征值为A 的一对共轭复特征值。 注：想加快迭代速度通常先将A 化为上Hessenberg 阵
3.
⎛***
***
∀A ∈R nn ⇒A 正交相似于一个n 阶上Hessenberg 矩阵 0**
00 ⎝
（i >j +1⇒a ij =0） 证明：见（P125）
*⎫⎪ *⎪ *⎪
⎪ ⎪**⎪⎭
§3.3 QR 方法
QR 方法即使用QR 分解构造迭代序列，是目前求一般矩阵全部特征值的最有效并广泛
使用的方法之一。
3.3.1 QR 方法的计算公式
思想：从A 1 = A 出发用正交相似变换得序列{A k }使当k → ∞时，A k 本质收敛到块上三角阵 方法：设A 1 = Q 1R 1（QR 分解），令A 2 = R 1Q 1，设A 2 = Q 2R 2，令A 3 = R 2Q 2，即
⎧A k =Q k R k ⎨
⎩A k +1=R k Q k
(QR 分解) (迭代定义)
k = 1，2，… (3.3-1)
{A k }的性质：
① A k ~ A ：A k +1 = R k Q k = (Q k -1A k ) Q k = R k = Q k -1A k =…= (Q 1…Q k ) -1A (Q 1…Q k )
记G k = Q 1…Q k —— 正交，故有A k ~ A ，且有A 1G k = G k A k +1 ② 记H k = R k …R 1，则A k = G k H k ——QR 分解
G k H k = (Q 1…Q k )(R k …R 1) = G k -1Q k R k Hk -1 = G k -1A k H k -1 = A 1G k -1H k -1 = …= A 1k = A k
注：为求得A 的特征值，只须{A k }能趋于块上三角阵。
定义：① 矩阵列{A k }，当k → ∞时，若其对角元均收敛，且严格下三角部分元素收敛到0，则称{A k }本质收敛到上三角阵。
② 矩阵列{A k }，当k → ∞时，若其对角子块收敛到1阶或2阶的方阵，其下部收敛到0，则称{A k }本质收敛到块上三角阵。
定理3.3-1设A 的特征值满足条件：|λ1| > |λ2| > … > |λn | > 0，v i 为λi 对应的特征向量，i = 1，2，…，n 。记X = (v 1，v 2，…，v n ) ，若有直接三角分解X -1 = LU （杜利特尔分解），则(3.3-1)序列{A k }本质收敛于上三角阵，其主对角元素均为A 的特征值。
⎛5-2-5-1⎫ ⎪10-32⎪ 例1 用QR 方法求的A 特征值，其中A =
022-3⎪ ⎪001-2⎝⎭
见（P322）
注：若A 不满足定理条件，{A k }不一定本质收敛于上三角矩阵。
3.3.2 上Hessenberg 矩阵的QR 方法及带原点平移的QR 方法
在使用QR 方法之前，先A 将作正交相似变换化为上Hessenberg 矩阵H ，然后对H 作QR 迭代，可大量节省运算量。 ① Givens 变换
记s = sinθ，c = cosθ，则J = -s
⎝
⎛c
s ⎫⎪为旋转变换——正交阵。推广到n 维： ⎪c ⎭
⎡1⎤⎢ ⎥⎢⎥
1⎢⎥
⎢⎥
c s ⎢⎥i
⎢⎥1⎢⎥
J (i , k , θ) =⎢ ⎥
⎢⎥1⎢⎥
-s c ⎢⎥k
⎢⎥1⎢⎥
⎥⎢
⎢1⎥⎣⎦
称为Givens 矩阵或Givens 变换（旋转变换）。易知J (i ，k ，θ) 为正交阵。
② 对上Hessenberg 矩阵用Givens 变换作QR 分解
⎡1⎢ ⎢
1⎢
⎢
c i
⎢J (i , i +1, θi ) H =⎢-s i ⎢⎢⎢⎢⎣
s i c i
⎤
⎥⎡h 11h 12 h 1i ⎥⎢h 21h 22 h 2i ⎥⎢
⎥⎢
⎥⎢h ii ⎥⎢
h i +1, i ⎥⎢
1⎥⎢ ⎥⎢
⎥⎢1⎦⎣
h 1, n -1h 2, n -1 h i , n -1 h n , n -1
h i +1, n -1
h 1n ⎤h 2n ⎥⎥ ⎥⎥h i , n ⎥h i +1, ⎥
⎥ ⎥h n , n ⎥⎦
令h i +1,i * = s i h ii + c i h i +1,i = 0，即选择θi 使右边第i +1行第i 列元素为0，而H 的第i 行与第i +1
行零元素位置上左乘J (i ，i +1，θi ) 后仍为0，其他行则不变。（可以证明）
定理
这样i = 1，2，„，n -1共n -1次左乘J 后H 变为上三角阵R 。即U = R ，其中U T = J (n -1，n ，θn -1) „J (1，2，θ1) 正交，且为下Hessenberg 阵，⇒ U 为上Hessenberg 阵 ⇒ H = UR (QR 分解)
③ 记H 1 = H ，设H 1= U 1R 1 令H 2= R 1U 1，则H 2为上Hessenberg 阵。此变换约需4n 2次乘除和加减运算。
⎧H k =U k R k
一般的有⎨得上Hessenberg 阵列{H k }。
H =R U k k ⎩k +1
注：对n 阶上H 阵施行几次迭代后，主对角线下的次对角线上会出现小元素近似为0，将上
H 阵分块，分别使用QR 算法这样可节省计算时间且便于并行处理。 ④ 带原点平移的QR 方法
若在上述QR 迭代中产生得H k 都不可约（即不能将H k 分块）则可仿照幂法中的原点平移法。
取平移量μ对H 1 –μI 作QR 分解：H 1 –μI = Q 1R 1，令H 2 = R 1Q 1 + μI ， 即⎨
⎧H k -μI =Q k R k
k = 1，2，…
⎩H k +1=R k Q k +μI
则∀k 有H k ~ H 1 ~ A 。 若{λi }为A 的特征值：|λ1 – μ | ≥ |λ2 – μ | ≥ … ≥ |λn – μ |，则H k 的第j 个次对角元收敛速度
λj +1-μ取决于。若μ接近于λn 收敛将会很快
λj -μ
k
第3章 矩阵特征值与特征向量的计算
一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。
3.1 特征值的估计
较粗估计ρ(A ) ≤ ||A ||
欲将复平面上的特征值一个个用圆盘围起来。
3.1.1 盖氏图
定义3.1-1 设A = [a ij ]n ⨯n ，称由不等式z -a ii ≤记为G i ，i = 1，2，…，n 。
∑a
j =1
j ≠i
n
ij
所确定的复区域为A 的第i 个盖氏图，
j =1j ≠i
n
n
定理3.1-1 若λ为A 的特征值，则λ∈
G
i =1
i
证明：设Ax = λx (x ≠ 0)，若k 使得x k =max x i =x
1≤i ≤n
∞
因为
∑a
j =1
n
kj
x j =λx k
⇒(λ-a kk ) x k =
∑a
j ≠k
n
kj
x j
⇒λ-a kk =
∑
j ≠k n
n
a kj x j x k
≤∑a kj
j =1j ≠k
n
x j x k
≤∑a kj
j =1j ≠k
n
⇒λ∈G k ⊂
G
i =1
i
0. 10. 2⎡1⎢0. 530. 1⎢例1 估计方阵A =⎢10. 3-1⎢
⎣0. 2-0. 3-0. 1
解：
0. 3⎤0. 2⎥⎥特征值的范围 0. 5⎥⎥-4⎦
G 1 = {z ：|z – 1|≤ 0.6}；G 2 = {z ：|z – 3|≤ 0.8}； G 3 = {z ：|z + 1|≤ 1.8}；G 4 = {z ：|z + 4|≤ 0.6}。
注：定理称A 的n 个特征值全落在n 个盖氏圆上，但未说明每个圆盘内都有一个特征值。
3.1.2 盖氏圆的连通部分
称相交盖氏圆之并构成的连通部分为连通部分。
孤立的盖氏圆本身也为一个连通部分。
定理3.1-2 若由A 的k 个盖氏圆组成的连通部分，含且仅含A 的k 个特征值。 证明: 令D = diag(a 11，a 12，…，a nn ) ，M = A – D ，记
⎛a 11⎫⎛0a 12 ⎪
a ⎪ a 21022
A (ε) =D +εM = ⎪+ε
⎪ ⎪a nn ⎭ ⎝⎝a n 1a n 2 a 1n ⎫
⎪
a 2n ⎪
(0≤ε≤1)
⎪
⎪
0⎪⎭
则显然有A (1) = A ，A (0) = D ，易知A (ε) 的特征多项式的系数是ε的多项式，从而A (ε) 的特征
值λ1(ε) ，λ2(ε) ，…，λn (ε) 为ε的连续函数。 A (ε) 的盖氏圆为：G i (ε) ={z :z -a ii ≤
∑|εa
j =1j ≠i
n
ij
|=ε∑|a ij |}⊂G i , (0≤ε≤1)
j =1j ≠i
n
因为A (0) = D 的n 个特征值a 11，a 12，…，a nn ，恰为A 的盖氏圆圆心，当ε由0增大到1时，λi (ε) 画出一条以λi (0) = a ii 为始点，λi (1) = λi 为终点的连续曲线，且始终不会越过G i ；
不失一般性，设A 开头的k 个圆盘是连通的，其并集为S ，它与后n – k 个圆盘严格分离，显然，A (ε) 的前k 个盖氏圆盘与后n – k 个圆盘严格分离。
当ε = 0时，A (0) = D 的前k 个特征值刚好落在前k 个圆盘G 1，…，G k 中，而另n – k 个特征值则在区域S 之外，ε从0变到1时，
。连续曲 G (ε) 与 G (ε) 始终分离（严格）
i
i
i =1
i =k +1
k
n
线始终在S 中，所以S 中有且仅有A 的k 个特征值。
注：1) 每个孤立圆中恰有一个特征值。
2) 例1中G 2，G 4为仅由一个盖氏圆构成的连通部分，故它们各有一个特征值，而G 1，G 3构成的连通部分应含有两个特征值。
3) 因为例1中A 为实方阵，所以若λ为A 的特征值，则λ也是A 的特征值，所以G 2，G 4
中各有一个实特征值。
3.1.3 盖氏圆与相似变换
由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换将盖氏圆的半径变小，以得到更好的估计。
原理：取对角阵作相似变换阵：P = diag(b 1，b 2，…，b n ) 其中b i > 0，i = 1，2，…，n 则B =P AP =diag (
-1
111
, , , ) A diag (b 1, b 2, , b n ) 与A 有相同特征值. b 1b 2b n
而B 的第i 个盖氏圆为：{z :z -a ii ≤
∑a ij
j =1
j ≠i
n
b j b i
，
⎛1 b 1 B =
⎝
1b 2
⎫⎪
⎪⎛a 11a 12⎪ a a 22⎪ 21
⎪ ⎪
1⎪⎝a n 1a n 2b n ⎪⎭
a 1n ⎫⎛b 1⎫
⎪ ⎪
a 2n ⎪ b 2⎪
⎪ ⎪
⎪ ⎪⎪ a nn ⎭⎝b n ⎪⎭
⎛
a 11
b 1 a = 21b
2 a b 1 n 1b
n ⎝
a 12
b 2
b 1
a 22 b a n 22
b n
b n b 1b
a 24n
b 2
a 1n
a nn
⎫
⎪⎪⎪⎪ ⎪⎪⎪⎪⎭
适当选取b 1，b 2，…，b n 就有可能使B 的某些盖氏圆的半径比A 的相应盖氏圆的半径小。 1) 欲缩小G i ，可取b i 最大。
2) 欲缩小除G i 外的圆，可取b i 最小。
⎛0. 90. 010. 12⎫ ⎪
例2，估计A = 0. 010. 80. 13⎪的特征值范围。
0. 010. 020. 4⎪⎝⎭
解：A 的三个盖氏圆分别为：
G 1 = {z ：| z – 0.9|≤ 0.13}；G 2 = {z ：| z – 0.8|≤ 0.14}； G 3 = {z ：| z – 0.4|≤ 0.03} λ3 ∈ G 3，较好。
为了更好地估计另外两个特征值可取b 3最小：
⎛1⎫⎛0. 90. 010. 012⎫ ⎪ ⎪-1
1取b 1 = b 2 = 1，b 3 = 0.1即P = ⎪，则B =P AP = 0. 010. 80. 013⎪ 0. 100. 20. 1⎪0. 4⎪⎝⎭⎝⎭
所以G 1' = {z ：| z – 0.9|≤ 0.022}；G 2' = {z ：| z – 0.8|≤ 0.023}；G 3' = {z ：| z – 0.4|≤ 0.3}
三个盖氏圆分离，故有λ1 ∈ G 1' ，λ2 ∈ G 2' ，λ3 ∈ G 3。
3.2 幂法与反幂法
幂法是求方阵的最大特征值及对应特征向量的一种迭代法。
3.2.1 幂法
设A n 有n 个线性相关的特征向量v 1，v 2，…，v n ，对应的特征值λ1，λ2，…，λn ，满足 |λ1| > |λ2| ≥ …≥ |λn | 1. 基本思想
因为{v 1，v 2，…，v n }为C 的一组基，所以任给x ≠ 0，x 所以有
n
(0)
(0)
(3.2.1)
=∑a i v i —— 线性表示
i =1
n
A x
k (0)
=A (∑a i v i ) =∑a i A k v i
k
i =1
i =1
n n
λk
=∑a i λi k v k =λ1[a 1v 1+∑(i ) k a i v i ]i =1i =2λ1
n
n
(3.2-2)
若a 1 ≠ 0，则因
λi
max(λ1a 1v 1) λ1max(a 1v 1) max(A k x (0) ) ≈==λ1
max(A k -1x (0) ) max(λ1k -1a 1v 1) λ1k -1max(a 1v 1)
k
k
另一方面，记max(x ) = x i ，其中|x i | = ||x ||∞，则当k 充分大时，
若a 1 = 0，则因舍入误差的影响，会有某次迭代向量在v 1方向上的分量不为0，迭代下去可求得λ1及对应特征向量的近似值。 2. 规范化
在实际计算中，若|λ1| > 1则|λ1k a 1| →∞，若|λ1|
⎧(k ) x (k ) ⎪y =
max(x (k ) ) ， k = 0，1，2，… (3.2-4) ⎨(k +1) ⎪=Ay (k ) ⎩x
(0)
定理3.2-1 任给初始向量x
v 1⎧(k )
lim y =⎪max(v 1) ≠0有，⎨k →∞
⎪lim max(x (k ) ) =λ1⎩k →∞
特征向量
(3.2-5)
特征值
证明：
y
(k )
x (k ) Ay (k -1)
==
max(x (k ) ) max(Ay (k -1) )
x (k -1)
A
A k x (0) max(x (k -1) )
== =(k -1) k (0)
x max(A x )
max(A )
max(x (k -1) )
λi k
) v i ](3. 2-2)
λi =21
=n
⎧k ⎫λ
max ⎨λ1[a 1v 1+∑a i (i ) k v i ]⎬
λ1i =2⎩⎭
λ1[a 1v 1+∑a i (
k
n
λi k
) v i ]
k →∞λa 1v 1v 1i =21
=→=
n
max(a v ) max(v ) ⎧⎫λ111
max ⎨[a 1v 1+∑a i (i ) k v i ]⎬
λ1i =2⎩⎭
[a 1v 1+∑a i (
n
而
max(x (k ) ) =max(Ay (k -1) )
k →∞
→max(A
v 1max(Av 1) ) =
max(v 1) max(v 1)
=
max(λv 1) =λ1
max(v 1)
注：若A 的特征值不满足条件(3.2.1)，幂法收敛性的分析较复杂，但若λ1 = λ2 = … = λ r且|λ1| > |λ r +1| ≥ …≥ |λn |则定理结论仍成立。
此时不同初始向量的迭代向量序列一般趋向于λ1的不同特征向量。 3. 算法
求maxa(x ) 的流程，设数组x (n ) 数向量x 的n 个分量
幂法流程：
⎛246⎫ ⎪
例1，用幂法求A = 3915⎪的最大模特征值及对应特征向量
41636⎪⎝⎭
见P312
function y = maxa(x) k=1;n=length(x); for i=2:n
if (abs(x(i))>abs(x(k)),k=i; end; end; y=x(k);
A=[2,4,6;3,9,15;4,16,36]; x0=[1;1;1]; y=x0/maxa(x0) x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) x1=A*y end; y
maxa(x1)
3.2.2 加速方法
⎧(k ) x (k ) ⎪y =
幂法的迭代公式：⎨max(x (k ) )
(k +1) ⎪=Ay (k ) ⎩x
当k →∞时y
(k )
→
v 1
，max(x (k ) ) → λ1，其中|λ1| > |λ2| ≥ …≥ |λn |
max(v 1)
注：幂法的收敛速度取决于比值|λ2| / |λ1|，考虑收敛加速
1. 特征值的Aitken 加速法
(1) 思想：由定理3.2.1的证明知
max(x (k ) ) =max(Ay (k -1) )
n
⎧⎫λ
max ⎨A [a 1v 1+∑a i (i ) k -1v i ]⎬
λ1i =2⎩⎭
=n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
i =2
λ1
λi k
) v i ]λi =21
=λ1
n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
max[a 1v 1+∑a i (
n
n
λi k λ
max[a 1v 1+∑a i () v i ]-max[a 1v 1+∑a i (i ) k -1v i ]
λ1λ1i =2i =2
max(x (k ) ) -λ1=λ1
n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
n
=λ1(
⇒
λ2k -1
) λ1
max[∑a i (
λi k -1λi
) (-1) v i ]λλ1i =22n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
n
(3.2.6)
=λ1(
λ2k -1
) λ1
n
λi λλ
max[a 2(-1) v 2+∑a i (i ) k -1(i -1) v i ]
λ1λ2λ1i =2
n
λ
max[a 1v 1+∑a i (i ) k -1v i ]
λ1i =2
≈λ1(
λ2k -1
) M →0(当k 充分大时) λ1
(≈
max(x (k +2) ) -λ1max(x (k +1) ) -λ1⇒≈(k +1)
max(x ) -λ1max(x (k ) ) -λ1
解之得
λ2
) λ1
max(x (k +2) ) -max(x (k ) ) -[max(x (k +1) )]2
λ1≈
max(x (k +2) ) -2max(x (k +1) ) +max(x (k ) ) =max(x (k +2) ) -
[max(x ) -max(x )](k +1)
=λ1
max(x (k +2) ) -2max(x (k +1) ) +max(x (k ) )
∆
(k +2)
(k +1)
2
(3.2.7)
使用λ1(k +2)作为λ1的近似值的算法称为Aitken 加速法。 (2) Aitken加速法
设{x k }线性收敛到x *，即存在c ，|c |
k →∞
2
x k x k +2-x k -x *+1
=0 令=则lim *k →∞x k +2-2x k +1+x k x k -x
算法：
x (0) →max(x (0) ) →y (0) →x (1) →max(x (1) ) →y (1)
→ →x
计算λ1
(k +2)
(k +2)
→max(x
(k +2)
(k +2)
) →y
(k +2)
=max(x
[max(x (k +2) ) -max(x (k +2) )]2
) -
max(x (k +2) ) -2max(x (k +1) ) +max(x (k ) )
流程图
例2 用幂法求方阵A 的最大模特征值，并用Aitkem 加速法
⎛-4140⎫ ⎪A = -5130⎪
-102⎪⎝⎭
解：见（P314）
x0=[1;1;1];y0=x0/maxa(x0) x1=A*y0;y1=x1/maxa(x1) x2=A*y1;y2=x2/maxa(x2)
l0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0)) while (abs(l1-l0))>0.01 x0=x1;x1=x2;l1=l0; x2=A*y2
maxk=maxa(x2) y2=x2/maxk
l0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0)) end;
2. 原点平移法
思想：由矩阵论知，若λ为A 的特征值则λ – a 为A – aI 的特征值，且特征向量相同。 若λ1 – a 为A – aI 的最大模特征值，且
λk -a λ2
λ1-a λ1
则对A – aI 计算λ1 – a 及对应的特征向量比对A 计算收敛得快，此即为原点平移法。
⎧(k ) x (k ) ⎪y =
计算λ1 – a 及特征向量的迭代公式⎨max(x (k ) )
(k +1) ⎪=(A -aI ) y (k ) ⎩x
特征向量：y
(k )
→
v 1
，max(x (k ) ) → λ1– a ，⇒ a + max(x (k ) ) → λ1。
max(v 1)
注：a 的选取较为困难。
⎛-310⎫⎛0⎫ ⎪ ⎪(0)
例3 设A = 1-3-3⎪，x = 0⎪，求最大模特征值及特征向量。
0-34⎪ 1⎪⎝⎭⎝⎭
解：（P315）
幂法：
A=[-3,1,0;1,-3,-3;0,-3,4]; x0=[0;0;1];k=1; y=x0/maxa(x0) x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) k=k+1 x1=A*y end; y
maxa(x1)
原点平移法：
A=[-3,1,0;1,-3,-3;0,-3,4]; x0=[0;0;1];k=1; y=x0/maxa(x0) x1=(A+4*eye(3))*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) k=k+1
x1=(A+4*eye(3))*y end; y
maxa(x1)-4
3. 对称矩阵的Rayleigh 商加速法
x T Ax
定义 设A 对称，x ≠ 0，则称R (x ) =T 为x 关于A 的Rayleigh 商
x x
思想：A 对称⇒特征值λ1，λ2，…，λn 均为实数，且存在特征向量v 1，v 2，…，v n 为标准正
交基。 设x
(0)
=∑a i v i ，a 1 ≠ 0，则x
i =1
n
(k )
=Ax
(k -1)
=... =A x
n
k (0)
=∑λi k a i v i
i =1
n
R (x (k ) ) =
(x ) A (x )
=(k ) T (k )
(x ) (x ) a
(k ) T (k )
(∑λa i v i ) A (∑λi k a i v i )
k i
T
n
(∑λi k a i v i ) T (∑λi k a i v i )
i =1
i =1
n
2
i =1n i =1n
=
∑λ
i =1n i =1
n
2k +12i i
λ
=
2k +12
11
a +∑λi 2k +1a i
i =2
n i =2
22k
λa ∑i i 2k +1211
n
222k
λ1a 1+∑λi 2k a i
λi 2k +1a i 2
) () ]λa 1i =21
=n
λa 22k
λ1a 1[1+∑(i ) 2k (i ) 2]
a 1i =2λ1
λ
a [1+∑(
λi 2k +1a i 2
) () λa 1i =21
=λ1
n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
1+∑(
n
λi 2k +1a i 2n λi 2k a i 2() () -∑() () ∑λa 1a 1i =2λ1
R (x (k ) ) -λ1=λ1i =21
n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
n
λi 2k a i 2λi
) () [-1]λa 1λ1
=λ1i =21n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
∑(
n
=λ1(
λ22k
) λ1
∑(
λi 2k a i 2λi
) () [-1]
a 1λ1i =2λ2
n
λa
1+∑(i ) 2k (i ) 2
a 1i =2λ1
n
≈λ1(
λ22k
) M ' λ1
当k 充分大时，M ' 与k 无关)
注；此比Aitken 加速中的(3.2-6)更快
⎧
y (k ) =⎪
⎪(k +1)
=公式⎨x
⎪(k )
⎪R (y ) =⎩
其中R (y
(k )
x (k )
max(x (k ) ) Ay (k ) 称为Rayleigh 商加速法。 (y (k ) ) T x (k +1) (y (k ) ) T (y (k ) )
v 1
max(v 1)
) →λ1, y (k ) →
注：有了R (x (k ) ) ，R (x (k +1)) ，R (x (k +2)) ，的值，可再用Aitken 加速法得到λ1的一个更好的近似值：
R (x (k +2) ) -λ1R (x (k +1) ) -λ1
≈因为
R (x (k +1) ) -λ1R (x (k ) ) -λ1
所以λ1≈R (x
(k +2)
[R (x (k +2) ) -R (x (k +1) )]2(k +2)
) -=λ 1(k +2) (k +1) (k )
R (x ) -2R (x ) +R (x )
⎛621⎫
⎪
例4 设A = 231⎪，用Rayleigh 商加速法求A 的最大模特征值及特征向量，并与幂法
111⎪⎝⎭
相比较。
解：（P317） 幂法：
A=[6,2,1;2,3,1;1,1,1]; x0=[1;1;1];k=1; y=x0/maxa(x0) x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) x1=A*y k=k+1 end; y
maxa(x1)
Rayleigh 商加速法： A=[6,2,1;2,3,1;1,1,1]; x0=[1;1;1];k=1;r=0; y=x0/maxa(x0) x1=A*y
while(abs(r1-r))>0.001 x0=x1;r1=r;
y=x0/maxa(x0) x1=A*y
r = y'*x1/(y'*y) k=k+1 end; y
maxa(x1) r
3.2.3 反幂法
——用A 代替A 作幂法，即反幂法 1. 求最小模特征值及相应的特征向量 若A 可逆，|λ1| > |λ2| ≥ …≥ |λn |为其特征值，则
-1
1
λn
为A -1的最大模特征值。
迭代公式：x (k +1) = A –1x (k ) ，k = 0，1，2，…，
但A –1不易求，通常可解方程组Ax (k +1) = x (k ) 来求x (k +1)
⎧x (k ) (k ) ⎪y =即有⎨max(x (k ) )
(k +1) ⎪=y (k ) ⎩Ax
k =0, 1, 2,... (3.2.12)
v n ⎧(k )
y →⎪⎪max(v n )
当k → ∞时有⎨
1
⎪max(x (k ) ) →⎪λn ⎩
1
或→λn
(k )
max(x )
注：为解(3.2-12)中的方程组。对A 作LR 分解（带行交模）P A = LR 则有
⎧x (k ) (k )
⎪y =max(x (k ) ) ⎪
⎪(k +1)
=Py (k ) ⎨Lz
⎪
⎪(k +1) (k +1) Rx =z ⎪⎩
k =0, 1, 2,...
2. 求任一特征值及相应特征向量 ——反幂法结合原点平移法
ˆ为λj 的近似值，则(A -λˆ) 的特征值是 思想：若已知λj j
-1
1λ1-λj
,...,
1λj -λj
,...,
1λn -λj
ˆλj -λj
而显然非常大（最大）比值很小
λi -λj λj -λj
1
迭代公式：
⎧(k ) x (k )
⎪y =
max(x (k ) ) ⎨
ˆI ) x (k +1) =y (k ) ⎪(A -λj ⎩
k =0, 1, 2,...
v j ⎧(k )
y →⎪max(v j ) ⎪
当k → ∞时有⎨
1
⎪max(x (k ) ) →
λj -λ⎪j ⎩
1ˆ+或λ→λn j
max(x (k ) )
ˆI ) =LR 注：(1) 若有分LR 解P (A -λj
则迭代公式
⎧x (k ) (k )
⎪y =max(x (k ) ) ⎪
⎪(k +1)
=Py (k ) ⎨Lz
⎪
⎪(k +1) (k +1) Rx =z ⎪⎩
k =0, 1, 2,... (3.2-16)
(2) 在(3.2-16)中直接取z (1) = (1，…，1) T 作初值开始迭代称为半次迭代法
1⎤⎡-12
⎢⎥ =-6.42，用带原点平移的反幂例5 设A =2-41的一个特征值的λ的近似值λ⎢⎥
1-6⎥⎢⎣1⎦
法求λ及相应的特征向量
见[P320] format long;
A=[-1,2,1;2,-4,1;1,1,-6]; x0=[1;1;1];
B=A+6.42*eye(3); C=lu(B); R = triu(C,0);
L =eye(3)+ tril(C,-1); y=x0/maxa(x0); z=[1,1,1]'; x1=inv(R)*z
while(abs(maxa(x1)-maxa(x0)))>0.001 x0=x1;
y=x0/maxa(x0) z=inv(L)*y
x1=inv(R)*z end;
-6.42+1/maxa(x1)
预备知识：矩阵论
1. 矩阵QR 分解
定理 设A ∈R n ⨯n 可逆，则存在正交阵Q 与上三角阵R 使A=QR 注：方法 1) 使用史密斯正交变换
2) 使用Householder 变换（反） 3) 使Givens 变换
2. 矩阵Schur 分解
⎛R 11
n ⨯n T
定理 设A ∈R ，则存在正交阵Q 使Q AQ =
⎝R 12 R 1m ⎫
⎪
R 22 R 2m ⎪
——实Schur 型 ⎪ ⎪R mm ⎭
其中R ii 至多2阶。若1阶，其元素即A 的特征值若2阶其特征值为A 的一对共轭复特征值。 注：想加快迭代速度通常先将A 化为上Hessenberg 阵
3.
⎛***
***
∀A ∈R nn ⇒A 正交相似于一个n 阶上Hessenberg 矩阵 0**
00 ⎝
（i >j +1⇒a ij =0） 证明：见（P125）
*⎫⎪ *⎪ *⎪
⎪ ⎪**⎪⎭
§3.3 QR 方法
QR 方法即使用QR 分解构造迭代序列，是目前求一般矩阵全部特征值的最有效并广泛
使用的方法之一。
3.3.1 QR 方法的计算公式
思想：从A 1 = A 出发用正交相似变换得序列{A k }使当k → ∞时，A k 本质收敛到块上三角阵 方法：设A 1 = Q 1R 1（QR 分解），令A 2 = R 1Q 1，设A 2 = Q 2R 2，令A 3 = R 2Q 2，即
⎧A k =Q k R k ⎨
⎩A k +1=R k Q k
(QR 分解) (迭代定义)
k = 1，2，… (3.3-1)
{A k }的性质：
① A k ~ A ：A k +1 = R k Q k = (Q k -1A k ) Q k = R k = Q k -1A k =…= (Q 1…Q k ) -1A (Q 1…Q k )
记G k = Q 1…Q k —— 正交，故有A k ~ A ，且有A 1G k = G k A k +1 ② 记H k = R k …R 1，则A k = G k H k ——QR 分解
G k H k = (Q 1…Q k )(R k …R 1) = G k -1Q k R k Hk -1 = G k -1A k H k -1 = A 1G k -1H k -1 = …= A 1k = A k
注：为求得A 的特征值，只须{A k }能趋于块上三角阵。
定义：① 矩阵列{A k }，当k → ∞时，若其对角元均收敛，且严格下三角部分元素收敛到0，则称{A k }本质收敛到上三角阵。
② 矩阵列{A k }，当k → ∞时，若其对角子块收敛到1阶或2阶的方阵，其下部收敛到0，则称{A k }本质收敛到块上三角阵。
定理3.3-1设A 的特征值满足条件：|λ1| > |λ2| > … > |λn | > 0，v i 为λi 对应的特征向量，i = 1，2，…，n 。记X = (v 1，v 2，…，v n ) ，若有直接三角分解X -1 = LU （杜利特尔分解），则(3.3-1)序列{A k }本质收敛于上三角阵，其主对角元素均为A 的特征值。
⎛5-2-5-1⎫ ⎪10-32⎪ 例1 用QR 方法求的A 特征值，其中A =
022-3⎪ ⎪001-2⎝⎭
见（P322）
注：若A 不满足定理条件，{A k }不一定本质收敛于上三角矩阵。
3.3.2 上Hessenberg 矩阵的QR 方法及带原点平移的QR 方法
在使用QR 方法之前，先A 将作正交相似变换化为上Hessenberg 矩阵H ，然后对H 作QR 迭代，可大量节省运算量。 ① Givens 变换
记s = sinθ，c = cosθ，则J = -s
⎝
⎛c
s ⎫⎪为旋转变换——正交阵。推广到n 维： ⎪c ⎭
⎡1⎤⎢ ⎥⎢⎥
1⎢⎥
⎢⎥
c s ⎢⎥i
⎢⎥1⎢⎥
J (i , k , θ) =⎢ ⎥
⎢⎥1⎢⎥
-s c ⎢⎥k
⎢⎥1⎢⎥
⎥⎢
⎢1⎥⎣⎦
称为Givens 矩阵或Givens 变换（旋转变换）。易知J (i ，k ，θ) 为正交阵。
② 对上Hessenberg 矩阵用Givens 变换作QR 分解
⎡1⎢ ⎢
1⎢
⎢
c i
⎢J (i , i +1, θi ) H =⎢-s i ⎢⎢⎢⎢⎣
s i c i
⎤
⎥⎡h 11h 12 h 1i ⎥⎢h 21h 22 h 2i ⎥⎢
⎥⎢
⎥⎢h ii ⎥⎢
h i +1, i ⎥⎢
1⎥⎢ ⎥⎢
⎥⎢1⎦⎣
h 1, n -1h 2, n -1 h i , n -1 h n , n -1
h i +1, n -1
h 1n ⎤h 2n ⎥⎥ ⎥⎥h i , n ⎥h i +1, ⎥
⎥ ⎥h n , n ⎥⎦
令h i +1,i * = s i h ii + c i h i +1,i = 0，即选择θi 使右边第i +1行第i 列元素为0，而H 的第i 行与第i +1
行零元素位置上左乘J (i ，i +1，θi ) 后仍为0，其他行则不变。（可以证明）
定理
这样i = 1，2，„，n -1共n -1次左乘J 后H 变为上三角阵R 。即U = R ，其中U T = J (n -1，n ，θn -1) „J (1，2，θ1) 正交，且为下Hessenberg 阵，⇒ U 为上Hessenberg 阵 ⇒ H = UR (QR 分解)
③ 记H 1 = H ，设H 1= U 1R 1 令H 2= R 1U 1，则H 2为上Hessenberg 阵。此变换约需4n 2次乘除和加减运算。
⎧H k =U k R k
一般的有⎨得上Hessenberg 阵列{H k }。
H =R U k k ⎩k +1
注：对n 阶上H 阵施行几次迭代后，主对角线下的次对角线上会出现小元素近似为0，将上
H 阵分块，分别使用QR 算法这样可节省计算时间且便于并行处理。 ④ 带原点平移的QR 方法
若在上述QR 迭代中产生得H k 都不可约（即不能将H k 分块）则可仿照幂法中的原点平移法。
取平移量μ对H 1 –μI 作QR 分解：H 1 –μI = Q 1R 1，令H 2 = R 1Q 1 + μI ， 即⎨
⎧H k -μI =Q k R k
k = 1，2，…
⎩H k +1=R k Q k +μI
则∀k 有H k ~ H 1 ~ A 。 若{λi }为A 的特征值：|λ1 – μ | ≥ |λ2 – μ | ≥ … ≥ |λn – μ |，则H k 的第j 个次对角元收敛速度
λj +1-μ取决于。若μ接近于λn 收敛将会很快
λj -μ
k