迭代改善法

212 2 1 22 2
9
定理25
设有初等反射阵 H I 2wwT , 其中
wT w 1, 则：
(1) H 是对称矩阵，即 H T H .
(2) H 是正交矩阵，即 H 1 H .
说明 x2 就是 Ax b 的精确解.
2
但在实际计算中，由于有舍入误差，x2 只是方程组的 近似解，重复(6.14)，(6.15)的过程，就产生一近似解序列 ，有时可能得到比较好的近似. { xk } 算法5 (迭代改善法) 设 Ax b，其中 A R nn 为非
奇异矩阵，且 Ax b为病态方程组(但不过分病态)，用选 主元分解法得到 PA 及计算解 x1 . LU
5.6.2
迭代改善法
设 Ax b，其中 A R nn 为非奇异矩阵，且为病态 方程组(但不过分病态).
本节研究的问题是，当求得方程组的近似解 x1后，
如何对其精度进行改善.
首先用选主元三角分解法实现分解计算
PA LU ,
其中 P为置换阵， L为单位下三角阵， U为上三角阵， 且求得计算解 x1 .
(1) 计算 rk b Axk (用原始 A 及双精度计算)
(2) 求解 LUdk P rk , (3) (4) 如果 d k

即
Ly P rk , (用单精度) Udk y.
/ xk

10t 则输出 k , xk , rk , 停机
改善 xk 1 xk d k (用单精度计算)

i1

n

1/ 2
在计算 时，为了避免上溢或下溢，将 x
d x

,
x
x ( 设d 0)， d
16
则有 H 使 H x e1 , 其中
H I ( ) 1 u u T , 2 / d , u u / d , / d , H H .
2
显然
15

1 u 2
2 2

1 2 2 (( x1 ) 2 x2 x2 ) ( x1 ). 2
如果 和 x1 异号，那么计算 x1 时有可能出现两相
近数相减的情况，有效数字可能损失. 取 和 x1 有相同的符号，即取
sgn( x1 ) x 2 sgn( x1 ) xi2 .
H I 1uu T ,
sgn( x1 ) x 2 ,

