第 章 矩阵特征值的数值解法

第9章 矩阵特征值的数值解法
9.1 引言
矩阵特征值问题有广泛的应用背景. 例如动力系统和结构系统中的振动问题、电力系统的静态稳定分析上、工程设计中的某些临界值的确定等，都归结为矩阵特征值问题. 数学中诸如方阵的对角化及解微分方程组等问题，都要用到特征值的理论. 本章介绍n 阶实矩阵n n ⨯∈R A 的特征值与特征向量的数值解法.
定义9.1.1 已知n 阶实矩阵()n n ij a ⨯=∈R A ，如果存在常数λ和非零向量x ，使
λ=Ax x 或 ()λ-=0A I x (9.1.1)
那么称λ为A 的特征值(eigenvalue)，x 为A 的相应于λ的特征向量(eigenvector). 多项式
11121212221
2
()det()n n n n n nn a a a a a a p a a a λ
λλλλ-⎡⎤⎢⎥-⎢
⎥=-=⎢⎥
⎢⎥-⎣⎦
L L M M O M L
A I (9.1.2) 称为特征多项式(characteristic polynomial)，
det()0λ-=A I (9.1.3)
称为特征方程(characteristic equation).
注 式(9.1.3)是以λ为未知量的一元n 次代数方程，
()det()n p λλ=-A I 是λ的n 次多项式. 显然，A 的特征值就是特征方程(9.1.3)的根. 特征方程(9.1.3)在复数范围内恒有解，其个数为方程的次数(重根按重数计算)，因此n 阶矩阵A 在复数范围内有n 个特征值. 除特殊情况 (如2,3n =或A 为上(下)三角矩阵)外，一般不通过直接求解特征方程(9.1.3)来求A 的特征值, 原因是这样的算法往往不稳定. 在计算上常用的方法是幂法与反幂法和相似变换方法. 本章只介绍求矩阵特征值与特征向量的这两种基本方法. 为此将一些特征值和特征向量的性质列在此处.
定理9.1.2 设n 阶方阵()ij n n a ⨯=A 的特征值为12,,,n λλλL ，那么 (1) 121122n nn a a a λλλ+++=+++L L ; (2) 12det n λλλ=L A .
定理9.1.3 如果λ是方阵A 的特征值，那么 (1) k λ是k A 的特征值，其中k 是正整数;
(2) 当A 是非奇异阵时，1λ
是1-A 的特征值. (3) ()n p λ是()n p A 的特征值，其中()n p x 是多项式 2012()n n n p x a a x a x a x =++++L .
定义9.1.4 设,A B 都是n 阶方阵. 若有n 阶非奇异阵P ，使得1-=P AP B ，则称矩阵A 与B 相似(similar)，1-P AP 称为对A 进行相似变换(similarity
transformation)，P 称为相似变换矩阵(similarity transformation matrix).
定理9.1.5 若矩阵A 与B 相似，则A 与B 的特征值相同. 定理9.1.6 如果A 是n 阶正交矩阵，那么 (1) 1T -=A A ，且det 1=A 或1-;
(2) 若=y Ax ，则22=y x , 即T T ⋅=⋅x x y y .
定理9.1.7 设A 是任意n 阶实对称矩阵，则 (1) A 的特征值都是实数；
(2) A 有n 个线性无关的特征向量.
定理9.1.8 设A 是任意n 阶实对称矩阵，则必存在n 阶正交矩阵P ，使得
1T -==P AP P AP Λ，
其中12diag(,,,
)n λλλ=L Λ是以A 的n 个特征值12,,,n λλλL 为对角元素的对角矩阵.
定理9.1.9 (圆盘定理) 矩阵()ij n n a ⨯=A 的任意一个特征值至少位于复平面上的几个圆盘
中的一个圆盘上。

