计算方法第七章(r)

(1e1,
a(2) 2
,
, an )
,
a(2) n
)
1
a(2) 12
a(2) 22
a(2) n2
a(2) 12
a(2) 2n
a(2) nn
1
a(2) 2 A1
对矩阵 A1 又再重复前面的过程，即求出Householder矩阵 H2
H2 (a2(22) ,
,
a(2) n2
)
e(2
21
mk 2mk1 mk2
M k 1
1.3 反幂法
如果 A 非奇异，用其逆矩阵代替 A 进行乘幂法，称为反幂 法。
逆矩阵的特征值与A 互为倒数。即为：
1/ n 1/ n1 1/ 1
用 A 的逆进行乘幂法，得到 1/ n
迭代格式为：
yk A1zk1 (k 1, 2, , z0任意给定) mk max( yk ) ( yk中绝对值最大的分量),
（尤其，对于对称矩阵可以化为三对角矩阵）
（2）利用带原点位移的QR 方法构造矩阵序列 Ak
（3）对矩阵
Ak
取加速因子
a(k nn
)
进行加速
（4）判断矩阵 Ak 的最后一行非对角元素（由于是拟上三
选取 ，使得 Ak (ai(jk) )nn
,
a(k ) pq
0
第 k 步迭代矩阵的元素为：
a(k) pj
a(k 1) pj
cos
a(k 1) qj
sin
a(k) jp
a(k) qj
a
(k pj
1)
sin
a(k 1) qj
cos
a(k) jq
( j p, q)
a(k) pp
a(k 1) pp
2 1
称为收敛率
由于 zk 线性收敛于 x1 ，于是可以对之进行埃特金加速，
Wi
(zk )i (zk2 )i (zk1)i2 (zk )i 2(zk1)i (zk2 )i
(i 1, 2, , n)
以上是计算特征向量的埃特金加速，同样可以得到关于计算特 征值的埃特金加速，
Mk
mk mk2
m2 k 1
A QR
3.2 基本 QR 方法 利用矩阵的 QR 分解，立即可以得到矩阵的一系列相似矩阵
A A1 Q1R1
A2 R1Q1 Q2R2
A3 R2Q2 Q3R3
Ak Rk1Qk1 Qk Rk 其中，Qk 为正交矩阵，Rk 为上三角矩阵，Ak 称为 QR 序列
Ak Rk 1Qk 1 Qk 1 Ak Qk R Q 1 k 1 k 1 Ak 1Qk 1
征值和特征向量。
第三节 QR 分解方法 3.1 QR 分解 设 u 为n维实单位向量，称下面矩阵为Householder矩阵：
H I 2uu
容易验证Householder矩阵为正交阵，同时又是对称阵：
HH I , H H
对任意的向量 x 以及单位向量 g，存在Householder矩阵，使
cos
a(k iq
1)
sin
p(k) iq
p(k ip
1)
sin
p(k 1) iq
cos
p p (k )
(k 1)
ij
ij
( j p, q)
i 1, 2,
,n
最后，雅可比方法的计算步骤可以归纳为： （1）确定非对角绝对值最大元位置（p，q），并计算sin和
cos的值； （2）计算迭代矩阵的元素； （3）计算特征向量； （4）与计算精度进行比较，以决定是否终止计算，并输出特
值，通常按如下步骤计算：
y
a( k 1) pp
a( k 1) qq
,
x
sign(a(pkp1)
a( k 1) qq
)
2a(pkq1)
cos 2
tg2 x / y
y x2 y2
cos
1 cos 2
2
sin 2 x
x2 y2
sin sin 2 2 cos
实际计算中，一般预先给一个计算精度 ，当第 m 步满足
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
1.2 加速技术： 显然，乘幂法的收敛速度依赖
2 1
，如此比值接近1，则收敛
速度会很慢。
用 A－ pI 代替A，进行乘幂法。迭代速度可能会大大加快。
cos2
2a
(k 1) pq
sin
cos
a(k 1) qq
sin2
a(k) qq
a(k 1) pp
sin2
2a(pkq1)
sin
cos
a(k 1) qq
cos2
a(k) pq
( aq( kq 1)
a
(k 1) pp
)
sin
cos
a(k 1) pq
(cos
2
sin2 )
a(k) qp
a a (k )
R1
R1
A(1) 22
R1
于是，取R为Householder矩阵，则 R1c1 化为除第一个分量外
其余分量为零的向量。
如此下去，可以将矩阵 A 化为拟上三角矩阵。
利用前面带原点位移的QR 方法的性质，可以看出，用拟上三 角矩阵进行原点位移的QR 方法计算是非常方便的。
QR 方法的总结： （1）利用Householder矩阵，将矩阵 A 相似于拟上三角矩阵
Lx zk1
Ryk
x
1/ n
lim
k
mk
lim
k
zk
xn max(xn )
反幂法的一个主要应用是已知矩阵的近似特征值后，求其特征 向量。
如果已求得矩阵某个特征值 m 的近似值 m ，则
0 m m i m (i m) 于是，用反幂法可以求出 A m I 的按模最小特征值及相应
的特征向量。此时，迭代为：
zk zk1 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
如此，产生一个新的阵，然后再重复上面的步骤，直到最后将
A 化为对角矩阵，则对角元就是所要求的特征值！
将上述过程数学化，首先，记 A A0 ，则
A1 R1 A0R1 , A2 R2 A1R2 ,
,
Ak
Rk
Ak
R
1 k
,
Ak 1 (ai(jk 1) )nn
,
ak 1 pq
max 1i jn
a(k 1) ij
Hx x g 2
x x g
u
2
x x g
2
特别，取 g = e = ( 1 , 0 , …… , 0)
n
Hx x 2 e1 ( x 2 , 0, , 0) ( xi2 , 0, , 0) i 1
将矩阵 A 记为
A (a1, a2, , an ) , aj (a1j , a2 j , , anj )
计算，并划掉矩阵的最后一行和最后一列，产生一个子矩
阵
（4）对子矩阵重复进行上面的加速计算
一个新的计算方案：
对于矩阵
a11 a12
A1
a21
a22
an1 an2
取变换
1
U1
R1
a1n
a2n
a11
c1
ann
A(1) 12
A(1) 22
A2
U1A1U1
a11
R1c1
A(1) 12
zk
yk
/ mk
1/ n
lim
k
mk
为避免矩阵的求逆运算，通常也采取如下的算法：
Ayk zk1 (k 1, 2, mk max( yk ) zk yk / mk
, z0任意给定)
每次迭代需要解 Ayk zk1 ，为此，可将 A 进行LR分解，
则每次迭代只需解两个三角方程组
最后求得：
第一节 乘幂法与反幂法 1.1 乘幂法：用于求矩阵的模（绝对值）最大的特征值。
记矩阵 A 的特征值为： 1 2 相应的特征向量为： 1,2 , ,n
任取非零向量 z0 c11 c22
n
cnn ，记 zk Ak z0
则有： zk 1k c11
1
lim
k
( zk 1 )i (zk )i
)
(2 , 0,
, 0)
于是，我们记
1 0
Q2
H
2
则
1
a(2) 12
2
a(3) 13
a(3) 23
Q2
A1
Q2Q1 A
a(3) 33
0
a(3) 3n
a(3) 12
a(3) 2n
a(3) 3n
1
a(2) 12
2
*
*
A2
a(3) nn
如此一直下去，最后得到
Qn1Qn2 Q1 A R 记 Qˆ Qn1Qn2 Q1 ，注意到这是一个正交矩阵，令 Qˆ Q
这叫原点移位法。
1 p 2 p i p (i 3, 4, , n)
2 p 2 1 p 1 求得A pI的按模最大特征值1, A的按模最大特征值1 1 p
埃特金加速： 可以证明：乘幂法线性收敛
mk1 1 2
mk 1
1
[zk1 10 ] i 2
[zk 10 ] i
1
Ak
Q k 1
Ak 1Qk
1
Ak相似于Ak－1 Ak相似于A1 A
最后，可以证明， Ak 的对角线下面的元素（不包括对角线）
收敛于零，由相似关系，不难推出其对角线上的元素收敛
到矩阵 A 的特征值！
最后，要指出的是，当用 QR 方法求出特征值后（准确讲，是 特征值的近似值），我们可以用反幂法求出更加精确的特征 值，更为重要的是可以求出相应的特征向量。
n
a(m) ij
i j1
停止计算，这时，
Am RmRm1 R1AR1 Rm1Rm Pm APm
则对角阵的对角元为特征值近似值，矩阵 P 的列向量为特征向
量近似值。实际计算中，矩阵 P 是按如下步骤计算：
P0 I
,
Pk
Pk
R
1 k
(k 1, 2,
, m)
p(k) ip
p(k 1) ip
于是，可以求得Householder矩阵，将 A 的第一个列向量化简。
x a1 , g s ign(a11)e1
u1
a1 1e1 a1 1e1 2
, 1 s ign(a11) a1 2
, Q1 I 2u1u1
1
a1 1e1
2
[(a1 1e1) (a1 1e1)
2(12 1a11)
Q1A Q1(a1, a2 ,
3.3 带原点位移的QR 方法
为加速收敛，每次选取位移 sk ，作
Ak Ak
1
sk
I Qk Rk RkQk sk
I
(k 1, 2,
)
该矩阵序列有如下性质：
（1） Ak1相似于Ak
（2）如 Ak 为拟上三角，则 Ak 也为拟上三角矩阵（拟上三角
矩阵指的是次对角线下的元素为零的矩阵）
（3）如取位移
( A mk
m ) yk
max(
yk
zk )
1
zk
yk
/ mk
lim
k
zk
xm max(xm )
(k 1, 2, , z0任意给定)
lim
k
mk
m
1
m
通常，初值选为：
z0 Le e (1,1, ,1)
这里，矩阵 L 为 角矩阵。
A m I LR 分解中的单位下三
第二节 对称矩阵的雅可比方法 两个重要的基本性质：
sk
a 为
(k ) nn
，则
Ak
最后一行非对角元二阶收
敛于零（特别对于对称矩阵，能达到三阶收敛），其余次对
角元收敛于零的速度会慢一些。
加速技术下的算法：
（1）确定计算精度10E-m
（2）对矩阵 Ak
取加速因子
a(k) nn
进行加速
（3）判断矩阵 Ak 的最后一行非对角元素是否小于要求的精
度，如果不小于，继续加速迭代，如已经小于精度，停止
（1）如 A 为实对称矩阵，则一定存在正交矩阵 Q ，使之相似 于一个对角矩阵，而该对角矩阵的对角元正是 A 的特征值。
Q Q1 , QQ I , QAQ
（2）一个矩阵左乘一个正交矩阵或右乘一个正交矩阵，其E范
数不变。
nn
A 2 E
ai2j
trace( A A) trace( AQQA)
(k 1)
ij
ij
(i, j p, q)
为使
a(k ) pq
0
，必须
tg 2
2a(pkq1)
a a (k1) pp
(k 1) qq
在这里，我们通常，限制
/4
，如果
a( k 1) pp
a( k 1) qq
，
当
a( k 1) pq
0
时，取 / 4
，当
a( k 1) pq
0
时， / 4
在具体计算第 k 步迭代矩阵的元素时，需要计算正弦值和余弦
这里，( z)i 表示向量的第 i 个分量
具体计算时，对于任意取的初始向量，按以下格式计算：
yk Azk1 (k 1, 2, , z0任意给定) mk max( yk ) ( yk中绝对值最大的分量), zk yk / mk
则有
1
lim
k
mk
if mk mk1 or
lim
k
zk
x1 / max(x1)
QA
2 E
i1 j1
A2
A 2
Q A 2
( AQ) 2
Байду номын сангаас
2
AQ
E
E
E
E
E
下面的矩阵是一个 n 阶正交矩阵：
1
cos
sin
R( p, q, )
1
sin
cos
1
(p)
(q)
2.1 雅可比算法 算法的思想：
设 A 为对称矩阵，选出 A 的除对角元外的所有元素中绝对值 最大的一个，然后用前一页中的正交矩阵将此元化为零。