2004-12-1 1
第八章矩阵特征值与特征向量的计算
§8.0 问题描述
§8.1 乘幂法与反幂法
§8.2 雅可比方法
§8.3 QR方法
§8.4 求实对称三对角矩阵特征值的二分法
2004-12-1 2
§8.0 问题描述设A为n×n矩阵，所谓A的特征问题是求数λ和非零向量x，使
Ax＝λx
成立。数λ叫做A的一个特征值，非零向量x叫做与特征值λ对应的特征向量。这个问题等价于求使方程组
(A- λI)x=0有非零解的数λ和相应的非零向量x。
线性代数理论中是通过求解特征多项式det(A- λI)=0的零点而得到λ，然后通过求解退化的方程组(A- λI)x=0
而得到非零向量x。当矩阵阶数很高时，这种方法极为困难。目前用数值方法计算矩阵的特征值以及特征向量比较有效的方法是迭代法和变换法。
Back
2004-12-1 3
§8.1 乘幂法与反幂法一、乘幂法通过求矩阵特征向量求出特征值的一种迭法方法，
它用以求按模最大的特征值和相应的特征向量。
n
λλλ ≥≥≥ L
21
∑
=
=+++==
n
j
jjjnnn
xaxaxaxaAVV
1
222111
0)1(
λλλλ L
）（
设实矩阵A的特征值为λ
1
，λ
2
，…,λ
n
，相应的特征向量线性无关。设A的特征值按模排序为：
n
xxx,,,
21
L
∑
=
=+++=
n
j
jjnn
xaxaxaxaV
1
2211
)0(
L
则对任一非零向量，可以得到：
n
RV ∈
)0(
令，可以构造一个向量序列，
L,2,1,0,
)()1(
==
+
kAVV
kk
2004-12-1 4
乘幂法（续）
∑
=
=+++==
n
j
jjjnnn
xaxaxaxaAVV
1
22
22
2
211
2
1
1)2(
λλλλ L
）（
M
∑
=
=+++==
n
j
jj
k
jnn
k
n
kkkk
xaxaxaxaAVV
1
222111
1)(
λλλλ L
）（
下面分四种情况讨论：
))((
2
1
111
)0()(
∑
=
+==
n
j
jj
k
j
kkk
xaxaVAV
λ
λ
λ
),,3,2(
1
nj
j
L=> λλ
若由于，故k充分大时，
)2(1
1
≥< i
i
λ
λ
,0
1
≠a
111111
)(
)( xaxaV
k
k
kk
λελ ≈+=
1.
2004-12-1 5
乘幂法（续）
)(k
V是相应于的近似特征向量
1
λ
设表示
)(k
l
V
个分量。的第lV
k )(
nl
xa
xa
V
V
lkl
k
lkl
k
k
l
k
l
,,2,1,
1
111
111
1
1
)(
)1(
L=≈
+
+
=
+
++
λ
ελ
ελ
）（）（
）（）（
2.主特征值是实数，但不是单根，如
r
λλλ === L
21 nr
λλλ ≥≥>
+
L
11
而，则
Remark:具体计算时，V
(0)
的选取很难保证一定有
α
1
≠0。但是，由于舍入误差的影响，只要迭代次数足够多，如，就会有
，因而最后结论是成立的。对于的情形，
由于对任意l均有上面的结论，故只要取另外的l使即可。
nn
m
xaxaxaV
′
++
′
+
′
= L
2211
)(
0
1
≠
′
a
0
)(
=
k
l
V
0
)(
≠
k
l
V
2004-12-1 6
乘幂法（续）
））（（
ii
k
n
ri
i
i
r
i
i
kk
xaxaV
∑∑
+==
+=
1
1
1
1
)(
λ
λ
λ
)(k
V
仍然是相对于的近似特征向量。
1
λ
nl
V
V
k
l
k
l
,,2,1,
)(
)1(
1
L=≈
+
λ
Remark:由于相应于λ
1
的特征向量子空间可能不是一维的，由上式得到的V
(k)
只是该子空间的一个特征向量。而且不同的V
(0)
可能得到线性无关的V
(k)
。
3.且为实数，
21
λλ?= ),,4,3(
21
nj
j
L=>= λλλ
n
k
n
n
kkkk
xaxaxaxa )()()1((
1
3
1
3
322111
)(
λ
λ
λ
λ
λ +++?+= LV
由于
( )
*
22111
)1(
k
kk
xaxa ελ +?+=
2004-12-1 7
乘幂法（续）
当k充分大以后，有或
2
1
2
λ≈
+
）（
）（
k
l
k
l
V
V
）（）（k
l
k
l
VV /
2
1
+
≈λ
又可推出：
）（）（kk
VV
1
1
λ+
+
)11
*
2211
1
1
*
122
1
11
1
1 k
kk
k
kk
xaxaxaxa ελελ +?+++?+=
+
+
++
）（（））（（
11
1
1
2 xa
k+
≈ λ
22
1
1
1
1
1
21 xaVV
kkkk +++
≈? λλ）（
）（）（
）（）（kk
VV
1
1
λ?
+
和
）（）（kk
VV
1
1
λ+
+
分别看作
21
λλ和的近似特征向量。
2004-12-1 8
乘幂法（续）
4.
n
λλλλλλ ≥≥>== L
32121
，且此处只限于讨论实矩阵的情形。若有
12121
xxee
ii
====
，而且，必有
θθ
ρλλρλ
∑
=
++=+++=
n
j
jjnn
xaxaxaxaxaxaV
3
11112211
)0(
L
当V
(0)
为实向量时，有
111111
3
111111
)(
xaxaxaxaxaV
kk
n
j
jj
k
j
kkk
λλλλλ +≈++=
∑
=
11
1
111
1
1
)1(
xaxaV
kkk +++
+≈ λλ
11
2
111
2
1
)2(
xaxaV
kkk +++
+≈ λλ
2004-12-1 9
乘幂法（续）
容易验证，
0
)(
11
)1(
11
)2(
≈++?
++ kkk
VVV λλλλ）（
若存在二实数p，q，使得
0
)()1()2(
=++
++ kkk
qVpVV
则可以得到
=++
=++
0)(
0)(
1
2
11
1
2
11
qp
qp
k
k
λλλ
λλλ
上式说明是方程λ
2
＋pλ＋q＝0的一对共轭复根。
这里的系数p、q可由（*）式的分量形式确定，即
11
,λλ
(*)
(**)
jlnjl
qVpVV
qVpVV
k
j
k
j
k
j
k
l
k
l
k
l
≠≤≤
=++
=++
++
++
,,1
0
0
)()1()2(
)()1()2(
2004-12-1 10
乘幂法（续）
解出p和q后，即可由（**）式求得
2
2,1
22
）（
p
qi
p
±?=λ
又因为
≈?
≈?
+
+
22122
)(
1
)1(
11211
)(
2
)1(
)(
)(
xaVV
xaVV
kkk
kkk
λλλλ
λλλλ
故与λ
1
，λ
2
对应的特征向量分别为V
(k+1)
- λ
2
V
(k)
和
V
(k+1)
- λ
1
V
(k)
。
2004-12-1 11
乘幂法（续）
于是得到以下规范化算法：
L,2,1
max/
0
)()()(
)1()(
)0()0(
=
=
=
≠=
k
VVU
AUV
VU
kkk
kk
）（
征向量的方法称为幂法。
k
A
)0(
V上述由已知非零向量及矩阵A的乘幂构造
{ }
）（k
V
1
λ
向量序列以计算A 的特征值及相应的特
)(
)(k
V
其中，max表示向量的最大分量。
)(k
V
以上讨论只是说明了乘幂法的基本原理。当太小或太大时，将会使过小或过大，以致运算无法继续进行。因此，实际计算时，每步要进行规范化。
1
λ
∞
)(k
V
2004-12-1 12
乘幂法（续）
对于规范化的幂法有如下收敛性定理：
证明：取）构造向量序列（00
1
1
)0(
≠≠=
∑
=
axaV
j
n
j
j
其特征值满足
n
λλλ ≥≥> L
21,则按规范化幂法定理：
nn
R
×
∈
设A 有n个线性无关的特征向量，
()
)max(
lim
1
1
x
x
U
k
k
=
∞→
构造的序列有
( )
1
)max(lim λ=
∞→
k
k
V
0)0()1(
AVAUV ==
）（
）（
）（
）（
）（）（
0
0
)1(11
max
max/
AV
AV
VVU ==
2004-12-1 13
==
)1()2(
AUV
）（
）（
）（
0
02
max AV
VA
乘幂法（续）
）（
）（）（）（222
max/ VVU =
)
)max(
)(
max(
1
.
max
)0(
)0(20
02
AV
VAAV
VA
）（
）（
）（
=
)max(
)0(2
)0(2
VA
VA
=
2004-12-1 14
乘幂法（续）
==
）（1)( kk
AUV
）（
）（
）（
01
0
max VA
VA
k
k
M
)max(
)0(
)0(
)(
VA
VA
U
k
k
k
=
][
1
2
111
1
0
j
k
j
n
j
j
k
i
n
i
k
ii
k
xaxaxaVA）（
）（
λ
λ
λλ
∑∑
==
+==
]max[
][max
][
1
1
2
1
111
2
1
111
)(
x
x
xaxa
xaxa
U
n
j
j
j
j
k
n
j
j
j
j
k
k
→
+
+
=
∑
∑
=
=
）（
）（
λ
λ
λ
λ
λ
λ
即
=
∞→
）（K
k
Ulim
]max[
1
1
x
x
而故
2004-12-1 15
=
）（
）（
）（
01
0
)(
max
max)max(
VA
VA
V
k
k
k
+
+
=
∑
∑
=
=
n
j
j
k
j
j
k
n
j
j
k
j
j
k
xaxa
xaxa
2
1
1
11
1
1
2
1
111
]max[
][
max
）（（
）（
λ
λ
λ
λ
λ
λ
1
2
λ
λ
1
maxlim λ=
∞→
）（
）（k
k
V
即，收敛速度依赖于。
]
)(
]max[
[
max
1
2
1
11
2
1
111
∞→→
+
+
=
∑
∑
=
=
k
xaxa
xaxa
n
j
j
kj
j
n
j
j
kj
j
λ
λ
λ
λ
λ
λ
）（（
）（
乘幂法（续）
证毕
2004-12-1 16
乘幂法（续）
Remark：上述定理仅适用于的情形。
n
λλλ ≥≥> L
21
对于规范化幂法，若当分别有两
{ } { }
）（）（
和
122 +kk
UU
个不同的极限，说明。此时若有m，使
21
λλ?=
)1()2()(1
,
+++
==
mmmm
AVVAUV
ε<?
∞
）（）（2mm
UU
，再做两次迭代
2
1
)(
)2(
λ≈
+
m
l
m
l
U
V
12
)(
)2(
1
,
)(
)(
λλλ?==
+
m
l
m
l
U
V
)(
11
)1(
1
)2(
λλ属于xVV
mm
≈+
++
特征向量分别为则有即
( )
)(
22
)1
1
)2(
λλ属于xVV
mm
≈?
++
2004-12-1 17
乘幂法（续）
若当规范化序列无规律可循时，这时说明，
。在m次迭代后，由
{ }
)(k
U
21
λλ =,
)()1( mm
AUV =
+
( ) )2()3()1(2
,
++++
==
mmmm
AVVAVV
做三次迭代，
然后按
=++
=++
+++
+++
0
0
1)2()3(
1)2()3(
）（
）（
m
j
m
j
m
j
m
l
m
l
m
l
qVpVV
qVpVV
计算q，p
再用所得之p，q计算
，
∞
+++
++
)1()2()3( mmm
qVpVV
2004-12-1 18
乘幂法（续）
再由
++
++
2
)2(
1
3
1
)2(
2
)3(
xVV
xVV
mm
mm
作为作为
）（
λ
λ
，计算特征向量。
,
21，
求λ
，）（，就由小于
，
2
21
22
p
qi
p
±?=λε若其值作为初始向量，
( )
）（
）（3
3
)3(
max
+
+
+
=
m
m
m
V
V
U否则以重复以上计算步骤,直至停止。
）（）（
ε<++
∞
+++ )1(23 mmm
qVpVV
2004-12-1 19
二、乘幂法的加速技术
1
2
λ
λ
=r
乘幂法收敛速度依赖于比值。表示满足的那个下标。当r<1但接近于1
时，收敛可能很慢。因此必须研究加速收敛的方法。
0
j
λ
njj
λλλλ ≥≥>==
LL
00
11
1、原点平移法用B=A-PI来代替A解乘幂法，P是一个可选参数。
如果A的特征值为，则B 的相应特征
n
λλλ,
21
L，，
PPP
n
λλλ,,,
21
L
征值为，而且A，B的特征向量相同，XPXPIABX )()(?=?= λ
2004-12-1 20
原点平移法（续）
若A的主特征值为λ
1
，则P的选取要求仍然是B 的主特征值，即
P?
1
λ
),,4,3,2(
1
njPP
j
L=?>? λλ
且使
1
2
1
2
max
λ
λ
λ
λ
<
≤≤
p
p
j
nj
这时对A-PI解乘幂法，就会有较快的收敛速度。
的按模最大的特征值之后，A的按模求得A-PI
1
μ
最大的特征值为,
而A-PI对应于的特征向量就是1
μ
1
λ
A对应于的特征向量。
参数P的选取依赖于对A的特征值分布的大体了解。这可通过盖尔圆定理得到矩阵特征值的分布。
1
λ
2004-12-1 21
原点平移法（续）
定理（Gerschgorin圆盘定理）：设A为n×n实矩阵，则
（1）A的每一特征值必属于下述几个圆盘（称为盖尔圆）
niara
n
ij
j
ijiii
,,2,1,
1
L==≤?
∑
≠
=
λ
的并集之中；
（2）由矩阵A的所有盖尔圆组成的连通部分中任取一个，
若它是由k个盖尔圆构成，则在这个连通部分中有且仅有k个特征值（盖尔圆相重时重复计算，特征值相同也重复计算）。
Remark：原点平移法由于p的选取困难，故很少单独使用。通常将原点平移法与反幂法或QR方法结合使用。
2004-12-1 22
2、Rayleigh商加速法（用于对称矩阵）
定义：设A 为n阶实对称矩阵，x为任意非零向量，则
称为A关于x的Rayleigh商。
），（
）（
xx
xAx
xR
,
)( =
定理：设A 为实对称矩阵，特征值满足
nn
R
×
∈
n
λλλ ≥≥> L
21
,应用规范化幂法公式计算A 的
）（K
U
1
λ
主特征值的规范化向量 的Rayleigh商可给出
1
λ
的较好的近似（对任意的）。0
)0(
≠V
))((
2
1
2
1
k
kk
kk
k
O
UU
UAU
UR
λ
λ
λ +==
），（
），（
）（
）（）（
）（）（
）（
2004-12-1 23
Rayleigh商加速法（续）
∑
=
=
n
j
jj
k
j
k
xaVA
1
)0(
λ
)max(
)0(
)0(
)(
VA
VA
U
k
k
k
=
证明：
由
20001
]/[max)(
）（）（）（）（）（
），（，VAVAVAUAU
kkkkk +
=
），（
），（
）（
）（）（
）（）（
）（
kk
kk
k
UU
UAU
UR =
2000
]/[max)(
）（）（）（）（）（
），（，VAVAVAUU
kkkkk
=
∑
∑
=
=
+
+
==
n
j
k
jj
n
j
k
jj
kk
kk
a
a
VAVA
VAVA
1
22
1
122
00
001
λ
λ
），（
），（
）（）（
）（）（
))((
1
1
2
1
2
1
2
2
1
2
1
2
12
1
2
1
1
k
n
j
k
jj
n
j
k
jj
O
a
a
a
a
λ
λ
λ
λ
λ
λ
λ
λ +=
+
+
∑
∑
=
=
+
）（）（
）（）（
注：此处取的A的特征向量是标准正交化的。
证毕
=
2004-12-1 24
三、反幂法对用乘幂法求解按模最大的特征值是，特
1?
A
n
λ
1
设矩阵A非奇异，用代替A作幂的方法就称为反
1?
A
n
λλλ
111
21
≤≤≤ L
的特征值满足，并且A对应于A
-1
的相应的特征向量是相同的。
n
λλλ ≥≥≥ L
21
1?
A
幂法。当A的特征值满足时，
n
x
征向量是，即是A的按模最小的特征值和特征向量。
计算矩阵按模最小的特征值及相应的特征向量。
反幂法的算法为：取，构造向量序列。
0
)0()0(
≠=UV
=
=
)max(
)(
)(
)(
)1(1)(
k
k
k
kk
V
V
U
UAV
），，（L21=k
2004-12-1 25
反幂法（续）
为避免求逆，可将上面的序列改变为
=
=
)max(
)(
)(
)(
)1()(
k
k
k
kk
V
V
U
UAV
），，（L21=k
当时，由反幂法构造的
0
21
>>≥≥
n
λλλ L
{ }{ }
)()(
,
kk
UV
向量序列满足
n
k
k
n
n
k
k
V
x
x
U
λ
1
)max(lim).2(
)max(
lim).1(
)(
)(
=
=
∞←
∞→
收敛速度的依赖于比值为。
1?n
n
λ
λ
2004-12-1 26
反幂法（续）
Remark：若已知矩阵A的某个特征值λ
i
的相对分离较好的近似值p。不要求p的近似程度有多好，只要求j≠i
时，，则便是的主特征值。
这样一来，结合原点平移法的反幂法的计算公式为：
pp
ji
<? λλ
p
i
λ
1
1
)(
pIA
L,2,1
)max(
)(
0
)(
)(
)(
)1(1)(
)0()0(
=
=
=
≠=
k
V
V
U
UpIAV
VU
k
k
k
kk
2004-12-1 27
反幂法（续）
则向量序列满足
{ }{ }
)()(
,
kk
UV
()
=
=
∞→
∞←
p
V
x
x
U
i
k
k
i
i
k
k
λ
1
maxlim).2(
)max(
lim.1
)(
)(
）（
）（
）（k
k
i
V
p
max
1
lim
∞→
+=λ
其收敛速度依赖于。
p
p
i
j
λ
λ
Back
2004-12-1 28
§8.2 雅可比方法用以求实对称矩阵全部特征值和特征向量的一个变换方法。其基本思想是通过一系列正交变换，化实对称矩阵为对角阵，这个对角阵的对角元素就是该矩阵的特征值，这些正交矩阵的乘积矩阵的各个列就是相应的特征向量。
有关代数知识：
A1.矩阵和相似矩阵的特征值相同。
1?
= PAPB
Q
2.若矩阵满足，则称为正交阵，当然
Q
IQQ
T
=
T
QQ =
1
k
QQQ,,,
21
L
k
QQQQ L
21
=。正交阵的乘积仍为正交阵。
2004-12-1 29
有关代数知识（续）
Q
3.若为实对称矩阵，则存在正交阵使A
∧ A∧=
T
QAQ
。其中为对角阵，的对角元是
∧
的特
i
q
T
Q ),,2,1( ni L= A
征值，而的列就是的特征向量。
4.设为n阶对称矩阵，，其中为正交A
T
QAQC =
Q
22
FF
AC =
。
矩阵，则
5.设为实对称矩阵，为任意正交矩阵，则
A
Q
∑
=
====
n
j
j
F
T
FFF
AQAQAQQAA
1
2
2
222
)(λ
2004-12-1 30
6.Givens旋转矩阵
p
=
1
1
cossin
1
1
sincos
1
1
),,(
O
LLL
MM
MOM
MM
LLL
O
θθ
θθ
θqpR
q
p
称为旋转阵。是一个正交阵，即
),,( θqpR
。
IRR
T
=
T
RR =
1
q
有关代数知识（续）
2004-12-1 31
考虑实矩阵
=
2221
1211
aa
aa
A
)0(
2112
≠=aa
一、古典Jacobi方法
=
θθ
θθ
cossin
sincos
R
令易知，对任何θ，R为正交矩阵。称R为平面旋转矩阵。
22
)1(
2221
1211
)1(
)(
cossin
sincos
cossin
sincos
×
=
==
ij
T
a
aa
aa
RARA
θθ
θθ
θθ
θθ
令由矩阵乘法可得
+?==
+?=
++=
)sin(coscossin)(
coscossin2sin
sincossin2cos
22
121122
)1(
21
)1(
12
2
2212
2
11
)1(
22
2
2212
2
11
)1(
11
θθθθ
θθθθ
θθθθ
aaaaa
aaaa
aaaa
2004-12-1 32
古典Jacobi方法（续）
，当时，可选。
2211
21
2
2
aa
a
tg
=θ
2211
aa =
4
π
θ =
由上式最后一式知，当θ满足时，有
0
)1(
21
)1(
12
== aa
这也就是说，当A为实二阶对称阵时，用适当的正交相似变换一次即可把A化为对角阵。当A为n阶实对称阵时，
要用到n阶旋转矩阵（Givens旋转矩阵），它与单位矩阵的区别仅在于(p,p),(p,q),(q,p),(q,q)四个位置的元素不一样。容易验证，R(p,q,θ)具有如下性质：
（1）R(p,q,θ)为正交阵；
（2）RA只改变A的第p行与第q行元素，AR
T
只改变
A的第p列与第q列元素，RAR
T
只改变A的第p、q行与
p、q列的元素。
2004-12-1 33
古典Jacobi方法（续）
Jacobi方法就是用一系列的旋转相似变换逐渐将A
化为对角阵的过程：
L,2,1,0,
111
0
=
=
=
+++
k
RARA
AA
T
kkkk
恰当地选取每个旋转阵R
k+1
，就可以使A
k
趋于对角化。
设R
k+1
＝R(p,q,θ) 。由R
k+1
正交可知，A
k+1
与A
k
相似，
且A
k
都为是对称阵。A
k+1
与A
k
的差别仅在于p、q行与
p、q列的元素。由矩阵乘法可得：
qpj
aaaa
aaaa
k
jq
k
qj
k
pj
k
qj
k
jp
k
qj
k
pj
k
pj
,,
cossin
sincos
)1()()()1(
)1()()()1(
≠
=+?=
=+=
++
++
θθ
θθ
2004-12-1 34
古典Jacobi方法（续）
=?+?=
+?=
++=
++
+
+
)1(22)()()()1(
2)()(2)()1(
2)()(2)()1(
)sin(coscossin)(
coscossin2sin
sincossin2cos
k
qp
k
pq
k
pp
k
qq
k
pq
k
qq
k
pq
k
pp
k
qq
k
qq
k
pq
k
pp
k
pp
aaaaa
aaaa
aaaa
θθθθ
θθθθ
θθθθ
qpjiaa
k
ij
k
ij
,,,
)()1(
≠=
+
由上面几式容易知道：
qpjiaa
qpjaaaa
k
ij
k
ij
k
qj
k
pj
k
qj
k
pj
,,,)()(
,,)()()()(
2)(2)1(
2)(2)(2)1(2)1(
≠=
≠+=+
+
++
再由正交矩阵的性质（代数知识5）有：
2)(2)(2)(2)1(2)1(2)1(
)(2)()()(2)()(
k
pq
k
qq
k
pp
k
pq
k
qq
k
pp
aaaaaa ++=++
+++
2004-12-1 35
古典Jacobi方法（续）
若，选取θ，使，只需θ满足：
0
)(
≠
k
pq
a 0
)1(
=
+k
pq
a
)()(
)(
2
)2tan(
k
qq
k
pp
k
pq
aa
a
=θ
4
π
θ ≤
这里假定。若，则取
，故有
)()( k
qq
k
pp
aa =
4
)(
)(
π
θ?=
k
pp
asign
=
=
θθ
θ
cos)(sin
2
2
cos
)(k
pq
asign
当时，令
)()( k
qq
k
pp
aa ≠
)()(
)(
2
1
)2tan(
k
qq
k
pp
k
pq
aa
a
d?
==θ
2004-12-1 36
古典Jacobi方法（续）
则
1tan,01tan2tan
22
+±?==?+ ddd θθθ
为避免相近的数相减，取
1
)(
tan
2
++
==
dd
dsign
t θ
θθθ cossin,
1
1
cos
2
t
t
=
+
=
最后得此时有
2)(2)(2)(2)1(2)1(
)(2)()()()(
k
pq
k
qq
k
pp
k
qq
k
pp
aaaaa ++=+
++
引入记号与
∑
≠
=
=
n
ji
ji
ij
aAS
1,
2
)(
∑
=
=
n
i
ii
aAD
1
2
)(
其中D(A)表示对角元的平方和，S(A)表示非对角元的平方和。
2004-12-1 37
古典Jacobi方法（续）
则
( )
()
=
+=
+
+
2
)(
1
2
)(
1
2)()(
2)()(
k
pqkk
k
pqkk
aASAS
aADAD
这说明，只要，则按上述方法构造的旋转矩阵R(p,q,θ)对A
k
变换后，就会使对角线元素平方和增加，非对角线元素平方和减少。
0
)(
≠
k
pq
a
Remark：使用旋转矩阵变换后，虽然，但可能又会变成非零。因此，并不能通过有限次旋转变换就可以将矩阵A化为对角阵。但是，
Jacobi方法确实是收敛的，我们可以通过反复变换将
A化为近似的对角阵。
0
)1(
=
+k
pq
a
)1(
)(
+> kma
m
pq
下面我们来证明古典Jacobi方法的收敛性。
2004-12-1 38
古典Jacobi方法（续）
证明：
我们在计算时是选取，故
)()(
max
k
ij
ji
k
pq
aa
<
=
定理：设A为实对称矩阵，则由A
0
=A出发，经古典Jacobi
方法计算所产生的矩阵序列{ A
k
}收敛于对角阵，且就是A的全部特征值。
),,,(
21 n
diag λλλ L=Λ
n
λλλ,,,
21
L
( )
2
)(
1
2)()(
k
pqkk
aASAS?=
+
由前面的推导过程知，
从而
2)(2
1,
2)(
1,
2)(
))(()()()(
k
pq
n
ji
ji
k
pq
n
ji
ji
k
ijk
annaaAS?=≤=
∑∑
≠
=
≠
=
)(
)1(
1
)(
2)(
k
k
pq
AS
nn
a
≥
)()
)1(
2
1()(
)1(
2
)()(2)()(
2)(
1 kkk
k
pqkk
AS
nn
AS
nn
ASaASAS
=
≤?=
+
)()
)1(
2
1()(
0
1
1
AS
nn
AS
k
k
+
+
≤
当n>2时，显然有，即非对角元的平方和趋于零，A
k
趋于对角阵，Jacobi方法收敛。
0)(
lim
=
∞→
k
k
AS
证毕
2004-12-1 39
特征向量的计算
T
kkkk
TT
RARARARAARRA
12122111
,
=== L，，
由
T
k
T
kkkk
T
kkkk
RRARRRARA
1211
==
T
k
T
k
TT
kk
RRRARRRRR
121121
= LL
IH
T
=
0
T
k
TTT
k
RRRH L
21
=
记
L,2,1,
0
1
=
=
=
k
HAHA
RHH
T
kkk
T
k
T
k
T
k
则有
2004-12-1 40
特征向量的计算（续）
若矩阵序列A
k
收敛于对角阵，
则的第j列就是A对应于的标准正交特征向量的近似值。
T
k
H
j
λ
),,,(
21 n
diag λλλ L=Λ
T
k
H的计算可以使用上述的递推关系式的第一式，
其元素间的关系为：
≠=
=+?=
+=
qpjhh
nihhh
hhh
k
ij
k
ij
k
iq
k
ip
k
iq
k
iq
k
ip
k
ip
,,
,,2,1,cossin
sincos
)1()(
)1()1()(
)1()1()(
Lθθ
θθ
2004-12-1 41
Jacobi方法的计算步骤
1.选主元，即确定p，q(p<q)，使；
)()(
max
k
ij
ji
k
pq
aa
<
=
2.计算cosθ，sinθ；
3.按照递推公式计算新的正交矩阵的元素；
T
k
H
1+
4.计算A
k+1
的元素，其中。0
)1(
=
+k
pq
a
反复执行以上各步，直至为止。这里的ε为给定的误差界。
ε<
+
)(
1k
AS
2004-12-1 42
二、Jacobi过关法使用古典Jacobi方法时，每次先要在非对角元中挑选主元，这要花费很多机时。实用中常采用Jacobi
过关法，方法如下：
1.给定控制误差限ε；
2.计算非对角元素的平方和；
∑∑
+=
=
=
n
ij
ij
n
i
a
1
2
1
1
0
2ν
3.设置一个阀值，如；
0
1
>ν
n
0
1
ν
ν =
4.对A的非对角元a
ij
(i<j)逐个扫描，若某个，则立即对A做一次旋转相似变换，之后对所得新矩阵继续扫描，只要有非对角元的绝对值大于，就用旋转相似变换将其变为零。如此多次扫描和变换，直到每个非对角元；
1
ν
1
ν>
ij
a
1
ν≤
ij
a
2004-12-1 43
Jacobi过关法（续）
5.若，则结束计算，得到特征值与特征向量；
否则转向6；
εν ≤
1
6.缩小阀值，比如用代替，重复4～5。
n
1
ν
1
ν
Remark：Jacobi算法的优点是数值稳定性比较好，精度高，求得的特征向量正交性好。
Jacobi算法的缺点是当Ａ为稀疏矩阵时，Jacobi
变换将破坏其稀疏性质。
Back
2004-12-1 44
§8.3 QR方法一、基本思想
QR方法是求一般矩阵的全部特征值与特征向量的最有效的方法之一。它是在1961年由J.G.F.Francis
首次提出的。该方法的基本原理如下：
设一般的实矩阵A＝(a
ij
)
n×n
可以分解为A＝QR的形式，其中Q为正交阵，R为上三角阵。
1111111112111
,QAQQRQQQRARQAA
TT
=====
令由此可以产生矩阵序列
,
11 kkkkk
RQQRA ==
L,2,1,
111
===
+++
kRQQRA
kkkkk
2004-12-1 45
基本思想（续）
kk
T
kkkk
T
kkkk
QAQQRQQQRA ===
+1由于，故与彼此相似。事实上，中所有矩阵彼此相似，从而它们有相同的特征值。
A
k 1+
A
k
{ }
A
k
k
Q
k
R
序列称为QR序列，其中为正交阵，为上三角阵。
{ }
A
k
{ }
A
k
若趋于一个对角阵或上（下）三角阵，则当k充分大时，的对角元素就可以做为Ａ的特征值。
{ }
A
k
现在的问题是：（1）Ａ的QR分解是否存在；若存在，分解是否唯一？（2）在何种条件下收敛于上述的对角阵或上（下）三角阵？
{ }
A
k
2004-12-1 46
QR方法的实施方法
（1）用反射变换先将A相似化为上Hessenberg阵B，
BAQQ
T
=
此时A与B的特征值相同，即。
（2）用旋转变换对上Hessenberg阵B进行QR分解
RQ
111
RQBB ==
并且有仍为上Hessenberg阵。
（3）形成QR序列
kkk
RQB =一般地，对上Hessenberg阵B作QR分解，
再形成QR序列，
kkk
QRB =
+1
直到B
k＋1
趋于一个上三角阵。
2004-12-1 47
二、反射矩阵（Householder矩阵）
设ν∈R
n
，，称矩阵H＝I－2νν
T
为反射阵。
1
2
=ν
由HIIIH
TTTTTT
=?=?=?= νννννν 2)(2)2(
TTTTTT
TTTTTT
HHIII
IIIIHH
==+?=+?=
==
νννννννννν
νννννννν
44)(44
)2)(2()2()2(
知H是对称阵，也是正交阵。
S
-y
y
x
0
w
w′
x
ν
考虑以ν为法向量，过原点x
0
的超平面S={x|(ν,x)=0,x∈R
n
}（如右图所示）。设w为R
n
中任一向量，则有w=x+y，其中x∈S，
y∈S
⊥
（S的正交补空间）。
2004-12-1 48
反射矩阵（Householder矩阵）（续）
显然，ν与x正交，从而有ν
T
x＝0。
SxxxxIHx
TT
∈=?=?= )(2)2( νννν
由于y∈S
⊥
，故y=cν（c为常数），从而有
cccy
TT
===
2
2
νννν
ycccyyyIHy
TT
=?=?=?=?= ννννννν 2)(2)2(
S
-y
y
x
0
w
w′
x
ν
wyxHyHxyxHHw
′
=?=+=+= )(
从而
Remark：可以看出，w′是w关于S的镜面反射，S
的法向量就是H中的单位向量ν。这就是反射阵的几何意义。
2004-12-1 49
反射矩阵（Householder矩阵）（续）
定理（反射阵的性质）：?x，y∈R
n
，，
总?反射阵H，s.t,Hx＝y。
22
yx =
若x=y，则取ν⊥y，且即可。
1
2
=ν
证明：
今设x≠y，为确定ν，s.t,Hx＝(I-2νν
T
)x=y
证毕只需-2ν(ν
T
x)=y-x即可，即取ν为平行于
y-x的单位向量即可。故可取
2
xy
xy
±=ν
Remark：利用该定理的结论知道，利用反射阵，
可将任一非零向量x变成与另一向量y平行的向量，
且该向量的长度与x的长度相同。
2004-12-1 50
反射矩阵（Householder矩阵）（续）
证明：
显然
2
2
1
1
2
2
1
2
1
)( xxee
n
i
i
==?=?
∑
=
σσ
推论：设x=(x
1
,x
2
,…,x
n
)
T
≠0,则?反射阵H，
s.t,Hx＝-σe
1
，其中
=+==
=+=
==
=
∑
=
xuxu
eexu
xxsignxxsign
uuIH
T
T
n
i
i
T
)(
2
1
)0,,0,1(,
)()(
1
2
2
11
2
1
1
2
1
2
1
1
σσβ
σ
σ
β
L
利用反射阵的性质定理，若记y＝-σe
1
，则必存在反射阵H＝I－2νν
T
，
2
yx
yx
=ν
2004-12-1 51
反射矩阵（Householder矩阵）（续）
使Hx=y＝-σe
1
，而
22
1
1
2
u
u
ex
ex
yx
yx
=
+
+
=
=
σ
σ
ν
T
n
xxxexu ),,,(
211
Lσσ +=+=
其中故
T
T
T
uuI
u
uu
IIH
1
2
2
22
=?=?= βνν
xuxxxxu
T
n
=+=++++== )(])[(
2
1
2
1
1
22
2
2
1
2
2
σσσβ L
证毕
Remark：由该推论知，反射阵H可将一个向量中的多个元素一次性变为零。
2004-12-1 52
反射矩阵（Householder矩阵）（续）
反射变换的具体做法：
1.设x= 计算
,0),,,(
21
≠
T
n
xxx L
2
1
2
1
2
)( xx
n
i
i
=
∑
=
2.计算
∑
=
=
n
i
i
xxsign
1
2
1
2
1
))((σ
3.计算
T
n
xxxexu ),,,(
211
Lσσ +=+=
xux
T
u
=+== )(
2
1
1
2
2
σσβ
4.计算
5.计算
T
uuIH
1?
= β
Remark：若x
1
为零，一般可取其符号函数值为1。
2004-12-1 53
三、平面旋转矩阵（Givens矩阵）
反射阵H对于大量引进零元素是方便的。然而在许多计算中，必须有选择消去一些元素，解决这一问题的有效工具是Jacobi方法中介绍过的Givens变换矩阵。
下面我们给出下文要用到的一个重要定理：
定理：
设已知向量,其中不全为零，则可选择平面旋转矩阵R(i,j,θ)，s.t.
T
n
xxxx ),,,(
21
L=
ji
xx,
T
nji
xxxxRx ),,,,,,(
1
LLL
′′
=
其中
′
=
′
=
=
′
≠+=
′
i
j
i
i
jjii
x
x
x
x
xxxx
θθ sin,cos
0,0)(
2
1
22
2004-12-1 54
平面旋转矩阵（Givens矩阵）（续）
=Rx
证明：令
1
1
cossin
1
1
sincos
1
1
O
LLL
MM
MOM
MM
LLL
O
θθ
θθ
=
x
x
x
x
x
x
x
x
n
j
i
n
j
i
'
'
'
'
1
1
M
M
M
M
M
M
M
M
M
M
M
M
M
M
2004-12-1 55
平面旋转矩阵（Givens矩阵）（续）
有
≠=
′
+?=
′
+=
′
jikxx
xxx
xxx
kk
jij
jii
,,
cossin
sincos
θθ
θθ
令
0cossin =+?=
′
jij
xxx θθ
得0,
cos
sin
tan ≠==
i
i
j
x
x
x
θ
θ
θ
又利用
1cossin
22
=+ θθ
可得
1cos)(cos
222
=+ θθ
i
j
x
x
22
2
2
cos
ji
i
xx
x
+
=θ
22
2
2
sin
ji
j
xx
x
+
=θ
2004-12-1 56
平面旋转矩阵（Givens矩阵）（续）
取
22
cos
ji
i
xx
x
+
=θ
22
sin
ji
j
xx
x
+
=θ
当时，可取
0=
i
x
。1sin,0cos == θθ
此时有
()
=
′
+=
+
+
=+=
′
0
sincos
2
1
22
22
22
j
ji
ji
ji
jii
x
xx
xx
xx
xxx θθ
证毕
2004-12-1 57
平面旋转矩阵（Givens矩阵）（续）
平面旋转变换的具体做法：
①计算
22
jii
xxx +=
′
0=
′
j
x0≠
j
x
对指定的，需变换后有。
②计算
22
cos
ji
i
xx
x
+
=θ
22
sin
ji
j
xx
x
+
=θ
③给出R(i,j,θ)与
T
nji
xxxxRx ),,,,,,(
1
LLL
′′
=
′
=
′
=
=
′
≠+=
′
i
j
i
i
jjii
x
x
x
x
xxxx
θθ sin,cos
0,0)(
2
1
22
其中
2004-12-1 58
四、矩阵的QR分解定理：设A= 为实非奇异矩阵，则A有正交
nnij
a
×
)(
分解A=QR，其中Q为正交阵，R为上三角阵。
且当R具有正对角元时，分解A=QR是唯一的。
下面用反射矩阵证明（也可用旋转矩阵证明）。
证明：
对非奇异矩阵A按列分块，记为
),,,()(
21
)1(
1 nnnij
aaaaAA L===
×
Step1：由于a
1
≠0，故?反射阵H
1
，s.t,H
1
a
1
=-σ
1
e
1
，
2
22
221
)2()2(
2
)2(
2
)2(
22
)2(
1
)2(
121
)1(
1
)1(
21
)1(
111
0
0
0
),,,( A
Dc
Br
aa
aa
aa
aHaHaHAH
nnn
n
n
T
n
=
=
==
σ
σ
L
MMMM
L
L
L
2004-12-1 59
矩阵的QR分解（续）
其中，H
1
的定义见反射阵推论部分。
()0,,
1)2(
2
)2(
222
≠∈=
n
T
n
RaaC L
Step2：由于C
2
≠0，故?（n-1）阶反射阵,s.t.
2
H
′
1222
eCH
′
=
′
σ
1
e
′
其中为n-1维欧式空间的第一个标准正交基向量。
定义，则H
2
仍然为反射阵，且
′
=
2
2
1
H
H
3
33
33
2
)2(
121
)()3(
3
)3(
3
)3(
33
)3(
2
)3(
232
)2(
1
)2(
13
)2(
121
12
A
DC
Br
a
aa
aa
aa
aaa
AHH
n
nnn
n
n
n
=
=
=
σ
σ
σ
σ
L
MOM
L
L
L
2004-12-1 60
矩阵的QR分解（续）
其中，H
2
的定义见反射阵推论部分。
0
2
3
≠∈
n
RC
Stepk：设已经完成对Ａ的第一步至第k步约化过程，即
反射阵H
1
，H
2
，…，H
k –1
，s.t.
==
+
+
+
+?
+?
aaa
aaa
aaa
aaaa
aaaaa
k
nn
k
kn
k
nk
k
kn
k
kk
k
kk
k
nk
k
kk
k
kkk
nkkk
nkkk
kkk
AAHHH
)()(
)1(,
)(
)()(
)1(,
)(
)(
),1(
)(
)1(),1(
)(
),1(1
)3(
2
)3(
)1(,2
)3(
2
)3(
1,22
)2(
1
)2(
)1(,1
)2(
1
)2(
1,1
)2(
121
121
L
MOMM
L
L
MMMMMO
LL
LL
L
σ
σ
σ
2004-12-1 61
矩阵的QR分解（续）
=
kk
kkk
DC
BrR
0
其中，H
i
(i=1,2,…,k-1)的定义见反射阵推论部分。
0
1
≠∈
+?kn
k
RC
由于C
k
≠0，故?（n-k+1）阶反射阵,s.t.
k
H
′
1
~
eCH
kkk
′
=
′
σ
1
~
e
′
其中为（n-k+1）维欧式空间的第一个标准正交基向量，计算反射矩阵的计算公式为：
k
H
′
+=
′
=
+=
′
=
′′
=
′
+
=
+?
∑
)(
2
1
),,,(
)()(
)(
2
2
)()(
,1
)(
2)()(
1
1
k
kkkkkk
Tk
nk
k
kkk
k
kkk
n
ki
k
ik
k
kkk
T
kkkknk
au
aaau
aasign
uuIH
σσβ
σ
σ
β
L
2004-12-1 62
矩阵的QR分解（续）
定义，则H
k
仍然为反射阵，且
′
=
k
k
k
H
I
H
1
′′
=
′
=
==
+?
kkk
kkk
kk
kkk
k
k
kkkkk
DHe
BrR
DC
BrR
H
I
AAHHHAH
1
1
111
~
00 σ
L
这样就使A的三角分解过程前进了一步。
当进行步后，有，且为上三角阵，可记。
)1(?n
nnn
AAHHH =
121
L
n
A
)(n
AR =
由于反射阵H
1
，H
2
，…，H
n –1
都是正交阵，正交阵的乘积、逆阵仍然是正交阵，故若记，则
T
n
TT
HHHQ
121?
= L
QRRHHHA
T
n
TT
==
121
L
2004-12-1 63
矩阵的QR分解（续）
下面证明分解的唯一性。
设A=Q
1
R
1
=Q
2
R
2
，Q
1
、Q
2
均为正交阵，R
1
、R
2
为非奇异上三角阵，且对角元均为正，则
222222
111111
RRRQQRAA
RRRQQRAA
TTTT
TTTT
==
==
由假设（R具有正的对角元）及Cholesky分解的唯一性可得，R
1
＝R
2
。
证毕
Remark：实际计算H
k
时，可按更简单的方法进行，即在公式中直接将写成，
从而得到。
k
u
′
Tk
nk
k
kkk
k
kkk
aaau ),,,,0,,0(
)()(
,1
)(
LL
+
+= σ
T
kkkk
uuIH
1?
= β
2004-12-1 64
矩阵的QR分解（续）
其中每个R
ii
是1×1或2×2 的矩阵。若是1×1的，
其元素就是A的特征值。若是2×2的，则R
ii
的特征值就是A的一对共轭复特征根。
nn
RA
×
∈
对于一般的实矩阵，有如下分解定理：
定理：（实Schur分解定理）设，则?正交阵，s.t.
nn
RA
×
∈
nn
RQ
×
∈
==
mm
m
m
T
R
RR
RRR
RAQQ
MO
L
L
222
11211
Remark：该定理表明，可以通过正交相似变换将A化为块上三角阵。此时可直接求得A的特征值，但很难直接求得定理中的Q和R。
2004-12-1 65
五、Householder方法
=
nnnn
n
n
bb
bbb
bbb
B
1,
22221
11211
MOO
L
L
定义：若时，的方阵B称为上Hessenberg阵，
即
1+> ji
0=
b
ij
如果，则称B为不可约上Hessenberg
阵。类似的也有下Hessenberg阵的定义。
)1,,2,1(0
,1
=≠
+
nib
ii
L
基本思想：Householder方法是先用反射阵将A正交相似变换为上Hessenberg阵，然后对所得之上
Hessenberg阵使用QR方法，这时的QR分解通常用平面旋转矩阵来实现。
2004-12-1 66
Householder方法（续）
定理：设A为n阶实矩阵，则存在反射阵H
1
，
H
2
，…，H
n-2
，s.t.
BHHAHHHH
nnn
=
221132
LL
其中B为上Hessenberg阵。
推论：当A为n阶实对称矩阵时，存在反射阵
H
1
，H
2
，…，H
n-2
，s.t.
=
)(
1
1
)3(
332
2
)2(
221
1
)1(
11
221132
n
nnn
n
nnn
a
a
a
a
HHAHHHH
σ
σ
σ
σσ
σ
OO
O
LL
2004-12-1 67
Householder方法（续）
QR方法的计算过程为：
1.用Householder变换化矩阵为上Hessenberg阵B；
2.用平面旋转变换对B进行QR分解：B＝QR
（可以证明，RQ仍为上Hessenberg阵）；
3.反复进行QR分解，直到矩阵序列B
k
趋于一个对角阵或上三角阵；
4.计算出3中所得之矩阵的对角元，即是原矩阵的特征值。
Remark：详细的计算实施办法还可参考关冶，陈景良编写的《数值计算方法》一书。
2004-12-1 68
六、QR方法的收敛性在A受相当限制的假设下，由QR方法产生的矩阵序列{A
k
}基本上收敛于上三角阵。
定理（QR方法的收敛性）：设且
nn
RA
×
∈
2.存在非奇异矩阵X，使Ａ＝XDX
-1
，D＝diag(λ
1
，
λ
2
，…，λ
n
)，且设X
-1
有三角分解X
-1
＝LU（L为单位下三角阵，U为上三角阵），则QR算法产生的矩阵序列{A
k
}基本上收敛于上三角阵，即
1.Ａ的特征值满足；
0
21
>>>>
n
λλλ L
==
∞→
n
k
k
RA
λ
λ
λ
MO
L
L
*
**
2
1
lim
2004-12-1 69
QR方法的收敛性（续）
这里基本收敛的意思是
)(,0),,,2,1(,
)()(
limlim
jiania
k
ij
k
i
k
ii
k
>===
∞→∞→
Lλ
但不一定存在。
)(,
)(
lim
jia
k
ij
k
<
∞→
推论：若A为n阶实对称矩阵且满足QR方法收敛性定理的条件，则由QR方法产生的矩阵序列{A
k
}
收敛于对角阵D＝diag(λ
1
，λ
2
，…，λ
n
) 。
Remark：对于一般的实矩阵，QR方法产生的矩阵序列{A
k
}的收敛情况很复杂。在某些条件下，{A
k
}收敛于实Schur型的矩阵。
2004-12-1 70
七、带原点平移的QR方法假定A已化为上Hessenberg阵（仍用A表示）。
为了使QR方法收敛的更快，通常使用带原点平移的QR方法，其一般形式为：
=+=
=?
=
+
L,2,1,
1
1
kQRISA
RQISA
AA
kkkk
kkkk
k
S其中称为平移量。这样产生的矩阵序列{A
k
}具有下列性质：
1.由于故有A
k+1
与A
k
相似；
kkkkkkkkk
QAQQIsAQISA
11
1
)(
+
=?+=
2004-12-1 71
带原点平移的QR方法（续）
2,当A
k
为上Hessenberg阵时，若QR分解采用平面旋转变换，则A
k+1
也是上Hessenberg阵。
迭代若干步后，若足够小，就可以认为是特征值的近似值，然后将A
k+1
的第n行、
第n列去掉（剩下的n-1阶矩阵仍为上Hessenberg
阵），对压缩后的矩阵重复上述分解过程。
)1(
1,
+
k
nn
a
)1(
,
+k
nn
a
Remark1：恰当选取平移量可以加速过程的收敛。
的选取通常有两种方法：
k
S
k
S
（1）直接选取；
)(k
nnk
as =
2004-12-1 72
带原点平移的QR方法（续）
（2）当初始矩阵A为实对称矩阵时，经反射变换后得到的上Hessenberg阵实际上是一个对称三对角阵。当再用一系列旋转变换进行分解时，对称三对角性保持不变。可以证明，这种情况下平移量s
k
取二阶矩阵的最接近的那一个特征值，能使很快趋于零。
aa
aa
k
nn
k
nn
k
nn
k
nn
)()(
1,
)(
,1
)(
1,1
a
k
nn
)(
a
k
nn
)(
1,?
Remark2：当A的特征值为复数时，计算过程限定在实数范围内时，应该使用双步QR方法。该方法可参阅聂铁军编写的《数值计算方法》一书。
Back
2004-12-1 73
§8.4 求实对称三对角矩阵特征值的二分法一、问题由上节知道，实对称矩阵通过正交相似变换，
可以化为如下的对称三对角矩阵。
设
=
nn
n
A
αβ
β
αβ
βαβ
βα
1
1
32
221
11
OO
O
且。否则，可将A表示成适当的分块对角阵，每一块的对角元均不为零。
1,,2,1,0?=≠ ni
i
Lβ
对于这类矩阵，可以用二分法求其全部的特征值。
2004-12-1 74
二、矩阵A的特征多项式序列及其性质定义：称实系数多项式序列?
0
(x),?
1
(x),…,?
m
(x)为?
m
(x)
在(a，b)区间内的一个施铎姆（Sturm）序列，如果它具有以下三个性质：
（1）?
0
(x)在(a，b)区间内无实根；
（2）序列中任意两个相邻的多项式在(a，b)内无公共根；
（3）设x
0
∈(a，b)且?
j
(x
0
)＝0，则?
j-1
(x
0
)与?
j+1
(x
0
)反号。
为了计算f
n
(λ)=det(A- λI)，引进多项式序列
nkf
kk
k
k
≤≤
=
1,)(
1
1
32
221
11
λαβ
β
λαβ
βλαβ
βλα
λ
OO
O
2004-12-1 75
矩阵A的特征多项式序列及其性质（续）
即f
k
(λ)表示A的k阶顺序主子矩阵的特征多项式。
规定f
0
(λ)＝1，直接计算有，
λαλ?=
11
)(f
)()()()(
0
2
1122
λβλλαλ fff=
一般地用归纳法可以证明：
nkfff
kkkkk
,,3,2),()()()(
2
2
11
L==
λβλλαλ
序列{f
k
(λ)}称为矩阵A的特征多项式序列。
定理：序列{f
k
(λ)}是f
n
(λ)在(-∞，∞)内的一个Sturm序列。
证明：只需证明{f
k
(λ)}满足定义中的三条即可。
（1）f
0
(λ)＝1，自然满足。
2004-12-1 76
矩阵A的特征多项式序列及其性质（续）
（2）用反证法证明。假设存在某个λ∈(-∞，∞),s.t.
f
j
(λ)= f
j-1
(λ)=0，则由递推关系式知，f
j-2
(λ)=0。同理可知，f
j-3
(λ)=…= f
0
(λ)=0，这与f
0
(λ)≡1矛盾。故不存在任何实数λ，使f
j
(λ)= f
j-1
(λ)=0。
（3）设f
j
(λ)=0，则由递推关系式知，
又因β
j
≠ 0，f
j-1
(λ)≠ 0，故f
j＋1
(λ)≠ 0，且与f
j-1
(λ)反号。
)()(
1
2
1
λβλ
+
=
jjj
ff
证毕定理：设A为对称三对角矩阵，，
{f
k
(λ)}为A的特征多项式序列，则f
k
(λ)＝0有k个单根，
且f
k
(λ)＝0与f
k-1
(λ)＝0的根相互交替。
1,,2,1,0?=≠ ni
i
Lβ
2004-12-1 77
矩阵A的特征多项式序列及其性质（续）
证明：用归纳法。首先容易知道，
λαλ?=
11
)(f
有一个单根λ
1
＝α
1
。
而
2
1122
))(()( βλαλαλ=f
故
0)(,)(
2
1122
<?=+∞=±∞ βαff
故f
2
(λ)=0在(-∞,α
1
)及(α
1
，+∞)内各有一根。这说明
f
1
(λ)与f
2
(λ)符合定理结论。
设对于f
1
(λ)，f
2
(λ)，…,f
k-1
(λ)定理成立，再证f
k
(λ)与
f
k-1
(λ)的零点相互交错，且f
k
(λ)有k个单根。
将f
k-1
(λ)与f
k-2
(λ)的零点分别记为λ
i
(i=1,2,…,k-1)和
μ
i
(i=1,2,…,k-2)，且顺序为：
23211221
,
<<<<<<<<
kkkk
μμμμλλλλ LL
2004-12-1 78
矩阵A的特征多项式序列及其性质（续）
再由归纳法假设，有
1222211
<<<<<<<
kkk
λμλμλμλ L
由于f
k
(λ)＝(-1)
k
λ
k
+λ的低次项，故
,...2,1,
)1()(
0)(
=
=+∞
>?∞
k
signsignf
f
k
k
k
再由归纳法假设有，
f
k-2
(λ
1
)>0，f
k-2
(λ
2
)<0，
f
k-2
(λ
3
)>0，…，
即f
k-2
(λ
j
)的符号为
(-1)
j+1
。
2004-12-1 79
矩阵A的特征多项式序列及其性质（续）
再由
)()(
2
2
1 jkkjk
ff λβλ
=
知f
k
(λ
j
)的符号为(-1)
j
，即证毕
L,0)(,0)(,0)(
21
><>?∞ λλ
kkk
fff
于是在(-∞,λ
1
)，(λ
1
,λ
2
)，…,(λ
k-1
,+∞)内都有f
k
(λ)=0
的根。这里共有k个区间，而f
k
(λ)＝0只有k个根，故在每个区间上有且仅有一个根。
定义：整值函数S
n
(λ)表示数列{f
0
(λ)，f
1
(λ)，…，
f
n
(λ)}中相邻两数符号相同的个数，将数S
n
(λ)称为数列{f
j
(λ)}的同号数。计算同号数时约定，当f
j
(λ)＝0时，
f
j
(λ)与f
j-1
(λ)同号。
2004-12-1 80
矩阵A的特征多项式序列及其性质（续）
定理：设A为实对称三对角矩阵，
{f
k
(λ)}是A的特征多项式序列，则
（1）S
n
(c)表示为f
n
(λ)＝0在区间[c,+∞)上根的个数；
（2）若c<d，则f
n
(λ)＝0在区间[c,d]上根的个数为：
S
n
(c) －S
n
(d)。
1,,2,1,0?=≠ ni
i
Lβ
证明：（2）的结论容易由（1）得出，故下面只证明结论（1）。使用数学归纳法来证明。
{ }{ }ccfcf?=
110
,1)(,)( α
当n＝1时，
当α
1
≥c时，由于f
1
(-∞)>0，f
1
(α
1
)＝0，f
1
(c)≥0，故S
1
(c)＝1，
S
1
(c)就是[c,+∞)上f
1
(λ)＝0的根的个数。
而当α
1
<c时，f
1
(c)<0，故S
1
(c)＝0，f
1
(λ)＝0在区间[c,+∞)
上无根，故S
1
(c)也是[c,+∞)上f
1
(λ)＝0的根的个数。
2004-12-1 81
矩阵A的特征多项式序列及其性质（续）
假设对于{f
0
(λ)，f
1
(λ)，…,f
k-1
(λ)}定理成立，并设
S
k-1
(c)=m，即f
k-1
(λ)＝0的不小于c的零点为m个。将
f
k-1
(λ)＝0与f
k
(λ)＝0的根分别记为λ
i
(i=1,2,…,k-1)和
μ
i
(i=1,2,…,k)，且由归纳法假定有：
<<<<<<<<<
<<<≤<<<<
++
+
111111
12121
μλλμλμλμ
λλλλλλ
LL
LL
mmmkkk
mmkk
c
由此可得到c∈(λ
m+1
，λ
m
]，f
k
(λ)不小于c的零点个数只可能是m或m＋1个。下面分四种情况分别讨论。
（1）c=λ
m
此时f
k
(λ)＝0的不小于c的根为m个，而f
k-1
(c)=0。根据Sturm序列的性质有，f
k
(c)≠0。按约定，f
k-1
(c)与
f
k
(c)不同号，故S
k
(c)＝S
k-1
(c)=m。
2004-12-1 82
矩阵A的特征多项式序列及其性质（续）
（2）c= μ
m＋1
此时f
k
(λ)＝0的不小于c的根为m+1个，而f
k
(c)=0。按约定，f
k-1
(c)与f
k
(c)同号，故S
k
(c)＝S
k-1
(c)+1=m+1。
（3）μ
m＋1
< c<λ
m
此时f
k
(λ)＝0的不小于c的根为m个。将f
k-1
(λ)与f
k
(λ)
写成
=
=
)())(()(
)())(()(
21
1211
λμλμλμλ
λλλλλλλ
kk
kk
f
f
L
L
容易看出，signf
k-1
(c)＝(-1)
k-1-m
，signf
k
(c)＝(-1)
k-m
，
即f
k-1
(c)与f
k
(c)反号，故S
k
(c)与S
k-1
(c)相等，均为m。
2004-12-1 83
矩阵A的特征多项式序列及其性质（续）
（4）λ
m＋1
< c<μ
m＋1
此时f
k
(λ)＝0的不小于c的根为m+1个。将f
k-1
(λ)与
f
k
(λ)写成
=
=
)())(()(
)())(()(
21
1211
λμλμλμλ
λλλλλλλ
kk
kk
f
f
L
L
容易看出，signf
k-1
(c)＝(-1)
k-1-m
，signf
k
(c)＝(-1)
k-(m-1)
，
即f
k-1
(c)与f
k
(c)同号，故S
k
(c)＝S
k-1
(c)+1=m+1。
综合以上四种情况，即可得到定理的结论。
证毕
2004-12-1 84
三、特征值的计算设A为实对称三对角矩阵，，
{f
k
(λ)}是A的特征多项式序列。由特征多项式序列的性质知，若a
0
<b
0
,且S
n
(a
0
)=m，S
n
(b
0
)=m+1，则A的第m个特征值λ
m
∈[a
0
,b
0
)（λ按从小到大的顺序编号）。此时，
令，若S
n
(c
0
)=m，则置a
1
=c
0
,b
1
=b
0
（否则，置a
1
=a
0
,b
1
=c
0
），可得新的区间[a
1
,b
1
)，且
λ
m
∈[a
1
,b
1
)。这样继续做下去，可得一区间序列{[a
k
,b
k
)}
且始终保持λ
m
∈[a
k
,b
k
)。当k充分大时，可用作为λ
m
的近似值。这就是二分法求实对称三对角矩阵特征值的方法。
1,,2,1,0?=≠ ni
i
Lβ
)(
2
1
000
bac +=
)(
2
1
kkk
bac +=
2004-12-1 85
特征值的计算（续）
Remark：二分法求特征值简便，稳定，精度可以很高，且由很大的灵活性，既可以求A的全部特征值，
也可以求A的部分特征值。不论哪一种情形，都是以计算同号数S
n
(λ)为基础的。当A的阶数不高时，一般来说容易实现这一算法。但当A的阶数较高时，计算
f
k
(λ)(k=0,1,2,…,n)的值时，可能产生“溢出”，从而使计算中断。
为了避免“溢出”发生，常采用一种修正的二分法，
其算法过程如下：
定义一个新的序列：
{ }
n
k
k
G
1
)(
=
λ
λα
λ
λ
λ?==
1
0
1
1
)(
)(
)(
f
f
G
2004-12-1 86
特征值的计算（续）
=∞?
=?
≠≠
=
==
0)(,
0)(,
0)(&0)(,
)(
)(
)()()(
)(
)(
)(
1
2
21
1
2
1
1
2
2
11
1
λ
λλα
λλ
λ
β
λα
λ
λβλλα
λ
λ
λ
k
kk
kk
k
k
k
k
kkkk
k
k
k
f
f
ff
G
f
ff
f
f
G
由于序列{f
k
(λ)}相邻两项不同时为零，故f
j
(λ)＝0等价于G
j
(λ)＝0(j=1,2,…,n)。这样，{G
k
(λ)}可以写成如下形式：
2004-12-1 87
特征值的计算（续）
nk
G
G
GG
G
G
k
kk
kk
k
k
k
k
,,3,2,
0)(,
0)(,
0)()(,
)(
)(
1
2
21
1
2
1
L=
=∞?
=?
≠
=
λ
λλα
λλ
λ
β
λα
λ
若f
k
(λ)≠0，f
k
(λ)与f
k-1
(λ)同号等价于G
k
(λ)>0；而当
f
k
(λ)＝0时，按约定f
k
(λ)与f
k-1
(λ)同号，这时G
k
(λ)＝0。
所以，序列{f
k
(λ)}的同号数S
n
(λ)等于序列{G
k
(λ)}中非负项的数目。这样，计算{f
k
(λ)}的同号数S
n
(λ)可以改为计算{G
k
(λ)}的非负项数目。
λαλ?=
11
)(G
2004-12-1 88
特征值的计算（续）
Remark1：当按递推关系式计算{G
k
(λ)}时，若某个| G
j
(λ) |
充分小，可以认为G
j
(λ)＝0；若出现某个G
j
(λ)＝-∞时（由
G
j
(λ)＝0来判断），只要给G
j
(λ)任一负值即可。
Remark2：利用二分法求A的全部特征值时，需要用盖尔圆定理确定A的特征值的范围。设A为实对称三对角矩阵，
，约定β
0
＝0，β
n
＝0，并令
1,,2,1,0?=≠ ni
i
Lβ
++=
=
≤≤
≤≤
}{
}{
1
1
1
1
max
min
iii
ni
iii
ni
d
c
ββα
ββα
则由盖尔圆定理知，所有的特征值λ满足c< λ<d。
这时就知道从那些λ处计算{f
k
(λ)}的同号数或{G
k
(λ)}
的非负项数目，进而进行根的隔离，用二分法进行特征值的计算。
Back
第八章矩阵特征值与特征向量的计算
§8.0 问题描述
§8.1 乘幂法与反幂法
§8.2 雅可比方法
§8.3 QR方法
§8.4 求实对称三对角矩阵特征值的二分法
2004-12-1 2
§8.0 问题描述设A为n×n矩阵，所谓A的特征问题是求数λ和非零向量x，使
Ax＝λx
成立。数λ叫做A的一个特征值，非零向量x叫做与特征值λ对应的特征向量。这个问题等价于求使方程组
(A- λI)x=0有非零解的数λ和相应的非零向量x。
线性代数理论中是通过求解特征多项式det(A- λI)=0的零点而得到λ，然后通过求解退化的方程组(A- λI)x=0
而得到非零向量x。当矩阵阶数很高时，这种方法极为困难。目前用数值方法计算矩阵的特征值以及特征向量比较有效的方法是迭代法和变换法。
Back
2004-12-1 3
§8.1 乘幂法与反幂法一、乘幂法通过求矩阵特征向量求出特征值的一种迭法方法，
它用以求按模最大的特征值和相应的特征向量。
n
λλλ ≥≥≥ L
21
∑
=
=+++==
n
j
jjjnnn
xaxaxaxaAVV
1
222111
0)1(
λλλλ L
）（
设实矩阵A的特征值为λ
1
，λ
2
，…,λ
n
，相应的特征向量线性无关。设A的特征值按模排序为：
n
xxx,,,
21
L
∑
=
=+++=
n
j
jjnn
xaxaxaxaV
1
2211
)0(
L
则对任一非零向量，可以得到：
n
RV ∈
)0(
令，可以构造一个向量序列，
L,2,1,0,
)()1(
==
+
kAVV
kk
2004-12-1 4
乘幂法（续）
∑
=
=+++==
n
j
jjjnnn
xaxaxaxaAVV
1
22
22
2
211
2
1
1)2(
λλλλ L
）（
M
∑
=
=+++==
n
j
jj
k
jnn
k
n
kkkk
xaxaxaxaAVV
1
222111
1)(
λλλλ L
）（
下面分四种情况讨论：
))((
2
1
111
)0()(
∑
=
+==
n
j
jj
k
j
kkk
xaxaVAV
λ
λ
λ
),,3,2(
1
nj
j
L=> λλ
若由于，故k充分大时，
)2(1
1
≥< i
i
λ
λ
,0
1
≠a
111111
)(
)( xaxaV
k
k
kk
λελ ≈+=
1.
2004-12-1 5
乘幂法（续）
)(k
V是相应于的近似特征向量
1
λ
设表示
)(k
l
V
个分量。的第lV
k )(
nl
xa
xa
V
V
lkl
k
lkl
k
k
l
k
l
,,2,1,
1
111
111
1
1
)(
)1(
L=≈
+
+
=
+
++
λ
ελ
ελ
）（）（
）（）（
2.主特征值是实数，但不是单根，如
r
λλλ === L
21 nr
λλλ ≥≥>
+
L
11
而，则
Remark:具体计算时，V
(0)
的选取很难保证一定有
α
1
≠0。但是，由于舍入误差的影响，只要迭代次数足够多，如，就会有
，因而最后结论是成立的。对于的情形，
由于对任意l均有上面的结论，故只要取另外的l使即可。
nn
m
xaxaxaV
′
++
′
+
′
= L
2211
)(
0
1
≠
′
a
0
)(
=
k
l
V
0
)(
≠
k
l
V
2004-12-1 6
乘幂法（续）
））（（
ii
k
n
ri
i
i
r
i
i
kk
xaxaV
∑∑
+==
+=
1
1
1
1
)(
λ
λ
λ
)(k
V
仍然是相对于的近似特征向量。
1
λ
nl
V
V
k
l
k
l
,,2,1,
)(
)1(
1
L=≈
+
λ
Remark:由于相应于λ
1
的特征向量子空间可能不是一维的，由上式得到的V
(k)
只是该子空间的一个特征向量。而且不同的V
(0)
可能得到线性无关的V
(k)
。
3.且为实数，
21
λλ?= ),,4,3(
21
nj
j
L=>= λλλ
n
k
n
n
kkkk
xaxaxaxa )()()1((
1
3
1
3
322111
)(
λ
λ
λ
λ
λ +++?+= LV
由于
( )
*
22111
)1(
k
kk
xaxa ελ +?+=
2004-12-1 7
乘幂法（续）
当k充分大以后，有或
2
1
2
λ≈
+
）（
）（
k
l
k
l
V
V
）（）（k
l
k
l
VV /
2
1
+
≈λ
又可推出：
）（）（kk
VV
1
1
λ+
+
)11
*
2211
1
1
*
122
1
11
1
1 k
kk
k
kk
xaxaxaxa ελελ +?+++?+=
+
+
++
）（（））（（
11
1
1
2 xa
k+
≈ λ
22
1
1
1
1
1
21 xaVV
kkkk +++
≈? λλ）（
）（）（
）（）（kk
VV
1
1
λ?
+
和
）（）（kk
VV
1
1
λ+
+
分别看作
21
λλ和的近似特征向量。
2004-12-1 8
乘幂法（续）
4.
n
λλλλλλ ≥≥>== L
32121
，且此处只限于讨论实矩阵的情形。若有
12121
xxee
ii
====
，而且，必有
θθ
ρλλρλ
∑
=
++=+++=
n
j
jjnn
xaxaxaxaxaxaV
3
11112211
)0(
L
当V
(0)
为实向量时，有
111111
3
111111
)(
xaxaxaxaxaV
kk
n
j
jj
k
j
kkk
λλλλλ +≈++=
∑
=
11
1
111
1
1
)1(
xaxaV
kkk +++
+≈ λλ
11
2
111
2
1
)2(
xaxaV
kkk +++
+≈ λλ
2004-12-1 9
乘幂法（续）
容易验证，
0
)(
11
)1(
11
)2(
≈++?
++ kkk
VVV λλλλ）（
若存在二实数p，q，使得
0
)()1()2(
=++
++ kkk
qVpVV
则可以得到
=++
=++
0)(
0)(
1
2
11
1
2
11
qp
qp
k
k
λλλ
λλλ
上式说明是方程λ
2
＋pλ＋q＝0的一对共轭复根。
这里的系数p、q可由（*）式的分量形式确定，即
11
,λλ
(*)
(**)
jlnjl
qVpVV
qVpVV
k
j
k
j
k
j
k
l
k
l
k
l
≠≤≤
=++
=++
++
++
,,1
0
0
)()1()2(
)()1()2(
2004-12-1 10
乘幂法（续）
解出p和q后，即可由（**）式求得
2
2,1
22
）（
p
qi
p
±?=λ
又因为
≈?
≈?
+
+
22122
)(
1
)1(
11211
)(
2
)1(
)(
)(
xaVV
xaVV
kkk
kkk
λλλλ
λλλλ
故与λ
1
，λ
2
对应的特征向量分别为V
(k+1)
- λ
2
V
(k)
和
V
(k+1)
- λ
1
V
(k)
。
2004-12-1 11
乘幂法（续）
于是得到以下规范化算法：
L,2,1
max/
0
)()()(
)1()(
)0()0(
=
=
=
≠=
k
VVU
AUV
VU
kkk
kk
）（
征向量的方法称为幂法。
k
A
)0(
V上述由已知非零向量及矩阵A的乘幂构造
{ }
）（k
V
1
λ
向量序列以计算A 的特征值及相应的特
)(
)(k
V
其中，max表示向量的最大分量。
)(k
V
以上讨论只是说明了乘幂法的基本原理。当太小或太大时，将会使过小或过大，以致运算无法继续进行。因此，实际计算时，每步要进行规范化。
1
λ
∞
)(k
V
2004-12-1 12
乘幂法（续）
对于规范化的幂法有如下收敛性定理：
证明：取）构造向量序列（00
1
1
)0(
≠≠=
∑
=
axaV
j
n
j
j
其特征值满足
n
λλλ ≥≥> L
21,则按规范化幂法定理：
nn
R
×
∈
设A 有n个线性无关的特征向量，
()
)max(
lim
1
1
x
x
U
k
k
=
∞→
构造的序列有
( )
1
)max(lim λ=
∞→
k
k
V
0)0()1(
AVAUV ==
）（
）（
）（
）（
）（）（
0
0
)1(11
max
max/
AV
AV
VVU ==
2004-12-1 13
==
)1()2(
AUV
）（
）（
）（
0
02
max AV
VA
乘幂法（续）
）（
）（）（）（222
max/ VVU =
)
)max(
)(
max(
1
.
max
)0(
)0(20
02
AV
VAAV
VA
）（
）（
）（
=
)max(
)0(2
)0(2
VA
VA
=
2004-12-1 14
乘幂法（续）
==
）（1)( kk
AUV
）（
）（
）（
01
0
max VA
VA
k
k
M
)max(
)0(
)0(
)(
VA
VA
U
k
k
k
=
][
1
2
111
1
0
j
k
j
n
j
j
k
i
n
i
k
ii
k
xaxaxaVA）（
）（
λ
λ
λλ
∑∑
==
+==
]max[
][max
][
1
1
2
1
111
2
1
111
)(
x
x
xaxa
xaxa
U
n
j
j
j
j
k
n
j
j
j
j
k
k
→
+
+
=
∑
∑
=
=
）（
）（
λ
λ
λ
λ
λ
λ
即
=
∞→
）（K
k
Ulim
]max[
1
1
x
x
而故
2004-12-1 15
=
）（
）（
）（
01
0
)(
max
max)max(
VA
VA
V
k
k
k
+
+
=
∑
∑
=
=
n
j
j
k
j
j
k
n
j
j
k
j
j
k
xaxa
xaxa
2
1
1
11
1
1
2
1
111
]max[
][
max
）（（
）（
λ
λ
λ
λ
λ
λ
1
2
λ
λ
1
maxlim λ=
∞→
）（
）（k
k
V
即，收敛速度依赖于。
]
)(
]max[
[
max
1
2
1
11
2
1
111
∞→→
+
+
=
∑
∑
=
=
k
xaxa
xaxa
n
j
j
kj
j
n
j
j
kj
j
λ
λ
λ
λ
λ
λ
）（（
）（
乘幂法（续）
证毕
2004-12-1 16
乘幂法（续）
Remark：上述定理仅适用于的情形。
n
λλλ ≥≥> L
21
对于规范化幂法，若当分别有两
{ } { }
）（）（
和
122 +kk
UU
个不同的极限，说明。此时若有m，使
21
λλ?=
)1()2()(1
,
+++
==
mmmm
AVVAUV
ε<?
∞
）（）（2mm
UU
，再做两次迭代
2
1
)(
)2(
λ≈
+
m
l
m
l
U
V
12
)(
)2(
1
,
)(
)(
λλλ?==
+
m
l
m
l
U
V
)(
11
)1(
1
)2(
λλ属于xVV
mm
≈+
++
特征向量分别为则有即
( )
)(
22
)1
1
)2(
λλ属于xVV
mm
≈?
++
2004-12-1 17
乘幂法（续）
若当规范化序列无规律可循时，这时说明，
。在m次迭代后，由
{ }
)(k
U
21
λλ =,
)()1( mm
AUV =
+
( ) )2()3()1(2
,
++++
==
mmmm
AVVAVV
做三次迭代，
然后按
=++
=++
+++
+++
0
0
1)2()3(
1)2()3(
）（
）（
m
j
m
j
m
j
m
l
m
l
m
l
qVpVV
qVpVV
计算q，p
再用所得之p，q计算
，
∞
+++
++
)1()2()3( mmm
qVpVV
2004-12-1 18
乘幂法（续）
再由
++
++
2
)2(
1
3
1
)2(
2
)3(
xVV
xVV
mm
mm
作为作为
）（
λ
λ
，计算特征向量。
,
21，
求λ
，）（，就由小于
，
2
21
22
p
qi
p
±?=λε若其值作为初始向量，
( )
）（
）（3
3
)3(
max
+
+
+
=
m
m
m
V
V
U否则以重复以上计算步骤,直至停止。
）（）（
ε<++
∞
+++ )1(23 mmm
qVpVV
2004-12-1 19
二、乘幂法的加速技术
1
2
λ
λ
=r
乘幂法收敛速度依赖于比值。表示满足的那个下标。当r<1但接近于1
时，收敛可能很慢。因此必须研究加速收敛的方法。
0
j
λ
njj
λλλλ ≥≥>==
LL
00
11
1、原点平移法用B=A-PI来代替A解乘幂法，P是一个可选参数。
如果A的特征值为，则B 的相应特征
n
λλλ,
21
L，，
PPP
n
λλλ,,,
21
L
征值为，而且A，B的特征向量相同，XPXPIABX )()(?=?= λ
2004-12-1 20
原点平移法（续）
若A的主特征值为λ
1
，则P的选取要求仍然是B 的主特征值，即
P?
1
λ
),,4,3,2(
1
njPP
j
L=?>? λλ
且使
1
2
1
2
max
λ
λ
λ
λ
<
≤≤
p
p
j
nj
这时对A-PI解乘幂法，就会有较快的收敛速度。
的按模最大的特征值之后，A的按模求得A-PI
1
μ
最大的特征值为,
而A-PI对应于的特征向量就是1
μ
1
λ
A对应于的特征向量。
参数P的选取依赖于对A的特征值分布的大体了解。这可通过盖尔圆定理得到矩阵特征值的分布。
1
λ
2004-12-1 21
原点平移法（续）
定理（Gerschgorin圆盘定理）：设A为n×n实矩阵，则
（1）A的每一特征值必属于下述几个圆盘（称为盖尔圆）
niara
n
ij
j
ijiii
,,2,1,
1
L==≤?
∑
≠
=
λ
的并集之中；
（2）由矩阵A的所有盖尔圆组成的连通部分中任取一个，
若它是由k个盖尔圆构成，则在这个连通部分中有且仅有k个特征值（盖尔圆相重时重复计算，特征值相同也重复计算）。
Remark：原点平移法由于p的选取困难，故很少单独使用。通常将原点平移法与反幂法或QR方法结合使用。
2004-12-1 22
2、Rayleigh商加速法（用于对称矩阵）
定义：设A 为n阶实对称矩阵，x为任意非零向量，则
称为A关于x的Rayleigh商。
），（
）（
xx
xAx
xR
,
)( =
定理：设A 为实对称矩阵，特征值满足
nn
R
×
∈
n
λλλ ≥≥> L
21
,应用规范化幂法公式计算A 的
）（K
U
1
λ
主特征值的规范化向量 的Rayleigh商可给出
1
λ
的较好的近似（对任意的）。0
)0(
≠V
))((
2
1
2
1
k
kk
kk
k
O
UU
UAU
UR
λ
λ
λ +==
），（
），（
）（
）（）（
）（）（
）（
2004-12-1 23
Rayleigh商加速法（续）
∑
=
=
n
j
jj
k
j
k
xaVA
1
)0(
λ
)max(
)0(
)0(
)(
VA
VA
U
k
k
k
=
证明：
由
20001
]/[max)(
）（）（）（）（）（
），（，VAVAVAUAU
kkkkk +
=
），（
），（
）（
）（）（
）（）（
）（
kk
kk
k
UU
UAU
UR =
2000
]/[max)(
）（）（）（）（）（
），（，VAVAVAUU
kkkkk
=
∑
∑
=
=
+
+
==
n
j
k
jj
n
j
k
jj
kk
kk
a
a
VAVA
VAVA
1
22
1
122
00
001
λ
λ
），（
），（
）（）（
）（）（
))((
1
1
2
1
2
1
2
2
1
2
1
2
12
1
2
1
1
k
n
j
k
jj
n
j
k
jj
O
a
a
a
a
λ
λ
λ
λ
λ
λ
λ
λ +=
+
+
∑
∑
=
=
+
）（）（
）（）（
注：此处取的A的特征向量是标准正交化的。
证毕
=
2004-12-1 24
三、反幂法对用乘幂法求解按模最大的特征值是，特
1?
A
n
λ
1
设矩阵A非奇异，用代替A作幂的方法就称为反
1?
A
n
λλλ
111
21
≤≤≤ L
的特征值满足，并且A对应于A
-1
的相应的特征向量是相同的。
n
λλλ ≥≥≥ L
21
1?
A
幂法。当A的特征值满足时，
n
x
征向量是，即是A的按模最小的特征值和特征向量。
计算矩阵按模最小的特征值及相应的特征向量。
反幂法的算法为：取，构造向量序列。
0
)0()0(
≠=UV
=
=
)max(
)(
)(
)(
)1(1)(
k
k
k
kk
V
V
U
UAV
），，（L21=k
2004-12-1 25
反幂法（续）
为避免求逆，可将上面的序列改变为
=
=
)max(
)(
)(
)(
)1()(
k
k
k
kk
V
V
U
UAV
），，（L21=k
当时，由反幂法构造的
0
21
>>≥≥
n
λλλ L
{ }{ }
)()(
,
kk
UV
向量序列满足
n
k
k
n
n
k
k
V
x
x
U
λ
1
)max(lim).2(
)max(
lim).1(
)(
)(
=
=
∞←
∞→
收敛速度的依赖于比值为。
1?n
n
λ
λ
2004-12-1 26
反幂法（续）
Remark：若已知矩阵A的某个特征值λ
i
的相对分离较好的近似值p。不要求p的近似程度有多好，只要求j≠i
时，，则便是的主特征值。
这样一来，结合原点平移法的反幂法的计算公式为：
pp
ji
<? λλ
p
i
λ
1
1
)(
pIA
L,2,1
)max(
)(
0
)(
)(
)(
)1(1)(
)0()0(
=
=
=
≠=
k
V
V
U
UpIAV
VU
k
k
k
kk
2004-12-1 27
反幂法（续）
则向量序列满足
{ }{ }
)()(
,
kk
UV
()
=
=
∞→
∞←
p
V
x
x
U
i
k
k
i
i
k
k
λ
1
maxlim).2(
)max(
lim.1
)(
)(
）（
）（
）（k
k
i
V
p
max
1
lim
∞→
+=λ
其收敛速度依赖于。
p
p
i
j
λ
λ
Back
2004-12-1 28
§8.2 雅可比方法用以求实对称矩阵全部特征值和特征向量的一个变换方法。其基本思想是通过一系列正交变换，化实对称矩阵为对角阵，这个对角阵的对角元素就是该矩阵的特征值，这些正交矩阵的乘积矩阵的各个列就是相应的特征向量。
有关代数知识：
A1.矩阵和相似矩阵的特征值相同。
1?
= PAPB
Q
2.若矩阵满足，则称为正交阵，当然
Q
IQQ
T
=
T
QQ =
1
k
QQQ,,,
21
L
k
QQQQ L
21
=。正交阵的乘积仍为正交阵。
2004-12-1 29
有关代数知识（续）
Q
3.若为实对称矩阵，则存在正交阵使A
∧ A∧=
T
QAQ
。其中为对角阵，的对角元是
∧
的特
i
q
T
Q ),,2,1( ni L= A
征值，而的列就是的特征向量。
4.设为n阶对称矩阵，，其中为正交A
T
QAQC =
Q
22
FF
AC =
。
矩阵，则
5.设为实对称矩阵，为任意正交矩阵，则
A
Q
∑
=
====
n
j
j
F
T
FFF
AQAQAQQAA
1
2
2
222
)(λ
2004-12-1 30
6.Givens旋转矩阵
p
=
1
1
cossin
1
1
sincos
1
1
),,(
O
LLL
MM
MOM
MM
LLL
O
θθ
θθ
θqpR
q
p
称为旋转阵。是一个正交阵，即
),,( θqpR
。
IRR
T
=
T
RR =
1
q
有关代数知识（续）
2004-12-1 31
考虑实矩阵
=
2221
1211
aa
aa
A
)0(
2112
≠=aa
一、古典Jacobi方法
=
θθ
θθ
cossin
sincos
R
令易知，对任何θ，R为正交矩阵。称R为平面旋转矩阵。
22
)1(
2221
1211
)1(
)(
cossin
sincos
cossin
sincos
×
=
==
ij
T
a
aa
aa
RARA
θθ
θθ
θθ
θθ
令由矩阵乘法可得
+?==
+?=
++=
)sin(coscossin)(
coscossin2sin
sincossin2cos
22
121122
)1(
21
)1(
12
2
2212
2
11
)1(
22
2
2212
2
11
)1(
11
θθθθ
θθθθ
θθθθ
aaaaa
aaaa
aaaa
2004-12-1 32
古典Jacobi方法（续）
，当时，可选。
2211
21
2
2
aa
a
tg
=θ
2211
aa =
4
π
θ =
由上式最后一式知，当θ满足时，有
0
)1(
21
)1(
12
== aa
这也就是说，当A为实二阶对称阵时，用适当的正交相似变换一次即可把A化为对角阵。当A为n阶实对称阵时，
要用到n阶旋转矩阵（Givens旋转矩阵），它与单位矩阵的区别仅在于(p,p),(p,q),(q,p),(q,q)四个位置的元素不一样。容易验证，R(p,q,θ)具有如下性质：
（1）R(p,q,θ)为正交阵；
（2）RA只改变A的第p行与第q行元素，AR
T
只改变
A的第p列与第q列元素，RAR
T
只改变A的第p、q行与
p、q列的元素。
2004-12-1 33
古典Jacobi方法（续）
Jacobi方法就是用一系列的旋转相似变换逐渐将A
化为对角阵的过程：
L,2,1,0,
111
0
=
=
=
+++
k
RARA
AA
T
kkkk
恰当地选取每个旋转阵R
k+1
，就可以使A
k
趋于对角化。
设R
k+1
＝R(p,q,θ) 。由R
k+1
正交可知，A
k+1
与A
k
相似，
且A
k
都为是对称阵。A
k+1
与A
k
的差别仅在于p、q行与
p、q列的元素。由矩阵乘法可得：
qpj
aaaa
aaaa
k
jq
k
qj
k
pj
k
qj
k
jp
k
qj
k
pj
k
pj
,,
cossin
sincos
)1()()()1(
)1()()()1(
≠
=+?=
=+=
++
++
θθ
θθ
2004-12-1 34
古典Jacobi方法（续）
=?+?=
+?=
++=
++
+
+
)1(22)()()()1(
2)()(2)()1(
2)()(2)()1(
)sin(coscossin)(
coscossin2sin
sincossin2cos
k
qp
k
pq
k
pp
k
k
pq
k
k
pq
k
pp
k
k
k
pq
k
pp
k
pp
aaaaa
aaaa
aaaa
θθθθ
θθθθ
θθθθ
qpjiaa
k
ij
k
ij
,,,
)()1(
≠=
+
由上面几式容易知道：
qpjiaa
qpjaaaa
k
ij
k
ij
k
qj
k
pj
k
qj
k
pj
,,,)()(
,,)()()()(
2)(2)1(
2)(2)(2)1(2)1(
≠=
≠+=+
+
++
再由正交矩阵的性质（代数知识5）有：
2)(2)(2)(2)1(2)1(2)1(
)(2)()()(2)()(
k
pq
k
k
pp
k
pq
k
k
pp
aaaaaa ++=++
+++
2004-12-1 35
古典Jacobi方法（续）
若，选取θ，使，只需θ满足：
0
)(
≠
k
pq
a 0
)1(
=
+k
pq
a
)()(
)(
2
)2tan(
k
k
pp
k
pq
aa
a
=θ
4
π
θ ≤
这里假定。若，则取
，故有
)()( k
k
pp
aa =
4
)(
)(
π
θ?=
k
pp
asign
=
=
θθ
θ
cos)(sin
2
2
cos
)(k
pq
asign
当时，令
)()( k
k
pp
aa ≠
)()(
)(
2
1
)2tan(
k
k
pp
k
pq
aa
a
d?
==θ
2004-12-1 36
古典Jacobi方法（续）
则
1tan,01tan2tan
22
+±?==?+ ddd θθθ
为避免相近的数相减，取
1
)(
tan
2
++
==
dd
dsign
t θ
θθθ cossin,
1
1
cos
2
t
t
=
+
=
最后得此时有
2)(2)(2)(2)1(2)1(
)(2)()()()(
k
pq
k
k
pp
k
k
pp
aaaaa ++=+
++
引入记号与
∑
≠
=
=
n
ji
ji
ij
aAS
1,
2
)(
∑
=
=
n
i
ii
aAD
1
2
)(
其中D(A)表示对角元的平方和，S(A)表示非对角元的平方和。
2004-12-1 37
古典Jacobi方法（续）
则
( )
()
=
+=
+
+
2
)(
1
2
)(
1
2)()(
2)()(
k
pqkk
k
pqkk
aASAS
aADAD
这说明，只要，则按上述方法构造的旋转矩阵R(p,q,θ)对A
k
变换后，就会使对角线元素平方和增加，非对角线元素平方和减少。
0
)(
≠
k
pq
a
Remark：使用旋转矩阵变换后，虽然，但可能又会变成非零。因此，并不能通过有限次旋转变换就可以将矩阵A化为对角阵。但是，
Jacobi方法确实是收敛的，我们可以通过反复变换将
A化为近似的对角阵。
0
)1(
=
+k
pq
a
)1(
)(
+> kma
m
pq
下面我们来证明古典Jacobi方法的收敛性。
2004-12-1 38
古典Jacobi方法（续）
证明：
我们在计算时是选取，故
)()(
max
k
ij
ji
k
pq
aa
<
=
定理：设A为实对称矩阵，则由A
0
=A出发，经古典Jacobi
方法计算所产生的矩阵序列{ A
k
}收敛于对角阵，且就是A的全部特征值。
),,,(
21 n
diag λλλ L=Λ
n
λλλ,,,
21
L
( )
2
)(
1
2)()(
k
pqkk
aASAS?=
+
由前面的推导过程知，
从而
2)(2
1,
2)(
1,
2)(
))(()()()(
k
pq
n
ji
ji
k
pq
n
ji
ji
k
ijk
annaaAS?=≤=
∑∑
≠
=
≠
=
)(
)1(
1
)(
2)(
k
k
pq
AS
nn
a
≥
)()
)1(
2
1()(
)1(
2
)()(2)()(
2)(
1 kkk
k
pqkk
AS
nn
AS
nn
ASaASAS
=
≤?=
+
)()
)1(
2
1()(
0
1
1
AS
nn
AS
k
k
+
+
≤
当n>2时，显然有，即非对角元的平方和趋于零，A
k
趋于对角阵，Jacobi方法收敛。
0)(
lim
=
∞→
k
k
AS
证毕
2004-12-1 39
特征向量的计算
T
kkkk
TT
RARARARAARRA
12122111
,
=== L，，
由
T
k
T
kkkk
T
kkkk
RRARRRARA
1211
==
T
k
T
k
TT
kk
RRRARRRRR
121121
= LL
IH
T
=
0
T
k
TTT
k
RRRH L
21
=
记
L,2,1,
0
1
=
=
=
k
HAHA
RHH
T
kkk
T
k
T
k
T
k
则有
2004-12-1 40
特征向量的计算（续）
若矩阵序列A
k
收敛于对角阵，
则的第j列就是A对应于的标准正交特征向量的近似值。
T
k
H
j
λ
),,,(
21 n
diag λλλ L=Λ
T
k
H的计算可以使用上述的递推关系式的第一式，
其元素间的关系为：
≠=
=+?=
+=
qpjhh
nihhh
hhh
k
ij
k
ij
k
iq
k
ip
k
iq
k
iq
k
ip
k
ip
,,
,,2,1,cossin
sincos
)1()(
)1()1()(
)1()1()(
Lθθ
θθ
2004-12-1 41
Jacobi方法的计算步骤
1.选主元，即确定p，q(p<q)，使；
)()(
max
k
ij
ji
k
pq
aa
<
=
2.计算cosθ，sinθ；
3.按照递推公式计算新的正交矩阵的元素；
T
k
H
1+
4.计算A
k+1
的元素，其中。0
)1(
=
+k
pq
a
反复执行以上各步，直至为止。这里的ε为给定的误差界。
ε<
+
)(
1k
AS
2004-12-1 42
二、Jacobi过关法使用古典Jacobi方法时，每次先要在非对角元中挑选主元，这要花费很多机时。实用中常采用Jacobi
过关法，方法如下：
1.给定控制误差限ε；
2.计算非对角元素的平方和；
∑∑
+=
=
=
n
ij
ij
n
i
a
1
2
1
1
0
2ν
3.设置一个阀值，如；
0
1
>ν
n
0
1
ν
ν =
4.对A的非对角元a
ij
(i<j)逐个扫描，若某个，则立即对A做一次旋转相似变换，之后对所得新矩阵继续扫描，只要有非对角元的绝对值大于，就用旋转相似变换将其变为零。如此多次扫描和变换，直到每个非对角元；
1
ν
1
ν>
ij
a
1
ν≤
ij
a
2004-12-1 43
Jacobi过关法（续）
5.若，则结束计算，得到特征值与特征向量；
否则转向6；
εν ≤
1
6.缩小阀值，比如用代替，重复4～5。
n
1
ν
1
ν
Remark：Jacobi算法的优点是数值稳定性比较好，精度高，求得的特征向量正交性好。
Jacobi算法的缺点是当Ａ为稀疏矩阵时，Jacobi
变换将破坏其稀疏性质。
Back
2004-12-1 44
§8.3 QR方法一、基本思想
QR方法是求一般矩阵的全部特征值与特征向量的最有效的方法之一。它是在1961年由J.G.F.Francis
首次提出的。该方法的基本原理如下：
设一般的实矩阵A＝(a
ij
)
n×n
可以分解为A＝QR的形式，其中Q为正交阵，R为上三角阵。
1111111112111
,QAQQRQQQRARQAA
TT
=====
令由此可以产生矩阵序列
,
11 kkkkk
RQQRA ==
L,2,1,
111
===
+++
kRQQRA
kkkkk
2004-12-1 45
基本思想（续）
kk
T
kkkk
T
kkkk
QAQQRQQQRA ===
+1由于，故与彼此相似。事实上，中所有矩阵彼此相似，从而它们有相同的特征值。
A
k 1+
A
k
{ }
A
k
k
Q
k
R
序列称为QR序列，其中为正交阵，为上三角阵。
{ }
A
k
{ }
A
k
若趋于一个对角阵或上（下）三角阵，则当k充分大时，的对角元素就可以做为Ａ的特征值。
{ }
A
k
现在的问题是：（1）Ａ的QR分解是否存在；若存在，分解是否唯一？（2）在何种条件下收敛于上述的对角阵或上（下）三角阵？
{ }
A
k
2004-12-1 46
QR方法的实施方法
（1）用反射变换先将A相似化为上Hessenberg阵B，
BAQQ
T
=
此时A与B的特征值相同，即。
（2）用旋转变换对上Hessenberg阵B进行QR分解
RQ
111
RQBB ==
并且有仍为上Hessenberg阵。
（3）形成QR序列
kkk
RQB =一般地，对上Hessenberg阵B作QR分解，
再形成QR序列，
kkk
QRB =
+1
直到B
k＋1
趋于一个上三角阵。
2004-12-1 47
二、反射矩阵（Householder矩阵）
设ν∈R
n
，，称矩阵H＝I－2νν
T
为反射阵。
1
2
=ν
由HIIIH
TTTTTT
=?=?=?= νννννν 2)(2)2(
TTTTTT
TTTTTT
HHIII
IIIIHH
==+?=+?=
==
νννννννννν
νννννννν
44)(44
)2)(2()2()2(
知H是对称阵，也是正交阵。
S
-y
y
x
0
w
w′
x
ν
考虑以ν为法向量，过原点x
0
的超平面S={x|(ν,x)=0,x∈R
n
}（如右图所示）。设w为R
n
中任一向量，则有w=x+y，其中x∈S，
y∈S
⊥
（S的正交补空间）。
2004-12-1 48
反射矩阵（Householder矩阵）（续）
显然，ν与x正交，从而有ν
T
x＝0。
SxxxxIHx
TT
∈=?=?= )(2)2( νννν
由于y∈S
⊥
，故y=cν（c为常数），从而有
cccy
TT
===
2
2
νννν
ycccyyyIHy
TT
=?=?=?=?= ννννννν 2)(2)2(
S
-y
y
x
0
w
w′
x
ν
wyxHyHxyxHHw
′
=?=+=+= )(
从而
Remark：可以看出，w′是w关于S的镜面反射，S
的法向量就是H中的单位向量ν。这就是反射阵的几何意义。
2004-12-1 49
反射矩阵（Householder矩阵）（续）
定理（反射阵的性质）：?x，y∈R
n
，，
总?反射阵H，s.t,Hx＝y。
22
yx =
若x=y，则取ν⊥y，且即可。
1
2
=ν
证明：
今设x≠y，为确定ν，s.t,Hx＝(I-2νν
T
)x=y
证毕只需-2ν(ν
T
x)=y-x即可，即取ν为平行于
y-x的单位向量即可。故可取
2
xy
xy
±=ν
Remark：利用该定理的结论知道，利用反射阵，
可将任一非零向量x变成与另一向量y平行的向量，
且该向量的长度与x的长度相同。
2004-12-1 50
反射矩阵（Householder矩阵）（续）
证明：
显然
2
2
1
1
2
2
1
2
1
)( xxee
n
i
i
==?=?
∑
=
σσ
推论：设x=(x
1
,x
2
,…,x
n
)
T
≠0,则?反射阵H，
s.t,Hx＝-σe
1
，其中
=+==
=+=
==
=
∑
=
xuxu
eexu
xxsignxxsign
uuIH
T
T
n
i
i
T
)(
2
1
)0,,0,1(,
)()(
1
2
2
11
2
1
1
2
1
2
1
1
σσβ
σ
σ
β
L
利用反射阵的性质定理，若记y＝-σe
1
，则必存在反射阵H＝I－2νν
T
，
2
yx
yx
=ν
2004-12-1 51
反射矩阵（Householder矩阵）（续）
使Hx=y＝-σe
1
，而
22
1
1
2
u
u
ex
ex
yx
yx
=
+
+
=
=
σ
σ
ν
T
n
xxxexu ),,,(
211
Lσσ +=+=
其中故
T
T
T
uuI
u
uu
IIH
1
2
2
22
=?=?= βνν
xuxxxxu
T
n
=+=++++== )(])[(
2
1
2
1
1
22
2
2
1
2
2
σσσβ L
证毕
Remark：由该推论知，反射阵H可将一个向量中的多个元素一次性变为零。
2004-12-1 52
反射矩阵（Householder矩阵）（续）
反射变换的具体做法：
1.设x= 计算
,0),,,(
21
≠
T
n
xxx L
2
1
2
1
2
)( xx
n
i
i
=
∑
=
2.计算
∑
=
=
n
i
i
xxsign
1
2
1
2
1
))((σ
3.计算
T
n
xxxexu ),,,(
211
Lσσ +=+=
xux
T
u
=+== )(
2
1
1
2
2
σσβ
4.计算
5.计算
T
uuIH
1?
= β
Remark：若x
1
为零，一般可取其符号函数值为1。
2004-12-1 53
三、平面旋转矩阵（Givens矩阵）
反射阵H对于大量引进零元素是方便的。然而在许多计算中，必须有选择消去一些元素，解决这一问题的有效工具是Jacobi方法中介绍过的Givens变换矩阵。
下面我们给出下文要用到的一个重要定理：
定理：
设已知向量,其中不全为零，则可选择平面旋转矩阵R(i,j,θ)，s.t.
T
n
xxxx ),,,(
21
L=
ji
xx,
T
nji
xxxxRx ),,,,,,(
1
LLL
′′
=
其中
′
=
′
=
=
′
≠+=
′
i
j
i
i
jjii
x
x
x
x
xxxx
θθ sin,cos
0,0)(
2
1
22
2004-12-1 54
平面旋转矩阵（Givens矩阵）（续）
=Rx
证明：令
1
1
cossin
1
1
sincos
1
1
O
LLL
MM
MOM
MM
LLL
O
θθ
θθ
=
x
x
x
x
x
x
x
x
n
j
i
n
j
i
'
'
'
'
1
1
M
M
M
M
M
M
M
M
M
M
M
M
M
M
2004-12-1 55
平面旋转矩阵（Givens矩阵）（续）
有
≠=
′
+?=
′
+=
′
jikxx
xxx
xxx
kk
jij
jii
,,
cossin
sincos
θθ
θθ
令
0cossin =+?=
′
jij
xxx θθ
得0,
cos
sin
tan ≠==
i
i
j
x
x
x
θ
θ
θ
又利用
1cossin
22
=+ θθ
可得
1cos)(cos
222
=+ θθ
i
j
x
x
22
2
2
cos
ji
i
xx
x
+
=θ
22
2
2
sin
ji
j
xx
x
+
=θ
2004-12-1 56
平面旋转矩阵（Givens矩阵）（续）
取
22
cos
ji
i
xx
x
+
=θ
22
sin
ji
j
xx
x
+
=θ
当时，可取
0=
i
x
。1sin,0cos == θθ
此时有
()
=
′
+=
+
+
=+=
′
0
sincos
2
1
22
22
22
j
ji
ji
ji
jii
x
xx
xx
xx
xxx θθ
证毕
2004-12-1 57
平面旋转矩阵（Givens矩阵）（续）
平面旋转变换的具体做法：
①计算
22
jii
xxx +=
′
0=
′
j
x0≠
j
x
对指定的，需变换后有。
②计算
22
cos
ji
i
xx
x
+
=θ
22
sin
ji
j
xx
x
+
=θ
③给出R(i,j,θ)与
T
nji
xxxxRx ),,,,,,(
1
LLL
′′
=
′
=
′
=
=
′
≠+=
′
i
j
i
i
jjii
x
x
x
x
xxxx
θθ sin,cos
0,0)(
2
1
22
其中
2004-12-1 58
四、矩阵的QR分解定理：设A= 为实非奇异矩阵，则A有正交
nnij
a
×
)(
分解A=QR，其中Q为正交阵，R为上三角阵。
且当R具有正对角元时，分解A=QR是唯一的。
下面用反射矩阵证明（也可用旋转矩阵证明）。
证明：
对非奇异矩阵A按列分块，记为
),,,()(
21
)1(
1 nnnij
aaaaAA L===
×
Step1：由于a
1
≠0，故?反射阵H
1
，s.t,H
1
a
1
=-σ
1
e
1
，
2
22
221
)2()2(
2
)2(
2
)2(
22
)2(
1
)2(
121
)1(
1
)1(
21
)1(
111
0
0
0
),,,( A
Dc
Br
aa
aa
aa
aHaHaHAH
nnn
n
n
T
n
=
=
==
σ
σ
L
MMMM
L
L
L
2004-12-1 59
矩阵的QR分解（续）
其中，H
1
的定义见反射阵推论部分。
()0,,
1)2(
2
)2(
222
≠∈=
n
T
n
RaaC L
Step2：由于C
2
≠0，故?（n-1）阶反射阵,s.t.
2
H
′
1222
eCH
′
=
′
σ
1
e
′
其中为n-1维欧式空间的第一个标准正交基向量。
定义，则H
2
仍然为反射阵，且
′
=
2
2
1
H
H
3
33
33
2
)2(
121
)()3(
3
)3(
3
)3(
33
)3(
2
)3(
232
)2(
1
)2(
13
)2(
121
12
A
DC
Br
a
aa
aa
aa
aaa
AHH
n
nnn
n
n
n
=
=
=
σ
σ
σ
σ
L
MOM
L
L
L
2004-12-1 60
矩阵的QR分解（续）
其中，H
2
的定义见反射阵推论部分。
0
2
3
≠∈
n
RC
Stepk：设已经完成对Ａ的第一步至第k步约化过程，即
反射阵H
1
，H
2
，…，H
k –1
，s.t.
==
+
+
+
+?
+?
aaa
aaa
aaa
aaaa
aaaaa
k
nn
k
kn
k
nk
k
kn
k
kk
k
kk
k
nk
k
kk
k
kkk
nkkk
nkkk
kkk
AAHHH
)()(
)1(,
)(
)()(
)1(,
)(
)(
),1(
)(
)1(),1(
)(
),1(1
)3(
2
)3(
)1(,2
)3(
2
)3(
1,22
)2(
1
)2(
)1(,1
)2(
1
)2(
1,1
)2(
121
121
L
MOMM
L
L
MMMMMO
LL
LL
L
σ
σ
σ
2004-12-1 61
矩阵的QR分解（续）
=
kk
kkk
DC
BrR
0
其中，H
i
(i=1,2,…,k-1)的定义见反射阵推论部分。
0
1
≠∈
+?kn
k
RC
由于C
k
≠0，故?（n-k+1）阶反射阵,s.t.
k
H
′
1
~
eCH
kkk
′
=
′
σ
1
~
e
′
其中为（n-k+1）维欧式空间的第一个标准正交基向量，计算反射矩阵的计算公式为：
k
H
′
+=
′
=
+=
′
=
′′
=
′
+
=
+?
∑
)(
2
1
),,,(
)()(
)(
2
2
)()(
,1
)(
2)()(
1
1
k
kkkkkk
Tk
nk
k
kkk
k
kkk
n
ki
k
ik
k
kkk
T
kkkknk
au
aaau
aasign
uuIH
σσβ
σ
σ
β
L
2004-12-1 62
矩阵的QR分解（续）
定义，则H
k
仍然为反射阵，且
′
=
k
k
k
H
I
H
1
′′
=
′
=
==
+?
kkk
kkk
kk
kkk
k
k
kkkkk
DHe
BrR
DC
BrR
H
I
AAHHHAH
1
1
111
~
00 σ
L
这样就使A的三角分解过程前进了一步。
当进行步后，有，且为上三角阵，可记。
)1(?n
nnn
AAHHH =
121
L
n
A
)(n
AR =
由于反射阵H
1
，H
2
，…，H
n –1
都是正交阵，正交阵的乘积、逆阵仍然是正交阵，故若记，则
T
n
TT
HHHQ
121?
= L
QRRHHHA
T
n
TT
==
121
L
2004-12-1 63
矩阵的QR分解（续）
下面证明分解的唯一性。
设A=Q
1
R
1
=Q
2
R
2
，Q
1
、Q
2
均为正交阵，R
1
、R
2
为非奇异上三角阵，且对角元均为正，则
222222
111111
RRRQQRAA
RRRQQRAA
TTTT
TTTT
==
==
由假设（R具有正的对角元）及Cholesky分解的唯一性可得，R
1
＝R
2
。
证毕
Remark：实际计算H
k
时，可按更简单的方法进行，即在公式中直接将写成，
从而得到。
k
u
′
Tk
nk
k
kkk
k
kkk
aaau ),,,,0,,0(
)()(
,1
)(
LL
+
+= σ
T
kkkk
uuIH
1?
= β
2004-12-1 64
矩阵的QR分解（续）
其中每个R
ii
是1×1或2×2 的矩阵。若是1×1的，
其元素就是A的特征值。若是2×2的，则R
ii
的特征值就是A的一对共轭复特征根。
nn
RA
×
∈
对于一般的实矩阵，有如下分解定理：
定理：（实Schur分解定理）设，则?正交阵，s.t.
nn
RA
×
∈
nn
RQ
×
∈
==
mm
m
m
T
R
RR
RRR
RAQQ
MO
L
L
222
11211
Remark：该定理表明，可以通过正交相似变换将A化为块上三角阵。此时可直接求得A的特征值，但很难直接求得定理中的Q和R。
2004-12-1 65
五、Householder方法
=
nnnn
n
n
bb
bbb
bbb
B
1,
22221
11211
MOO
L
L
定义：若时，的方阵B称为上Hessenberg阵，
即
1+> ji
0=
b
ij
如果，则称B为不可约上Hessenberg
阵。类似的也有下Hessenberg阵的定义。
)1,,2,1(0
,1
=≠
+
nib
ii
L
基本思想：Householder方法是先用反射阵将A正交相似变换为上Hessenberg阵，然后对所得之上
Hessenberg阵使用QR方法，这时的QR分解通常用平面旋转矩阵来实现。
2004-12-1 66
Householder方法（续）
定理：设A为n阶实矩阵，则存在反射阵H
1
，
H
2
，…，H
n-2
，s.t.
BHHAHHHH
nnn
=
221132
LL
其中B为上Hessenberg阵。
推论：当A为n阶实对称矩阵时，存在反射阵
H
1
，H
2
，…，H
n-2
，s.t.
=
)(
1
1
)3(
332
2
)2(
221
1
)1(
11
221132
n
nnn
n
nnn
a
a
a
a
HHAHHHH
σ
σ
σ
σσ
σ
OO
O
LL
2004-12-1 67
Householder方法（续）
QR方法的计算过程为：
1.用Householder变换化矩阵为上Hessenberg阵B；
2.用平面旋转变换对B进行QR分解：B＝QR
（可以证明，RQ仍为上Hessenberg阵）；
3.反复进行QR分解，直到矩阵序列B
k
趋于一个对角阵或上三角阵；
4.计算出3中所得之矩阵的对角元，即是原矩阵的特征值。
Remark：详细的计算实施办法还可参考关冶，陈景良编写的《数值计算方法》一书。
2004-12-1 68
六、QR方法的收敛性在A受相当限制的假设下，由QR方法产生的矩阵序列{A
k
}基本上收敛于上三角阵。
定理（QR方法的收敛性）：设且
nn
RA
×
∈
2.存在非奇异矩阵X，使Ａ＝XDX
-1
，D＝diag(λ
1
，
λ
2
，…，λ
n
)，且设X
-1
有三角分解X
-1
＝LU（L为单位下三角阵，U为上三角阵），则QR算法产生的矩阵序列{A
k
}基本上收敛于上三角阵，即
1.Ａ的特征值满足；
0
21
>>>>
n
λλλ L
==
∞→
n
k
k
RA
λ
λ
λ
MO
L
L
*
**
2
1
lim
2004-12-1 69
QR方法的收敛性（续）
这里基本收敛的意思是
)(,0),,,2,1(,
)()(
limlim
jiania
k
ij
k
i
k
ii
k
>===
∞→∞→
Lλ
但不一定存在。
)(,
)(
lim
jia
k
ij
k
<
∞→
推论：若A为n阶实对称矩阵且满足QR方法收敛性定理的条件，则由QR方法产生的矩阵序列{A
k
}
收敛于对角阵D＝diag(λ
1
，λ
2
，…，λ
n
) 。
Remark：对于一般的实矩阵，QR方法产生的矩阵序列{A
k
}的收敛情况很复杂。在某些条件下，{A
k
}收敛于实Schur型的矩阵。
2004-12-1 70
七、带原点平移的QR方法假定A已化为上Hessenberg阵（仍用A表示）。
为了使QR方法收敛的更快，通常使用带原点平移的QR方法，其一般形式为：
=+=
=?
=
+
L,2,1,
1
1
kQRISA
RQISA
AA
kkkk
kkkk
k
S其中称为平移量。这样产生的矩阵序列{A
k
}具有下列性质：
1.由于故有A
k+1
与A
k
相似；
kkkkkkkkk
QAQQIsAQISA
11
1
)(
+
=?+=
2004-12-1 71
带原点平移的QR方法（续）
2,当A
k
为上Hessenberg阵时，若QR分解采用平面旋转变换，则A
k+1
也是上Hessenberg阵。
迭代若干步后，若足够小，就可以认为是特征值的近似值，然后将A
k+1
的第n行、
第n列去掉（剩下的n-1阶矩阵仍为上Hessenberg
阵），对压缩后的矩阵重复上述分解过程。
)1(
1,
+
k
nn
a
)1(
,
+k
nn
a
Remark1：恰当选取平移量可以加速过程的收敛。
的选取通常有两种方法：
k
S
k
S
（1）直接选取；
)(k
nnk
as =
2004-12-1 72
带原点平移的QR方法（续）
（2）当初始矩阵A为实对称矩阵时，经反射变换后得到的上Hessenberg阵实际上是一个对称三对角阵。当再用一系列旋转变换进行分解时，对称三对角性保持不变。可以证明，这种情况下平移量s
k
取二阶矩阵的最接近的那一个特征值，能使很快趋于零。
aa
aa
k
nn
k
nn
k
nn
k
nn
)()(
1,
)(
,1
)(
1,1
a
k
nn
)(
a
k
nn
)(
1,?
Remark2：当A的特征值为复数时，计算过程限定在实数范围内时，应该使用双步QR方法。该方法可参阅聂铁军编写的《数值计算方法》一书。
Back
2004-12-1 73
§8.4 求实对称三对角矩阵特征值的二分法一、问题由上节知道，实对称矩阵通过正交相似变换，
可以化为如下的对称三对角矩阵。
设
=
nn
n
A
αβ
β
αβ
βαβ
βα
1
1
32
221
11
OO
O
且。否则，可将A表示成适当的分块对角阵，每一块的对角元均不为零。
1,,2,1,0?=≠ ni
i
Lβ
对于这类矩阵，可以用二分法求其全部的特征值。
2004-12-1 74
二、矩阵A的特征多项式序列及其性质定义：称实系数多项式序列?
0
(x),?
1
(x),…,?
m
(x)为?
m
(x)
在(a，b)区间内的一个施铎姆（Sturm）序列，如果它具有以下三个性质：
（1）?
0
(x)在(a，b)区间内无实根；
（2）序列中任意两个相邻的多项式在(a，b)内无公共根；
（3）设x
0
∈(a，b)且?
j
(x
0
)＝0，则?
j-1
(x
0
)与?
j+1
(x
0
)反号。
为了计算f
n
(λ)=det(A- λI)，引进多项式序列
nkf
kk
k
k
≤≤
=
1,)(
1
1
32
221
11
λαβ
β
λαβ
βλαβ
βλα
λ
OO
O
2004-12-1 75
矩阵A的特征多项式序列及其性质（续）
即f
k
(λ)表示A的k阶顺序主子矩阵的特征多项式。
规定f
0
(λ)＝1，直接计算有，
λαλ?=
11
)(f
)()()()(
0
2
1122
λβλλαλ fff=
一般地用归纳法可以证明：
nkfff
kkkkk
,,3,2),()()()(
2
2
11
L==
λβλλαλ
序列{f
k
(λ)}称为矩阵A的特征多项式序列。
定理：序列{f
k
(λ)}是f
n
(λ)在(-∞，∞)内的一个Sturm序列。
证明：只需证明{f
k
(λ)}满足定义中的三条即可。
（1）f
0
(λ)＝1，自然满足。
2004-12-1 76
矩阵A的特征多项式序列及其性质（续）
（2）用反证法证明。假设存在某个λ∈(-∞，∞),s.t.
f
j
(λ)= f
j-1
(λ)=0，则由递推关系式知，f
j-2
(λ)=0。同理可知，f
j-3
(λ)=…= f
0
(λ)=0，这与f
0
(λ)≡1矛盾。故不存在任何实数λ，使f
j
(λ)= f
j-1
(λ)=0。
（3）设f
j
(λ)=0，则由递推关系式知，
又因β
j
≠ 0，f
j-1
(λ)≠ 0，故f
j＋1
(λ)≠ 0，且与f
j-1
(λ)反号。
)()(
1
2
1
λβλ
+
=
jjj
ff
证毕定理：设A为对称三对角矩阵，，
{f
k
(λ)}为A的特征多项式序列，则f
k
(λ)＝0有k个单根，
且f
k
(λ)＝0与f
k-1
(λ)＝0的根相互交替。
1,,2,1,0?=≠ ni
i
Lβ
2004-12-1 77
矩阵A的特征多项式序列及其性质（续）
证明：用归纳法。首先容易知道，
λαλ?=
11
)(f
有一个单根λ
1
＝α
1
。
而
2
1122
))(()( βλαλαλ=f
故
0)(,)(
2
1122
<?=+∞=±∞ βαff
故f
2
(λ)=0在(-∞,α
1
)及(α
1
，+∞)内各有一根。这说明
f
1
(λ)与f
2
(λ)符合定理结论。
设对于f
1
(λ)，f
2
(λ)，…,f
k-1
(λ)定理成立，再证f
k
(λ)与
f
k-1
(λ)的零点相互交错，且f
k
(λ)有k个单根。
将f
k-1
(λ)与f
k-2
(λ)的零点分别记为λ
i
(i=1,2,…,k-1)和
μ
i
(i=1,2,…,k-2)，且顺序为：
23211221
,
<<<<<<<<
kkkk
μμμμλλλλ LL
2004-12-1 78
矩阵A的特征多项式序列及其性质（续）
再由归纳法假设，有
1222211
<<<<<<<
kkk
λμλμλμλ L
由于f
k
(λ)＝(-1)
k
λ
k
+λ的低次项，故
,...2,1,
)1()(
0)(
=
=+∞
>?∞
k
signsignf
f
k
k
k
再由归纳法假设有，
f
k-2
(λ
1
)>0，f
k-2
(λ
2
)<0，
f
k-2
(λ
3
)>0，…，
即f
k-2
(λ
j
)的符号为
(-1)
j+1
。
2004-12-1 79
矩阵A的特征多项式序列及其性质（续）
再由
)()(
2
2
1 jkkjk
ff λβλ
=
知f
k
(λ
j
)的符号为(-1)
j
，即证毕
L,0)(,0)(,0)(
21
><>?∞ λλ
kkk
fff
于是在(-∞,λ
1
)，(λ
1
,λ
2
)，…,(λ
k-1
,+∞)内都有f
k
(λ)=0
的根。这里共有k个区间，而f
k
(λ)＝0只有k个根，故在每个区间上有且仅有一个根。
定义：整值函数S
n
(λ)表示数列{f
0
(λ)，f
1
(λ)，…，
f
n
(λ)}中相邻两数符号相同的个数，将数S
n
(λ)称为数列{f
j
(λ)}的同号数。计算同号数时约定，当f
j
(λ)＝0时，
f
j
(λ)与f
j-1
(λ)同号。
2004-12-1 80
矩阵A的特征多项式序列及其性质（续）
定理：设A为实对称三对角矩阵，
{f
k
(λ)}是A的特征多项式序列，则
（1）S
n
(c)表示为f
n
(λ)＝0在区间[c,+∞)上根的个数；
（2）若c<d，则f
n
(λ)＝0在区间[c,d]上根的个数为：
S
n
(c) －S
n
(d)。
1,,2,1,0?=≠ ni
i
Lβ
证明：（2）的结论容易由（1）得出，故下面只证明结论（1）。使用数学归纳法来证明。
{ }{ }ccfcf?=
110
,1)(,)( α
当n＝1时，
当α
1
≥c时，由于f
1
(-∞)>0，f
1
(α
1
)＝0，f
1
(c)≥0，故S
1
(c)＝1，
S
1
(c)就是[c,+∞)上f
1
(λ)＝0的根的个数。
而当α
1
<c时，f
1
(c)<0，故S
1
(c)＝0，f
1
(λ)＝0在区间[c,+∞)
上无根，故S
1
(c)也是[c,+∞)上f
1
(λ)＝0的根的个数。
2004-12-1 81
矩阵A的特征多项式序列及其性质（续）
假设对于{f
0
(λ)，f
1
(λ)，…,f
k-1
(λ)}定理成立，并设
S
k-1
(c)=m，即f
k-1
(λ)＝0的不小于c的零点为m个。将
f
k-1
(λ)＝0与f
k
(λ)＝0的根分别记为λ
i
(i=1,2,…,k-1)和
μ
i
(i=1,2,…,k)，且由归纳法假定有：
<<<<<<<<<
<<<≤<<<<
++
+
111111
12121
μλλμλμλμ
λλλλλλ
LL
LL
mmmkkk
mmkk
c
由此可得到c∈(λ
m+1
，λ
m
]，f
k
(λ)不小于c的零点个数只可能是m或m＋1个。下面分四种情况分别讨论。
（1）c=λ
m
此时f
k
(λ)＝0的不小于c的根为m个，而f
k-1
(c)=0。根据Sturm序列的性质有，f
k
(c)≠0。按约定，f
k-1
(c)与
f
k
(c)不同号，故S
k
(c)＝S
k-1
(c)=m。
2004-12-1 82
矩阵A的特征多项式序列及其性质（续）
（2）c= μ
m＋1
此时f
k
(λ)＝0的不小于c的根为m+1个，而f
k
(c)=0。按约定，f
k-1
(c)与f
k
(c)同号，故S
k
(c)＝S
k-1
(c)+1=m+1。
（3）μ
m＋1
< c<λ
m
此时f
k
(λ)＝0的不小于c的根为m个。将f
k-1
(λ)与f
k
(λ)
写成
=
=
)())(()(
)())(()(
21
1211
λμλμλμλ
λλλλλλλ
kk
kk
f
f
L
L
容易看出，signf
k-1
(c)＝(-1)
k-1-m
，signf
k
(c)＝(-1)
k-m
，
即f
k-1
(c)与f
k
(c)反号，故S
k
(c)与S
k-1
(c)相等，均为m。
2004-12-1 83
矩阵A的特征多项式序列及其性质（续）
（4）λ
m＋1
< c<μ
m＋1
此时f
k
(λ)＝0的不小于c的根为m+1个。将f
k-1
(λ)与
f
k
(λ)写成
=
=
)())(()(
)())(()(
21
1211
λμλμλμλ
λλλλλλλ
kk
kk
f
f
L
L
容易看出，signf
k-1
(c)＝(-1)
k-1-m
，signf
k
(c)＝(-1)
k-(m-1)
，
即f
k-1
(c)与f
k
(c)同号，故S
k
(c)＝S
k-1
(c)+1=m+1。
综合以上四种情况，即可得到定理的结论。
证毕
2004-12-1 84
三、特征值的计算设A为实对称三对角矩阵，，
{f
k
(λ)}是A的特征多项式序列。由特征多项式序列的性质知，若a
0
<b
0
,且S
n
(a
0
)=m，S
n
(b
0
)=m+1，则A的第m个特征值λ
m
∈[a
0
,b
0
)（λ按从小到大的顺序编号）。此时，
令，若S
n
(c
0
)=m，则置a
1
=c
0
,b
1
=b
0
（否则，置a
1
=a
0
,b
1
=c
0
），可得新的区间[a
1
,b
1
)，且
λ
m
∈[a
1
,b
1
)。这样继续做下去，可得一区间序列{[a
k
,b
k
)}
且始终保持λ
m
∈[a
k
,b
k
)。当k充分大时，可用作为λ
m
的近似值。这就是二分法求实对称三对角矩阵特征值的方法。
1,,2,1,0?=≠ ni
i
Lβ
)(
2
1
000
bac +=
)(
2
1
kkk
bac +=
2004-12-1 85
特征值的计算（续）
Remark：二分法求特征值简便，稳定，精度可以很高，且由很大的灵活性，既可以求A的全部特征值，
也可以求A的部分特征值。不论哪一种情形，都是以计算同号数S
n
(λ)为基础的。当A的阶数不高时，一般来说容易实现这一算法。但当A的阶数较高时，计算
f
k
(λ)(k=0,1,2,…,n)的值时，可能产生“溢出”，从而使计算中断。
为了避免“溢出”发生，常采用一种修正的二分法，
其算法过程如下：
定义一个新的序列：
{ }
n
k
k
G
1
)(
=
λ
λα
λ
λ
λ?==
1
0
1
1
)(
)(
)(
f
f
G
2004-12-1 86
特征值的计算（续）
=∞?
=?
≠≠
=
==
0)(,
0)(,
0)(&0)(,
)(
)(
)()()(
)(
)(
)(
1
2
21
1
2
1
1
2
2
11
1
λ
λλα
λλ
λ
β
λα
λ
λβλλα
λ
λ
λ
k
kk
kk
k
k
k
k
kkkk
k
k
k
f
f
ff
G
f
ff
f
f
G
由于序列{f
k
(λ)}相邻两项不同时为零，故f
j
(λ)＝0等价于G
j
(λ)＝0(j=1,2,…,n)。这样，{G
k
(λ)}可以写成如下形式：
2004-12-1 87
特征值的计算（续）
nk
G
G
GG
G
G
k
kk
kk
k
k
k
k
,,3,2,
0)(,
0)(,
0)()(,
)(
)(
1
2
21
1
2
1
L=
=∞?
=?
≠
=
λ
λλα
λλ
λ
β
λα
λ
若f
k
(λ)≠0，f
k
(λ)与f
k-1
(λ)同号等价于G
k
(λ)>0；而当
f
k
(λ)＝0时，按约定f
k
(λ)与f
k-1
(λ)同号，这时G
k
(λ)＝0。
所以，序列{f
k
(λ)}的同号数S
n
(λ)等于序列{G
k
(λ)}中非负项的数目。这样，计算{f
k
(λ)}的同号数S
n
(λ)可以改为计算{G
k
(λ)}的非负项数目。
λαλ?=
11
)(G
2004-12-1 88
特征值的计算（续）
Remark1：当按递推关系式计算{G
k
(λ)}时，若某个| G
j
(λ) |
充分小，可以认为G
j
(λ)＝0；若出现某个G
j
(λ)＝-∞时（由
G
j
(λ)＝0来判断），只要给G
j
(λ)任一负值即可。
Remark2：利用二分法求A的全部特征值时，需要用盖尔圆定理确定A的特征值的范围。设A为实对称三对角矩阵，
，约定β
0
＝0，β
n
＝0，并令
1,,2,1,0?=≠ ni
i
Lβ
++=
=
≤≤
≤≤
}{
}{
1
1
1
1
max
min
iii
ni
iii
ni
d
c
ββα
ββα
则由盖尔圆定理知，所有的特征值λ满足c< λ<d。
这时就知道从那些λ处计算{f
k
(λ)}的同号数或{G
k
(λ)}
的非负项数目，进而进行根的隔离，用二分法进行特征值的计算。
Back
