<<

. 36
( 83 .)



>>

Now, if W k; n and Lk; n denote the subspaces W k; n = span{ 2 sn ; : : : ; 2 sn+k’1 } and Lk; n =
(n) (n)
span{y1 ; : : : ; yk }, then from (2.3) and (2.4), the generalized residuals satisÿes
˜
˜ (n)
r(tk ) ’ sn ∈ W k; n (2.5)

and
˜
˜ (n)
r(tk )⊥ Lk; n : (2.6)

˜ (n)
Conditions (2.5) and (2.6) show that the generalized residual r(tk ) is obtained by projecting, the
˜ ˜
vector sn onto the subspace W k; n , orthogonally to Lk; n .
(n)
In a matrix form, r(tk ) can be written as
˜

˜ (n) 2
Sk; n (LT n 2
Sk; n )’1 LT n sn ;
r(tk ) = sn ’ (2.7)
k; k;

(n) (n)
where Lk; n , Sk; n and 2 Sk; n are the k —k matrices whose columns are y1 ; : : : ; yk ; sn ; : : : ; sn+k’1
˜ (n)
and 2 sn ; : : : ; 2 sn+k’1 respectively. Note that r(tk ) is well deÿned if and only if the k — k matrix
152 K. Jbilou, H. Sadok / Journal of Computational and Applied Mathematics 122 (2000) 149“165


LT n 2 Sk; n is nonsingular; a necessary condition for this is that the matrices Lk; n and 2
Sk; n are full
k;
(n)
rank. In this case, tk exists and is uniquely given by
(n)
Sk; n (LT n 2
Sk; n )’1 LT n sn :
tk = sn ’ (2.8)
k; k;
(n)
The approximation tk can also be expressed as
k
(n) (n)
tk = ÿj sn+j
j=0

with
k
(n)
ÿj = 1
i=0

and
k
(n) [n)
i; j ÿj = 0; j = 0; : : : ; k ’ 1;
j=0

(n)
where the coe cients are deÿned by
i; j
(n)
= ( sn+i ; sn+j ) for the MPE method;
i; j
(n) 2
=( sn+i ; sn+j ) for the RRE method;
i; j
(n)
= (yi+1 ; sn+j ) for the MPE method; i = 0; : : : ; k ’ 1 and j = 0; : : : ; k:
i; j
(n)
From these relations it is not di cult to see that tk can also be written as a ratio of two
determinants as follows:
sn sn+1 : : : sn+k 1 1 ::: 1
(n) (n) (n) (n) (n) (n)
::: :::
0; 0 0; 1 0; 0 0; 1
0; k 0; k
(n)
tk = : (2.9)
. . . . . .
. . . . . .
. . . . . .
(n) (n) (n) (n) (n) (n)
::: :::
k’1; 0 k’1; 1 k’1; k k’1; 0 k’1; 1 k’1; k

The determinant in the numerator of (2.9) is the vector obtained by expanding this determinant
with respect to its ÿrst row by the classical rule.
Note that the determinant in the denominator of (2.9) is equal to det(LT n 2 Sk; n ) which is as-
k;
(n)
sumed to be nonzero. The computation of the approximation tk needs the values of the terms
sn ; sn+1 ; : : : ; sn+k+1 .

2.2. Application to linear systems

Consider the system of linear equations
Cx = f; (2.10)
where C is a real nonsingular N — N matrix, f is a vector of RN and x— denotes the unique solution.
K. Jbilou, H. Sadok / Journal of Computational and Applied Mathematics 122 (2000) 149“165 153


Instead of applying the extrapolation methods for solving (2.10), we will use them for the pre-
conditioned linear system
M ’1 Cx = M ’1 f; (2.11)
where M is a nonsingular matrix.
Starting from an initial vector s0 , we construct the sequence (sj )j by
sj+1 = Bsj + b; j = 0; 1; : : : (2.12)
with B = I ’ A; A = M ’1 C and b = M ’1 f.
Note that if the sequence (sj ) is convergent, its limit s = x— is the solution of the linear system
(2.10).
From (2.12) we have sj = b ’ Asj = r(sj ), the residual of the vector sj . Therefore using (2.3)
(n)
and (2.12), it follows that the generalized residual of the approximation tk is the true residual
˜ (n) (n) (n)
r(tk ) = r(tk ) = b ’ Atk : (2.13)
Note also that, since 2 sn = ’A sn , we have 2 Sk; n = ’A Sk; n .
(0)
For simplicity and unless speciÿed otherwise, we set n = 0, we denote tk = tk and we drop the
index n in our notations. Let d be the degree of the minimal polynomial P d of B for the vector
s0 ’ x— and, as A = I ’ B is nonsingular, Pd is also the minimal polynomial of B for r0 = s0 .
Therefore, the matrices Sk = [ s0 ; : : : ; sk’1 ] and 2 Sk = [ 2 s0 ; : : : ; 2 sk’1 ] have full rank for
k6d. We also note that the approximation td exits and is equal to the solution of the linear system
(2.10).
The three extrapolation methods make use implicitly of the polynomial P d and since this poly-
nomial is not known in practice, the aim of these methods is to approximate it.
When applied to the sequence generated by (2.12), the vector extrapolation methods above produce
approximations tk such that the corresponding residuals rk = b ’ Atk satisfy the relations
˜ ˜
rk ∈ W k = AV k (2.14)
and
˜
rk ⊥ Lk ; (2.15)
˜ ˜ ˜ ˜ ˜ ˜ ˜
where V k = span{ s0 ; : : : ; sk’1 } and Lk ≡ W k for RRE, Lk ≡ V k for MPE and Lk ≡ Y k = span
{y1 ; : : : ; yk } for MMPE where y1 ; : : : ; yk are linearly independent vectors.
˜
Note that, since W k ≡ Kk (A; Ar0 ), the extrapolation methods above are Krylov subspace methods.
RRE is an orthogonal projection and is theoretically equivalent to GMRES while MPE and MMPE
are oblique projection methods and are equivalent to the method of Arnoldi and to the Hessenberg
method [38], respectively. From this observation, we conclude that for k6d, the approximation tk
exists and is unique, unconditionally for RRE, and this is not always the case for MPE and MMPE.
In fact, for the last two methods the approximation tk (k ¡ d) exists if and only if det( Sk 2 Sk ) = 0
T

