迭代改善法

例11 用迭代改善法解
1.0303 0.99030
0.99030 0.95285
x1 x2
2.4944 2.3988
(记为Ax b) (这里 10,，t 用55位浮点数运算).
5
解 精确解 x* (1.2240,1(.2舍4入54到)T小数后 第4位).
cond ( A)
A
A1
22000 4000.
如果 需xk要更多的数位，迭代可以继续.
7
5.7 矩阵的正交三角化及应用
8
5.7.1 初等反射阵
定义9
设向量 w R且n w，T w 1
H (w) I 2wwT
为初等反射阵(或称为豪斯霍尔德（Householder）变换).
如果记 w (1,,2则, ,n )T
1 212 212
r1 5.7 107 3.3715 105
d1
x2
r2
0.03220 1.2238 1.18106
0.033502 1.2456 9107
d2 2.285 104 2.365104
x3 (1.2240,1.2454)T , r3 (0.682105 , 0.659105 )T , d3 (0.2717104 , 0.3515104 )T .
i
j
Px (x1, xi, 0, xn )T
其中 xi
xi2
x
2 j
,
arctan(x j / xi ).
证明 取 c cos xi / xi, s .sin x j / xi 由 P(i, j, )x x (x1, x2 , ,, xi, , xj , , xn )T
利用矩阵乘法，显然有
H
I
2
uu u
T 2
I 1uuT ,
2
其中
u (x1 , x2 ,
, xn )T ,
1 2
u
2.
2
显然
则有
15
1
2
u
2 2
1 2
((
x1
)
2
x22
x22 )
( x1).
如果 和 异x号1 ，那么计算 时x1有可能出现两相
近数相减的情况，有效数字可能损失.
取 和 有x1 相同的符号，即取
1
0 0
a(2) 12
a(2) 22
a(2) m2
a(2) 1n
a(2) 2n
0
1
r2 c2
a(2) mn
B2 D2
,
其中 c2
(a2(22) ,
,
a(2) m2
)T
R m1, D2 R . (m1)(n2)
26
(2) 第 步k约化：
设已完成对 上A述第1步…第 k步的1约化，再进行第 k步约化. 即存在初等反射阵 H1, H2 , , Hk1 使
设任意向量 v R，n 则 v x, y 其中 x S, y. S
于是
Hx (I 2wwT )x
x 2wwT x x. (wT x 0)
11
对于 y ，S Hy (I 2wwT ) y y 2wwT y y. 从而对任意向量 v ，Rn 总有
Hv x y v,
其中 v为 关v 于平面S的镜面反射(见图5-1).
i
j
1
1
cos
sin
i
1
P P(i, j, )
sin
1 cos
j
1
ห้องสมุดไป่ตู้
1
20
称为 R中n 平面 {xi的, x旋j}转变换(或称为吉文斯（Givens） 变换)，
P P(i, j, ) P(i, j) 称为平面旋转矩阵.
显然， P(i具, j有,性) 质：
(1) 与P单位阵 只是I 在 元素不一样，其他相同.
证明 只证 H的正交性，其他都可通过验证得到.
H T H H 2 (I 2wwT )(I 2wwT ) I 4wwT 4w(wT w)wT I.
10
设向量
u, 0
uu T H I 2 u 2
2
是一个初等反射阵.
初等反射阵的几何意义.
考虑以 为w法向量且过原点 的O超平面 S : w. T x 0
首先实现分解计算 PA,L且U求
x1
1
01.0303 0.99030
A 0.9118 1 0
0.00099 LU ,
计算解 x1 (1.2560,1.2121)T . 应用迭代改善法需要用原始矩阵 且A用双倍字长精度
计算剩余向量r b ，A其x他计算用单精度.
6
计算结果如下表.
x1 1.2560 1.2121
sgn(
x1 )
x
sgn(
2
x1 )
n i1
xi2
1/
2
.
在计算 时，为了避免上溢或下溢，将 x
d x ,
x x (设d 0)，
d
16
则有 H使 H x , e1 其中
H I ( )1uuT , / d , u u / d , / d 2 ,
H H .
图5-1
12
定理26 设 x,为y 两个不相等的 维n向量， x y ,
2
2
则存在一个初等反射阵 H， 使 Hx y.
证明
令 w x， y 则得到一个初等反射阵
x y 2
而且
H I 2wwT
I 2
x y x y 2
( xT
yT
),
2
Hx x 2
x y x y 2
(xT
yT
)x
2
x
2
(
x
y)( x
18
5.7.2 平面旋转矩阵
设 x, y, 则R变2 换
y1 y2
cos sin
sin cos
x1 x2
,
或 y Px
是平面上向量的一个旋转变换，其中
P(
)
cos sin
sin cos
为正交矩阵.
19
中R变n 换： y Px
其中 x (x1, x2 , , xn )T , y ( y1, y2 , , yn )T , 而
H k 1 , H 2 H1 A A(k ) ,
其中
1
a(2) 12
2
A(k )
a(2) 1k 1
a (3) 2 k 1
a(2) 1k
a (3) 2k
k1
a (k1) k 1k
a(k) kk
a(k) mk
a(2) 1n
a (3) 2n
a (k1) k 1n a(k) kn
(i, i), (i, j位), (置j, i), ( j, j)
(2) 为P正交矩阵 (P1 PT ).
(3) P((i左, j乘) A)只需计算第 行与第i 行元素，j
即对 A (aij )mn
21
ail ajl
c s
s c
ail a jl
其中 c cos , s sin .
(l 1,2, , n).
C(n, n)保存三角矩阵 L及 ，U用 Ip记(n录) 行交换信息，
3
x(n)存贮 x及1 ,xk r保(n存) 或 r.k d k 1. 用选主元三角分解实行分解计算 PA且求LU计
算解 x(1用单精度)
2. 对于
k 1,2, , N0
(1) 计算 rk (用b 原 A始xk及双精度A计算)
{x，k }有时可能得到比较好的近似. 算法5 (迭代改善法) 设 Ax ，b其中
A为 非R nn
奇异矩阵，且 Ax 为b病态方程组(但不过分病态)，用选 主元分解法得到 PA 及L计U算解 . x1
本算法用迭代改善法提高近似解 的x1精度. 设计算机字长为 ，t用数组 A保(n存, n)元素，A数组
算法6 设 x (x1, x2 , ,，xn )T 0 本算法计算 u,
及 , 使 (I 1uuT )x e1, u的分量冲掉 x的分量.
(1) 计算d
d max 1in
xi
(2) xi ui xi /d (i 1,2, ,n)
n
(3) sgn( u1)
ui2
i1
17
(4) u1 u1 (5) u1 (6) d
关于 的HA计算， 设 A (a1, a2 , , a, n ) R mn 其中 ai为 A的第 列i 向量， 则
HA (Ha1 Ha2
因此计算 H就A 是要计算
Han ),
Hai (I 1uuT )ai ai ( 1uT ai )u (i 1,2, ,n). 于是计算 H只a需i 计算两向量的数量积和两向量的加法. 计算 H只A 需作 2次nm乘法运算.
H k I k 1
k 1
H k
m
k
, 1
28
第 k步约化为
H k A(k ) H k H k1 H 2 H1 A A(k1)
I
k 1
0
0 H k
Rk 0
rk ck
Bk Dk
Rk 0
rk H kck
Bk H k Dk
,
其中 A(k左1) 上角的 阶k子矩阵为 阶上k三角阵
这就使 A的三角化过程前进了一步.
(2) 求解 LUd, k即Prk
(用单LU精yd度k )P
rk , y.
(3) 如果
dk / 则xk输出 10t , 停机 k, xk , rk
(4) 改善
xk1 (用x单k 精d度k 计算)
4
3. 输出迭代改善方法迭代 次失N0败信息 当 Ax不是b过分病态时，迭代改善法是比较好的改进 近似解精度的一种方法，当 Ax 非b常病态时， 可{x能k } 不收敛.
(4) AP(右(i乘, j))只需计算第 列与第i 列元素 j
ali ,
alj ali ,
alj
c s
s c
(l 1,2, , m).
利用平面旋转变换，可使向量 x中的指定元素变为零.
22
定理28 (约化定理) 设 x (x1, xi ,, x j , xn )T
其中 xi ,不x j全为零， 则可选择平面旋转阵 P(i, j,, ) 使
求解 Ad , 得r1 到的解记为 . d1 然后改善
x2 x1 d1.
（6.14） （6.15）
如果(6.14)，(6.15)及解 A的d 计 算r1 没有误差，则
Ax2 A(x1 d1) Ax1 Ad1
Ax1 r1 b,
说明 x就2 是 Ax的精b 确解.
2
但在实际计算中，由于有舍入误差， 只x2是方程组的 近似解，重复(6.14)，(6.15)的过程，就产生一近似解序列
a(k) mn
27
Rk 0
rk ck
Bk Dk
mkk11,
其中 R为k k阶1上三角阵， ck R mk1, Dk R (mk1)(nk ) . 不妨设 ck， 0 否则这一步不需要约化 (如果 A列满秩，
则 ck ).0 于是，可选取初等反射阵 Hk I k1ukukT , 使
Hkck k e1. 令
xT x y
2
yT
x)
.
2
13
因为 所以
x y
2 2
(x y)T
(x y)
2(xT x yT x),
Hx x (x y) y.
是w使 H成x 立 的y 唯一长度等于1的向量(不计符号). 定理27 (约化定理) 设 x (x1, x2 , ,，xn )T 0 则存在初等反射阵 H, 使 Hx ，e1 其中
xi
cxi
sx
j
,
xj sxi cx j ,
xk xk (k i, j).
23
由 c,的s 取法得
xi
xi2
x
2 j
,
xj 0.
24
5.7.3 矩阵的QR分解
设 A且R为m非n 零矩阵，则存在初等反射阵 H1,H2 , ,H s , 使
H s H2 H1 A A(s1)（上梯形）.
设
a11 a12
A
A(1)
a21 am1
a22
am2
a1n
a2n a mn
(a1
a2
an )
(按列分块).
25
(1) 第1步约化：
如果 a1, 取 0 ，H1 I 即这一步不需要约化， 不妨设 a1 ，0 于是可选取初等反射阵 H1 I ,11u1u1T 使
H1a1 1e1. 于是 H1 A(1) (H1a1 H1a2 H1an )
H I 1uuT ,
sgn(
x1 )
x
,
2
u xe1,
1 2
u
2 (
2
x1).
14
证明
记 y ,e1
设 x , y
取 ，x 2
x 2
y
,
2
于是由定理22存在H变换
H I 2wwT ,
其中 w x ,使e1 x e1 2
Hx y e1.
记 u x e1 (u1, u2 , , un )T . 于是
. Rk 1
令 s min( m，继1续, n上) 述过程，
H
(w)
221
1 222
2n1 2n2
21n
22n
.
1 2n2
9
定理25 设有初等反射阵 H I 2, wwT 其中
wT w1, 则：
(1) 是H对称矩阵，即 H T H .
(2) 是H正交矩阵，即 H 1 H .
(3) 设 为A对称矩阵，那么 称矩阵.
A1 H 1亦AH是对 HAH
5.6.2 迭代改善法
设 Ax，其b中
A为非R奇n异n 矩阵，且为病态
方程组(但不过分病态).
本节研究的问题是，当求得方程组的近似解 后x1， 如何对其精度进行改善.
首先用选主元三角分解法实现分解计算
PA LU ,
其中 P为置换阵， 为L 单位下三角阵， 为U上三角阵， 且求得计算解 x1.
1
现利用 的x1剩余向量来提高 的精x1 度. 计算剩余向量 r1 b Ax1,