矩阵特征值的计算

λi = 15 − i , i = 1, 2, 3, 4
比值
r
=
| λ2 | λ1
| |
=
0.9
，令
p
= 12
作变换
B = A − pI
则 B 的特征值为 µ1 = 2 , µ2 = 1 , µ3 = 0 , µ4 = −1
应用幂法计算 B 的按模最大的特征值 µ1 时，确定收敛速度的比
值为：
µ2 = λ2 − p = 0.5 < | λ2 | ≈ 0.9
如果 p 是矩阵 A 的特征值 λi 的一个近似值，且
| λi − p |<| λ j − p | , i ≠ j
1 则 λ i − p 是矩阵 ( A − pI )−1 的按模最大的特征值。因此，当给
10
出特征值 λ i 的一个近似值 p 时，可对矩阵 A − pI 应用反幂法， 求出对应于 λ i 的特征向量。反幂法迭代公式中的 yi 通过方程组
机上计算时就可能溢出停机。为了避免这一点，在计算过程中，
常采用把每步迭代的向量 xk 进行规范化，即用 xk 乘以一个常
数，使得其分量的模最大为1。
� 幂法迭代公式：
⎧ yk = Ax k −1
⎪ ⎨
m
k
=
max(
yk )
⎪ ⎩
x
k
=
yk
/ mk
( k = 1, 2, ⋅ ⋅⋅)
其中 mk 是 yk 模最大的第一个分量。相应地取
1
1
2
2
-2
2
2
1
3
3
-4
3
-4 -0.75
4 -2.5 3.5 -2.5 3.5 -0.714
5 -2.428 3.428 -2.428 3.428 -0.708
6 -2.416 3.416 -2.416 3.416 -0.707
7 -2.414 3.414 -2.414 3.414 -0.707
意矩阵全部特征值的 QR 方法。
1
§1 幂法与反幂法
一、幂法
� 幂 法 ：是一种求任意矩阵 A 的按模最大特征值及其对应特
征向量的迭代算法。
该方法最大的优点是计算简单，容易在计算机上实现，对稀
疏矩阵较为合适，但有时收敛速度很慢。
为了讨论简单，我们假设：
(1) n 阶方阵 A 的特征值 λ1 , λ2 ,⋅ ⋅ ⋅, λn 按模的大小排列
方法为反幂法。
反幂法也是一种迭代算法，每一步都要解一个系数矩阵相同
的线性方程组。
� 原点平移的反幂法
设 p 为任一实数，如果矩阵 A − pI 可逆，则 ( A − pI )−1 的特
征值为
11
1
,
,⋅ ⋅ ⋅,
λ1 − p λ2 − p λn − p
对应的特征向量仍为 v1 , v2 ,⋅ ⋅ ⋅, vn 。
λ2 λ1
− −
p p
,
λn λ1
− −
p p
⎫ ⎬ ⎭
=
min
当 λ2
−
p
=
−(λn
−
p) ，即
p
=
λ2
+ λn 2
时， r
为最小。这时用
幂法计算 λ 1 及 v1 时得到加速。
如果需要计算 λn 及 vn 时，应选择 p 使
| λn − p |>| λ1 − p |
且确定收敛速度的比值：
r
=
max
称为幂法。
由(4)式知，幂法的收敛速度取决于比值
| |
λ2 λ1
| |
的大小。比
3
值越小，收敛越快，但当比值
| |
λ2 λ1
| |
接近于1时，收敛十分缓慢。
用幂法进行计算时，如果 | λ 1 |> 1 ，则迭代向量 xk 的各个
不为零的分量将随着 k 无限增大而趋于无穷。反之，如果
| λ 1 |< 1 ，则 xk 的各分量将趋于零。这样在有限字长的计算
v
1
(7)
这说明当 k 充分大时，两个相邻迭代向量 xk +1 与 xk 近似地相差
一个倍数，这个倍数便是矩阵 A 的按模最大的特征值 λ1 。若用
( xk )i 表示向量 xk 的第 i 个分量，则
λ1
≈
( xk +1 )i (xk )i
(8)
即两个相邻迭代向量对应分量的比值近似地作为矩阵 A 的按模
(A − λI)x = 0
的非零解，即是对应于 λ 的特征向量。这对于阶数较小的矩阵是
可以的，但对于阶数较大的矩阵来说，求解是十分困难，所以用 这种方法求矩阵的特征值是不切实际的。
� 如果矩阵 A 与 B 相似，则 A 与 B 有相同的特征值。 因此希望在相似变换下，把 A 化为最简单的形式。一般矩阵
矩阵特征值的计算
物理、力学和工程技术中的许多问题在数学上都归结为求矩 阵的特征值和特征向量问题。
� 计算方阵 A 的特征值，就是求特征多项式方程：
| A − λI |= 0 即 λn + p1λn−1 + ⋅ ⋅ ⋅ + pn−1λ + pn = 0
的根。求出特征值 λ 后，再求相应的齐次线性方程组：
⎧ ⎨ ⎩
λ1 λn
− −
p p
,
λn−1 − p λn − p
⎫ ⎬ ⎭
=
min
当 λ1
−
p
=
− ( λ n −1
−
p)
，即
p
=
λ1
+ λn−1 2
时， r
为最小。这时
用幂法计算 λ1 及 v1 时得到加速。
原点平移的加速方法，是一种矩阵变换方法。这种变换容易
计算，又不破坏 A 的稀疏性，但参数 p 的选择依赖于对 A 的特
( A − pI ) yk = xk −1 , k = 1, 2, ⋅ ⋅ ⋅
求得。
例3 用反幂法求矩阵
⎡ 2 −1 0 0⎤
A = ⎢⎢−1
2
−1
0
⎥ ⎥
⎢ 0 −1 2 −1⎥
矩阵 A 进行三角分解：
A = LU
(16)
再解三角形方程组：
⎧Lzk = xk −1 ⎩⎨Uyk = zk
k = 1, 2, ⋅⋅⋅
(17)
当 A 是三对角方阵，或是非零元素较少且分布规律的方阵
时，无论存储或计算都比较便。根据幂法的讨论，我们知道，在
一定条件下，可求得 A−1 的按模最大的特征值和相应的特征向 量，从而得到 A 的按模最小的特征值和对应的特征向量，称这种
| |
λ2 λ1
| |
的大
小。当比值接近于1时，收敛可能很慢，一个补救的方法是采用
原点平移法。
设矩阵
B = A − pI
(11)
其中 p 为要选择的常数。
A 与 B 除了对角线元素外，其它元素都相同，而 A 的特征
值 λi 与 B 的特征值 µ i 之间有关系 µ i = λi − p ，并且相应的特
当 k = 7 时， xk 已经稳定，于是得到： λ1 ≈ m 7 = 3 .414
及其相应的特征向量 v1 为：
v1 ≈ x7 = (−0.707 , 1, − 0.707 )T
xk
0
1
-1
1
1 -0.75
1 -0.714
1 -0.708
1 -0.707
1 -0.707
� 应用幂法时，应注意以下两点： (1)应用幂法时，困难在于事先不知道特征值是否满足(1)
⎧λ1 ≈ m k
⎨ ⎩
v1
≈
xk
( 或 yk )
例1 设
(9) (10)
⎡ 2 −1 0⎤
A = ⎢⎢−1
2
−1
⎥ ⎥
⎢⎣ 0 −1 2⎥⎦
用幂法求其模为最大的特征值及其相应的特征向量(精确到小数
点后三位)。
解 取 x0 = (1, 1, 1)T ，计算结果如表1所示。
4
表1
k
ykT
mk
1
1
0
1
(13)
为了防止溢出，计算公式为
⎧ Ay k = xk −1
⎪ ⎨
m
k
=
max(
yk )
( k = 1, 2, ⋅ ⋅⋅)
⎪ ⎩
x
k
=
yk
/ mk
(14)
相应地取
⎧ ⎪
λ
n
⎨
≈
1 mk
⎪⎩ v n ≈ y k ( 或 x k )
(15)
9
(13)式中方程组有相同的系数矩阵 A ，为了节省工作量，可先对
(2)
称为迭代向量。
由于 v1, v2 ,⋅⋅⋅, vn 线性无关，构成 n 维向量空间的一组基，所
以，初始向量 x0 可唯一表示成：
x0 = a1v1 + a2v2 + ⋅⋅⋅ + anvn
(3)
于是
xk = Axk −1 = A2 xk −2 = ⋅ ⋅ ⋅ = Ak x0
= a1λ1k v1 + a2λk2v2 + ⋅ ⋅ ⋅ + anλknvn
最大的特征值。
因为 xk +1 ≈ λ1 xk ，又 xk +1 = Ax k ，所以有 Ax k ≈ λ1 x k ，因
此向量 xk 可近似地作为对应于 λ1 的特征向量。
这种由已知的非零向量 x 0 和矩阵 A 的乘幂构造向量序列
{xk } 以计算矩阵 A 的按模最大特征值及其相应特征向量的方法
征值的分布有大致了解。
三、反幂法 � 反 幂 法 ：用于求矩阵 A 的按模最小的特征值和对应的特征
向量，及其求对应于一个给定的近似特征值的特征向量。 设 n 阶方阵 A 的特征值按模的大小排列为：
8
| λ1 |≥| λ2 |≥ ⋅ ⋅ ⋅ ≥| λn−1 |>| λn |> 0
相应的特征向量为 v1, v2 ,⋅⋅⋅, vn 。则 A−1 的特征值为：
xk = A k x0 。按照基向量 v1 , v2 ,⋅ ⋅ ⋅, vn 展开时， v1 的系数可能不 等于零。把这一向量 xk 看作初始向量，用幂法继续求向量序列 xk +1 , xk + 2 ,⋅ ⋅ ⋅ ，仍然会得出预期的结果，不过收敛速度较慢。
如果收敛很慢，可改换初始向量。
二、原点平移法
由前面讨论知道，幂法的收敛速度取决于比值
λ1 > λ2 ≥ λ3 ≥ ⋅ ⋅ ⋅ ≥ λn−1 > λn 则对于任意实数 p ， B = A − pI 的按模最大的特征值 λ1 − p 或λn − p 。
如果需要计算 λ 1 及 v1 时，应选择 p 使
7
| λ1 − p |>| λn − p |
且确定的收敛速度的比值：
r
=
max
⎧ ⎨ ⎩
征向量相同。
这样，要计算 A 的按模最大的特征值，就是适当选择参数
p ，使得 λ1 − p 仍然是 B 的按模最大的特征值，且使
λ2 − p < | λ2 | λ1 − p | λ1 |
对 B 应用幂法，使得在计算 B 的按模最大的特征值 λ1 − p
6
的过程中得到加速，这种方法称为原点平移法。
例2 设4阶方阵 A 有特征值
µ1 λ1 − p
| λ1 |
所以对 B 应用幂法时，可使幂法得到加速。
虽然选择适当的 p 值，可以使得幂法得到加速，但由于矩阵
的特征值的分布情况事先并不知道，所以在计算时，用原点平移
法有一定的困难。
下面考虑当 A 的特征值为实数时，如何选择参数 p ，以使
得用幂法计算 λ 1 时得到加速的方法。 设 A 的特征值满足：
式，以及方阵 A 是否有 n 个线性无关的特征向量。克服上述困难
的方法是：先用幂法进行计算，在计算过程中检查是否出现了预 期的结果。如果出现了预期的结果，就得到特征值及其相应特征 向量的近似值；否则，只能用其它方法来求特征值及其相应的特 征向量。
5
Fra Baidu bibliotek
(2)如果初始向量 x0 选择不当，将导致公式(3)中 v1 的系数 a 1 等于零。但是，由于舍入误差的影响，经若干步迭代后，
| λ1 |>| λ2 |≥ ⋅ ⋅ ⋅ ≥| λn |
(1)
(2) vi 是对应于特征值 λi 的特征向量 (i = 1,2,⋅ ⋅ ⋅, n) ；
(3) v1, v2 ,⋅⋅⋅, vn 线性无关。
任取一个非零的初始向量 x0 ，由矩阵 A 构造一个向量序列
xk = Axk −1 , k = 1, 2, ⋅ ⋅ ⋅
(12)
用(12)式计算向量序列{xk }时，首先要计算逆矩阵 A−1 。由于计
算 A−1 时，一方面计算麻烦，另一方面当 A 为稀疏阵时， A−1 不
一定是稀疏阵，所以利用 A−1 进行计算会造成困难。在实际计算
时，常采用解线性方程组的方法求 xk 。（12）式等价于：
Axk = xk −1 , k = 1, 2, ⋅ ⋅ ⋅
11
11
≤ ≤ ⋅⋅⋅ ≤
<
λ1 λ2
λn −1
λn
对应的特征向量仍然为 v1, v2 ,⋅⋅⋅, vn 。因此，计算矩阵 A 的按模
最小的特征值，就是计算 A−1 的按模最大的特征值。
� 反幂法的基本思想：把幂法用到 A−1 上。
任取一个非零的初始向量 x0 ，由矩阵 A−1 构造向量序列：
xk = A−1xk−1 , k = 1, 2, ⋅⋅⋅
的最简单的形式是约当（Jordan）标准形。由于在一般情况下，
用相似变换把矩阵 A 化为约当标准形是很困难的，所以设法对矩 阵 A 依次进行相似变换，使其逐步趋向于一个约当标准形，从而 求出 A 的特征值。
下面主要介绍求部分特征值和特征向量的幂法、反幂法，求 实对称矩阵全部特征值和特征向量的雅 可 比（Jacobi)方法，求任
=
λ1k [a1v1
+
a2
( λ2 λ1
)k
v2
+
⋅⋅⋅
+
an
( λn λ1
)k
vn
]
(4)
2
因为比值
| |
λi λ1
| |
<1
(i
= 1,2,⋅ ⋅ ⋅, n) ，所以
lim
k→∞
xk λ1k
= a1v1
(5)
当 k 充分大时，有
x k ≈ λ1k a1v1
(6)
从而
xk +1
≈
λ
k 1
+
1
a
1