数值分析第8章——矩阵特征值问题计算

1k
a1x1
a2
2 1
k
x2
an
n 1
k
xn
1k a1x1 k
16
其中
k
a2
2 1
k
x2
an
n 1
k
xn
由假设条件 从而
j 1
1 j 2,,n, 所以
lim
k
vk
1k
a1x1
lim
k
k
0
所以当k充分大时，有
vk 1k a1x1
17
vk 1k a1x1
即为矩阵 A 的对应特征值 1 的近似特征向量。
0.5, 1, 0.7501 0.5, 1, 0.7500 0.5, 1, 0.7500
从而
24
二、幂法的加速
1、原点平移法
如果是矩阵 A 的特征值,则对任意的实数p, 矩阵
A-pE 的特征值为 -p,且 A 与 A-pE 的特征向量相同.
据此, 如果要计算 A 的主特征值 1 , 只要选择合适的 数 p,使 1-p 为矩阵 A-pE 的主特征值,且
n
( Ax, x) (x, x)
1
（2）
n
min x0
( Ax, x) (x, x)
（3）
1
max x0
( Ax, x) (x, x)
6
定理6 （Gerschgorin圆盘定理） 设AR nn，则
（1）A的每一个特征值必属于下述某个圆盘之中，
n
z aii aij j 1 ji
, i 1, 2,, n
i xi
A1xi
1 i
xi
其相应的特征向量 x1 , x2, …, xn 线性无关,则 A-1 的特
征值为1/ i ，对应的特征向量仍为 xi (i=1,2, …,n).
33
此时,A-1 的特征值满足 1 1 1
n n1
1
因此,对 A-1 应用幂法,可求出其主特征值
1/ n k
和特征向量 xn uk ,
对应的特征向量。 解 取p=-2.5, 做平移变换B=A-pE，
则
2 3 2
4.5 3 2
A 10 3 4 3 6 1
B
10
5.5
4
3 6 3.5
对B应用幂法，仍 取 x0=(0,0,1)T , 则
v0 u0 0 vk Auk1,
k max(vk ),
(k 1,2,)
v1
2,4,3.5T
解Lzk uk1, 求出zk
解Uvk zk ,
k
maxvk
求出vk k 1,2,
uk vk / k
n
1
k
,
xn uk
35
对给定的误差 ，当 | k – k1 | < 时，得
n
1
k
,
xn vk
显然，反幂法的收敛速度取决于比值
n ，比值越小，
n1
收敛越快。
例3 用反幂法求矩阵 2 1 0 0
39
若A为奇异方阵,则零为A的特征值.任取一数 p 不是A 的特征值,则 A-pI 为非奇异方阵.只要求 出 A-pI 的特征值,就很容易求出A的特征值, 所以 假设A 为非奇异方阵,并不防碍讨论的一般性.
设A 为非奇异方阵,令 A1 A , 对A1 进行QR 分解, 即把 A1 分解为正交矩阵 Q1与上三角形矩阵 R1 的乘积 A1 Q1R1
特征向量仍为v1 =(0.5,1,0.7500)T。
28
29
30
2、瑞利商加速 设A为n阶实对称矩阵，称
R x
Ax, x x, x
为向量 x 的瑞利商，其中 ( x, x)= xT x 为内积。不难
证明，对实对称矩阵A，如果其特征值满足
1 2 n1 n
由幂法公式生成的 xk 的瑞利商满足
38
8.4 QR算法
QR 算法也是一种迭代算法,是目前计算任 意实的非奇异矩阵全部特征值问题的最有效的
方法之一.该方法的基础是构造矩阵序列 Ak,并
对它进行QR 分解. 由线性代数知识知道,若A 为非奇异方阵,则
A 可以分解为正交矩阵Q 与三角形矩阵R 的乘 积,即A=QR ,而且当R的对角线元素符号取定时, 分解式是唯一的.
性质1 所有都相似,它们具有相同的特征值.
41
证明 因为 Ak1 RkQk QkT AkTk QkTQkT1 Ak1Qk1Qk QkTQkT1Q1T AQ1Q2 Qk
征向量 x1 , x2, …, xn ，即
Axi i xi (i 1,2,, n)
任取非零向量 v0 , 则可
唯一表示为 v0 a1x1 a2 x2
an xn
15
v0 a1x1 a2 x2 an xn
则
vk Akv0 Ak a1x1 a2 x2 an xn
a11k x1 a2k2 x2 ankn xn
（1）A与B的特征值完全相同； （2）若x是B的特征向量，则Px便为A的特征向量。
2
定理2： 设AR nn具有完全的特征向量系，即存在n个线性无关
的特征向量构成Rn的一组基底，则经相似变换可化A为 对角阵，即有可逆阵P，使
1
P 1 AP D
2
n
其中i为A的特征值，P的各列为相应于i的特征向量。
1 mk , v1 xk .
32
三、反幂法
反幂法是用于求非奇异矩阵A的按模最小的特征值 和对应特征向量的方法. 而结合原点平移法的反幂法 则可以求矩阵A的任何一个具有先验了解的特征值和对 应的特征向量。
设矩阵A非奇异,其特征值i (i=1,2,…,n) ,满足
1 2 n1 n 0
Axi
从而求得A的按模最小特征值
n 1/ k
和对应的特征向量 xn uk ，
这种方法称为反幂法。
34
v0 u0 0
vkk
A1uk 1 , max(vk ),
uk vk k
(k 1,2,)
为了避免求 A-1 ,可通过解线性方程组A vk= uk-1 得到yk ， 采用LU分解,即先对 A 进行LU分解 A=LU , 此时反幂法 的迭代公式为
第八章 矩阵特征值问题计算
对n 阶方阵A求数 和非零向量x ,使其满足Ax=x 这样的 值称为矩阵A的特征值，非零向量 x 称为矩 阵A的与特征值 相对应的一个特征向量。
1
8.1 预备知识
定义1 设矩阵A, BR nn，若有可逆阵P，使 B P1AP
则称A与B相似。
定理1 若矩阵A, BR nn且相似，则
13.5000
uk
0.5, 1, 0.875
0.5, 1, 0.7545 0.5, 1, 0.7507 0.5, 1, 0.7500 0.5, 1, 0.7500
可得到B的主特征值 113.5000
特征向量
v1 (0.5 ,1.0, 0.7500)T
因此,A的主特征值为
1 = 1 +p 11.0000,
n
表示以aii为中心，以 aij 半径为的复平面上的n个圆盘。 j 1 ji
（2）如果矩阵A的m个圆盘组成的并集S（连通的）与其余
n – m个圆盘不连接，则S内恰包含m个A的特征值。
7
8
9
10
11
12
定理7
13
14
8.2 幂法和反幂法
一 幂法
1 基本思想
设n 阶矩阵A 的特征值 i (i 1,2,, n) , 满足 | 1 || 2 | | n | ,且其对应有n个线性无关的特
1 p n p
且使
max
2 1
p p
,
n 1
p p
min
显然,当 2 - p = - (n- p ),即 P= ( 2+ n ) 2 时,上式取 最小值;如果希望计算 n , 类似的讨论可知应选取 p= (1 + n-1)2 。
26
例2 用原点平移加速法求例1中矩阵A的主特征值与其
,
1
4, u1uk1v1k
k
v1
0.5,
1,0.875
T
27
迭代5步的计算结果见下表
k
vk
k
1 2, 4,
3.5
4
2 7, 14, 10.5625 14
3 6.76, 13.5179, 10.1406
13.5179
4 6.7503, 13.5007, 10.1256
13.5007
5 6.7500, 13.5000, 10.1250
lim
k
vk
1k
a1x1
lim (vk1)i k (vk )i
1
19
迭代公式 (1)实质上是由矩阵A的乘幂 Ak与非零 向量 v0 相乘来构造向量序列 { xk }, 从而计算主特征 值及其对应的主特征向量,故称这种方法为幂法。
20
定理 设 A Rnn 有 n 个线性无关的特征向量，
主特征值 1 满足 1 2 3 n
的对应于特征值 的特征向量
A 1
2
1
0
0 1 2 1
0 0 1 2
0.4
36
解 取 x0 1,1,1,1T , 解方程组 A 0.4I y1 x0 得 y1 40,65,65,40T ,
m1 max y1 65,
x1
1 m1
y1
8
13 ,1,1 , 8
13
T
.
k max(vk ), uk vk k
(k 1,2,)
(1)
lim
k
uk
x1 max(
x1)
(2)
lim
k
k
1
22
例1 求矩阵
2 3 2
A 10 3 4 3 6 1
的主特征值与其对应的特征向量。
解 取 v0=(0,0,1)T , 则 v0 u0 0
vk Auk1,
k max(vk ),
11.4444 10.9223 11.0142
0.5, 1, 0.7360 0.5, 1, 0.7536 0.5, 1, 0.7494
6 5.4987, 10.9974, 8.2494 7 5.5002, 11.0005, 8.2501 8 5.5000, 11.0000, 8.2500
10.9974 11.0005 11.0000
且
vk 1 Avk 1k 1a1x1 1vk
用 (vk)i 表示 vk 的第 i 个分量,则当k充分大时,有
vk1 i
vk
i
1
即为主特征值的近似值。
18
定理 设 A Rnn 有 n 个线性无关的特征向量，
主特征值 1 满足
1 2 3 n
则对任意非零初始向量 v0 0 ，下面的式子成立
(k 1,2,)
uk vk k
v1
2,4,1T
,
1
4, u1
1
1
v1
0.5, 1,0.25T
23
直到k=8 时的计算结果见下表
k
vk
k
1 2, 4, 1,
4
uk
0.5, 1, 0.25
2 4.5, 9, 7.75
9
0.5, 1, 0.8611
3 5.7222, 11.4444, 8.361 4 5.4621, 10.9223, 8.2306 5 5.5075, 11.0142, 8.2576
解方程组 A 0.4I y2 x1
得 y2 445 13, 720 13, 720 13, 445 13T ,
m2 max y2 720 13,
37
来自百度文库
x2
1 m2
y2
89
144 ,1 ,1 , 89
144T
.
x1 与x2 的对应分向量大体上成正比, 所以对应于
0.4 的特征向量为 v 89 144,1,1,89 144T .
max 2in
i 1
p p
2 1
那么,对矩阵 A-pE 应用幂法求其主特征值 1-p ,收敛
速度将会加快.这种通过求 A-pE 的主特征值和特征
向量,进而得到A的主特征值和特征向量的方法叫原点
平移法。
25
则对任意实数p,矩阵 A-pE 的主特征值为1-p或 n-p , 要求 1 , 则选 p 使
3
定理3 ：AR nn，1, …, n为A的特征值，则
（1）A的迹数等于特征值之和，即
n
n
tr( A) aii i
i 1
i 1
（2）A的行列式值等于全体特征值之积，即
det( A) 12 n
4
定理4
5
定理5 设AR nn为对称矩阵，其特征值1≥2≥…≥n，则
（1）对任意AR n，x≠0，
Rxk
Axk xk ,
, xk xk
1
2 1
2k
由此可见，R(xk) 比 mk 更快的收敛于1 。
31
幂法的瑞利商加速迭代公式为
yk Axk1
mk
yk , xk1 xk1, xk1
xk yk / mk
k 1,2,
其中A为n阶实对称矩阵。
对给定的误差限 ，当 | mk – mk-1 | < 时，取
作矩阵 A2 R1Q1 Q1T A1Q1 继续对 A2进行QR分 解 A2 Q2R2 并定义 A3 Q2R2 Q2T A2Q2
40
一般地,递推公式为
A1 A Q1R1
Ak
1
RkQk
QkT AkQk
k 2,3,
QR 算法就是利用矩阵的QR 分解,按上述递推公
式构造矩阵序列 Ak .只要A 为非奇异方阵,则由QR算 法就完全确定 Ak 这个矩阵序列 Ak 具有如下性质.
则对任意非零初始向量 v0 0
按照下述方法构造的向量序列 uk , vk
v0 u0 0
vkk
A uk 1, max(vk
),
(k 1,2,)
uk vk k
则有
(1)
lim
k
uk
x1 max( x1)
(2)
lim
k
k
1
21
2. 幂法实用计算公式
v0 u0 0 vk Auk1,