9.2 幂法与反幂法
9.2.1 幂法及其加速
9.2.1.1 幂法
幂法是计算矩阵按模最大特征值(largest eigenvalue in magnitude)及相应特征向量的迭代法. 该方法稍加修改，也可用来确定其他特征值. 幂法的一个很有用的特性是：它不仅可以求特征值，而且可以求相应的特征向量. 实际上，幂法经常用来求通过其他方法确定的特征值的特征向量. 下面探讨幂法的具体过程.
设矩阵n n ⨯∈R A 的n 个特征值满足
1230n λλλλ>≥≥≥≥L , (9.2.1) 且有相应的n 个线性无关的特征向量12,,,n L x x x ，则12,,,n L x x x 构成n 维向量空间
n
R 的一组基, 因此1,,1,2,,n
n
i i i i i n αα=⎧⎫
==∈=⎨⎬⎩⎭∑R R L z z x .
在n R 中选取某个满足10α≠的非零向量
01n
i i i α==∑z x .
用矩阵A 同时左乘上式两边，得
01
1
n
n
i i i i i i i ααλ====∑∑Az Ax x .
再用矩阵A 左乘上式两边，得
2
201n
i i i i αλ==∑A z x .
这样继续下去，一般地有
01,1,2,.n
k
k i i i i k αλ===∑L A z x (9.2.2)
记10,1,2,k k k k -===L z Az A z ，则由式(9.2.2)得
0111121, 1,2,.k
n
n k k k i
k i i i i i i i k λαλλααλ==⎡⎤⎛⎫⎢⎥===+= ⎪⎢⎥⎝⎭⎣⎦
∑∑L z A z x x x (9.2.3)
由假设(9.2.1)，结合式(9.2.3)，得
111
lim
,k
k k αλ→∞
=z x (9.2.4)
于是对充分大的k 有
111.k k λα≈z x (9.2.5)
式(9.2.4)表明随着k 的增大，序列{}1k
k λz 越来越接近A 的对应于特征值1λ的特征
向量1x 的1α倍, 由此可确定对应于1λ的特征向量1x . 当k 充分大时，可得1x 的近似值.
上述收敛速度取决于比值21λλ. 事实上，由式(9.2.3)知，
32
112
23
311
1
1
k
k
k
k
n n
n k
λλλααααλλλλ≤++++L z x x x x . (9.2.6)
再由式(9.2.1)得
3
2111
1n λλλλλλ>
≥≥≥L . (9.2.7)
结合式(9.2.6)和式(9.2.7)知，序列{}1k
k λz 收敛速度取决于比值21λλ.
下面计算1λ. 由式(9.2.3)知
当k 充分大时, 11111k k λα++≈z x . 结合式(9.2.5)，得
11k k λ+≈z z .
这表明两个相邻向量大体上只差一个常数倍，这个倍数就是A 的按模最大特征值1λ. 记()
T (1)
(2)(),,,n k k
k
k
z z z
=L z , 则有
()11()
lim ,1,2,,j k j k k
z j n z λ+→∞==L ， (9.2.8)
即两个相邻的迭代向量所有对应分量的比值收敛到1λ.
定义9.2.1 上述由已知非零向量0z 及矩阵A 的乘幂k A 构造向量序列{}k z 来计算A 的按模最大特征值1λ及相应特征向量的方法称为幂法(power method)，其收敛速度由比值21γλλ=来确定，γ越小，收敛越快.
注 由幂法的迭代过程(9.2.3)容易看出，如果11λ>（或11λ<），那么迭代向量k z 的各个非零的分量将随着k →∞而趋于无穷（或趋于零），这样在计算机上实现时就可能上溢(或下溢). 为了克服这个缺点，需将每步迭代向量进行规范化: 记
()T
(1)(2)()
1,,,n k k k k k y y y -==L y Az .
若存在k y 的某个分量0
()j k y ，满足0()()1max j j k k j n
y y ≤≤=，则记0
()max()j k k y =y . 将k y 规范化
max()k k k =z y y ,
这样就把k z 的分量全部控制在[1,1]-中.
例如，设T (2,3,0,5,1)k =--y ，因为k y 的所有分量中，绝对值最大的的是5-，所以max()5k =-y ，故T max()(0.4,0.6,0,1,0.2)k k k ==--z y y .
综上所述，得到下列算法：
算法9.2.1 (幂法) 设A 是n 阶实矩阵，取初始向量0n ∈R z ，通常取T 0(1,1,,1)=L z ，其迭代过程是：对1,2,k =L ，有
1,max(),.
k k k k k
k k m m -=⎧⎪
=⎨⎪=⎩y Az y z y (9.2.9) 定理9.2.1 对式{}k z 和{}k m 有
()
1
1lim max k k →∞
=
x z x , 1lim k k m λ→∞=， ()
其收敛速度由21γλλ=确定. 证明 由迭代过程(9.2.9)知
001
max()k
k
k k i
i m
====∏L A z A z A z , (9.2.11)
其中01n
i i i α==∑z x . 若10α≠，则由(9.2.3)知：011121k
n k k
i i i i λλααλ=⎡⎤⎛⎫⎢⎥=+ ⎪⎢⎥⎝⎭⎣⎦
∑A z x x , 代入
式(9.2.11)得
11211121max k
n
i i i
i k k
n i i i i λααλλααλ==⎛⎫
+ ⎪⎝⎭=⎛⎫⎛⎫
⎪
+ ⎪ ⎪⎝⎭⎝⎭
∑∑x x z x x ， (9.2.12) 故
()
1
1lim max k k →∞
=
x z x . (9.2.13)
而
111211
1
11121max k
n k i
i i i k n k i i i i a a a a λλλλλλ=--=⎡⎤⎛⎫+⎢⎥
⎪⎢⎥⎝⎭⎣⎦=⎛⎫⎛⎫ ⎪+ ⎪ ⎪⎝⎭⎝⎭
∑∑x x x x , (9.2.14) 于是
11211
1
1121max max()max k
n i i i i k k k n i i i i a a m a a λλλλλ=-=⎛⎫⎛⎫
⎪
+ ⎪ ⎪⎝⎭⎝⎭==⎛⎫⎛⎫
⎪
+ ⎪ ⎪⎝⎭⎝⎭
∑∑x x y x x , (9.2.15) 故
1lim k k m λ→∞
=. (9.2.16)
由式(9.2.6)和式(9.2.7)知：上述收敛速度由21γλλ=确定. 证毕！ 例9.2.1 用幂法求方阵
按模最大特征值及相应的特征向量，要求3110k k m m ---<.
解 选取初始向量T 0(1,1,1)=z ，按式(9.2.9)迭代，结果见表
表9.2.1
因此，所求按模最大特征值18.3589λ≈，相应特征向量T 1(0.559817,0.559817,1)≈x . 事实上，A 的按模最大特征值148.3588989λ=≈L ，相应特征向量T 1(0.55981649,0.55981649,1)=L L x , 故所得结果具有较高的精度.
9.2.1.2 幂法的加速 从上面的讨论可知，由幂法求按模最大特征值，可归结为求数列{}k m 的极限值，其收敛速度由21γλλ=确定. 当21γλλ=接近1时, 收敛速度相当缓慢. 为了提高收敛速度, 可以利用外推法进行加速.
因为序列{}k m 的收敛速度由21γλλ=确定，所以若{}k m 收敛，当k 充分大时，
则有()
121k k m O λλλ⎡⎤
-=⎢⎥⎣⎦
，或改写为
()1
2
1
k
k m c λλλ-≈,
其中c 是与k 无关的常数. 由此可得
112
11
k k m m λλλλ+-≈-, (9.2.17)
这表明幂法是线性收敛的. 由式(9.2.17)知
1121
111
k k k k m m m m λλλλ+++--≈--.
由上式解出1λ，并记为2k m +%，即
2221
122121()22k k k k k k k k k k k k k
m m m m m m m m m m m m m ++++++++--=
=+-+-+%， (9.2.18) 这就是计算按模最大特征值的加速公式.
将上面的分析归结为如下算法：
算法9.2.2 (幂法的加速) 设A 是n 阶实矩阵，给定非零初始向量0n ∈R z ，通常取T 0(1,1,,1)=L z . 对1,2k =，用迭代式 求出12,m m 及12,z z . 再对3,4,k =L ，迭代过程为
1212212,max(),(),2.k k k k k k k k k k k k k k m m m m m m m m m
------=⎧⎪=⎪⎪
-⎨
=-⎪-+⎪
⎪=⎩%%y Az y z y () 当1k k m m ε--<%%(0ε>是预先给定的精度) 时，迭代结束，并计算k z ；否则继续迭代，直至满足迭代停止条件1k k m m ε--<%%.
有关加速收敛技术，读者请参考文献[11].
8.2.2 反幂法及其加速
反幂法是计算矩阵按模最小特征值及相应特征向量的迭代法, 其基本思想
是：设矩阵n n ⨯∈R A 非奇异，用其逆矩阵1-A 代替A , 矩阵A 的按模最小特征值就是矩阵1-A 的按模最大特征值. 这样用1-A 代替A 做幂法，即可求出1-A 的按模最大特征值，也就是矩阵A 的按模最小特征值；这种方法称为反幂法(inverse power method).
因为矩阵A 非奇异，所以由i i i λ=Ax x 可知：11
i i i
λ-=A x x . 这说明：如果A 的
特征值满足
1210n-n λλλλ≥≥≥>>L , (9.2.20)
那么1-A 的特征值满足
1
2
1
1
1
1
1
n
n-λλλλ>
≥≥
≥
L , (9.2.21)
且对A 的应于特征值i λ的特征向量i x 也是1-A 的对应于特征值1i λ的特征向量.
由上述分析知：对1-A 应用幂法求按模最大的特征值1n λ及相应的特征向量
n x ，就是求A 的按模最小的特征值n λ及相应的特征向量n x .
算法9.2.3 (反幂法) 任取初始非零向量0n ∈R z ，通常取T 0(1,1,,1)=L z . 为了避免求1-A ，对1,2,k =L
1,max(),.
k k k k k
k k m m -=⎧⎪
=⎨⎪=⎩Ay z y z y (9.2.22) 仿定理9.2.1的证明，可得： 定理9.2.2 对式{}k z 和{}k m 有
1
lim ,
lim max()
n
k k k k n n
m λ→∞
→∞
=
=
x z x , ()
其收敛速度由1n n γλλ-=%确定.
注 按(9.2.22)进行计算，每次迭代都需要解一个方程组1k k -=Ay z . 若利用三角分解法求解方程组，即=A LU ，其中L 是下三角矩阵，U 是上三角矩阵，这样每次迭代只需解两个三角方程组
9.2.3 原点平移法
为了提高收敛速度，下面介绍加速收敛的原点平移法.
设矩阵p =-B A I ，其中p 是一个待定的常数，A 与B 除主对角线上的元素外，其他元素相同. 设A 的特征值为12,,,n λλλL ，则B 的特征值为
12,,,n p p p λλλ---L ,且A 与B 的特征向量相同.
9.2.3.1 原点平移法的幂法
设1λ是A 的按模最大的特征值，选择p ，使
12,3,4,,i p p p i n λλλ->-≥-=L ，
及
22
11
p p λλλλ-<-. (9.2.24) 对B 应用幂法，得
算法9.2.4 对1,2,k =L ，有
1(),max(),
,
k k k k k
k k p m m -=-⎧⎪
=⎨⎪=⎩y A I z y z y (9.2.25) 且
()
1
11lim ,lim ,max k k k k m p λ→∞
→∞
=-=
x z x ()
其收敛速度由21()()p p λλ--确定.
由式(9.2.24)知：在计算B 的按模最大特征值1p λ-的过程(9.2.25)中，收敛速度得到加速；算法9.2.4又称为原点平移下的幂法(power method with shift).
9.2.3.2 原点平移下的反幂法
设n λ是A 的按模最小的特征值，选择p ，使
1,1,2,,2n n i p p p i n λλλ--<-≤-=-L . (9.2.27)
及
11
n n n n p p λλ
λλ---<-. (9.2.28)
若矩阵p =-B A I 可逆，则1-B 的特征值为 且有
1111
,1,2,,2n n i i n p p p
λλλ->≥=----L . (9.2.29) 对B 应用反幂法，得： 算法9.2.5 对1,2,k =L ，有
1(),max(),
,
k k k k k
k k p m m --=⎧⎪
=⎨⎪=⎩A I y z y z y () 且
()
111
lim ,lim ,max k k k k n m p λ→∞
→∞=
=-x z x ()
其收敛速度由1())n n p p λλ---确定.
由式(9.2.28)知：在计算1-B 的按模最大特征值
1
n p
λ-的过程(9.2.30)中，收敛速度得到加速. 算法9.2.5又称为原点平移下的反幂法 (inverse power method with shift).
定义9.2.8 原点平移下的幂法与原点平移下的反幂法统称为原点平移法. 注 有的资料上的原点平移法专指原点平移下的反幂法；而有的资料上的反幂法指的就是原点平移下的反幂法.
原点平移下的反幂法(算法9.2.7)的主要应用是：已知矩阵的近似特征值后，
求矩阵的特征向量. 其主要思想是：若已知A 的特征值m λ的近似特征值为m λ%，则m A I λ-%的特征值就是(1,2,,)i m
i m λλ-=%L ，其中(1,2,,)i i m λ=L 是A 的特征值. 而按模最小的特征值是m
m
λλ
-%，相应的特征向量与A 的特征向量相同. 利用公式(9.2.30)可求出{}{},k k m z ，并且有
1
lim ,lim max()
m
k k k k m m
m m λλ→∞
→∞
=
=
-%x z x ，
其收敛速度由0
11(min )m m i m i m i m i m
λλλλλλλλ≤≤---=--%%%%确定. 于是当k 充分大时，可取 m k ≈x z ， 1m m
k
m λλ≈+%, 只要近似值m
λ%适当好，收敛速度就很快，往往迭代几次就能得到满意的结果. 例9.2.2 已知特征值λ的近似值0.3589λ
=-%，用原点平移下的反幂法求方阵 得对应特征值λ的特征向量.
解 取0.3589p λ
==-%，对矩阵 1.3589
230.35892 1.3589333 5.3589A pI A I ⎡⎤⎢⎥-=+=⎢⎥
⎢⎥⎣⎦
. 迭代公式(9.2.30)中的k y 是通过解方程组
求得的. 为了节省工作量，可先将p -A I 进行LU 分解.
在LU 分解中尽量避免较小的rr u 当除数，通常可以先对矩阵p -A I 的行进行调
换后再分解. 为此，我们可用001010100⎡⎤
⎢⎥=⎢⎥
⎢⎥⎣⎦
P 乘p -A I 后再进行LU 分解，即 6
001 1.3589
2333 5.3589()0102 1.358932 1.3589310033 5.3589 1.3589231
0033 5.35890.66671000.64110.5726,0.45301100 3.0710p -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-===⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
⎡⎤⎡⎤⎢⎥⎢⎥=--⎢⎥⎢⎥
⎢⎥⎢⎥--⨯⎣⎦⎣⎦
P A I LU 1( 1.2679)k k --=P A I y Pz ， 即 1k k -=LUy Pz .
令11,.k k k k --==Lv Pz Uy v 选取0z ，使1T 10(1,1,1)-==Uy L Pz ，得
T 1(290929.45,290927.56,325732.90)=-y ，
11max()325732.90m ==-y ， T 111(0.8932,0.8931,1)m ==--z y .
由121-=Uy L Pz 得
T 2(845418.49,845418.49,946558.42)=--y ，
22max()946558.42m ==y ， T 222(0.8932,0.8932,1)m ==--z y .
由于1z 与2z 的对应分量几乎相等，故A 的特征值为
211
0.35890.35889894354117946558.42
m λλ≈+
=-+=-%, 相应的特征向量为T 2(0.8932,0.8932,1)==--x z .
而矩阵A
的一个特征值为40.358898943540674λ==L ，相应的特征向量为T (0.89315,0.89315,1)--，由此可见得到的结果具有较高的精度.
9.3 QR 算法
上一节我们介绍了求矩阵特征值的幂法和反幂法. 幂法主要用来求矩阵的按
模最大特征值，而反幂法主要用于求特征向量. 本节将介绍幂法的推广和变形——QR 算法，它是求一般中小型矩阵全部特征值的最有效的方法之一，其基本思想就是利用矩阵的QR 分解. 矩阵n n ⨯∈R A 的QR 分解就是：用Householder 变换将矩阵A 分解成正交矩阵Q 与上三角矩阵R 的乘积，即=A QR . 下面首先介绍Householder 变换.
9.3.1 豪斯霍尔德变换
定义9.3.1 设()ij n n b ⨯=B 是n 阶方阵，若当1i j >+时，0ij b =，则称矩阵B 为上Hessenberg 矩阵(Hessenberg matrix )，又称为准上三角矩阵，它的一般形式为
111212122
2323,1
n n n n n nn b b b b b b b b b b -⎡⎤⎢⎥⎢⎥
⎢⎥=⎢⎥⎢⎥⎢⎥⎣
⎦
L L L L L L O
M B . (9.3.1) 下面讨论如何将矩阵A 用正交相似变换化成(9.3.1)的形式. 为此先介绍一个
对称正交矩阵——Householder 矩阵.
定义9.3.2 设向量n ∈R u 的欧氏长度21=u ，I 为n 阶单位矩阵，则称n 阶方阵
T ()2==-H H u I uu (9.3.2)
为Householder 矩阵(Householder matrix). 对任何n ∈R x ，称由()=H H u 确定的变换
=y Hx (9.3.3)
为镜面反射变换(specular reflection transformation)，或Householder (Householder transformation).
注 Householder 变换，最初由A.C Aitken 在1932年提出. Alston Scott Householder 在1958年指出了这一变换在数值线性代数上的意义. 这一变换将一个向量变换为由一个超平面反射的镜像，是一种线性变换. 运用线性代数的知识，很容易证明：
定理9.3.3 (9.3.2)式定义的矩阵H 是对称正交矩阵；对任何n ∈R x ，由线性变换H =y x 得到y 的欧氏长度满足22=y x .
反之，有下列结论：
定理9.3.4 设,n ∈R x y ，≠x y . 若22=x y ，则一定存在由单位向量确定的镜面反射矩阵()H u ，使得H =x y .
证 令2
-=
-x y
u x y ，显然21=u . 构造单位向量u 确定的镜面反射矩阵 T ()2==-H H u I uu ，
()T T T T
22
22()()2()()
22⎡⎤----⎢⎥=-=-=---⎢⎥⎣⎦
x y x y x y x x y x Hx I uu x I x x x y x y . 又因为22=x y ，即T T =x x y y ，所以
————————————————————
海森伯格(Karl Adolf Hessenberg 1904年9月8日–1959年2月22日) 是德国数学家和工程师.
豪斯霍尔德(Alston Scott Householder 1904年5月5日–1993年7月4日) 是美国数学家，在数学生物学和数值
分析等领域卓有建树.
于是
T T 22
2()()
()--=-
=--=-x y x x y x Hx x x x y y x y
.
证毕！
由定理9.3.4得：
算法9.3.1 若T 12(,,,)n x x x =L x ，其中2,,n x x L 不全为零，则由
12T 112
12T
1T
22
sgn(),,(1,0,,0),(),
()2n x x σσρσσρ-⎧=⎪=+=∈⎪⎪=+⎨⎪⎪==-=-⎪⎩R L 12其中=x u x e e u uu H H u I I uu u (9.3.4) 确定的镜面反射矩阵H ，使得1σ=Hx e ，其中1,0,
sgn()0,0,1,0.a a a a >⎧⎪
==⎨⎪-<⎩
例9.3.1 设T (1,2,2)=--x ，按式(9.3.4)的方法构造镜面反射矩阵H ，使得
T (,0,0)*=Hx (*表示某非零元素).
解
12sgn()(3x σ==-=-x ；
T T T 1(1,2,2)(3,0,0)(4,2,2)σ=-=---=--u x e , 其中T 1(1,0,0)=e ；
2
12()3[3(1)]12x ρσσ==+=--+-=12
u ， 则所求镜面反射矩阵为
1T 100413232310102(4,2,2)2323131200122123ρ----⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥=-=---=⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦
H I uu ， 且
1232313232312023132320---⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦
Hx . 可以证明：
定理9.3.6 对任意n 阶方阵()ij n n a ⨯=A ，存在正交矩阵Q ，使得 为形如式(9.3.1)的上Hessenberg 阵.
证 记
(1)(1)
(1)1112111121(1)(1)
(1)21
222212221(1)(1)(1)121
2n n
n n n n nn n n nn a a a a a a a a a a a a a a a a a a ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥==
=⎢⎥
⎢⎥⎢⎥⎢
⎥⎢⎥⎣⎦⎣⎦L
L
L L M M M M M M L
L A A ，(1)
21(1)311(1)1n a a a ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦
M x . 由式(9.3.4)可构造1n -阶对称正交矩阵1H ：
12121121122T 11
111121111212
1T 1111
sgn()sgn(),
,(1,0,,0),
(),,n i i n a a a a σσρσσρ=--⎧⎛⎫
==-⎪ ⎪⎝⎭⎪
⎪=+=∈⎨⎪=+⎪⎪=-⎩∑R L 12其中=x u x e e u H I u u (9.3.5) 使得1111σ=H x e .
记1
11⎡⎤=⎢
⎥⎣
⎦
I Q H ，且1n n ⨯∈R Q ，1I 是1阶单位矩阵. 显然1Q 是对称正交矩阵. 用1Q 对A 作相似变换，得
(1)
(2)(2)11121(2)(2)1
2221(2)(2)111112323(2)(2)2
0n n n
n nn a a a a a a a a a σ-⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣
⎦
L
L L M M M L
记
Q AQ Q A Q A . (9.3.6)
记(2)(2)(2)T 2
232422(,,,)n n a a a -=∈R
L x ，同理可构造2n -阶对称正交矩阵2H ，使得2221σ=H x e (T 21(1,0,,0)n -=∈R L e ).
记2
22⎡⎤
=⎢
⎥⎣⎦
I Q H ，其中2I 为2阶单位矩阵，则2Q 仍是对称正交矩阵，用2Q 对2A 作相似变换，得
(1)(2)(3)
(3)1112131(2)(3)(3)12223
2(3)(3)1
233
32222223(3)
(3)434(3)
(3)300
00
0n
n n n n nn a a a a a a a a a a a a a σσ-⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
L L L L M M M M L
记
Q A Q Q A Q A . (9.3.7)
依此类推，经过k 步对称正交相似变换，得
(1)(2)(3)
(1)()
()
()
1112131,111,11(2)(3)(1)()
()()12223
2,122,12(3)(1)()()()233
3,133,133(1)()()()1,1
1,1,11,00
0000k k k k k k k n k k k k k k k n k k k k k k k n
k k k k k k k k k k k n
a a a a a a a a a a a a a a a a a a a a a a σσσ--+--+--+-----+-=L L L
L L L
O M M M
M O L ()()()
1
,1()()()1,1,1
1,()
()(),1
0000000
000
k k k k k kk k k kn k k k k k
k k k n k k k nk
n k nn
a a a a a a a a a σ-++++++⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢
⎥⎣⎦
O L L L M M M M M
M
M L
L
记
A . (9.3.8)
重复上述过程，则有
(1)(2)
(3)(1)
()
1112131,11(2)(3)(1)()122
232,12(3)(1)()1
2
333,1
3222
222
13
(1)()1,1
1,()
1
n n n n n n n n n n n n
n n n n n n n n n n n n n
n n nn a a a a a a a a a a a a a a a σσσσ-------------------⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦
L L L
O M
M O
记
Q A Q
Q A Q A . (9.3.9)
由式(9.3.6)-(9.3.9)，可得
122223332231132n n n n n n n n n n n n n -------------===L L A Q A Q Q Q A Q Q Q Q Q AQ Q Q .
若记1n -=B A ，122n -=L Q Q Q Q ，则Q 为正交矩阵，且有T =B Q AQ . 证毕！
注1 由定理9.3.6知：因为任意n 阶方阵A 与n 阶上Hessenberg 矩阵B 相似，所以求A 的阵特征值的问题，就可转化为求上Hessenberg 矩阵B 的特征值的问题.
注2 若A 是对称矩阵，则B 也是对称矩阵. 再由B 的形式(9.3.1)知，此时B
一定是对称三对角阵. 于是，求对称矩阵A 的阵特征值的问题，便可转化为求对称三对角阵B 的特征值问题.
例9.3.2 设矩阵
试用镜面反射变换求正交矩阵Q ，使T Q AQ 为上Hessenberg 矩阵.
解 第1步 记1=A A ，(1)(1)(1)T T 121
3141(,,)(2,1,2)a a a ==x ，利用式(9.3.4)构造三阶镜面反射阵1H ：
112sgn(2)3σ===x ；
T T T 1111(2,1,2)(3,0,0)(5,1,2)σ=+=+=u x e , 其中T 1(1,0,0)=e ；
2
(1)
1111212()3(32)15a ρσσ==+=+=12
u ； 则所求镜面反射矩阵为
T 2111111
3003 2.33330.46670.066700.4667 1.5733 1.346700.0667 1.34670.0933-⎡⎤⎢⎥
--⎢
⎥===⎢⎥⎢⎥-⎣⎦
A Q A Q Q AQ . 第2步 记(2)(2)T
T 232
42(,)(0.4667,0.0667)a a ==-x , 利用式(9.3.4)构造2阶镜面反射阵2H ：
222sgn(0.4667)0.4714σ===x ；
T T T 2221(0.4667,0.0667)(0.4714,0)(0.9381,0.0667)σ=+=-+=-u x e , 其中
T 1(1,0)=e ；
2
(2)
2222322
()0.4714(0.47140.4667)0.4422a ρσσ==+=+=12u ； 则所求镜面反射矩阵为
T 32222221
3003 2.33330.4714000.4714 1.5733 1.500000 1.50000.5000-⎡⎤⎢⎥--⎢
⎥===⎢⎥--⎢⎥-⎣⎦
A Q A Q Q A Q . 由此得正交矩阵
121
00000.66670.23570.707100.33330.94290.000100.66670.23570.7070⎡⎤⎢⎥
--⎢
⎥==⎢⎥--⎢⎥-⎣⎦
Q Q Q ， 使得
为上Hessenberg 矩阵。

9.3.2 QR 算法
QR 算法的基本思想是：利用QR 分解得到一系列与A 相似的矩阵{}k A ，在一定
的条件下，当k →∞时，{}k A 收敛到一个以A 的特征值(1,2,,)i i n λ=L 为主对角线元素的上三角矩阵. 首先介绍用QR 分解 (QR decomposition 或 QR factorization); 即用Householder 变换将矩阵A 分解成正交矩阵Q 与上三角矩阵R 的乘积，即=A QR .
9.3.2.1 QR 分解 算法9.3.2 (QR 分解)
第1步 记A 的第1列为(1)(1)(1)T 111211(,,,)n a a a =L x ，(1)1()ij n n a ⨯==A A . 利用式
(9.3.4)：
构造出的1H 是n 阶对称正交矩阵，使得1111σ=H x e ，从而有
第2步 记(2)(2)(2)T 222322(,,,)n a a a =L x ，同理可构造出1n -阶对称正交矩阵2%H ，使
得
2222σ=%H x e ，其中1(2)222222sgn()n i i a a σ=⎛⎫= ⎪⎝⎭
∑，T 12(1,0,,0)n -=∈R L e . 若记221
,⎛⎫
=
⎪⎝
⎭
%H H 它仍是对称正交矩阵，于是有 如此继续下去，直到完成第1n -步后，得到上三角矩阵 于是有
令n =R A ，121n -=L Q H H H ，其中Q 是对称正交矩阵，则=R QA . 因为Q 是对称正
交矩阵，所以得
=A QR .
注1 若A 非奇异，则上三角矩阵R 也非奇异，从而R 的主对角线元素不为零.
注2 若要求R 的主对角线元素取正数，则A 的QR 分解是唯一的. 例9.3.3 求矩阵
的QR 分解=A QR ，并使矩阵R 的主对角线上的元素都是正数。

解 对A 运用算法9.3.2.
第1步 记1=A A , T 1(1,2,2)=-x ，则
13σ==, 11111()3(31)12a ρσσ+=+==,
T T T T 11111(1,2,2)(3,0,0)(4,2,2),(1,0,0)σ=+=-+=-=u x e e ,
1T 11111004122110102(4,2,2)2211230012212ρ---⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥=-=--=-⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦
H I u u , 2113473053103072--⎡⎤
⎢⎥==⎢⎥
⎢⎥⎣⎦
A H A . 第2步 记T 2(53,73)=x
，2sign(5 2.86744σ==,
(2)
22222() 2.86744(2.8674453)13.0013a ρσσ+=+==,
T T T 2221(573)(2.86744,0)(4.53411,2.33333)σ=+=+=u x e ,
0.581240.813730.813730.58124--⎡⎤
=⎢⎥
-⎣⎦
, 记
221001000.581240.81373000.813730.58124H ⎡⎤
⎡⎤⎢⎥==--⎢⎥⎢⎥⎣⎦⎢⎥
-⎣⎦
%H ，
则
3223 1.33333 2.333330 2.86744 2.4799500 2.32494--⎡⎤
⎢⎥==--⎢⎥
⎢⎥-⎣⎦
A H A . 为了使R 的主对角线上的元素都是正数，取3100010001-⎡⎤
⎢⎥=-⎢⎥
⎢⎥-⎣⎦
H ，显然3H 是正交阵，且
4333 1.33333 2.333330 2.86744 2.4799500 2.32494-⎡⎤⎢⎥==⎢⎥
⎢⎥⎣⎦
A H A . 令
43 1.33333 2.333330 2.86744 2.4799500 2.32494-⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦R A ，1230.333330.154990.929980.666670.658740.348740.666670.736230.11625-⎡⎤⎢⎥==⎢⎥
⎢⎥--⎣⎦
Q H H H ， 且=A QR .
9.3.2.2 QR 算法
了解了QR 分解后，下面介绍QR 算法(QR algorithm ). 算法9.3.3 (QR 算法)
第1步 令1=A A ，利用算法9.3.2将1A 进行QR 分解，得111=A Q R ，其中1Q 为
正交矩阵，1R 为上三角矩阵；然后将1Q 与1R 逆序相乘，得
211=A R Q .
因为1111-=R Q A ，所以有1211111-==A R Q Q A Q ，即2A 与1A 相似.
第2步 以2A 代替1A ，再作QR 分解，得222=A Q R ，其中2Q 为正交矩阵，2
R 为上三角矩阵；再将2Q 与2R 逆序相乘，并记
322=A R Q .
因为1222-=R Q A ，所以1322222-==A R Q Q A Q ，即3A 与2A 相似.
依此类推，可得QR 算法公式：对1,2,k =L ，
111,
,
k k k k k k k k +++=⎧⎨
==⎩A Q R A R Q Q R () 因为1111k k k ----=R Q A , 所以111111k k k k k k ------==A R Q Q A Q ，即k A 与1k -A 相似. 故序列{}
k A 相似于1=A A .
这里，我们不加证明地给出QR 算法收敛的充分条件：
定理9.3.9 (QR 算法的收敛性) 设()n n ij a ⨯=∈R A , {}k A 是由QR 算法产生的矩
阵序列，其中()
()k k ij a =A . 若
(1) 1=A A 的特征值(1,2,,)i i n λ=L 满足
120n λλλ>>>>L ；
(2) 1-=A P DP , 其中12diag(,,,)n λλλ=L D ，且P 有三角分解=P LU (L 是单位下三角矩阵，U 是上三角矩阵)，则
()
0,,
lim ,,k ij
k i
i j a
i j λ→∞
>⎧=⎨
=⎩ () 即{}k A 收敛到一个以A 的特征值(1,2,,)i i n λ=L 为主对角线元素的上三角矩阵. 推论 若矩阵n n ⨯∈R A 是对称矩阵，且满足定理9.3.9中的条件，则由QR 算法产生的矩阵序列{}k A 收敛到对角矩阵12diag(,,,)n λλλ=L D .
9.3.3 带原点位移的QR 算法
QR 算法收敛速度是线性的，需要提高收敛速度. 经分析知：定理9.3.9中
()
lim k nn n k a λ→∞
= 的速度依赖于比值11n n n γλλ-=<，当n γ越小，收敛速度越快. 这启发我们把9.2节介绍的原点平移法加速幂法和反幂法的思想运用到QR 算法的收敛加速，这就是下面将要介绍的带原点位移的QR 算法(QR algorithm with shift). 因为求上Hessenberg 矩阵的特征值比一般的方阵简单，所以在下面的算法中首先将所考察的矩阵化成上Hessenberg 矩阵，然后再用位移加速方法进行加速.
算法9.3.4 (带原点位移的QR 算法)
第1步 将矩阵A 化成上Hessenberg 矩阵1A .
第2步 对1,2,k =L ，给定原点位移数列{}k s ，进行迭代加速
1
11QR (),
,k k k k k k k k k k s s +++-=⎧⎨
==+⎩分解A I Q R A R Q Q R I () 并称由k A 到1k +A 的变换为带原点位移的QR 变换 (QR transformation with
shift)，且记()
()k n n k ij a ⨯=∈R A . 其中{}k s 的选取方法主要有
(a) 取()
k k nn s a =；
或
(b) 取k s 为k A 的右下角主子矩阵
的2个特征值中最靠近()
k nn a 的那一个.
注1 由算法9.3.4产生的1k +A 与k A 相似，即1T 1k k k k k k k -+==A Q A Q Q A Q ；且它们都是上Hessenberg 矩阵.
注2 计算实践表明，取方法(b)中k s 的算法比取方法(a)中k s 的算法收敛速度更快. 特别是对称矩阵，带方法(b)中k s 位移的QR 算法是无条件收敛的，且收敛阶为3.
注3 在迭代过程中，当()
,10k n n a -≈时，
由定理9.3.4知：()k nn n a λ≈，因此可取()k nn a 作为A 的近似特征值. 而()lim k nn n k a λ→∞
=的速度依赖于比值()
1()()k n n k n k s s γλλ-=--，()
k n γ越小，收敛速度越快；因此平移值取()
k k nn s a =显然是一个很好选择.
注4 在迭代过程中，当()
1,20k n n a --≈时，取
的2个特征值作为A 的特征值1,n n λλ-的近似值.
注5 判别(),1k n n a -和()1,2k n n a --约等于0的标准可取为
其中0ε>是预先给定的精度.
注6 当求得矩阵A 的特征值n λ或1,n n λλ-后，可以将k A 作降阶处理：对k A 左上角的1n -阶或2n -阶主子矩阵继续运用带原点位移的QR 算法, 以求A 的其他特征值，这样可以大大节省运算量.
9.4 Jacobi 方法
上一节我们介绍了Householder 变换将矩阵化为上Hessenberg 矩阵的方法. 如果矩阵是实对称矩阵，用平面旋转变换将矩阵化为上Hessenberg 矩阵比用
Householder 变换更好. 本节介绍的Jacobi 方法就是这种方法. 它主要用于求实对称矩阵的全部特征值和特征向量.
Jacobi 方法(Jacobi method)的基本思想是：对实对称矩阵()ij n n a ⨯=A ，一定存在正交矩阵Q ，使
1T -==Q AQ Q AQ D , (9.4.1)
其中12diag(,,,)n λλλ=L D ，(1,2,,)j j n λ=L 就是矩阵A 的特征值，而正交矩阵Q 的第j 列就是对应于j λ的特征向量.
由此可见，Jacobi 方法的实质和关键就是找一个正交矩阵Q ，将A 化为对角矩阵.
9.4.1 Jacobi 方法
定义9.4.1 设()ij n n a ⨯=A 是n 阶实对称矩阵，称n 阶矩阵
()
()
1cos sin ()1
(,,)1sin cos ()1i j i i j j θ
θθθθ
⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥⎢
⎥=⎢⎥⎢⎥
⎢⎥-⎢⎥⎢⎥⎢⎥⎣
⎦
O L L L
M M M O
M M
M L
L
L
O
G (9.4.2) 为旋转矩阵(rotation matrix)，或Givens 矩阵(Givens matrix)，简记为ij G . 对
A 进行的变换
T
ij ij
G AG (9.4.3) 称为Givens 旋转变换(Givens rotation).
注1 Givens 矩阵是在n 阶单位矩阵I 的第i 行第i 列、第i 行第j 列、第j 行第
i 列、第j 行第j 列的交叉的位置上分别换上cos ii r θ=、sin ij r θ=、sin ji r θ=-、
cos jj r θ=而
成的.
注2 Givens 矩阵是正交矩阵，变换(9.4.3)是正交相似变换.
注3 Jacobi 方法就是通过一系列Givens 旋转变换，把A 化为对角矩阵，从而求得特征值及相应的 特征向量的方法. 因此，Jacobi 方法也称为平面旋转法(plane rotation method).
下面介绍具体介绍将n 阶实对称矩阵化为对角矩阵的Jacobi 方法. 设()ij n n a ⨯=A 是n 阶实对称矩阵，记
(1)T
1()ij n n ij ij
a ⨯==A G AG . (9.4.4) 因为
()T
T T T
11ij ij ij ij
===A G AG G AG A , 所以1A 仍是对称矩阵. 通过直接计算可得 ——————————————————————
基文斯(J. W. Givens 1910 年12月 14日–1993 年3月 5日) 是美国数学家，计算机领域的先驱之一.
(1)
22(1)22(1)(1)(1)(1)
(1)(1)(1)(1)12cos sin 2cos sin ,sin cos 2cos sin ,
cos sin ,
,,sin cos ,
,,,
,,,
()ii ii jj ij jj ii jj ij il li il jl jl lj il jl lm ml ml ij
ji jj
ii a a a a a a a a a a a a l i j a a
a a l i j a a a m l i j a a a a θθθθθθθθθθθθ=++=+-==+≠==-+≠==≠==-22sin 2(cos sin ).ij a θθθ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪+-⎩ (9.4.5)
不难看出，A 经过ij G 作用后，与A 相比，只有1A 的第i 行、第i 列、第j 行、第j 列的元素发生了变化，而其他元素与A 的相同.
特别地，若0ij a ≠，由(9.4.5)式中最后的式子知：若取θ满足关系式
21tan cot 2,,22tan 44
ii jj ij
a a a θππ
θθθ--=
=-<≤ (9.4.6) 可使(1)
(1)0ij ji a a ==，也就是说，用ij G 对A 进行变换，可将A 的2个非主对角线元素ij a 和ji a 化为零.
Jacobi 方法的一般过程是：记0=A A ，选择A 的一对最大的非零非主对角线元素ij a 和ji a ，使用Givens 矩阵ij G 对A 作正交相似变换得1A ，可使1A 的这对非
零非主对角线元素(1)
(1)0ij ji a a ==.
再选择1A 的一对最大的非零非主对角线元素作上述正交相似变换得2A ，可使2A 的这对非零非主对角线元素化为零.
如此不断地做下去，可产生一个矩阵序列01,,,,k =L L A A A A .
注 虽然A 至多只有(1)2n n -对非零非主对角线元素，但是不能期望通过(1)2n n -次变换使A 对角化. 因为每次变换能使一对非零非主对角线元素化为零，例如，ij a 和ji a 化为零. 但在下一次变换时，它们又可能由零变为非零.
不过可以证明，如此产生的矩阵序列01,,,,k =L L A A A A 将趋向于对角矩阵，即Jacobi 方法是收敛的. 而这个对角矩阵的主对角线元素就是矩阵A 的特征值.
用Jacobi 方法求矩阵A 的特征值的步骤为：
(1) 记0=A A ，在矩阵A 中找出按模最大的非主对角线元素ij a ，取相应的Givens 矩阵ij G ，记为1ij =G G ；
(2) 由条件22()sin 22(cos sin )0jj ii ij a a a θθθ-+-=，定出sin θ, cos θ. 为避免使用三角函数，令
21,2tan sgn(),cos (1),sin cos .
ii jj ij
a a d a t d d t t θθθθ--⎧=⎪
⎪⎪
==+
⎨⎪
⎪=+⎪
=⎩ (9.4.7)
(3) 按公式(9.4.5)计算T T
111
ij ij ==A G AG G AG 的元素； (4) 以1A 代替0A ，重复步骤(1)， (2) ，(3)，求出T 2212=A G A G ；依此类推，得
T 1,
1,2,3,k k k k k -==L A G A G . (9.4.8) 令0=Q I ，记T 1k k k -=Q Q G ，则k Q 是正交矩阵，且
T ,
1,2,3,k k k k ==L A Q AQ . (9.4.9) 若经过N 步旋转变换，N A 的所有非主对角线元素都小于允许误差ε时，停止计算. 此时N A 的主对角线元素就是A 的特征值的近似值. N Q 的列元素就是A 的对应于上述特征值的全部特征向量.
9.4.2 Jacobi 方法的收敛性
由矩阵理论知：
定理9.4.2 设()ij n n a ⨯=A ，P 是正交矩阵，T ()ij n n b ⨯==B P AP ，则
2
211
11
n n n n
ij
ij
i j i j b a
=====∑∑∑∑. (9.4.10)
定理9.4.3(收敛性) 设{}k A 是由Jacobi 方法产生的矩阵序列，其中
T k k k =A Q AQ 由式
12lim diag(,,,),lim k n k x x λλλ→∞
→∞
===L A D Q Q , (9.4.11)
其中(1,2,,)j j n λ=L 就是矩阵A 的特征值，而正交矩阵Q 的第j 列就是对应于j λ的特征向量，1,2,,j n =L .
分析：要证明Jacobi 方法的收敛性，只要能证明每次变换，总是使主对角线元素的平方和增大，而非主对角线元素的平方和减小即可.
证明 由定理9.4.2知：在正交相似变换下，矩阵元素的平方和不变. 设矩阵
A 的非主对角线元素平方和为
1
2,1()n ij
i j i j E a
-=≠=
∑A ， (9.4.12)
因此，要证明k A 收敛于对角矩阵D ，只要证明()E A 收敛于零即可. 假设k A 变为
1k +A 时，把k A 的绝对值最大的非主对角线元素()
k ij a 和()k ji a 化为零，由式(9.4.5)直
接计算，得
()()()()()()222
222(1)(1)()()()()()()cos sin sin cos k k k k k k k k il
jl
il
jl
il
jl
il jl
a a a
a a a
a
a θθθθ+++=++-+=+,
且有(1)
(1)0k k ij ji a a ++==，所以
()2
()
1()()2.k k k ij E E a +=-A A (9.4.13)
这表明经过变换后，1k +A 的非主对角线元素的平方和减小了()
2
()2k ij a ，由定理9.4.2
知：k A 与1k +A 的元素的平方和不变，故1k +A 的主对角线元素的平方和应增加
()2
()
2k ij a .
因为()
()k k ij ji a a =是k A 的绝对值最大的非主对角线元素化为零，所以
()
1
22(),1
()(1)n k k st
ij s t s t
E a
n n a -=≠=
≤-∑A ，
即
()2
()1
().(1)
k ij k a E n n ≥
-A (9.4.14) 于是
因为2
011(1)
n n <-
<-，所以当k →∞时，1()0k E +→A . 故 12lim diag(,,,)k n x λλλ→∞
==L A D . (9.4.15)
因为相似矩阵的特征值相同，所以(1,2,,)j j n λ=L 就是矩阵A 的特征值. 由式(9.4.9)和式(9.4.15)得
T T
12diag(,,,)lim n k k k λλλ→∞
===L D Q AQ Q AQ .
因为Q 是正交矩阵，所以T 1-=Q Q ，若记(1)(2)()(,,,)n Q Q Q =L Q ，其中
()(1,2,,)
j Q j n =L 是Q 的第j 列，则有=AQ QD ，即
(1)(2)()(1)(2)()(,,,)(,,,)n n Q Q Q Q Q Q =L L A D ，得
(1)(2)()(1)(2)()12(,,,)(,,,)n n n Q Q Q Q Q Q λλλ=L L A A A ，
亦即
()(),
1,2,,j j j Q Q j n λ==L A .
这说明正交矩阵Q 的第j 列就是对应于j λ的特征向量，1,2,,j n =L . 证毕！
例9.4.1 利用Jacobi 方法求矩阵
的全部特征值和特征向量，要求k A 的所有非主对角线元素的绝对值0.1ε<=.
解 记(0)
033()ij a ⨯==A A ，因为(0)122a =-是A 中所有非主对角线元素中绝对值最大的元素，取相应的Givens 矩阵12G . 因为(0)(0)
11
221,1a a ==-，所以 (0)(0)(0)
112212()20.5d a a a =-=-，
(tan sgn()
0.618034t d d θ===-，
于是
121cos sin 00.8506510.5257310sin cos 00.5257310.85065100
1001θ
θθ
θ-⎡⎤⎡⎤
⎢⎥⎢⎥=-=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
记
G G ,
T (1)
1011332.236068
00.5257310 2.2360680.850651()0.5257310.8506513ij a ⨯-⎡⎤⎢⎥=-=⎢⎥
⎢⎥-⎣⎦
记
G A G A .
记
则T 111=A Q AQ .
用1A 代替0A ，重复上述过程： 因为(1)
23
0.850651a =是1A 中所有非主对角线元素中绝对值最大的元素，取相应的Givens 矩阵23G .因为(1)(1)
22332.236068,3a a =-=，所以(1)(1)(1)
22
3323()2 3.077683d a a a =-=-
；(
tan sgn()0.158384t d d θ==+=-，
于是
T
(2)
2122332.2360680.0822410.5192580.082241 2.3707980()0.5192580 3.134730ij a ⨯-⎡⎤
⎢⎥=-=⎢⎥
⎢⎥-⎣⎦
记
G AG A ,
记
则T 222=A Q AQ .
用2A 代替1A ，重复上述过程： 因为(2)
13
0.519258a =-是2A 中所有非主对角线元素中绝对值最大的元素，取相应的Givens 矩阵13G .因为
(2)(2)
11332.236068, 3.134730a a ==，所以
(2)(2)(2)113313()20.865333d a a a =-=
；tan sgn()
0.457089t d d θ===，
于是
T
(3)
3233331.9987210.07479900.074799 2.3707880.034190()00.034190 3.372078ij a ⨯⎡⎤
⎢⎥=--=⎢⎥
⎢⎥-⎣⎦
记
G A G A .
记
则T 333=A Q AQ .
因为3A 的所有非主对角线元素的绝对值0.1ε<=，所以迭代停止. 此时A 的
特征值的近似值分别为3A 的对角元素：
1231.998721, 2.370788, 3.37208λλλ≈≈-≈.
相应的特征向量的近似值分别为
其中1230.410602,0.156434,0.898295k k k ==-=.
A 的特征值的精确值为
123112, 2.37228132, 3.372281322
222
λλλ==-
=-=+=L L .
相应的特征向量为
从例9.4.1可以看出，即使迭代的次数不是很多，迭代矩阵k A 的所有非对角元素的绝对值并不是很小时，用Jacobi 方法求得的结果精度都比较高，因此它是求实对称矩阵的全部特征值和特征向量的一个较好的方法.
9.4.3 改进的Jacobi 方法
由于每次旋转变换前选非零非主对角线元素的最大值很费时间，为此介绍如下改进方法.
第一种方法：把非主对角线元素按照行的次序12131232421,,,,,,,,,,n n n n a a a a a a a -L L L 依次化为零，称为一次扫描. 一次扫描后，前面已化为零的元素可能成为非零元素，需要再次扫描. 这一方法称为循环Jacobi 方法，这种方法的缺点是: 一些已经足够小的元素也作化零处理，浪费了时间.
第二种方法：首先对实对称矩阵A 计算
12
120112n n ij i j i v a -==+⎛⎫
= ⎪⎝⎭
∑∑， (9.4.16)
设置阀值10v v n =，按12131232421,,,,,,,,,,n n n n a a a a a a a -L L L 的顺序进行扫描.
若1ij a v ≥，则选取旋转矩阵ij G 作旋转相似变换将ij a 和ji a 化为零；否则让ij a 过关，即不进行旋转相似变换将其化为零.
因为某些绝对值小于1v 的元素的绝对值可能在后面的旋转变换中增长，所以应进行多次扫描，直到1A 的所有非零非主对角线元素的绝对值都小于1v .
再设置阀值21v v n =，重复上述过程，直到达到精度要求，即k v ε<为止(其中0ε>是指定的精度). 这种方法称为限值Jacobi 方法.
9.5 算法程序
9.5.1 幂法
%幂法
% A 是方阵，e 是精度，输出向量V 是A 的最大特征值MaxEig 对应的特征向量 function PowerMethod(A, e) V = ones(size(A,1), 1); for i=1:2 Y = A * V;
M(i) = norm(Y, Inf); V = Y / M(i); end
while abs(M(i)-M(i-1)) > e i = i + 1; Y = A * V;
M(i) = norm(Y, Inf); V = Y / M(i); end
MaxEig = M(i);
disp(sprintf('最大特征值MaxEig %.6f\n相应的特征向量', MaxEig))
V=V’,
end
例9.5.1用幂法求方阵
解在MATLAB命令窗口中输入：
A=[1 2 3;2 1 3; 3 3 5]; e=10^-3; PowerMethod(A,e) 回车，输出结果：
最大特征值MaxEig =8.358906
相应的特征向量
V =
0.5598 0.5598 1.0000
9.5.2 幂法的加速
%幂法的加速
%A是方阵，e是精度，输出向量V是A的最大特征值MaxEig对应的特征向量function Acc_ PowerMethod (A, e)
V = ones(size(A,1), 1);
for i=1:2
Y = A * V;
M(i) = norm(Y, Inf);
V = Y / M(i);
end
for i=3:4
Y = A * V;
M(i) = norm(Y, Inf);
M2(i) = M(i-2) - (M(i-1) - M(i-2))^2 / ( M(i)-2*M(i-1)+M(i-2) ); V=Y/M2(i);
end
while abs(M2(i)-M2(i-1)) > e
i = i + 1;
Y = A * V;
M(i) = norm(Y, Inf);
M2(i) = M(i-2) - (M(i-1) - M(i-2))^2 / ( M(i)-2*M(i-1)+M(i-2) ); V = Y / M2(i);
end
MaxEig = M2(i);
disp(sprintf('最大特征值MaxEig %.6f\n相应的特征向量', MaxEig))
V,
end
9.5.3 反幂法
%反幂法
% A是方阵，e是精度，输出向量V是A的最大特征值MaxEig对应的特征向量。