矩阵特征值的计算

矩阵特征值的计算
物理、力学和工程技术中的许多问题在数学上都归结为求矩 阵的特征值和特征向量问题。
� 计算方阵 A 的特征值，就是求特征多项式方程：
| A − λI |= 0 即 λn + p1λn−1 + ⋅ ⋅ ⋅ + pn−1λ + pn = 0
的根。求出特征值 λ 后，再求相应的齐次线性方程组：
(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 ，为了节省工作量，可先对
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, ⋅⋅⋅
如果 p 是矩阵 A 的特征值 λi 的一个近似值，且
| λi − p |<| λ j − p | , i ≠ j
1 则 λ i − p 是矩阵 ( A − pI )−1 的按模最大的特征值。因此，当给
10
出特征值 λ i 的一个近似值 p 时，可对矩阵 A − pI 应用反幂法， 求出对应于 λ i 的特征向量。反幂法迭代公式中的 yi 通过方程组
(12)
用(12)式计算向量序列{xk }时，首先要计算逆矩阵 A−1 。由于计
算 A−1 时，一方面计算麻烦，另一方面当 A 为稀疏阵时， A−1 不
一定是稀疏阵，所以利用 A−1 进行计算会造成困难。在实际计算
时，常采用解线性方程组的方法求 xk 。（12）式等价于：
Axk = xk −1 , k = 1, 2, ⋅ ⋅ ⋅
征向量相同。
这样，要计算 A 的按模最大的特征值，就是适当选择参数
p ，使得 λ1 − p 仍然是 B 的按模最大的特征值，且使
λ2 − p < | λ2 | λ1 − p | λ1 |
对 B 应用幂法，使得在计算 B 的按模最大的特征值 λ1 − p
6
的过程中得到加速，这种方法称为原点平移法。
例2 设4阶方阵 A 有特征值
λ2 λ1
− −
p p
,
λn λ1
− −
p p
⎫ ⎬ ⎭
=
min
−(λn
−
p) ，即
p
=
λ2
+ λn 2
时， r
为最小。这时用
幂法计算 λ 1 及 v1 时得到加速。
如果需要计算 λn 及 vn 时，应选择 p 使
| λn − p |>| λ1 − p |
且确定收敛速度的比值：
r
=
max
意矩阵全部特征值的 QR 方法。
1
§1 幂法与反幂法
一、幂法
� 幂 法 ：是一种求任意矩阵 A 的按模最大特征值及其对应特
征向量的迭代算法。
该方法最大的优点是计算简单，容易在计算机上实现，对稀
疏矩阵较为合适，但有时收敛速度很慢。
为了讨论简单，我们假设：
(1) n 阶方阵 A 的特征值 λ1 , λ2 ,⋅ ⋅ ⋅, λn 按模的大小排列
机上计算时就可能溢出停机。为了避免这一点，在计算过程中，
常采用把每步迭代的向量 xk 进行规范化，即用 xk 乘以一个常
数，使得其分量的模最大为1。
� 幂法迭代公式：
⎧ yk = Ax k −1
⎪ ⎨
m
k
=
max(
yk )
⎪ ⎩
x
k
=
yk
/ mk
( k = 1, 2, ⋅ ⋅⋅)
其中 mk 是 yk 模最大的第一个分量。相应地取
⎧ ⎨ ⎩
λ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 的特
的最简单的形式是约当（Jordan）标准形。由于在一般情况下，
用相似变换把矩阵 A 化为约当标准形是很困难的，所以设法对矩 阵 A 依次进行相似变换，使其逐步趋向于一个约当标准形，从而 求出 A 的特征值。
下面主要介绍求部分特征值和特征向量的幂法、反幂法，求 实对称矩阵全部特征值和特征向量的雅 可 比（Jacobi)方法，求任
矩阵 A 进行三角分解：
A = LU
(16)
再解三角形方程组：
⎧Lzk = xk −1 ⎩⎨Uyk = zk
k = 1, 2, ⋅⋅⋅
(17)
当 A 是三对角方阵，或是非零元素较少且分布规律的方阵
时，无论存储或计算都比较便。根据幂法的讨论，我们知道，在
一定条件下，可求得 A−1 的按模最大的特征值和相应的特征向 量，从而得到 A 的按模最小的特征值和对应的特征向量，称这种
(A − λI)x = 0
的非零解，即是对应于 λ 的特征向量。这对于阶数较小的矩阵是
可以的，但对于阶数较大的矩阵来说，求解是十分困难，所以用 这种方法求矩阵的特征值是不切实际的。
� 如果矩阵 A 与 B 相似，则 A 与 B 有相同的特征值。 因此希望在相似变换下，把 A 化为最简单的形式。一般矩阵
| λ1 |>| λ2 |≥ ⋅ ⋅ ⋅ ≥| λn |
(1)
(2) vi 是对应于特征值 λi 的特征向量 (i = 1,2,⋅ ⋅ ⋅, n) ；
(3) v1, v2 ,⋅⋅⋅, vn 线性无关。
任取一个非零的初始向量 x0 ，由矩阵 A 构造一个向量序列
xk = Axk −1 , k = 1, 2, ⋅ ⋅ ⋅
(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
µ1 λ1 − p
| λ1 |
所以对 B 应用幂法时，可使幂法得到加速。
虽然选择适当的 p 值，可以使得幂法得到加速，但由于矩阵
的特征值的分布情况事先并不知道，所以在计算时，用原点平移
法有一定的困难。
下面考虑当 A 的特征值为实数时，如何选择参数 p ，以使
得用幂法计算 λ 1 时得到加速的方法。 设 A 的特征值满足：
⎧λ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
征值的分布有大致了解。
三、反幂法 � 反 幂 法 ：用于求矩阵 A 的按模最小的特征值和对应的特征
向量，及其求对应于一个给定的近似特征值的特征向量。 设 n 阶方阵 A 的特征值按模的大小排列为：
8
| λ1 |≥| λ2 |≥ ⋅ ⋅ ⋅ ≥| λn−1 |>| λn |> 0
相应的特征向量为 v1, v2 ,⋅⋅⋅, vn 。则 A−1 的特征值为：
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
式，以及方阵 A 是否有 n 个线性无关的特征向量。克服上述困难
的方法是：先用幂法进行计算，在计算过程中检查是否出现了预 期的结果。如果出现了预期的结果，就得到特征值及其相应特征 向量的近似值；否则，只能用其它方法来求特征值及其相应的特 征向量。
5
(2)如果初始向量 x0 选择不当，将导致公式(3)中 v1 的系数 a 1 等于零。但是，由于舍入误差的影响，经若干步迭代后，
最大的特征值。
因为 xk +1 ≈ λ1 xk ，又 xk +1 = Ax k ，所以有 Ax k ≈ λ1 x k ，因