for MPE and det(YkT 2 Sk ) = 0 for MMPE where Yk = [y1 ; : : : ; yk ].
˜
Let Pk be the orthogonal projector onto W k . Then from (2.14) and (2.15), the residual generated
by RRE can be expressed as
rre
rk = r0 ’ Pk r0 : (2.16)
154 K. Jbilou, H. Sadok / Journal of Computational and Applied Mathematics 122 (2000) 149“165

˜ ˜ ˜
We also consider the oblique projectors Qk and Rk onto W k and orthogonally to V k and Y k respec-
tively. It follows that the residuals produced by MPE and MMPE can be written as
mpe
rk = r0 ’ Qk r0 (2.17)
and
mmpe
rk = r 0 ’ R k r0 : (2.18)
˜
The acute angle ‚k between r0 and the subspace W k is deÿned by
|(r0 ; z)|
cos ‚k = max : (2.19)
||r0 ||||z||
˜
z∈W k ’{0}

Note that ‚k is the acute angle between the vector r0 and Pk r0 .
In the sequel we give some relations satisÿed by the residual norms of the three extrapolation
methods.

Theorem 1. Let k be the acute angle between r0 and Qk r0 and let denote the acute angle
k
between r0 and Rk r0 . Then we have the following relations:
(1) ||rk ||2 = (sin2 ‚k )||r0 ||2 .
rre
mpe
(2) ||rk ||2 = (tan2 k )||r0 ||2 .
mpe
rre
(3) ||rk ||6(cos k )||rk ||.
Moreover if for MMPE yj = r0 for some j = 1; : : : ; k; then we also have
mmpe
(4) ||rk ||2 = (tan2 k )||r0 ||2 .
mmpe
rre
(5) ||rk ||6(cos k )||rk ||.

Proof. Parts (1) “ (3) have been proved in [18]
(4) From (2.18), we get
mmpe mmpe mmpe
(rk ; rk ) = (rk ; r0 ’ Rk r0 ):
mmpe
Since (rk ; r0 ) = 0, it follows that
mmpe mmpe mmpe
(rk ; rk ) = (rk ; ’Rk r0 )
mmpe mmpe
= ’||rk ||||Rk r0 ||cos(rk ; Rk r0 )
mmpe
= ||rk ||||Rk r0 ||sin k:

On the other hand,
||r0 || = ||Rk r0 ||cos k;

hence
mmpe
||rk || = ||r0 ||tan k:

(5) Using statements (1) and (4), we get
mmpe
||rk ||2 1 ’ cos2 k
(cos2 ’1
= k) :
rre 2
||rk || 1 ’ cos2 ‚k
K. Jbilou, H. Sadok / Journal of Computational and Applied Mathematics 122 (2000) 149“165 155


But cos k 6cos ‚k , therefore
mmpe
rre
||rk ||6||rk ||cos k:


Remark.

• From relations (1), (2) and (4) of Theorem 1, we see that the residuals of the RRE are always
deÿned while those produced by MPE and MMPE may not exist.
rre
• We also observe that if a stagnation occurs in RRE (||rk || = ||r0 || for some k¡d), then cos ‚k = 0
and, from (2.19), this implies that cos k = cos k = 0 and hence the approximations produced by
MPE and MMPE are not deÿned.

When the linear process (2.12) is convergent, it is more useful in practice to apply the extrapolation
methods after a ÿxed number p of basic iterations. We note also that, when these methods are used
in their complete form, the required work and storage grow linearly with the iteration step. To
overcome this drawback we use them in a cycling mode and this means that we have to restart the
algorithms after a chosen number m of iterations.
The algorithm is summarized as follows:
1. k = 0, choose x0 and the numbers p and m.
2. Basic iteration
set t0 = x0
z0 = t 0
zj+1 = B zj + b, j = 0; : : : ; p ’ 1.
3. Extrapolation scheme
s 0 = zp
sj+1 = B sj + b, j = 0; : : : ; m,
compute the approximation tm by RRE, MPE or MMPE.
4. Set x0 = tm , k = k + 1 and go to 2.

Stable schemes for the computation of the approximation tk are given in [32, 19]. In [32], Sidi gave
an e cient implementation of the MPE and RRE methods which is based on the QR decomposition
of the matrix Sk . In [19], we used an LU decomposition of Sk with a pivoting strategy. These
implementations require low work and storage and are more stable numerically.

2.3. Application to nonlinear systems

Consider the system of nonlinear equations
G(x) = x; (2.20)
where G : RN ’ RN and let x— be a solution of (2.20).
For any arbitrary vector x, the residual is deÿned by
r(x) = G(x) ’ x:
Let (sj )j be the sequence of vectors generated from an initial guess s0 as follows:

<<

. 36
( 83 .)



>>