矩阵计算与分析幂迭代法和逆幂迭代法

4.3.2 收缩方法
用幂法可求A的主特征值λ1，但不能求其他特征值. 现设|λ1| > |λ2| > |λ3| ≥ … ≥ |λn|．在求出λ1之后，用收 缩的方法来求λ2．
收缩的方法有多种，这里讨论一种以相似变换为基 础的方法．设λ1和x1已知，可找到一个Householder矩 阵H，使H x1 = k1e1 ，且有k1≠0．由Ax1=λ1x1，有
vk+1 ≈ α1λk1+1 x1 + α2λk2+1 x2
vk+2 ≈ α1λk1+2 x1 + α2λk2+2 x2
于是 vk+2 −（λ1 + λ 2）vk+1 + λ1λ 2vk
≈ [λk1+2 −（λ1 + λ 2）λk1+1 + λ1λ 2λk1 ]α1 x1
+ [λk2+2 −（λ1 + λ 2）λk2+1 + λ1λ 2λk2 ]α2 x2 = 0
=
λ1[1
+
O(
λ2 λ1
k
)]
→
λ1
k→∞
于是有
定理2 设矩阵 A∈Rn×n 有n个线性无关的特征向量， 其按模最大的特征值 λ1 满足
|λ1| > |λ2| ≥ |λ3| ≥ … ≥ |λn| 则对任意的非零向量 v0 (α1 ≠ 0)，按下列方法构造向量序列
⎧u0 = v0
⎪⎪⎨umkk
= =
4.3 幂迭代法和逆幂迭代法
4.3.1 幂迭代法
在实际问题中，具体要求也有不同．某些应用场合通常不 需要知道矩阵的全部特征值和特征向量，而仅要求得到矩阵 的按模最大的特征值（或称为矩阵的主特征值）和相应的特 征向量，如线性方程组迭代解法的收敛性问题．有的则要求 全部特征值和特征向量．根据这2种不同要求，矩阵的特征值 和特征向量的计算方法也大体上分成2种类型．本章就此2种 类型介绍其中最常用的2 种方法 ⎯ 幂法和QR法．
lim (vk )i k→∞ (vk−1 )i
=
λ1
lim
k→∞
α1( x1 )i + (εk )i α1( x1 )i + (εk−1 )i
= λ1
这种由已知非零向量 v0 及矩阵的乘幂 Ak 构造向量序列{vk} 以计算 A的按模最大的特征值λ1的方法称为乘幂法，简称为幂 法．而定理的证明过程事实上已给出了幂法的实施步骤．
∴ vk = Akv0 = α1 Ak x1 + α2 Ak x2 + + αn Ak xn
= α1λk1 x1 + α2λk2 x2 + + αnλkn xn
∑ =
λk1(α1 x1
+
n i=2
αi
⎜⎛ ⎝
λi λ1
⎟⎞k ⎠
xi
)
(4)
∑ 同理
vk −1
=
λk1−1(α1 x1
+
n i=2
α
i
⎜⎛ ⎝
幂法主要用于求矩阵的按模最大的特征值和相应的特征向 量．它是通过迭代产生向量序列，由此计算特征值和特征向 量．其算法的思想基于如下结论：
定理1 设矩阵 A∈Rn×n 有n个线性无关的特征向量 xi (i =1, 2,…,n)，其对应的特征值 λi (i =1, 2,…, n)满足
|λ1| > |λ2| ≥ |λ3| ≥ … ≥ |λn| (1)
HAH −1e1 = λ1e1 即HAH −1的第一列为λ1e1，其中e1=(1,0, … ,0)T∈C n 记
A2
=
HAH −1
=
⎜⎜⎝⎛
λ1 0
b1T B2
⎟⎟⎠⎞
其中B2为n −1阶方阵，它显然有特征值λ2 , … , λn．
在|λ2| > |λ3|条件下，可用幂法求B2的主特征值λ2和主
特征向量y2 ，其中B2 y2 = λ2 y2．设A2 z2 = λ2 z2，为求
说明vk+2 , vk+1 , vk 三个向量大体上线性相关，令 p = − (λ1 +λ2)， q = λ1λ2，则有 vk+2 + pvk+1 + q vk ≈ 0 ．
选定2个分量代入上述关系，构成二阶线性方程组，解出p ， q，然后再验证 p , q 对vk+2 , vk+1 , vk的其余分量是否也满足线性 关系方程，若满足，则再将 p , q 代入方程 λ2 + pλ + q = 0 就 可求得矩阵 A的按模最大的一对复特征值λ1 , λ2．
所以对任意的非零向量 v0都可用xi (i =1, 2,…,n)线性表示，即
v0 = α1 x1 + α2 x2 + + αn xn 假定 α1 ≠ 0
⎧v1 = Av0
用A构造向量序列
⎪⎪v2 = Av1 = A2v0 ⎨
⎪⎪vk = Avk−1 = Akv0 ⎩
∵ Axi = λ i xi ⇒ Ak xi = λki xi i = 1,2, , n
A的对应于λ1 , λ2的近似特征向量为vk+1 − λ2vk ，vk+1 − λ1vk ． 根据以上的讨论，在用乘幂法进行计算时，当k充分大时, 检 查是否出现下列三种情况之一：
(Ⅰ) (vk+1 )i /( vk )i 趋近于某一常数； (Ⅱ) (vk+2 )i /( vk )i 趋近于某一常数； (Ⅲ) vk 的波动不规律，但相继的三个向量vk+2 , vk+1 , vk 之 间满足vk+2 + pvk+1 + q vk ≈ 0 ． 若是，就可求出 A的按模最大的特征值和相应的特征向量．
出z2，设α为待定常数，y为n −1维向量，则
z2
=
⎜⎛ ⎝
α ⎞⎟ y⎠
,
⎩⎨⎧λB12αy
+ b1T y = λ2y
(1) 任取一个初始非零向量 v0 .
⎧u0 = v0
(2) 构造迭代序列
⎪⎪⎨mukk
= =
Avk −1 max(uk
)
⎩vk = uk / mk
k = 1,2,
(4.3)
式中用记号max(u)表示向量u 的绝对值最大的分量，从而有
迭代
u1 = Au0 = Av0
u2
=
Av1
=
A2v0 max( Av0 )
于是当 k → ∞ 时 , 有 : (1)当 λ1 > 1时 , (vk )i → ∞ ;
(2)当 λ1 < 1时 , (vk )i → 0 . 分量的变化太大，使得计算过程可能产生溢出．为了克服这
一不足，通常采用规范化迭代方式．具体做法为：把迭代向量 uk的最大分量化为1，于是得乘幂法的计算步骤如下：
(3) |λ1| = |λ2|，但 λ1=ρe iθ，λ2 =ρe −iθ，即λ1 , λ2为一对共轭复特 征值，且|λ1| > |λ3| ≥ … ≥ |λn|，则对任意初始向量
v0 = α1 x1 + α2 x2 + + αn xn 假定 α1 ≠ 0 , α2 ≠ 0
当 k 充分大时有 vk ≈ α1λk1 x1 + α2λk2 x2
n i=2
αi
(λi λ1
)k
xi
+
n i=2
α
i
(
λi λ1
)k
) xi
))
∑∑ =
α1 x1 + max(α1 x1
i
n =2
α
i
(
λi λ1
)k
+
n i=2
α
i
(
λi λ1
xi )k
xi
)
→
x1 max(
x1
)
k→∞
这说明规范迭代向量序列收敛到按模最大特征值对应的特征向 量．
同理，可得到
∑∑ ∑∑ mk
规范化
v1
=
u1 max(u1 )
=
Av0 max( Av0
)
v2
=
u2 max(u2 )
=
A2v0 max( A2v0 )
uk
=
Avk −1
=
Akv0 max( Ak−1v0 )
vk
=
uk max(uk )
=
Akv0 max( Akv0 )
∑∑ 则有
vk
=
Akv0 max( Akv0 )
=
λk1(α1 x1 + max(λk1(α1 x1
λi λ1
⎟⎞ k −1 ⎠
xi )
由于 λi < 1 (i = 2,3, , n) , 故对足够大的k , 有 λ1
vk = λk1(α1 x1 + εk ) vk−1 = λk1−1(α1 x1 + εk−1 ) (5)
式中 εk−1 , εk ∈ Rn , 且当 k → ∞ 时 , εk−1 , εk → 0 .当( x1 )i ≠ 0 时有
=
max(uk ) =
=λ1mmaaxx⎢⎢⎢⎢⎣⎡(mαa1 xx1(λλ+k1k1(−iα1=n(21αxα11ix(+1λλi+1i=n)2kiα=nx2iiα()λλi (1i λλ)k1i
max(α1 x1
+
n i=2
αi
(
λi λ1
)k
−1
xi
)
xi ) )k −1
xi
⎤ ⎥ ⎥ ))⎥⎥⎦
Avk −1 max(uk
)
⎩vk = uk / mk
k = 1,2,
(4.3)
则有
lim
k→∞
max(uk
)
=
λ1
lim
k→∞
vk
=
x1 max( x1 )
例4.4
⎜⎛ − 4 A = ⎜− 5
14 13
0 ⎟⎞ 0⎟
有特征值6,3,2。选
⎜⎝ − 1 0 2⎟⎠
v0 =(1,1,1)T，按(4.3)计算，有 u1 = Av0 = (10,8,1)T , m1 = max(u1 ) = 10,
$
1) 任取的非零向量 v0 应使α1 ≠ 0，通常可取v0 =(1,1,…,1)T．如 果初始向量v0 选择不当，以致使α1 = 0，上述结果理论不成立, 但由于计算中的舍入误差，经过若干步以后，可有
vk = Akv0 = b1 x1 + b2 x2 + + bn xn 其中b1 ≠ 0，但由于b1x1一项是由舍入误差得到， |b1|较小，而 vk 的第一项本在(4)式右端各项中占优势, 因此迭代收敛很慢,有 时，选择的v0虽没有使α1 = 0，但α1很小，接近于零时，迭代收 敛也会很慢．具体计算时，如果发现收敛很慢，不妨另取初始 向量v0再算，或者在计算方案中一开始就选择几个不同的v0来 进行试算．
v1
=
u1 10
=
(1,0.8,0.1)T
.
计算结果列成下表，可见到迭代12次主特征值有4 位有效数字的近似值．
k
(vk )T
mk
0
(1,1,1)
−
1
(1,0.8,0.1)
10
2
(1,0.75,−0.111)
7.2
11 (1,0.714346,−0.249790) 6.001675 12 (1,0.714316,−0.249895) 6.000837
vk+1 = λk1+1(α1 x1 + (−1)k+1 α2 x2 + εk+1 )
vk+2 = λk1+2 (α1 x1 + (−1)k+2 α2 x2 + εk+2 )
⇒
(vk+2 )i (vk )i
→
λ21Leabharlann (k→∞)
或
(vk+1 )i (vk−1 )i
→ λ21
(k
→ ∞)
可证明对应于λ1的A的近似特征向量为vk+1 +λ1vk ，对应于λ2的 为vk+1 − λ1vk ．
但上述三种情况也有可能均不出现，此时就需要考虑其它数 值解法，因此乘幂法在某种意义上来说只能用来试算．由于乘 幂法的程序简单，方法简便，在大多数场合下仍然是求按模最 大的特征值和相应特征向量的一个较好的方法．
另外由 vk = λk1(α1 x1 + εk ) 知当 k 充分大时 , 有 vk ≈ α1λk1 x1 ,
对任意的非零向量 v0∈Rn，用A构造向量序列 vk+1 = Avk k = 0,1,2, … (2)
则有
lim (vk )i k→∞ (vk−1 )i
= λ1
i = 1,2,
,n
(3)
式中，(vk)i 表示向量vk 的第i 个分量．
证明 因为 A有n个线性无关的特征向量 xi (i =1, 2,…,n)，
2) 由于 vk+1= Avk ，而vk+1≈ λ1vk ，故有Avk ≈ λ1vk，所以，vk可以 作为与λ1相应的特征向量的近似．
3) 定理的结论式(3)是假定按模最大的特征值是单个的．若 |λ1| > |λ2|的要求不能满足时，应视下列不同情形具体分析：
(1) 当λ1= λ2 = … = λr (即按模最大的特征值是实r 重的)时，且 |λ1| > |λr+1| ≥ … ≥ |λn| 则仿定理的证明，对任意初始向量
n
∑ v0 = αi xi (且 α1, , αr 不全为零)
∑ ∑ 则由幂法有
i =1
vk
=
Akv0
=
r
λk1( αi xi
i =r1
+
i
n = r +1
α
i
⎜⎛ ⎝
λi λ1
⎟⎞k ⎠
xi )
∑ = λk1( αi xi + εk )
∑ ⇒
(vk )i (vk−1 )i
→
λ1
(k
→
∞)
i =1
且有
lim
k→∞
vk λk1
=
r
αi xi
i =1
所以，vk可以作为与λ1相应的特征向量的近似．
(2) |λ1| = |λ2|，但 λ1= −λ2时，且|λ1| > |λ3| ≥ … ≥ |λn| 则对任意初 始向量 v0 = α1 x1 + α2 x2 + + αn xn 假定 α1 ≠ 0 , α2 ≠ 0 则由幂法有 vk = Akv0 = λk1(α1 x1 + (−1)k α2 x2 + εk )