1 u 2
2 2
u x e1 ,
( x1 ).
14
证明
x
2
记 y e1 , 设 x y , 取 x
H I 2wwT ,
2
，则有
y 2 , 于是由定理22存在 H变换
2 设 x, y R ,
y1 cos sin y 2
是平面上向量的一个旋转变换，其中
cos P( ) sin sin cos
为正交矩阵.
19
R n中变换： y Px
其中 x ( x1 , x2 ,, xn )T , y ( y1 , y2 , , yn )T , 而
本算法用迭代改善法提高近似解 x1 的精度.
设计算机字长为 t ，用数组 A(n, n)保存 A元素，数组
C (n, n)保存三角矩阵 L 及 U，用 Ip(n) 记录行交换信息，
3
x(n) 存贮 x1 及 xk , r ( n) 保存 rk 或 d k
.
1. 用选主元三角分解实行分解计算 PA LU 且求计 算解 x1 (用单精度) 2. 对于 k 1,2, , N 0
a ,
li
ali , alj
alj
c s
s c
(l 1,2, , m).
利用平面旋转变换，可使向量 x 中的指定元素变为零.
22
设 x ( x1 , xi , x j , xn )T , 其中 xi , x j不全为零，则可选择平面旋转阵 P(i, j , ), 使
( wT x 0)
11
对于 y S ，Hy ( I 2wwT ) y y 2wwT y y. 从而对任意向量 v R n ，总有
Hv x y v,
其中 v为 v 关于平面S的镜面反射(见图5-1).
图5-1
12
定理26
设 x, y为两个不相等的 n维向量， x 2 y 2 ,
P P(i, j , ) P(i, j ) 称为平面旋转矩阵.
显然， P(i, j , )具有性质： (1) P 与单位阵 I 只是在 (i, i ), (i, j ), ( j , i ), ( j , j ) 位置 元素不一样，其他相同. (2) P 为正交矩阵 ( P 1 PT ). (3) P(i, j ) A (左乘)只需计算第 i 行与第 j 行元素， 即对 A (aij ) mn
i 1 1 cos P P(i, j , ) sin j 1 1
sin 1 1 cos
i
j
20
称为 R n中平面 {xi , x j }的旋转变换(或称为吉文斯（Givens） 变换)，
其中 w x e1 , 使 Hx y e1. x e1 2 记 u x e1 (u1 , u2 ,, un )T . 于是
H I 2 uu T u
2 2
I 1uu T ,
其中 u ( x1 , x2 , , xn )T , 1 u 2 . 2
1
现利用 x1 的剩余向量来提高 x1 的精度.
计算剩余向量
r1 b Ax1 ,
（6.14）
求解 Ad r1,
得到的解记为 d1 . 然后改善
x2 x1 d1.
（6.15）
如果(6.14)，(6.15)及解 Ad r1 的计算没有误差，则
Ax2 A( x1 d1 ) Ax1 Ad1 Ax1 r 1 b,
计算解 x1 (1.2560,1.2121)T .
应用迭代改善法需要用原始矩阵 A且用双倍字长精度 计算剩余向量 r b Ax ，其他计算用单精度.
6
计算结果如下表.
x1 r1 1.2560 5.7 10 7 1.2121 3.3715 10 5 d1 x2 r2 0.03220 1.2238 1.18 10 6 0.033502 1.2456 9 10 7 d2 2.285 10 4 2.365 10 4
(3) sgn( u1 )
u
i1
n
2 i
17
(4) u1 u1 (5) u1
(6) d
关于 HA 的计算，设 A (a1 , a2 , , an ) R mn , 其中
ai 为 A的第 i 列向量，则 HA ( Ha1 Ha 2 Ha n ),
x3 (1.2240,1.2454)T ,
r3 (0.68210 5 , 0.65910 5 )T , d 3 (0.271710 4 , 0.351510 4 )T .
如果 xk 需要更多的数位，迭代可以继续.
7
5.7
矩阵的正交三角化及应用
8
5.7.1
定义9
21
c ail a jl s
s ail c a jl
(l 1,2, , n).
其中 c cos , s sin .
(4) AP(i, j )(右乘)只需计算第 i 列与第 j 列元素
4
3.
输出迭代改善方法迭代 N 0次失败信息
当 Ax b 不是过分病态时，迭代改善法是比较好的改进
近似解精度的一种方法，当 Ax b 非常病态时，{xk }可能
不收敛. 例11 用迭代改善法解
1.0303 0.99030 0.99030 x1 2.4944 0.95285 x2 2.3988
x y 则得到一个初等反射阵 ， x y 2
T
则存在一个初等反射阵 H， 使 Hx y. 证明 令 w
H I 2ww I 2
x y x y
T T ( x y ), 2 2
而且
Hx x 2 x y x y
T T ( x y )x 2 2
x2
( x y )( xT x y T x) x y
因此计算 HA 就是要计算
Hai ( I 1uu T )ai ai ( 1u T ai )u
(i 1,2,, n).
于是计算 Hai 只需计算两向量的数量积和两向量的加法.
计算 HA只需作 2nm 次乘法运算.
18
5.7.2
平面旋转矩阵
则变换
sin x1 , 或 y Px cos x2
(3) 设 A为对称矩阵，那么 A1 H 1 AH HAH 亦是对 称矩阵. 证明 只证 H的正交性，其他都可通过验证得到.
H T H H 2 ( I 2wwT )( I 2wwT )
I 4wwT 4w( wT w) wT I .
10
设向量 u 0 ,
初等反射阵
设向量 w R n 且 wT w 1 ，
H ( w) I 2wwT
为初等反射阵(或称为豪斯霍尔德（Householder）变换). 如果记 w (1 , 2 , , n )T , 则
2 1 21 221 H ( w) 2 n 1
i j Px ( x1 , xi, 0, xn )T
定理28 (约化定理)
2 2
.
13
因为
x y
2 T T T ( x y ) ( x y ) 2 ( x x y x), 2
所以
Hx x ( x y ) y.
w是使 Hx y成立的唯一长度等于1的向量(不计符号).
定理27 (约化定理) 设 x ( x1 , x2 , , xn )T 0 ， 则存在初等反射阵 H , 使 Hx e1， 其中
H I 2 uu T u
2 2
是一个初等反射阵. 初等反射阵的几何意义. 考虑以 w为法向量且过原点 O的超平面 S : wT x 0 . 设任意向量 v R n ， 则 v x y , 其中 x S , y S . 于是
Hx ( I 2wwT ) x
x 2wwT x x.
算法6
设 x ( x1 , x2 , , xn )T 0 ，本算法计算 u ,
及 , 使 ( I 1uuT ) x e1 , u的分量冲掉 x的分量.
(1) 计算d d max xi
1in
(2) xi ui xi / d (i 1,2,,n)
(记为Ax b) (这里 10, t 5 ，用5位浮点数运算).
5
解 第4位).
精确解 x* (1.2240,1.2454)T (舍入到小数后
cond ( A) A

A1

22000 4000.
首先实现分解计算 PA LU , 且求 x1
1 A 0.9118 0 1.0303 1 0 0.99030 LU , 0.00099