高阶对称矩阵特征值的计算毕业论文
高阶对称矩阵特征值的计
算毕业论文
目录
摘要.................................... 错误！未定义书签。Abstract ................................ 错误！未定义书签。目录.. (1)
引言.................................... 错误！未定义书签。
1 关于矩阵特征值的概念 (2)
1.1 矩阵特征值和特征向量的定义 (2)
1.2矩阵特征值的计算方法 (2)
1.3对称矩阵特征值的一些性质 (3)
2高阶对称矩阵特征值的计算方法 (4)
2.1Jacobi旋转法 (4)
2.1.1Jacobi旋转法的概念 (5)
2.2幂法 (7)
2.2.1幂法的概念 (7)
2.2.2反幂法的概念 (10)
2.3 QR方法 (12)
2.3.1豪斯霍尔德变换 (12)
高阶对称矩阵特征值的计算
2.3.2用正交相似变换约化一般矩阵为上海森伯格矩阵 ......... 12 2.3.3豪斯霍尔德约化对称矩阵为对称三对角矩阵 ............. 15 2.3.4QR 方法的算法及实例 ................................. 15 3结束语 ................................................ 17 参考文献 (19)
1 关于矩阵特征值的概念
本论文在第一部分首先介绍矩阵特征值的定义，接着介绍本文讨论的关于对称矩阵在特征值问题上的不同于一般矩阵的性质。这样在下面介绍其他矩阵特征值计算方法之前，可以对特征值有一个大致的了解。
1.1 矩阵特征值和特征向量的定义
定义1[3] 设矩阵()n n ij A a ?=∈ ，若存在实数或者复数λ和非零向量
12(,,,)T n n x x x x =∈ ，使
Ax x λ= （1） 则称λ为A 的特征值，x 为A 对应λ的特征向量。
1.2矩阵特征值的计算方法
上面说明了矩阵特征值的定义，下面我们来看矩阵特征值的计算方法。 由（1）式可知λ可使其次线性方程组
()0I A x λ-=
有非零解，故系数行列式det()0I A λ-=，记
（2）
()p λ称为矩阵A 的特征多项式，方程（2）称为矩阵A 的特征方程。因为n 次代
数方程()p λ在有n 个根12,,,n λλλ （存在于复数域中），故
12()()()()n p λλλλλλλ=---
把（2）式的行列式进行展开可以得到
12(1)(1)det n n n n c A λλλ=-=-
故矩阵()n n ij A a ?=∈ 的n 个特征值12,,,n λλλ 是它的特征方程（2）的n 个根[4]。并有
12det n A λλλ=
例1 求
的特征值
解 A 的特征方程为
故A 的特征值为122λλ==，37λ=。
1.3对称矩阵特征值的一些性质
11
12121
2221
2111()det()0n n
n n nn
n n n n a a a a a a p I A a a a c c c λλλλλλλλ--------=-=
---=++++= 1121n
ii
n i c a λλλ=-=+++=∑ 122224242A -?? ?
=-- ? ?-??32
21
22det()22432428242(2)(7)0
I A λλλλλλλλλ--?? ?-=+-=+-+ ?
?--+??=-+=
高阶对称矩阵特征值的计算
上面我们给出了矩阵特征值的定义及其计算方法，下面我们来看一看对称矩阵的一些独特的关于特征值的性质
定理1 设n n A ?∈ 为对称矩阵，则 （1）A 的特征值均为实数 （2）A 有n 个线性无关的特征向量 （3）存在一个正交矩阵P 使
且(1,2,,)i i n λ= 是A 的特征值，而12(,,,)n P u u u = ，它的列向量j u 为A 和j λ的特征向量是对应的。
定理2 设n n A ?∈ 为对称矩阵。其特征值按顺序写为12n λλλ≥≥≥ ，则
（1） （对任何非零向量n
x ∈ ）
（2） ，
2高阶对称矩阵特征值的计算方法 2.1Jacobi 旋转法
上文中提到的初等矩阵分解算法在遇到高阶矩阵的时候就显得不那么便于计算，而Jacobi 旋转法很好的解决了这一问题。Jacobi 通过他的著作《函数行列式》让矩阵计算进入新的纪元[5]Jacobi 旋转法可以应用于本文所讨论的高阶对称矩阵特征值计算中。Jacobi 算法主要的目的就是将对称矩阵A 化为对角矩阵，然后进行矩阵的特征值的求解。而要实现这种变换就必须使用旋转变换或者说是正交相似变换来达到目的[6]。另外，Jacobi 旋转法的计算步骤较多量，人为计算通常比较繁琐，这是因为其核心方法是迭代计算，所以一般使用计算机进行计算。但是它的优点是拥有较高的精确度，所以在计
12T n P AP λλλ?? ?
?= ? ?
??
1(,)(,)
n Ax x x x λλ≤≤10(,)
max
(,)
n x x Ax x x x λ∈≠= 0(,)min (,)n n x x Ax x x x λ∈≠=
算矩阵特征值时，它也是一种常用方法。 2.1.1Jacobi 旋转法的概念
令0A A =,依次构造矩阵序列k A ，使1(,)(,)(1,2)T
k k k k A R p q A R p q k -== ，其中
(,)R p q 是在(,)p q 平面上转过角度θ的一个旋转。其中 (,)1k R i i =， ,i p q ≠
(,)cos k R i i θ= ,i p q = (,)(,)sin k k R i j R j i θ=-=，i p =，j q =， (,)0R i j = ,i j =其他。
如何确定上面所提到的平面和旋转角θ就是我们下面需要考虑的问题了，这里我们可以观察S A 位于主对角线以上的全部元素，通过它来确定最大模的项()s pq a （这里由于对称性的原因，使得我们仅仅只要在A 的上部的元素中来求就可以了）[7]。剩下的问题就是确定旋转角θ，使得10s pq a +=，我们知道平面旋转R(p ，q)有只会对一部分行和列产生影响的相似变换。即第p ，q 行和列。p ，q ，θ由以下原则确定：
记 (1()1(),()k k k ij k ij A a A a --==，则有()(1)(,,)k k ij ij a a i j p q -=≠
()()(1)(1)cos sin k k k k pi ip ip iq a a a a θθ--==- (,)i p q ≠ ()(1)(1)()sin cos k k k k ip ip iq qi a a a a θθ--=-= (,)i p q ≠
()(1)2(1)(1)2cos 2cos sin sin k k k k pp pp pp pp a a a a θθθθ---=-+， ()(1)2(1)(1)2sin 2cos sin cos k k k k pp pp pp pp a a a a θθθθ---=-+， ()()(1)(1)(1)22()cos sin (cos sin )k k k k k pp qp pp pp pp a a a a a θθθθ---==-+-
为了达成
()0K pq
a =的条件，即消去(,)p q 位的元素，对于旋转角，必须满足下面的条件
（1） ；
（2）选择角θ满足： ；此外，若(1)(1)0k k pp qq a a ---=
那么选择 注意(1)0k pq a ->时 ，(1)0k pq a -<。 通过上述的方法得到的旋转角让它对应的旋转矩阵1T
k k k k A R A R -=中的
()0k pq a =，但是即使出现(1)0k pj a -=或(1)0k qj a -=的情况，在通过一系列的旋转变换后它们一般不会等于0，所以我们说这种对对称矩阵的非对角元素的旋转变换的处理方法没有结束的时候，这是因为使用了迭代的方法，所以我
(1)(1)(1)1
()sin 2cos 22k k k pp qq qq a a a θθ
---=-+(1)(1)(1)
22k qp k k pp qq a tg a a θ----=-44
ππ
θ-≤≤
(1)(1)||4
k qp k qp a a πθ--==4πθ=-,1,2
()i j n
ij
i j
E A a
=≠=
∑
高阶对称矩阵特征值的计算
们在进行这一过程之前都要先有指定的精度，方便我们的计算。Jacobi 方法的收敛性界定条件如下：矩阵中对角之外的元的平方和 在每一次正交相似变换后减少。
我们最后会算得一个矩阵序列。它是一个对角形矩阵。这个矩阵序列是确定的，我们观察这个对角矩阵，它和初始矩阵几乎没有不同的地方，同时，这也是稳定的一个过程。最终得到一个具有一定精度的对角矩阵，并有
T RAR D =，其中T R 各列就是矩阵A 的特征向量。同时其对稀疏带状矩阵的
平面旋转变换会将原矩阵的稀疏带状破坏。 例2 求矩阵特征值，要求使用Jacobi 方法
的特征值和特征向量。 解 首先取1,3i j ==由
得
sin 0.45969,cos 0.88808θθ≈≈
所以
不停的重复下去，当非对角元素趋近零的时候，九次变换之后，得到
210121012A -??
?=-- ?
?-?
?
tan 2θ=
=20.8880800.459690100.4596900.88808P -?? ?= ?
?-??
2220.633980.3250500.3250530.6279700.62797 2.36603T A P AP -??
?
==-- ?
?-??
9A 处在对角线上的元素就是我们要的结果，那么特征值就是
1230.58758, 2.00000, 3.41421λλλ≈≈≈
相应的特征向量为
相对应的特征值的精确值12322,2λλλ=== 相对应的特征向量为
2.2幂法
迭代算法中另外一个算法就是幂法。幂法有其不同于前一种算法的特殊用途，那就是算矩阵的主特征值，而主特征值就是矩阵按模最大的特征值[8]。该方法在大型稀疏矩阵中应用较多[9]。乘幂法也有一些缺点，譬如收敛速度慢，尤其是在出现了两个以上连续的绝对值相近的特征值的情况下。 2.2.1幂法的概念
设实矩阵()ij n n A a ?=有一个特征向量组。且这个向量组是完备的。矩阵A 有n 个线性无关的特征向量，其特征值为12,,,n λλλ ，对应的特征向量为
12,,,n x x x ，已知A 的按模最大的特征值是实根，并且满足
123n λλλλ>≥≥≥
90.587580.000000.000000.00000 2.000000.000000.000000.00000 3.41421A ?? ?= ?
???
1230.500000.707100.500000.70710,0.00000,0.707100.500000.707100.50000v v v ??????
? ? ?===- ? ? ?
? ? ?-?????
?
1231122,0,1122v v v ?
??? ? ? ? ?
=== ? ? ? ? ? ??
???
高阶对称矩阵特征值的计算
下面说明求1λ和1x 的方法。
幂法的基本思路是通过任意取一个非零初始向量0v ，由矩阵A 来构造一个向量序列
叫做迭代向量。由上面的假设，0v 可以写成
01122n n v x x x ααα=+++ （10α≠）
于是
其中12
(/)n k k i i i i x εαλλ==∑。由假设1/1(2,3,,)i i n λλ<= ，故， 从而
这就说明序列 逐渐的接近A 的对应1λ的特征向量，
也能说当k 充分大的时候
说明迭代向量k v 可以看成是1λ的特征向量，因为它们是相似的。这里还有个附加条件就是除一个因子外。
接着阐述如何计算主特征值1λ，假设()k i v 是k v 的第i 个分量，则
10
221011
0,,
k k k v Av v Av A v v Av A v ++=??==????==???
101112221
1111112[(/)](),
k k k k k k n n n
n
k
k i i i k i v Av A v x x x x x x αλαλαλλααλλλαε-====+++=+≡+∑ 11
1lim
k
k
k v x αλ→∞
=1k
k v λ1111111()()()()()()k i
i k i k i i
k i v x v x αελαε++??+=??+??111
k k v x αλ≈lim 0
k k ε→∞
=
故
这说明了两个相邻的迭代向量分量它们的比值在主特征值处收敛。 我们可以总结出来幂法先要构造一个向量序列{}k v ，为此就需要通过非零向量0v 和矩阵A 的乘幂k A ，完成之后就可以计算A 的主特征值1λ和该特征值所对应的特征向量了。
例3 用幂法计算
的主特征值和对应的特征向量
表 1 计算结果
11()lim
,
()k i
k k i
v v λ+→∞=1.0 1.00.51.0 1.00.250.50.25 2.0A ??
?= ?
???
高阶对称矩阵特征值的计算
下述结果有个限定条件，这里限定为8位浮点数字，k u 的分量值是四舍五入得到的。那么就可以得到
2.5365323λ≈
及相应的特征向量10.74820.64971.λT
（，，）和对应的特征向量的真值（8位数字）为
1 2.5365258λ=
1(0.74822116,0.64966116,1)T x
= 2.2.2反幂法的概念
设n n A ?∈ 这个矩阵是非奇异的，A 的特征值按非别从大到小写为
120n λλλ≥≥≥≥
则它所对应的特征向量是12,,,n x x x ，那么1A -的特征值如下
对应的特征向量为11,,,n n x x x -
因此无论是计算A 按模最小的特征值n λ，还是计算1A -的按模最大的特征值问题，两者没有本质区别[9]。
对于1A -使用幂法迭代，这种方法叫做反幂法。可计算出矩阵1A -的主特征值1/λ，进一步计算出矩阵A 按模最小的特征值n λ。
反幂法迭代公式为：设初始向量000v u =≠，向量序列构造如下
k=1,2，… 迭代向量k v 可以通过解线性方程组
1k k Av u -=
求得。
例4 用反幂法求
1
1
1
1
1
n
n λλλ-≥
≥≥
{}11
max k k k k k v A u v u v --?=??=??
的和计算特征值 1.2679λ=
对应（精确特征值是33λ=）的特征向量（限定运算为5位浮点数）。
解 用部分选主元的三角分解将A pI -（其中 1.2679p =）分解为
()P A pI LU -=
其中
由1(1,1,1)T Uv =，得
1(12692,9290.3,3400.8)T v =-
1(1,0.73198,0.26795)T u =-
由21LUv Pu =，得
2(20404,14937,5467.4)T v =-
2(1,0.73206,0.26796)T u =-
3λ对应的特征向量是
3(1,1(1,0.73205,0.26795)T T x =≈-
特征值321.26791/ 1.26794901u λ≈+=
210131014A ?? ?= ?
???
1
000100.73210.268071L ?? ?= ?
?-??
31 1.7321101 2.7321000.2940510U -??
?= ?
????010001100P ?? ?= ?
???
高阶对称矩阵特征值的计算
3λ
的真值为33 1.26794912λ==
通过上面这道例题，我们可以总结一下反幂法的大致步骤。 首先分解计算()P A pI LU -=，且保存L ，U 和P 的信息； 接着用反幂法迭代
（1）解1(1,1,,1)T Uv = 求1v
{}11max v μ=，111/u v μ=
（2）2,3,k =
①解1k k Ly Pu -=求k y 解k k Uv y =求k v ②{}max k k v μ= ③计算/k k k u v μ=
2.3 QR 方法
对于本文所讨论的高阶对称矩阵n n A ?∈ ，如果想用QR 方法计算其特征值，需要将其首先变换成上海森伯格矩阵或者对称三对角矩阵，然后再进行，此时就需要使用豪斯霍尔德变换和正交相似变换，通过这两种变换，将矩阵转化成上海森伯格矩阵，另一种情况是把对称矩阵转化成对称三角矩阵。 2.3.1豪斯霍尔德变换
设向量n w ∈ ，且1T w w =，称矩阵
()2T H w I ww =-
叫做豪斯霍尔德变换。它还有另外一种叫法，叫做初等反射矩阵。
如果，记12(,,,),n w T ωωω= 则
2.3.2用正交相似变换约化一般矩阵为上海森伯格矩阵
设()n n ij A a ?=∈ 。为了使A 经正交相似变换化成上海森伯格矩阵可选择初
2
11212
21
2221212222122()2212n n n n n H w ωωωωωωωωωωωω
ωωω??
---
?---
?
= ? ? ?---?
?