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

规范化
A2 v0 u2 A2 v0 u2 = Av1 = v2 = = max( Av0 ) max(u2 ) max( A2v0 ) uk = Av k −1 Ak v0 uk Ak v0 = = vk = k −1 k max(uk ) max( A v0 ) max( A v0 )
λi k λ (α1 x1 + ∑ α i ( ) xi ) k A v0 λ1 i =2 = 则有 vk = n k λi k max( A v0 ) k max(λ 1 (α1 x1 + ∑ α i ( ) xi )) λ1 i =2 n λi k α1 x1 + ∑ α i ( ) xi x1 λ1 i =2 = → k→∞ n λi k max( x1 ) max(α1 x1 + ∑ α i ( ) xi ) λ1 i =2
$
1) 任取的非零向量 v0 应使α1 ≠ 0，通常可取v0 =(1,1,…,1)T．如 果初始向量v0 选择不当，以致使α1 = 0，上述结果理论不成立, 但由于计算中的舍入误差，经过若干步以后，可有 v k = Ak v0 = b1 x1 + b2 x2 + + bn xn 其中b1 ≠ 0，但由于b1x1一项是由舍入误差得到， |b1|较小，而 vk 的第一项本在(4)式右端各项中占优势, 因此迭代收敛很慢,有 时，选择的v0虽没有使α1 = 0，但α1很小，接近于零时，迭代收 敛也会很慢．具体计算时，如果发现收敛很慢，不妨另取初始 向量v0再算，或者在计算方案中一开始就选择几个不同的v0来 进行试算． 2) 由于 vk+1= Avk ，而vk+1≈ λ1vk ，故有Avk ≈ λ1vk，所以，vk可以 作为与λ1相应的特征向量的近似．
定理1 设矩阵 A∈Rn×n 有n个线性无关的特征向量 xi (i =1, 2,…,n)，其对应的特征值 λi (i =1, 2,…, n)满足 |λ1| > |λ2| ≥ |λ3| ≥ … ≥ |λn| (1) 对任意的非零向量 v0∈Rn，用A构造向量序列 vk+1 = Avk k = 0,1,2, … (2) 则有
3) 定理的结论式(3)是假定按模最大的特征值是单个的．若 |λ1| > |λ2|的要求不能满足时，应视下列不同情形具体分析： (1) 当λ1= λ2 = … = λr (即按模最大的特征值是实r 重的)时，且 |λ1| > |λr+1| ≥ … ≥ |λn| 则仿定理的证明，对任意初始向量
v0 = ∑ α i x i ( 且 α1 ,
λ2 = λ 1[1 + O( )] → λ 1 λ1
于是有
k
k→∞
定理2 设矩阵 A∈Rn×n 有n个线性无关的特征向量， 其按模最大的特征值 λ1 满足 |λ1| > |λ2| ≥ |λ3| ≥ … ≥ |λn| 则对任意的非零向量 v0 (α1 ≠ 0)，按下列方法构造向量序列
⎧ u0 = v0 ⎪ uk = Av k −1 (4.3) ⎨ m = max(u ) k ⎪ k k = 1,2, ⎩vk = uk / mk
4.3 幂迭代法和逆幂迭代法
4.3.1 幂迭代法
在实际问题中，具体要求也有不同．某些应用场合通常不 需要知道矩阵的全部特征值和特征向量，而仅要求得到矩阵 的按模最大的特征值（或称为矩阵的主特征值）和相应的特 征向量，如线性方程组迭代解法的收敛性问题．有的则要求 全部特征值和特征向量．根据这2种不同要求，矩阵的特征值 和特征向量的计算方法也大体上分成2种类型．本章就此2种 类型介绍其中最常用的2 种方法 ⎯ 幂法和QR法． 幂法主要用于求矩阵的按模最大的特征值和相应的特征向 量．它是通过迭代产生向量序列，由此计算特征值和特征向 量．其算法的思想基于如下结论：
k 另外由 v k = λk ( α x + ε ) 知当 k 充分大时 , 有 v ≈ α λ 1 1 1 1 1 x1 , k k
于是当 k → ∞ 时 , 有 : (1) 当 λ 1 > 1 时 , (v k ) i → ∞ ;
( 2) 当 λ 1 < 1 时 , ( v k ) i → 0 . 分量的变化太大，使得计算过程可能产生溢出．为了克服这 一不足，通常采用规范化迭代方式．具体做法为：把迭代向量 uk的最大分量化为1，于是得乘幂法的计算步骤如下： (1) 任取一个初始非零向量 v0 .
k 则由幂法有 v k = Ak v0 = λk ( α x + ( − 1 ) α 2 x2 + ε k ) 1 1 1 k +1 +1 v k +1 = λk α x + − α 2 x 2 + ε k +1 ) ( ( 1 ) 1 1 1 k+2 +2 v k + 2 = λk α x + − α 2 x2 + ε k + 2 ) ( ( 1 ) 1 1 1
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的按模最大的特征值和相应的特征向量． 但上述三种情况也有可能均不出现，此时就需要考虑其它数 值解法，因此乘幂法在某种意义上来说只能用来试算．由于乘 幂法的程序简单，方法简便，在大多数场合下仍然是求按模最 大的特征值和相应特征向量的一个较好的方法．
k +1 k +2 ≈ [λk − （ λ + λ ） λ + λ λ λ 1 1 2 1 1 2 1 ]α1 x1 +2 k +1 k + [λk − （ λ + λ ） λ + λ λ λ 2 1 2 2 1 2 2 ]α 2 x 2 = 0
说明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．
i =1 n
, α r 不全为零 )
则由幂法有 v k = i + ε k )
i =1
⎛ λi ⎞ ⎟ xi ) ∑ α i xi + ∑ α i ⎜ ⎝ λ1 ⎠ i =1 i = r +1
r
r
n
k
所以，vk可以作为与λ1相应的特征向量的近似．
⎧ u0 = v0 ⎪ uk = Av k −1 (4.3) (2) 构造迭代序列 ⎨ mk = max(uk ) ⎪ k = 1,2, ⎩vk = uk / mk
式中用记号max(u)表示向量u 的绝对值最大的分量，从而有
迭代
u1 = Au0 = Av0
u1 Av0 v1 = = max(u1 ) max( Av0 )
则有
k →∞
lim max(uk ) = λ 1
x1 lim vk = k →∞ max( x1 )
例4.4
⎛ − 4 14 0 ⎞ ⎜ ⎟ A = ⎜ − 5 13 0 ⎟ ⎜ − 1 0 2⎟ ⎠ ⎝
有特征值6,3,2。选
v0 =(1,1,1)T，按(4.3)计算，有
u1 = Av0 = (10,8,1)T , m1 = max(u1 ) = 10, u1 v1 = = (1,0.8,0.1)T . 10
k 1 n
由于
λi < 1 ( i = 2, 3, λ1
, n) , 故对足够大的 k , 有
k −1 vk = λk ( α x + ε ) v = λ ( 5) 1 1 1 k k −1 1 ( α1 x1 + ε k −1 )
式中 ε k −1 , ε k ∈ R n , 且当 k → ∞ 时 , ε k −1 , ε k → 0 . 当( x1 ) i ≠ 0 时有 (v k )i α1 ( x1 ) i + (ε k )i lim = λ1 = λ 1 lim k → ∞ ( v k −1 ) i k → ∞ α1 ( x1 ) i + ( ε k −1 ) i 这种由已知非零向量 v0 及矩阵的乘幂 Ak 构造向量序列{vk} 以计算 A的按模最大的特征值λ1的方法称为乘幂法，简称为幂 法．而定理的证明过程事实上已给出了幂法的实施步骤．
(v k + 2 )i ( v k +1 ) i 2 ⇒ → λ1 (k → ∞) 或 → λ2 1 (k → ∞ ) (vk )i ( v k −1 ) i 可证明对应于λ1的A的近似特征向量为vk+1 +λ1vk ，对应于λ2的 为vk+1 − λ1vk ． (3) |λ1| = |λ2|，但 λ1=ρe iθ，λ2 =ρe −iθ，即λ1 , λ2为一对共轭复特 征值，且|λ1| > |λ3| ≥ … ≥ |λn|，则对任意初始向量
k 1 n
这说明规范迭代向量序列收敛到按模最大特征值对应的特征向 量． 同理，可得到
n ⎡ ⎤ λi k k λ 1 (α1 x1 + ∑ α i ( ) xi ) ⎥ ⎢ λ i=2 1 ⎥ mk = max(uk ) = max ⎢ n ⎢ max(λk −1 (α x + α ( λ i )k −1 x )) ⎥ ∑ 1 1 1 i i ⎥ ⎢ λ ⎣ ⎦ i=2 1 n λi k λ 1 max(α1 x1 + ∑ α i ( ) xi ) λ1 i=2 = n λ i k −1 max(α1 x1 + ∑ α i ( ) xi ) λ1 i=2