7 矩阵特征值与特征向量的计算
设 A 为 n 阶方阵,所谓 A 的特征值问题是求数和非零向量 x ,使
Ax 成立。数
称作 A 的一个特征值,非零向量 x 称作与特征值对应的特征向量。求给定方阵的特征值
与特征向量是先求解特征方程
x
( )
|
E A
| 0
然后对应于每一个特征值 i,再求解退化的齐次线性方程组
从而得到 A 的特征值 i及对应的特征向量 x 。
(
)
iE A x
0
但是这种方法计算机很大,计算过程复杂,因此有必要研究相对简单的数值解法。本章
主要介绍三类计算特征值的方法:计算大型(稀疏)矩阵主特征的幂法与反幂法,计算中小
型(实对称)矩阵全部特征值的 Jacobi 法,计算中小型矩阵全部特征值的 QR 法。
7.1 特征值估计 吗
在矩阵特征值计算中,有时需要对特征值所在范围给出一个估计。这里介绍一种从矩阵
的元素出发,运用较简便的运算估计特征值的方法。
定义 7-1 设
A
(
a
ij
)
n m
C
,称由不等式
|
z a
ii
|
R
i
在复平面上确定的区域为矩阵 A 的第 i 个盖尔圆(Gerschgorin 圆),并用 iG 表示。
其中
R
i
称为盖尔圆 iG 的半径(
i
a
ij
|
|
n
j
j
1
i
1,2,
。
, )
n
定理 7-1 矩阵
A
(
a
ij
)
n m
C
n
G
j
j
1
i
的一切特征值均落在它的 n 个盖尔圆的并集中,即
(
i
1,2,
, )
n
。
证明 设是 A 的任一特征值,
x
(
,
x x
1
2
,
,
x
n
)T
令
|
x
i
0
| max |
1
i n
x
i
|
,则
ix 。由 Ax
0
0
x ,可得
是对应的特征向量。
n
a x
i
j
0
x
i
0
。即
)
j
(
j
1
(
a
)
x
i
0
ii
00
xa
i
j
0
j
n
j
j
1
i
0
于是有
a
ii
00
a
i
0
j
n
j
j
1
i
0
j
x
x
i
0
a
i
0
j
n
j
j
1
i
0
x
x
j
i
0
R
i
0
这表明任一特征值
0iG ,从而也在 A 的第 n 个盖尔圆的并集中。
A
例 7-1 估计矩阵
0.1
3
0.3
0.3
1
0.5
1
0.2
解 A 的 4 个盖尔圆为:
2 :|
G z
1 :|
G z
1| 0.6
0.2
0.1
1
0.1
0.3
0.2
0.5
4
的特征值范围。
3| 0.8
3 :|
G z
1| 1.8
4 :|
G z
4 | 0.6
画在复平面上其区域如图 7-1 所示。
图 7-1 例 7-1 盖尔圆分布图
于是 A 的全部特征值就在这 4 个盖尔圆的并集中。为了更确切地知道某个特征值落在
哪个或哪几个盖尔圆的并集中,给出如下第二盖尔圆盘定理。
定理 7-2 若 A 的 n 个盖尔圆中,有 m 个盖尔圆构成的一个连通域(所谓连通域,是指
其中的任意两点都可以用位于该区域内的一条折线连接起来),且该连通域与其余 n m 个
盖尔圆严格分离,则在该连通域中恰好有 A 的 m 个特征值(重特征值按重数重复计算)。特
别地,每个孤立的盖尔圆恰有 A 的一个特征值(证明从略)。
由定理 2 可知,在例 1 中 2G 与 4G 中各有 A 的一个特征值,而 1G 与 3G 构成的连通部分
中有两个特征值,但不能确定这两个特征值具体落在哪个盖尔圆中。
例 7-2 估计矩阵
1
0.5
解 A 的两个盖尔圆为:
A
0.8
0
的特征值范围。
1 :|
G z , 2 :|
在复平面上的区域如图 7-2 所示。
1| 0.8
G z
0 | 0.5
图 7-2 例 7-2 盖尔圆分布图
此时只能判断 A 的两个特征值落在 1G 与 2G 的并集中,至于是每个盖尔圆中各有一个特
征 值 还 是 两 个 特 征 值 都 落 在 其 中 一 个 盖 尔 圆 上 则 无 法 确 定 。 实 际 上 , 由 于
1 (1
2
j
1,2
在盖尔圆 1G 中。
0.6)
, 1,2
|
|
0.4
,所以两个特征值都不会在盖尔圆 2G 中,而是落
0.5
对于某些矩阵,可利用相似变换矩阵具有相同特征值的性质得到更确切的特征值范围。
,对 A 作相似
d 构成对角阵
,取正数 1
diag(
)ij n m
,
d d
1
,
d d
)n
D
A
设
a
d
(
,
,
,
,
2
2
n
变换,令
B DAD
1
(
d
d
i
j
a
)
ij n n
,由于 B 相似于 A ,所以 B 与 A 的特征值完全相同,又
,
,
d d
由于 B 与 A 的主对角线元素对应相等,所以 B 与 A 的盖尔圆圆心相同。这表明,若适当选
d ,可以改变盖尔圆的半径,从而有可能将相交的盖尔圆分离得到仅含一
取正数 1
d 的一般方法是:欲使 A 的第i 个盖尔圆 iG 的半径
个特征值的孤立盖尔圆。选取 1
大而其余盖尔圆变小,就取
,
2
n
1(
j
。
)
i
n
jd
2
,
,
d d
1
,
id ,其余
0.8
1
10
5
10
2
j
例 7-3 求矩阵
A
20
4
1
的特征值范围。
解 A 的 3 个盖尔圆为:
1 :|
G z
20 | 5.8
, 2 :|
G z
10 | 5
, 3 :|
G z
10 | 3
j
其中 1G 与 2G 相交,而 3G 孤立。记 3G 中所含的一个特征值为 3,如图 7-3 所示。
为分离 2G 与 1G ,可以让 A 的第 3 行元素绝对值变大,第 3 列元素绝对值变小。
现取
,则
diag(1,1,2)
D
B DAD
1
20
4
2
5
10
4
0.4
0.5
10
j
图 7-3 例 3 盖尔圆分布图
其 3 个盖尔圆分别是: 1 :|
G z
显然, B 的盖尔圆是 3 个孤立的盖尔圆,如图 7-4,注意,此情况下, 3G 的半径变大了。
10 | 4.5
20 | 5.4
10 | 6
图 7-4 例 7-3 分离后盖尔圆分布图
, 2 :|
j
G z
, 3 :|
G z
例 7-4 设矩阵
证明 由线性代数知, A 可逆的充分条件是|
按行严格对角占优,则 A 可逆。
A
| 0A ,而
)ij n n
a
A
(
|
|
(其中 j 是 A 的
j
n
j
1
特征值),所以只要证明
j 即可 (
0
j
1,2,
。
, )
n
设是 A 的任一特征值,则必存在某个盖尔圆 iG 使
a
ii
R
i
a
ij
。
j
i
若
j ,则有
0
a
ii
的任意性,得|
| 0A 。
j
i
a
ij
,而这与 A 按行严格对角占优矛盾,故应有
0 ,由
7.2 幂法与反幂法
在线性代数中,设 A 是 n 阶方阵,若 A 存在 n 个线性无关的特征向量,则称这 n 个特
征向量构成 A 的一个完全的特征向量组。
例如,对矩阵
A
3
2
0
1 1 0
4 3 0
1
0 2
通过求解特征方程,不难求出 A 的三个特征值为 1
2 0
3
0
5
0
为 1
性无关的特征向量。我们称方阵 A 可对角化,而 B 不可对角化。
3
2,
1,
B
,
1
2
, B 的三个特征值
。方阵 A 可以找到三个线性无关的特征向量,而方阵 B 找不到三个线
3
5
2
7.2.1 幂法
幂法的基本思想是构造一个向量序列使之逼近主特征值(矩阵的按模最大的特征值)对
应的特征向量,然后求出主特征值。该方法简单易行,但收敛速度较慢。
A
,
。已知 A 的主特征值是单根 1,即特征值满足条件
1
有 一 个 完 全 的 特 征 向 量 组 1
现 设
,
,
)ij n n
a
,
x x
2
(
,
,
2
n
x , 其 对 应 的 特 征 值 是
n
|
1
2
|
|
|
|
|n
任取一个非零初始向量 0u ,由矩阵 A 构造向量序列
Au
0
Au
1
2
A u
0
Au
k
k
1
A u
0
u
1
u
2
u
1
k
由于 A 的完全特征向量组可以作为向量空间 nR 的一组基,因此 0u 可由 1
,
x x
2
,
性表示,即有
u
0
a x
1 1
a x
2 2
a x
n n
(设 1
a )
0
x 线
,
n
于是
其中
k
n
i
2
1
a
i
(
i
k
)
x
i
k
1
k
A u
0
i
1
k
1
ku
a x
1 1
u
k
a x
1 1
k
k
x
a
a
1
1 1
2 2
n
i
2
1
k
)
a
(
i
i
x
x
2
i
a
k
n
x
n
n
k
1
(
a x
1 1
k
)
。注意到
(1
i
,2
),
n
,故当 k 时,
k ,因此有
0
由于 1x 是主特征值 1对应的特征向量,其乘上常数因子 1
分大时,迭代向量 ku 是 1的特征向量的近似向量。
ka 仍为 1的特征向量,故当 k 充
1
为了利用迭代向量求出主特征值 1的近似值,设 (
)k iu 表示 ku 的第 i 个分量,则
(
)
1
k
(
)
k
(
a x
1
1
(
a x
1
1
(
u
k
(
u
)
1
)
[
1
)
i
)
]
k
i
i
i
i
i
于是
lim
k
(
u
k
(
u
)
1
)
i
k
i
1
这表明两相邻迭代向量对应分量的比值收敛于主特征值,亦即当 k 充分大时,可用两相邻迭
代向量的分量比作为主特征值的近似值,即
1
(
u
k
(
u
i
)
1
)
i
k
若主特征值是 A 的 r 重实特征值,即 1
2
r
(1
r
n
)
,对应的 r 个线性无
关特征向量为 1
,
x x
2
,
当 k 充分大时,
x 。则有
,
n
u
k
k
A u
0
k
1
r
i
1
a x
i
i
n
1
i r
a
i
(
i
1
k
)
x
i
u
k
k
1
a x
i
i
r
i
1
即 ku 仍为主特征值对应的特征向量的近似向量,相邻两迭代向量的分量比仍为主特征
值的近似值。综上所述,有
定理 7-3 设 A 是 n 阶实矩阵,具有完全的特征向量组,主特征值是 r 重根,即
|
1
r
|
|
|
|
1
r
|
|
n
|
2
(1
r
)
n
则对任意非零初始向量 0u ,迭代向量
ku
k
A u
0
满足
或
r
(
a
1
a x
i
i
u
lim
k
k
k
1
(
)
u
, 1
k
k
1
(
)
u
k i
a x
i
i
u
1
1
k
r
i
i
0)
,
lim
k
(
u
k
(
u
)
1
)
k i
i
1
i
1
这样用非零初始向量 0u 及矩阵 A 构造向量序列{ }ku 以计算 A 的主特征值 1及相应的
1 ,迭代向
特征向量的方法称为幂法。不过从上面的讨论中可以看到,如果 1
量 ku 当 k 时,其不为零的分量就会趋于无穷大或趋于零。为克服这个缺点,可以在每
步迭代中加上对向量规范化的步骤,使迭代向量的数量级保持在一个稳定的量级上,归纳起
| 1 或
1
|
来,幂法的计算步骤为:
步骤 1 给定非零初始向量 0u ,精度 1
max(
步骤 2 迭代
, ,令 0
v
,其中
2
)
u
k
)
0
u ;令 (0)
, 1k ;
1
max( ku 表示 ku 绝对值最大的分量;
max(
)v
0
)
,则 kv 即为 1的近似特征向量, (
k 即
|
1
)
k
k
v
k
ku
Av
, (
1
1
u
k
)
max(
(
且 (
|
1
1
1
1 k
,转步骤 2 继续迭代。
1)
2
;
u
k
k
k
)
步骤 3 规范化
v
k
步骤 4 若
v
1
k
k
为 1的近似值;否则,
例 7-5 用幂法计算
A
1.0
1.0
1.0
1.0
0.5 0.25
0.5
0.25
2.0
的主特征值和相应的特征向量,结果见表 7-1。
表 7-1
K
T
Kv (规范化向量)
0
1
5
10
15
16
17
18
19
20
而此题的准确值为
(1,1,1)
(0.9091,0.8182,1)
(0.7651,0.6674,1)
(0.7494,0.6508,1)
(0.7483,0.6497,1)
(0.7483,0.6497,1)
(0.7482,0.6497,1)
(0.7482,0.6497,1)
(0.7482,0.6497,1)
(0.7482,0.6497,1)
max(
)ku
1
2.7500000
2.558791
2.558791
2.5366256
2.5365840
2.5365598
2.5365598
2.5365374
2.5365323
1
2.5365258
x
1
(0.748221,0.649661,1.000000) T
7.2.2 幂法的加速
幂法的收敛速度由比值
r
2
1
来确定, r 越小收敛越快,而当 1r 时收敛可能很慢。
为了克服这一缺点,常采用原点平移法对幂法进行加速。
,其中 p 是待定参数。显然,若 A 的特征值为 1
设 B A pE
1,2,
应特征值 (
ik i
p
为 1
2
,而对 B 有特征方程|
| 0
为对 A 有特征方程|
i
,则 B 的相
,且 A 、 B 的特征向量相同。这是因
| 0
B k E
,所以
p
, )
n
E
,
p k k
i
)
p k E
n
i
|
|
A
A
p
p
(
,
,
,
2
n
i
i
,
i
,
,
i
另一方面,若 ix 是 A 的对应 i的特征向量,即
则
)
p x
i
i
原点平移法的思想是引入矩阵 B ,恰当的选择参数 p ,使 1
p
k
)
A pE x
i
Bx
i
px
i
(
(
1
是 B 的主特征值,
Ax
i
x
i
i
Ax
i
且其速比
r
B
max
2
1
p
p
r
A
2
1
,这样用幂法求 B 的主特征值 1k 的收敛速度就快于用
幂法求 A 的主特征值 1,而一旦 1k 求出,由 1
k
p 可得 A 的主特征值,达到了加速的
1
目的。但是为了选取恰当的选择参数 p ,需要对 A 的特征值的分布的大致了解。
例 7-6 设 4 阶方阵 A 有特征值
j
15
j
(
j
1,2,3,4)
其速比
Ar
2
1
。作变换
0.9
B A pE
(
p
12)
则 B 的特征值为 1
k , 2
k , 3
k , 4
2
1
0
k ,其速比
1
r
B
k
2
k
1
。
r
A
1
2
设 A 的实特征值满足
1
2
1n
n
若 2,
n 的值可大致估计出,若要求 1,考察待定参数 p 的选取。
在原点平移法通过变换
AB
pE
后,不论 p 如何选取,矩阵的 B 主特征值也只能是
在 n
p 或 1
p 。若希望求 1,就应选择 p ,使 1
p 称为 B 的主特征值,即
|
1
这时 B 的收敛速比 Br 是比值 2
p
|
r
B
|
|
|
p
p
n
| / |
|
|
p
和
1
n
|
|
|
p
2
|
|
|
p
1
max
,
n
1
p
| / |
|
p
|
p
1
中的较大者,即
p
|
显然 Br 依赖于 p 的选取,记做 (
快,我们希望选择最佳参数 *p ,使
Br p 。为了使应用幂法求 B 的主特征值的收敛速度尽可能
)
*
) min (
)
(
r p
r p
B
B
*prB
( prB 的最小化要求,只有当
(
)
由 Br 的表示式(求二者之间的较大值)和
对
|
|
|
p
p
2
n
p
会有得到矛盾的结果( 2
n
p
2
Br p 达到最小。由于 2
时, (
(
n
p
p
)
)
)
|
n ),所以只能是
即
*
p
类似地,若用反幂法求最小特征值 n ,若 1n , 1可大致估计,取最佳平移参数
2
n
2
*
p
1
1
n
2
p
例 7-7 取
解 作变换 B A pE
0.75
,则
,用原点平移法,计算例 7-7 中矩阵 A 的主特征值。
B
0.25
1
0.5
0.5
1
0.25 0.25
0.25 1.25
对 B 应用幂法,计算结果见表 7-2。即 1 1.7865914
,则 A 的主特征值 1为
k
1 0.75 2.5365914
k
1
与例 7-5 比较,上述结果比例 7-5 迭代 15 次还好。
表 7-2
k
0
5
6
7
8
9
10
T
Kv (规范化向量)
(1,1,1)
(0.7516,0.6522,1)
(0.7491,0.6511,1)
(0.7488,0.6501,1)
(0.7484,0.6499,1)
(0.7483,0.6497,1)
(0.7482,0.6497,1)
max(
)ku
1
1.7914011
1.7888443
1.7873300
1.7869152
1.7866587
1.7865914
7.2.3 反幂法
设方阵 A 按模最小的特征值是 n ,且
n ,则 A 可逆。于是,由 n
Ax
0
x
n n
,可得
1
A x
n
1
n
x
n
,这表明
1
n
是 1A 的主特征值。反幂法就是将幂法应用于 1A ,通过求出 1A
的主特征值得到 A 的按模最小的特征值及其对应的特征向量。
定理 7-4 设 A 是 n 阶实矩阵,具有完全的特征向量组,其特征值满足
|
1
2
|
|
|
|
n
| 0
u
则对任意非零初始向量 0
满足
v ,按下述方式构造的迭代向量
0
u
k
1
A v
k
1
,
v
k
u
k
max(
u
k
)
lim
k
v
k
x
n
max(
x
n
)
,
lim max(
k
u
)k
1
n
v
k
x
n
/ max(
x
n
)
,
max(
u
)k
1
n
在实际计算中,可先对 A 进行 LU 分解,通过求解
y
Ly
Uu
, k
v
1
k
k
k