计算方法第七章(r)

实际计算中，一般预先给一个计算精度 实际计算中，一般预先给一个计算精度 ε ，当第 m 步满足
τ=∑a
停止计算，这时， 停止计算，这时，
i≠ j =1
n
(m) ij
<ε
Τ Τ Am = RmRm−1 LR ARΤ LRm−1Rm = P APΤ ≈ Λ 1 1 m m
则对角阵的对角元为特征值近似值，矩阵 P 的列向量为特征向 则对角阵的对角元为特征值近似值， 量近似值。实际计算中， 是按如下步骤计算： 量近似值。实际计算中，矩阵 P 是按如下步骤计算：
第三节 QR 分解方法 3.1 QR 分解 维实单位向量， 矩阵： 设 u 为n维实单位向量，称下面矩阵为 维实单位向量 称下面矩阵为Householder矩阵： 矩阵
H = I − 2uu
Τ
容易验证Householder矩阵为正交阵，同时又是对称阵： 矩阵为正交阵，同时又是对称阵： 容易验证 矩阵为正交阵
λ1 − p > λ2 − p ≥ λ i − p
(i = 3, 4, L , n)
λ2 − p λ2 < λ1 − p λ1
求得A − pI的按模最大特征值µ1 , A的按模最大特征值λ1 = µ1 + p
埃特金加速： 埃特金加速： 可以证明： 可以证明：乘幂法线性收敛
mk +1 − λ1 λ2 → mk − λ1 λ1 [ zk +1 − ξ10 ] i [ zk − ξ10 ] i
A = A0
，则
Τ A1 = R1 A0 R1Τ , A2 = R2 A1 R2 , L , Ak = Rk Ak −1 RkΤ , L
( Ak −1 = (aijk −1) ) n×n
,
( a k −1 = max aijk −1) pq 1≤ i ≠ j ≤ n
k , a (pq) = 0
( Ak = (aijk ) ) n×n 选取θ ，使得
xm lim zk = k →∞ max( xm )
1 lim mk = % k →∞ λm − λm
通常，初值选为： 通常，初值选为：
z0 = Le e = (1,1, L ,1)
这里， 这里，矩阵 L 为 角矩阵。 角矩阵。
% A − λm I = LR 分解中的单位下三
第二节 对称矩阵的雅可比方法 两个重要的基本性质： 两个重要的基本性质： 为实对称矩阵， （1）如 A 为实对称矩阵，则一定存在正交矩阵 Q ，使之相似 ） 于一个对角矩阵， 的特征值。 于一个对角矩阵，而该对角矩阵的对角元正是 A 的特征值。
λ2 显然，乘幂法的收敛速度依赖 λ ，如此比值接近1，则收敛 显然， 如此比值接近 ， 1
速度会很慢。 速度会很慢。 代替A，进行乘幂法 迭代速度可能会大大加快。 乘幂法。 用 A－ pI 代替 ，进行乘幂法。迭代速度可能会大大加快。 这叫原点移位法。 这叫原点移位法。
1.2 加速技术： 加速技术：
yk = Azk −1 ( k = 1, 2,L , z0任意给定) mk = max( yk ) ( = yk中绝对值最大的分量), z = y / m k k k
则有
λ1 = lim mk
k →∞
lim zk = x1 / max(x1)
k →∞
if
mk ≈ mk −1 or
zk ≈ zk −1
迭代矩阵的元素为： 第 k 步迭代矩阵的元素为：
( k a(k ) = a(k−1) cosθ + aqjk−1) sinθ = a(jp) pj pj ( ( k aqjk ) = −a(k−1) sinθ + aqjk−1) cosθ = a(jq) ( j ≠ p, q) pj (k a(k ) = a(k−1) cos2 θ + 2a(k−1) sinθ cosθ + aqq−1) sin2 θ pp pp pq (k (k aqq) = a(k−1) sin2 θ − 2a(k−1) sinθ cosθ + aqq−1) cos2 θ pp pq ( (k a(k ) = (aqk−1) − a(k−1) )sinθ cosθ + a(k−1) (cos2 θ −sin2 θ ) = aqp) pq q pp pq ( ( aijk ) = aijk−1) (i, j ≠ p, q)
Τ P = I , P = P −1Rk (k =1,2,L, m) 0 k k
p
(k ) ip
=p
(k −1) ip
cosθ + a
(k −1) iq
sinθ
( (k (k piqk ) = − pip −1) sinθ + piq −1) cosθ ( ( pijk ) = pijk−1) ( j ≠ p, q) i =1,2,L, n
第七章
矩阵的
特征值与特征向量
第一节 乘幂法与反幂法 1.1 乘幂法：用于求矩阵的模（绝对值）最大的特征值。 乘幂法：用于求矩阵的模（绝对值）最大的特征值。 特征值为： 记矩阵 A 的特征值为： λ ≥ λ ≥ L ≥ λ 1 2 n 相应的特征向量为： 相应的特征向量为： ξ1 , ξ 2 , L , ξ n 任取非零向量 任取非零向量 则有： 则有：
每次迭代需要解 Ayk = zk −1 ，为此，可将 A 进行 分解， 为此， 进行LR分解 分解， 则每次迭代只需解两个三角方程组
Ayk = zk −1 ( k = 1, 2,L , z0任意给定) mk = max( yk ) z = y / m k k k
最后求得： 最后求得：
最后，雅可比方法的计算步骤可以归纳为： 最后，雅可比方法的计算步骤可以归纳为： ），并 （1）确定非对角绝对值最大元位置（p，q），并计算 和 ）确定非对角绝对值最大元位置（ ， ）， 计算sin和 cos的值； 的值； 的值 （2）计算迭代矩阵的元素； ）计算迭代矩阵的元素； （3）计算特征向量； ）计算特征向量； （4）与计算精度进行比较，以决定是否终止计算，并输出特 ）与计算精度进行比较，以决定是否终止计算， 征值和特征向量。 征值和特征向量。
λ2 → λ1
λ2 λ1
称为收敛率
由于
zk
于是可以对之进行埃特金加速， 线性收敛于 x1 ，于是可以对之进行埃特金加速，
2 k +1 i
( zk )i ( zk + 2 )i − ( z ) Wi = ( zk )i − 2( zk +1 )i + ( zk + 2 )i
(i = 1, 2, L , n)
m
% % 0 < λm − λm << λi − λm
的特征向量。此时，迭代为： 的特征向量。此时，迭代为：
(i ≠ m)
于是， 于是，用反幂法可以求出 A − λ I 的按模最小特征值及相应 % m
% ( A − λm ) yk = zk −1 ( k = 1, 2,L , z0任意给定) mk = max( yk ) z = y / m k k k
Q Τ = Q −1 , QQ Τ = I
数不变。 数不变。
, QAQ Τ = Λ
（2）一个矩阵左乘一个正交矩阵或右乘一个正交矩阵，其E范 ）一个矩阵左乘一个正交矩阵或右乘一个正交矩阵， 范
A
2 E
= ∑∑ a = trace( A A) = trace( A Q QA) = QA
i =1 j =1 2 ij Τ Τ Τ
HHΤ = I , H = HΤ
对任意的向量 x 以及单位向量 g，存在Householder矩阵，使 ，存在 矩阵， 矩阵
Hx = x 2 g
u=
x− x 2 g x− x 2 g
xi2 ,0,L,0) ∑
i=1 n
特别， 特别，取 g = e = ( 1 , 0 , …… , 0)
(p) (q)
2.1 雅可比算法 算法的思想： 算法的思想： 设 A 为对称矩阵，选出 A 的除对角元外的所有元素中绝对值 为对称矩阵， 最大的一个，然后用前一页中的正交矩阵将此元化为零 此元化为零。 最大的一个，然后用前一页中的正交矩阵将此元化为零。 如此，产生一个新的阵，然后再重复上面的步骤， 如此，产生一个新的阵，然后再重复上面的步骤，直到最后将 A 化为对角矩阵，则对角元就是所要求的特征值！ 化为对角矩阵，则对角元就是所要求的特征值！ 将上述过程数学化，首先， 将上述过程数学化，首先，记
为使 a ( k ) = 0 ，必须 pq
tg2θ =
2a(k−1) pq a
(k −1) pp
−a
(k −1) qq
在这里，我们通常， 在这里，我们通常，限制
k (k θ ≤ π / 4 ，如果 a (pp−1) = aqq −1) ，
当 a ( k −1) > 0 时，取 θ = π / 4 ，当 a ( k −1) < 0 时， = −π / 4 θ pq pq 迭代矩阵的元素时， 在具体计算第 k 步迭代矩阵的元素时，需要计算正弦值和余弦 通常按如下步骤计算： 值，通常按如下步骤计算：
stop
例子： 例子：求矩阵的模最大特征值及其特征向量
−1 2 1 A = 2 −4 1 1 1 −6
计算结果 程序
%用乘幂法计算矩阵模最大的特征值和相应的特征向量 A=[-1,2,1;2,-4,1;1,1,-6] z0=[1,1,1]'; errtel=1e-6; er=1; k=0; while er>errtel k=k+1; yk=A*z0; [c,p]=max(abs(yk)); mk=yk(p) zk=yk/mk; er=norm(zk-z0); z0=zk; end k,mk,zk
k (k k (k k y = a (pp−1) − aqq −1) , x = sign(a (pp−1) − aqq −1) ) ⋅ 2a (pq−1)
tg2θ = x / y
cos2θ = sin2θ = y x2 + y2 x x2 + y2
1+ cos2θ ⇒ cosθ = 2
⇒
sin2θ sinθ = 2cosθ
非奇异， 进行乘幂法 称为反幂 如果 A 非奇异，用其逆矩阵代替 A 进行乘幂法，称为反幂 法。 逆矩阵的特征值与A 互为倒数。即为： 逆矩阵的特征值与 互为倒数。即为：
1/ λn ≥ 1/ λn −1 ≥ L ≥ 1/ λ1
用 A 的逆进行乘幂法，得到 1/ λn 的逆进行乘幂法， 迭代格式为： 迭代格式为：
n
n
2 E
A
2 E
= A
Τ 2 E
% = Q A
Τ
Τ
2 E
% = ( AQ)
Τ
2 E
% = AQ
2 E
下面的矩阵是一个 阶正交矩阵： 下面的矩阵是一个 n 阶正交矩阵：
1 O cos θ sin θ O R ( p, q, θ ) = 1 O − sin θ cos θ O 1
以上是计算特征向量的埃特金加速，同样可以得到关于计算特 以上是计算特征向量的埃特金加速，同样可以得到关于计算特 计算特征向量的埃特金加速 征值的埃特金加速， 征值的埃特金加速， 的埃特金加速
mk mk + 2 − mk2+1 Mk = mk − 2mk +1 + mBiblioteka Baidu + 2
M k → λ1
1.3 反幂法
yk = A−1 zk −1 (k = 1, 2,L , z0任意给定) mk = max( yk ) ( = yk中绝对值最大的分量), z = y / m k k k
1/ λn = lim mk
k →∞
为避免矩阵的求逆运算，通常也采取如下的算法： 为避免矩阵的求逆运算，通常也采取如下的算法：
Lx = zk −1 Ryk = x
1/ λn = lim mk
k →∞
xn lim zk = k →∞ max( xn )
反幂法的一个主要应用是已知矩阵的近似特征值后， 反幂法的一个主要应用是已知矩阵的近似特征值后，求其特征 向量。 向量。 % 如果已求得矩阵某个特征值 λm 的近似值 λ ，则
z0 = c1ξ1 + c2ξ 2 + L + cnξ n ，记 zk = Ak z0
zk ≈ λ1k c1ξ1
这里， ( 表示向量的第 这里， z )i 表示向量的第 i 个分量
(zk+1)i λ1 = lim k →∞ (z ) k i
具体计算时，对于任意取的初始向量，按以下格式计算： 具体计算时，对于任意取的初始向量，按以下格式计算：