# Iterative Refinement for Ill-conditioned Linear by yyf11174

VIEWS: 19 PAGES: 4

• pg 1
```									                               2008 International Symposium on Nonlinear Theory and its Applications
NOLTA'08, Budapest, Hungary, September 7-10, 2008

Iterative Reﬁnement for Ill-conditioned Linear Equations.
Shin’ichi Oishi1,4 , Takeshi Ogita2,4 and Siegfried M. Rump3,1
1
Department of Applied Mathematics, Faculty of Science and Engineering,
Waseda University,
3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555 Japan
2
Department of Mathematics, College of Arts and Sciences,
Tokyo Woman’s Christian University,
2-6-1 Zempukuji, Suginami-ku, Tokyo 167-8585, Japan
3
Institute for Reliable Computing,
Hamburg University of Technology
Schwarzenbergstr. 95 21071 Hamburg Germany
4
CREST, JST, Japan

Abstract—This paper treats a linear equation                    provided that the residuals are evaluated by ex-
tended precision, in which the unit round oﬀ u is,
Av = b,                               for example, the order of u2 , before rounding back
where A ∈ Fn×n and b ∈ Fn . Here, F is a set of                   to the working precision. In this paper, we will
ﬂoating point numbers. Let u be the unit round-oﬀ                 treat ill-conditioned problems with
of the working precision and κ(A) = ∥A∥∞ ∥A−1 ∥∞
1 < uκ(A) < ∞.                   (4)
be the condition number of the problem. In this
paper, ill-conditioned problems with
We can assume without loss of generality that for a
1 < uκ(A) < ∞                               certain positive integer k the following is satisﬁed:

are considered and an iterative reﬁnement algo-                                       uk κ(A)     β < 1.               (5)
rithm for the problems is proposed. In this paper,
the forward and backward stability will be shown                     In Ref. [8], Rump has shown that for arbitrary
for this iterative reﬁnement algorithm.                           ill-conditioned matrices A, we can have good ap-
proximate inverses R1:k satisfying ∥R1:k A − I∥∞
1. Introduction                                                   α < 1. Here, R1:k is obtained as

In this paper, we will consider the convergence                              R1:k = R1 + R2 + · · · + Rk            (6)
of an iterative reﬁnement for a linear equation
with Ri ∈ Fn×n and I is the n-dimensional unit
Av = b,                        (1)    matrix. In Ref. [6], we have partially clariﬁed the
mechanism of the convergence of Rump’s method.
where A ∈ Fn×n and b ∈ Fn . Here, F is a set of                     Let A, B, C ∈ Fn×n . We assume that we have
ﬂoating point numbers. Let u be the unit round-oﬀ                 an accurate matrix product calculation algorithm
of the working precision and κ(A) = ∥A∥∞ ∥A−1 ∥∞                  [AB − C]k such that
be the condition number of the problem. Here,
∥ · ∥∞ is the maximum norm deﬁned by                                                 D1:k = [AB − C]k                  (7)
∥v∥∞ = max |vi |, for v = (v1 , v2 , · · · , vn ) ∈ R T      n
1 i n                                                    satisfying
(2)
and                                                                   ∑
k
Di − (AB − C)              cuk ∥AB − C∥∞ . (8)
∑
n
i=1                     ∞
∥A∥∞ = max                |Aij | for A = (Aij ) ∈ Rn×n .
1 i n
j=1                                          Here, D1:k is deﬁned as
(3)
The superscript T denotes the transpose and R is                                D1:k = D1 + D2 + · · · + Dk            (9)
the set of real numbers. For well posed problems,
i.e., in case of uκ(A) < 1, it has been shown [1]-                with Di ∈ Fn×n . Such algorithms have been pro-
[5] that the iterative reﬁnement improves the for-                posed by the present authors with c < 2.1 (See [7],
ward and backward errors of computed solutions                    [9] and [10]).

- 516 -
Now we propose the following iterative reﬁne-          Theorem 1 Let vn be generated from (20) with
ment algorithm:                                         any starting vector v0 ∈ Fn . We assume the as-
sumptions (17) and (19). If
v ′ = [v − R1:k [Av − b]k ]1 .      (10)
γ = (α + cβ + cαβ)(1 + cu) < 1,       (21)
Put rk = [Av − b]k and let Φ(v) = [v − R1:k rk ]1 .
Then, we can write                                      the relative forward error ∥vn − v ∗ ∥∞ /∥v ∗ ∥∞ re-
duces until
v ′ = Φ(v).                (11)
∥vn − v ∗ ∥∞        cu
∗∥
≈u+         .         (22)
The following holds:                                                    ∥v ∞           1−γ
v ′ = v − R1:k [(Av − b) + er ] + em ,   (12)   Here, for real numbers a and b, a ≈ b means that
a is approximately equal to b.
where er = rk − (Av − b) and em ∈ Rn satisfying           This implies the forward stability of the iterative
∥er ∥∞     cuk ∥Av − b∥∞            (13)   reﬁnement algorithm (20).

and                                                     3. Backward Stability
∥em ∥∞      cu∥v − R1:k rk ∥∞ .       (14)
In this paper, we will show the forward and back-       In this section, we will show the backward sta-
ward stability of the iterative algorithm (10). Fur-    bility of the iterative reﬁnement algorithm (20).
thermore, numerical examples are also given for il-        A normwise backward error of an approximation
lustrating the forward and backward stability of the    v is deﬁned by
iterative reﬁnement algorithm (10). The forward               η(v)   =   min{ε : (A + ∆A)v = b + ∆b,
stability of the algorithm guarantees that approx-
∥∆A∥∞ ε∥A∥∞ ,
imate solutions generated by the algorithm con-
verge, while the backward stability means the sta-                       ∥∆b∥∞       ε∥b∥∞ }.           (23)
bility of the algorithm against the rounding errors.    It is known [13] that
∥r∥∞
2. Convergence Theorem: Forward Stability                            η(v) =                    .        (24)
∥A∥∞ ∥v∥∞ + ∥b∥∞
Let us consider                                       Here, r = Av − b.
The next theorem shows the backward stability
Av = b,                   (15)
of the iterative reﬁnement algorithm (20):
where A ∈ Fn×n and b ∈ Fn . Let                         Theorem 2 Let vn be generated from (20) with
1 < uκ(A) < ∞.                 (16)   any starting vector v0 ∈ Fn . We assume the as-
sumptions (17) and (19). If
We assume that we have a good approximate in-
verses R1:k satisfying                                            γ = (α + cβ + cαβ)(1 + cu) < 1,       (25)

∥R1:k A − I∥∞          α < 1.       (17)   the backward error η(v) reduces until
η(v)    c2 u,             (26)
Here, R1:k is deﬁned as
where c2 is a certain constant. Here, for real num-
R1:k = R1 + R2 + · · · + Rk          (18)
bers a and b, a b means that a is approximately
with Ri ∈ Fn×n . As mentioned in the previous sec-      equal to b or a is less than b.
tion in Ref.[8], Rump has proposed a method of cal-        This implies the backward stability of the itera-
culating such approximate inverses and in Ref.[6],      tive reﬁnement algorithm (20).
we have partially clariﬁed the mechanism of the
convergence of Rump’s method. Further, we as-           4. Numerical Examples Illustrating For-
sume also that the following is satisﬁed:                  ward and Backward Stability
uk κ(A)     β < 1.             (19)      In this section, we will present numerical exam-
We propose the following iterative reﬁnement algo-      ples illustrating the forward and the backward sta-
rithm:                                                  bility of the iterative reﬁnement algorithm (20).
We have used the IEEE 754 double precision
vn = Φ(vn−1 ), Φ(v) = [v − R1:k rk ]1 ,        ﬂoating point number system in these numerical
rk = [Av − b]k (n = 1, 2, · · · )       (20)   calculations. Thus, in the following calculations,
the unit round-oﬀ u is given as
with any starting vector v0 ∈ Fn . The aim of this
section is to show the following theorem:                   1.11 × 10−16 < u = 2−53 < 1.12 × 10−16 .    (27)

- 517 -
4.1. Hilbert Matrix                                        where
Let H be the n × n Hilbert matrix. Let further                          R = R1 + R2 + · · · + R8            (34)
A = sH. Here, s is the minimum common multi-
with suitable R1 , R2 , · · · , R8 ∈ F. The iterative
plier of 1, 2, 3, · · · , n − 1. Furthermore,
reﬁnement algorithm (20) converges with 3 itera-
b = Az,                     (28)   tions.

where, z ∈ Fn and zi = 1. We have solved Ax = b
for n = 20. In this example, 1.92 × 1016 < ∥A∥∞ <                      Table 2: Rump’s matrix (n=100)
1.93 × 1016 , 1.92 × 1016 < ∥b∥∞ < 1.93 × 1016 and
2.44 × 1028 < κ(A) < 2.45 × 1028 .                                 i    ∥v ∗ − vi ∥∞ /∥v ∗ ∥∞       η(vi )
In this case, a good approximate inverse can be                  0        7.51 × 10−6         3.98 × 10−14
constructed with k = 2 such that                                   1        5.98 × 10−11        5.61 × 10−19
2        4.88 × 10−16        2.46 × 10−19
∥RA − I∥∞ < α = 4.16 × 10−4 ,           (29)           3        3.18 × 10−19        6.58 × 10−19
4        3.18 × 10−19        6.58 × 10−19
where
R = R1 + R2                    (30)
with suitable R1 , R2 ∈ F. The iterative reﬁnement
Moreover, it is seen that
algorithm (20) converges with 3 iterations. We ﬁ-
nally have an approximate solution with the rela-
β    = u8 κ(A) < (1.12 × 10−16 )8 × 1.75 × 10107
tive maximum error about 1.92 × 10−16 . Further-
more, it is seen that                                            < 4.34 × 10−21 .                        (35)

β   = u2 κ(A) < (1.2 × 10−16 )2 × 2.45 × 1028            Table 3 shows the relative errors and the backward
< 3.08 × 10−4 .                               (31)   errors of approximate solutions obtained by the it-
erative reﬁnement calculations (20). The calcula-
Table 1 shows the relative errors                          tions are done by the same computational environ-
∥v ∗ − vi ∥∞                         ment as that for the previous example.
(32)
∥v ∗ ∥∞
4.3. Rump’s Matrix (n=300)
and the backward errors η(vi ) of approximate so-
lutions obtained by the iterative reﬁnement calcu-           Let A be n × n matrix generated by Rump’s
lations (20). These calculations are done by MAT-          algorithm [12]. We choose n = 300 and b =
LAB on Intel core 2 duo CPU.                               (1, 1, · · · , 1)T ∈ Fn . In this example, 3.10 × 108 <
∥A∥∞ < 3.11 × 108 , ∥b∥∞ = 1 and 6.28 × 1059 <
Table 1: Hilbert matrix (n=20)                 κ(A) < 6.29 × 1059 .
In this case, a good approximate inverse can be
i   ∥v ∗ − vi ∥∞ /∥v ∗ ∥∞        η(vi )            constructed with k = 5 such that
0       3.50 × 10−4          4.55 × 10−6
1       4.03 × 10−9          4.12 × 10−11                      ∥RA − I∥∞ < α = 1.16 × 10−9 ,           (36)
2       5.10 × 10−14         5.04 × 10−16
3       1.91 × 10−16         1.77 × 10−18          where
4       1.91 × 10−16         1.77 × 10−18                          R = R1 + R2 + · · · + R5            (37)

with suitable R1 , R2 , · · · , R5 ∈ F. The iterative re-
ﬁnement algorithm converges (20) with 1 iteration.
4.2. Rump’s Matrix (n=100)                                   Moreover, it is seen that
Let A be n × n matrix generated by Rump’s
β     = u5 κ(A) < (1.12 × 10−16 )5 × 6.29 × 1059
algorithm [12]. We choose n = 100 and b =
(1, 1, · · · , 1)T ∈ Fn . In this example, 1.04 × 1016 <          < 1.11 × 10−20 .                       (38)
∥A∥∞ < 1.05 × 1016 , ∥b∥∞ = 1 and 1.74 × 10107 <
κ(A) < 1.75 × 10107 .                                      Table 3 shows the relative errors and the backward
In this case, a good approximate inverse can be          errors of approximate solutions obtained by the it-
constructed with k = 8 such that                           erative reﬁnement calculations (20). The calcula-
tions are done by the same computational environ-
∥RA − I∥∞ < α = 1.86 × 10−4 ,           (33)   ment as that for the previous example.

- 518 -
http://www.ti3.tu-harburg.de
Table 3: Rump’s matrix (n=300)                    /publications/rump.
i    ∥v ∗ − vi ∥∞ /∥v ∗ ∥∞       η(vi )          [11] T. Ohta, T. Ogita, S. M. Rump and S. Oishi:
0        8.02 × 10−12        1.42 × 10−17             ”A Method of Veriﬁed Numerical Computa-
1        8.10 × 10−23        4.07 × 10−19             tion for Ill-conditioned Linear System of Equa-
2        8.10 × 10−23        4.07 × 10−19             tions”, Journal of JSIAM,15:3 (2005), pp. 269–
287 in Japanese.
[12] S. M. Rump: A class of arbitrarily ill-
References                             conditioned ﬂoating-point matrices, SIAM J.
[1] G.B.Moler, ”iterative reﬁnement in ﬂoating             Matrix Anal. Appl., 12:4 (1991), 645–653.
point”, J. Assoc. Comput. Mach., 14 316-321       [13] J. D. Rigal and J. Gaches:”On the compati-
(1967)                                                 vility of a given solution with the data of a
[2] R. D. Skeel, ”Iterative reﬁnement implies nu-          linear equation”, J. Assoc. Comput. Mach., 14
merical stability for Gaussian elimination”,           (1967), 543-548.
Math. Comp., 35 817-832 (1980)
[3] N.J.Higham, ”Iterative reﬁnement for linear
systems and LAPACK”, IMA J. Numer. Anal.
17 495-509 (1997).
[4] M.Jankowsky and H.Woznlakowski, ”Iterative
reﬁnement implies numerical stability”, BIT
17 303-311 (1997).
[5] F.Tisseur, ”Newton’s method in ﬂoating point
arithmetic and iterative reﬁnement of gener-
alized eigenvalue problems”, SIAM J. Matrix
Anal. Appl. 22 No.4 1038-1057 (2001).
[6] S. Oishi, K. Tanabe, T.Ogita, and S.M. Rump,
”Convergence of Rump’s method for inverting
arbitrary ill-conditioned matrices”, J. Comp.
and Appl. Math, 205 533-544 (2007).
[7] T. Ogita, S. M. Rump and S. Oishi: ”Accu-
rate Sum and Dot Product”, SIAM Journal on
Scientiﬁc Computing, 26/6,1955-1988,(2005)
[8] S.M.Rump:”Approximate inverses of almost
singular matrices still contain useful informa-
tion”, Technical Report 90.1, Faculty of In-
formation and Communication Sciences, Ham-
burg University of Technology (1990).
[9] S.M. Rump, T. Ogita, and S. Oishi, ”Accurate
Floating-point Summation I: Faithful Round-
ing”, accepted for publication in SIAM Jour-
nal on Scientiﬁc Computing. Preprint is avail-
able from

http://www.ti3.tu-harburg.de
/publications/rump.

[10] S.M. Rump, T. Ogita, and S. Oishi, ”Accu-
rate Floating-point Summation II: Sign, K-
fold Faithful and Rounding to Nearest”, sub-
mitted for publication in SIAM Journal on Sci-
entiﬁc Computing. Preprint is available from

- 519 -

```
To top