- Email: [email protected]
Chapter 7 Look-ahead methods
Chapter 7
Look-ahead methods As we have already mentioned, under certain circumstances (see Causes of breakdown, page 24) both GODir and GOMin break down because P f GP^, whose inverse is required by the algorithm, is either singular or too badly conditioned for accurate inversion. It is possible in certain circumstances to overcome this problem by taking the appropriate remedial action, the nature of this remedial action being dictated by the precise cause of the singularity or near-singularity. If it is due to the columns of P^ becoming linearly dependent or nearly so then we have one type of instability which is dealt with in the next chapter. If, on the other hand, P^ has full rank and the singularity of P ^ G P ^ is due solely to the indefiniteness of G, then we have serious breakdown, and a different course of action is required. This action is referred to as "look-ahead" and it is this approach that we now describe. We assume in this chapter, therefore, that P^ always has full rank except on termination. It is convenient to use another of our formal algorithms to give an overview of the look-ahead technique but before doing so we first, for purposes of comparison, give a formal definition of the basic algorithm, i.e. without look-ahead. This basic algorithm is just GODir (or GOMin) with Step 2(e) replaced by compute
P^^i =
QfWi^i
where W^+i is either KGP^ or KFj^i as appropriate, and is as follows: The basic algorithm (1) Select G, H, Pi = KFi (2) while PjGPi (a) compute (b) compute (c) compute (d) compute (e) compute
K and Xi. Set Qo = I and i = 1. Compute Fi = G X i - H and is nonsingular D^ = PJGP^ Xi^i = Xi - P^D^^Pf F^ F^+i = F^ - GP^D^-^Pf F^ Q^ = Qi_i - G P ^ D ^ l p f Pi+i = Qf W^^i for the appropriate Wi4.i 133
134
Chapter
7. Look-ahead
methods
(f) set 2 = z -f 1 (3) end while (4) end t h e basic algorithm. This algorithm terminates when PfGPi becomes singular, something which must inevitably occur sooner or later (see Chapter 2). When this happens prematurely it is not possible t o compute t h e new projection matrix Q i , b u t there is no reason why t h e matrices P i + j , j = 1,2,..., should not be computed using t h e old oblique projector Q ^ - i instead of t h e updated one. This at least ensures t h a t the matrices P i + j , j = 0 , 1 , . . . , are conjugate t o t h e matrices P^;, A: = 1,2, ...,z — 1, even though they will not be conjugate t o each other. However this process cannot be allowed t o continue indefinitely since t h e algorithm functions by constructing a sequence of conjugate matrices, t h e linear-algebraic costs of not doing so being prohibitive. T h e obvious way forward is t o keep adjoining t h e matrices Pi-^j t o t h e matrix P i j u s t as in t h e basic algorithm in t h e hope t h a t , sooner or later, Qi as defined by equation (1.46) will become computable. If it does then it may be calculated a n d t h e algorithm can proceed happily on its way. If it does not then we have t h e problem of incurable hreakdovm for which, as t h e name implies, there is no known remedy apart from restarting from a different initial approximate solution. T h e formal representation of t h e look-ahead version is The look-ahead version (1) Select G , H , K and X i . Set Qo = I , z = 1 a n d A: - 1. Compute F i = G X i - H and J » i = K F i (2) set Pk = "Pi (3) while P^GPk is singular (a) if Pfc is rank-deficient, stop, else (b) compute P i ^ i = Q ^ . ^ K C P ^ (c) s e t P A : = \Pk Pi-hi (d) set z = 2-h 1 (4) endwhile (5) c o m p u t e hk
= P^GPA:
(6) compute X ^ + i = XA. - PfcD'^Pl'F^t (7) compute F ^ + i = F ^ - G P ^ D ^ ^ P ^ F ^ (8) c o m p u t e QA^ = QAT-I -
(9) (10) (11) (12)
GPkDj^^Pl
compute P i + i = Q ^ W ^ + i for t h e appropriate W j ^ . ! set z = 2 + 1, A: = A: 4- 1 return t o Step 2 end t h e look-ahead version.
R e m a r k 2 . This algorithm terminates when t h e columns of P ^ become linearly dependent, as they must eventually do. If thej*ank of Pk is zero then t h e solution has been computed. If it is equal t o t h a t of D ^ t h e retrieval of a partial solution
135 should be possible. If it exceeds t h a t of Djt then incurable breakdown has occurred. See Chapter 2, and Chapter 8 for further details of obtaining a partial solution. Note also t h a t in Step 3(b) it is not possible t o replace K G P i by KFfc^-l since Ffc_|_i can only be calculated after the while-loop has terminated. It follows from t h e above description of the look-ahead version t h a t if P f G P ^ is nonsingular for all i t h e while-loop is never entered and the algorithm reduces t o t h e basic algorithm. If, on the other hand, P f G P ^ is singular for some i then the whileloop will be activated and the resulting matrix computed, P ^ , will have dimensions determined by t h e number of times the while-loop is traversed. Suppose t h a t the loop is traversed tk — I times. Then P ^ will consist of tk blocks each consisting of one of t h e matrices P j for appropriate values of i. In fact, if we define ^o = 0 ^^^ Sk, k = 1,2,..., by k
sk = ^tj
(7.1)
it is readily seen t h a t P f c = [ P . . _ . +lP»._>+2... P a J .
(7.2)
It follows from this and equation (1.31) t h a t an alternative expression for P^^ when look-ahead is used is
(7.3)
p . . = P i P2 Now a glance at Step 3 and Step 9 of the algorithm indicates t h a t Pk =
Q^^i^k
for some matrix W ^ while from Step 8, Q, = I - ^ G P , D 7 i p J .
(7.4)
Comparison of these two equations with Steps 2a and 2b of the generalised GramSchmidt algorithm (see page 23) indicates t h a t the look-ahead versions of both GODir and GOMin are simply more sophisticated examples of this algorithm. They thus generate sequences of conjugate matrices \Pi\ and their projectors, by analogy with equations (1.49) and (1.50), satisfy QlPs. = 0
(7.5)
and
T h e matrices Fk^i and X^j+i are computed by Steps 6 and 7 of t h e algorithm only when a new nonsingular D ^ has been computed, and are given by Xfc+i = Xjfc — PjfcD^
^k^k
136
Chapter 7. Look-ahead methods
and F;,^i = Fk- GPkB-'PlFk.
(7.7)
or F^^i = ( l - G P f c D - i p ^ ) F f c . Thus k
F.+i = n ( l - G P , D - ' P j ) F i or, from the conjugacy of the matrices P j , F^^i = Q^Fi
(7.8)
where Qk is defined by equation (7.4) so that PjFk^i
= 0.
l
(7.9)
Equation (7.8) is essentially the same as equation (1.43) except that the rank of Qk is less than that of Qi if any look-ahead steps have been taken. Finally, it is convenient at this point to define T^ by Tfc = P , D - > = [ T . , _ , + i T . , , , + 2 . . . T . J
(7.10)
where the partitioning of T^^ echoes of that of P^; in equation (7.2), and to express Qk in terms of TA:. With TA; thus defined equation (7.7) may be written
F;,+i = Ffc - G I
^
P,Tj I Fk
(7.11)
(note: since D^ is a full matrix rather than a block-diagonal one there is no simple relationship between the matrices Pj and Tj). The projector Qk becomes, since T>k is symmetric,
Qk = I-GY,^,fJ
= l-Gf2^jTj.
(7.12)
7.0.1. Note on notation From equations (7.2) - (7.4) and the fact that that
Q.=i-GP.,(pr^GP.,)"'pr,
PJGPA;
=
O,
j ^ k, it follows
7.1.
The computational
versions
137
and so might be designated Qs^. However in this context Q j only exists for j = Si, i = 1, 2,..., /c, and even for these values of j is not in general given by equation (1.45). For these reasons the designation Q^ is preferred.
7.1. T h e computational versions In just the same way t h a t , if G and K are both symmetric, GODir and GOMin reduce t o t h e block Lanczos and block HS algorithms respectively, similar simplifications may be made t o the look-ahead versions of the basic algorithm in the case when these two matrices are symmetric. Since the matrices P j are mutually conjugate it follows from equation (7.4) t h a t Q^-lGPfc+l = GPfc+1,
1 < i < A: + 1.
(7.13)
Now for 5i_i + 1 < j < 5i the matrices P j are computed by P,=Qr_iW,
(7.14)
SO t h a t premultiplying equation (7.13) by W ^ gives P
GP^+i=W
GP^+i.
This equation holds for b o t h Si^i -\- I < j < Si and i = 1,2, ...,/c + 1, and since SQ = 0 by definition (see text immediately preceding equation (7.1)) it follows t h a t P j G P ; , ^ ! = WjGPfc+i,
1 < J < Sk^i.
(7.15)
Hence, since P f GPjt-t-i = O for 1 < i < k, W j G P ^ + i = O,
1 < J < Sk.
(7.16)
Similarly, by analogy with equation (7.13), equations (7.4) and (7.9) give a _ i F f c + i = Ffc+i,
l
(7.17)
Premultiplying equation (7.17) by W j then yields, from equation (7.14), PjFfc+i = WjFfc+i,
5,_i + 1 < j < 5^.
A parallel argument t o t h a t used to establish equation (7.15) then gives P j F f c + i = WjFk^i,
1 < j < Sk^i
so t h a t , from equations (7.2) and (7.9), WjF;,.^i=.0,
l
(7.18)
Equations (7.16) and (7.18) permeate the whole of conjugate gradient theory. T h e y were used implicitly in Chapter 3 t o derive the short recurrences and will
138
Chapter
7. Look-ahead
methods
be called upon t o perform the same office again in Chapter 8. It is the properties of these equations t h a t , if the matrices Wjt are chosen appropriately, enable us t o obtain short recurrences when computing the matrix sequences { P i } for t h e look-ahead algorithms. J u s t as in the case without look-ahead we have essentially two choices for W/^ leading again t o t h e Lanczos and to the HS versions of t h e algorithms, and as before t h e Lanczos choice leads t o a slightly more complicated algorithm but by a simpler route. T h e analysis for t h e look-ahead methods is inherently more complicated t h a n t h a t for the standard ones as not only do we have to consider the existence of the while-loop (Step 3 of the algorithm) but also the fact t h a t this loop is not executed a fixed number of times. As a result of this the notation required for the lookahead methods is significantly more complicated t h a n t h a t for the standard ones. Moreover, for the HS versions of the algorithms, we are forced t o use two different values of W^^ since the one t h a t essentially defines t h e HS method is not available until after t h e while-loop is terminated. T h e approach here will be similar t o t h a t of Chapter 3. We first obtain look-ahead variants of b o t h main versions, HS and Lanczos, of the block-CG algorithm in terms of G and K . We then substitute in the equations so obtained t h e particular values of these two matrices appropriate t o the specific algorithms. 7.LL
The look-ahead block Lanczos
algorithm
For this version we choose W j by Wj = KGPj-i
(7.19)
for all values of j so equation (7.16) becomes, on replacing j by j -f 1 for clarity, PjGKGPk+i
= O,
l
(7.20)
Suppose now t h a t the matrices P j , j = 1, 2 , . . . , /c, have already been computed, a n d consider the calculation of PAJ+I. Equations (7.14) and (7.19) yield Pi = Q^KGPi_i,
Sk + l
sk^i.
(7.21)
P u t t i n g i = Sk-\- I then gives, from equation (7.12),
or, from equation (7.20) and since P^^ is a block column of P^;,
P , , + i = [ l - f;
T.PjGJKGP,,.
(7.23)
7.1. The computational versions
139
An appeal to equation (7.10) then yields P.,+i= (l-T,,_,Pj;_^G-TfcP^G)KGP,,.
(7.24)
For 5/e. + 2 < 2 < Sk^i equation (7.21) still holds but the matrices Pi_i are now block-columns of Pk-^i. This equation and equation (7.12) give
P^=
JI-X^T,PJG|KGP,_I
and this, from equation (7.20), becomes P^ = (I -
T,,PJ;G)
KGPi_i.
(7.25)
Note that, if the while-loop (Step 3) of the look-ahead version is never invoked, equation (7.24) reverts to the standard block Lanczos algorithm. 7.1.2. Uie Hestenes-Stiefel version In this version, based on GOMin, the matrices W^ are chosen by W^ = KFi
(7.26)
in the absence of breakdown but if, for some i, P f GP^ is singular then it is not possible to compute Fi4.i and we thus cannot choose Wi+i using this equation. We can, though, make the Lanczos choice, i.e. Wi^i = KGP^, so that if we want to apply look-ahead to the HS version we arrive at a hybrid algorithm, owing something to both the Lanczos and the Hestenes-Stiefel ideas. The derivation of the computational form of the look-ahead HS version of the algorithm is, notwithstanding the greater simplicity of the final result, considerably more complicated than that of the Lanczos form. An inspection of the formal lookahead version of the algorithm shows that, just as in the Lanczos form, it generates sequences of matrices \Pi\ and ^ Qi [, defined by equations (7.2) and (7.4), where the matrices P^ are mutually conjugate with respect to G. If the integers Si are as previously defined it follows from the algorithm and equation (7.2) that the block columns of P^ are computed by P j = Q ^ ^ W j where W^- = KFi
(7.27)
for j = Si_i -h 1 and W^ = KGP^-_i.
(7.28)
for 5i_i 4- 2 < j < ^i, where i = 1,2,..., /c. It now follows from these two equations and from equations (7.18) and (7.27) that F f KFfc^i = O
(7.29)
140
Chapter 7. Look-ahead
methods
for 2 = 1,2,..., k, and from equations (7.18) and (7.28) that PjGKF;t+i = O
(7.30)
for all values of j between 1 and Sk — I inclusive except si,S2,^..,Sk-i- We shall now show^ that equation (7.30) holds even for these values of j if certain conditions are satisfied. Substituting i for k in equation (7.11) and transposing gives
Ff+i - Ff = - F f which on post multiplying by KF^^+i gives, from equation (7.29), f2
T , P j I GKF^+i = O,
1 < z < A: - 1.
l
(7.31)
If we now make the assumption (see below) that F f T^. is nonsingular for 1 < i < k — 1 this leads immediately to the conclusion that Pj;GKFfc+i = 0 ,
l
and combining this equation with equation (7.30) gives P j G K F ; t + i = O,
l
(7.32)
This is one of the key results that enable us to obtain a simple computational implementation of the look-ahead version of the block HS algorithm. The other key result is equation (7.20), which is proved for the HS methods by noting that equation (7.16) is simply equation (7.18) with GPA;_^I substituted for F/c+i. Making the same substitution in equations (7.29) - (7.32) and following the same argument used to derive equation (7.32) gives the desired result. With these two equations established we can now obtain the computational form of the HS version. Suppose then that the matrices P^, i = 1,2, .../c, together with Qk and FA;+I have already been computed. The first (and possibly only) block column of P^-f-i is then given by P^.^i^Q^KF;,^!.
(7.33)
or, from equation (7.12),
Sk-^l
I-5^T,PjG
KF.^i.
7.i.
The computational
versions
141
It then follows immediately from equation (7.32) t h a t P,,+i = (l-T,,Pj;G)KFfc+i.
(7.34)
T h e remaining blocks of Pfc+i, if any, are computed by P^ = Q j K G P i _ i ,
sk-^2
Sk^i,
where the matrices P i - i are themselves submatrices of Pfc+i- Equations (7.12) and (7.20) then give Pi = {l-
T , , P J ; G ) KGPi_i,
sk + 2
6'fc+i,
(7.35)
and this equation, together with equation (7.34), defines the final computational form of the algorithm. Note t h a t unlike t h e Lanczos case, the projection matrices of equations (7.34) and (7.35) are identical. We note t h a t in order to obtain these two equations it is necessary t o assume t h a t F ^ T s ^ is nonsingular. For t h e basic version (no look-ahead) a similar assumption is necessary in order t o derive the corresponding formula but in t h a t case we have been able t o relate the nonsingularity of P j F j t t o t h e original Krylov sequence (see Theorem 19 above). Although we conjecture t h a t a similar result holds in t h e look-ahead case we know of no comparable theorem t o Theorem 19. To this extent our analysis of look-ahead methods is incomplete. T h e above analysis completes the derivation of the basic computational forms of the two look-ahead algorithms but before we go on to examine particular examples we prove the following theorem. T h e o r e m 2 7 . Let G be symmetric and nonsingular and K be symmetric positive definite. Assume t h a t P ^ , X ^ + i and F^^+i, /c = L 2 , . . . , i , have been successfully computed by one or other version of t h e look-ahead algorithm and let the first block column of P i + i , Ps^+i say, be computed either by P,,+i = Q f K G F , . (Lanczos) or P^i+i = Qi KFi_^i (HS). Then if P s , + i has full rank and Pj]^_iGPs,4-i = O (serious breakdown), only one look-ahead step will be necessary to restore the algorithm t o normality. Proof. Since Ps^+i is a submatrix of P^+i it follows from conjugacy and equation (7.4) t h a t QiGP,,+i = GP,,+i. Now Psi+2 is computed by Psi+2 = Q f K G P g . + i so t h a t D,+i =
P^
G [ P , , + i P,,.+2]
(7.36)
142
Chapter
7. Look-ahead
methods
or
(7.37)
Di4-1 = and it follows from hypothesis and equation (7.36) t h a t O B B C
where B = P f . _ ^ i G K G P s .^_i. From t h e hypotheses of the theorem, B is nonsingular and hence D i + i is nonsingular so t h a t Xi+2, Fi^2 and Qi+i may be duly computed. • This theorem underlines yet again the crucial importance of the properties of the matrices G and K t o the functioning of conjugate-gradient algorithms and indicates why it is desirable (at least in theory) for at least one of these two matrices to be positive definite. It may be thought t h a t the requirement t h a t Pj]_^2GPs._|_i = O (rather t h a n being merely singular) is a weakness of t h e theorem but in most practical cases, as we shall see, Pj^_^iGP5.+i is either itself a scalar or has the [O a l form r. for some scalar a . It is thus either nonsingular (in which case there is
Ia UJ
no breakdown) or null (when the theorem appUes).
7.2. Particular algorithms We turn now t o the look-ahead versions of some common algorithms, which we obtain from equations (7.24) etc. by making the appropriate substitutions for G and K (see C h a p t e r 3 for the precise values of these two matrices). We begin with the b£isic CG algorithm generating vector sequences and applied t o a symmetric b u t indefinite system A x = b . Thus G = A and is indefinite and p j G p ^ is either zero or nonzero. K is t h e identity matrix and is thus positive definite so t h a t , from Theorem 27, a t most only one step of look-ahead is ever needed. _ _ _ _ If ts^ is the last column (first or second as appropriate) of T ^ where Tjt = ^kJ^k ' equations (7.34) and (7.35) give p,,^i = ( I - t , , p i ; A ) r , + i
(7.38)
and, if ps^_,_2 needs to be computed, P«.+-2 = P . . , , = (I - t . ^ p j . A ) A p . ^ ^ i .
(7.39)
A method similar t o this for overcoming serious breakdown was first proposed by Luenberger [180] who called it the method of hyperbolic pairs. T h e method given above is equivalent to the improved version of Luenberger's algorithm given by Fletcher [110].
7.2. Particular
143
algorithms
Now for all the other look-ahead versions outlined here (BiCG, MRZ, t h e HegediisGalerkin algorithm and Q M R ) , P j has the form of equation (3.26) and G is given either by
G =
O A^ A O
0 I
or
1
O
0 a for some a 0 scalar a so t h a t P ^ G P j is either nonsingular or null, making possible an appeal t o Theorem 27. T h e second is t h a t This impUes two things. The first is t h a t P ^ G P j is equal t o
Us^_i+l 0
0 U, fc-i+2 v,,_j+i 0
0 V,^_^+2
0
V,
Now the matrix P ^ D ^ ^ P ^ is unchanged if P A ^ M is substituted for Pj^, where M is any nonsingular matrix, and since column permutation is equivalent^to postmultipUcation by such a matrix we may assume when evaluating PfcD^ P;, t h a t Ufc 0 ^ 0 V;, where Ufc -= [ U s ^ _ j + i , Us^_^+2, . • - , Us J
with a similar expression for V^. Galerkin algorithms, where G =
Ql
I-TJCZ^VIA O
Thus for the BiCG, MRZ and the HegediisO A T-[ , we have A O
O
(7.40)
I-VfcC-^U^A^
where Ck - V [ A U f c . Let now (7.41) and denote its last column by Zs.. The look-ahead version of BiCG becomes, from equation (7.34), Usfc + l = ( I - ZsfcVj^ A ) Tk+l
(7.42)
while t h e corresponding equation for MRZ is, from equation (7.24), u.,+1 = ( l - ZfcVjA - z,,_,vf^_^ A ) A U , , .
(7.43)
144
Chapter
7. Look-ahead
methods
For both these algorithms, for Sk -\- I < i < Sk+i — 1, Ui+i is given by
(7.44)
Ui+i = ( I - 2 - . v . . A)Au^.
with corresponding equations for v^^-^i and Vi-^i. For the BiCG and MRZ algorithms K is indefinite so t h a t Theorem 27 does not apply, and for these algorithms incurable breakdown is a remote possibility. For the Hegedtis algorithms, on the other hand, K is equal t o the unit matrix so t h a t only one step of look-ahead is necessary [150]. Making the substitutions appropriate to the Hegediis-Galerkin algorithm in equations (7.24) and (7.34) then yields, for the HS version (HG),
and for the Lanczos version (HGL), u , , + i = (I - ZkVlA
- z,,_,vj;_^ A ) A ^ v , ,
while for both versions, if needed, equation (7.35) yields Usf^+2 =
Us
(I-z«fcvJ^A)A^Vs^+i.
We finally derive the equations for t h e look-ahead version of QMR. For this algorithm equation (7.40) becomes
Ql =
i-Ufc(v^Ufc)"'v:
o
o
Vk{uiVky'ui
and if Zfc is now defined by Zfc = U f c ( v r U f c ) * '
(7.45)
with, as before, its last column denoted by Zs^ we obtain u , , + i = ( l - ZfcVf - Zsfc_,vf^_^) A u , ,
(7.46)
and, for Sfe + 1 < j < Sk+i - 1,
Uj+i = (I - z^^vJJ Auj
(7.47)
with corresponding equations for Vs^_^i and Vj_|_i. These equations are somewhat simpler t h a n those given in [117], which take into account certain numerical considerations designed t o make the matrices V ^ U ^ less ill-conditioned. For further details see the original reference. However, since G and K are both indefinite it is still remotely possible t h a t Q M R can suffer from incurable breakdown even if look-ahead is fully implemented.
7.3. More Krylov aspects'''
145
7.3. M o r e Krylov aspects* We now consider the look-ahead algorithms from the point of view of Chapter 4, a n d see how Theorems 16 and 17 have t o be modified for the more sophisticated algorithms of this chapter. We then discuss the implications of the analysis for the look-ahead algorithms. T h e o r e m 2 8 . Let the matrices Sj be defined by equation (4.2) and be nonsingular for j = 5i, 3 2 , . . . , Sk-^i where si < 52 < . . . < s^+i. Then the sequence of matrices { P j } is computable by the look-ahead version of GODir for I < j < Sk-\-i and the matrices P j satisfy
p.- = P . U , -
(7.48)
for some block unit upper triangular, and hence nonsingular, matrix U j . Proof. By induction. We show t h a t if the theorem is true for 1 < j < i, where i is any integer such t h a t s^
Qk^I-GPsA'K
(7.49)
Denote now t h e submatrices of JJj by Urs-, etc. Since equation (7.48) is assumed t o hold for j = i we have, from equations (1.31) and (4.1), i
Pi = ^ ( K G ) ^ - ' P i U j , j=i
so that KGPi = Y, (KG)^ PlUj^•
(7.50)
For the look-ahead version of GODir and for Sk
(7.51)
SO t h a t since, from equations (4.1) and (7.49), Q ^ {KGY
P i = O,
0
(7.52)
it follows t h a t Pi.^i = Q r ^ ( K G ) ^ P i U , , .
(7.53)
146
Chapter
7. Look-ahead
methods
Hence since Ua = I by hypothesis it follows from equation (7.49) t h a t P i + i = Pi,
Yi4.l I
for some Y i + i G W^^^.
But since equation (7.48) is t r u e for j = i this gives
Pi-hi = P i ^ i U i + i where U 14-1
O
(7.54)
I
so t h a t equation (7.48) is t r u e for 2 + 1 , establishing t h e induction. Since t h e theorem is trivially true for j = I ( U i = U n = I) it is true for 1 < j < ^^-^1 — 1, completing t h e proof. • T h e corresponding theorem for GOMin is: T h e o r e m 2 9 . Let Sj be defined by equation (4.2) a n d be nonsingular for j = si^S2,... ,Sk-\-i where ^i < 52 < . . • < ^Ar+i- If t h e matrices Tg. are defined by equation (7.10) a n d t h e matrices Tj.Fi computed by the look-ahead version of GOMin are nonsingular for I < i < k then the sequences of matrices { P j } and {F^} are computable by GOMin for 1 < j < ^^.^i, and t h e matrices P^ satisfy P,- = P , U ,
(7.55)
for some nonsingular block upper triangular matrix U^. Proof. By induction. We show t h a t if the theorem is true for 1 < j < z, where i is any integer such t h a t s^ < i < Sk-^-i — 1, then it is also true for j = i -{-1. Since S^^ is nonsingular from t h e hypotheses of t h e theorem it follows from Lemma 6 and equations (4.2), (7.3) a n d (7.4) t h a t Qk is computable and is given by equation (7.49). Let first i = s^. For this case P.,+1 = QlKFk^i
(7.56)
so t h a t , from equation (7.7), P . , + 1 = Q ^ K ( F ^ - GPfcD^-ip^Ffc) . Now from equations (7.6) a n d (7.56),
QjKFfc = QlQl_,KFk
= QJP.,_,+I = O
since Q ^ P ^ = O and Psfe_i-Hi is a submatrix of P ^ . Hence, from equation (7.10), P..+1 = -QlKGPkTlFk. Now equation (7.55) is assumed to be satisfied for j — Sk so t h a t QjKGP,,
=Q^KGP,,U,,.
(7.57)
7.3. More Krylov
147
aspects*
Hence, from equation (7.52), Q^KGP,, =
O . . . O Q ^ (KG)''^ P i
U Sk
a n d if
u 11
u
u.,= o
u..
l,Sfc
u.
o we have Q^KGP,, =
O...OQnKGrPiU.„.,
where the number of null block columns in t h e right-hand side is Sk — 1. Since Pk is a submatrix of P^^ (see equation (7.3)) this implies t h a t Q ^ K G P i t = [ o . . O Q j ( K G ) ^ ^ P i U Ste,Ste where t h e right-hand side now has Sk — Sk~i — 1 null block columns. Postmultiplying this equation by T ^ then gives, from equation (7.10), QlKGPkTl
= Q ^ (KG)^^ P j U Sk,Sk
-^ Sk
SO t h a t , from equation (7.57),
This equation may be written, from equations (4.1) and (7.49), as Sk+l
Sk + l
Ysfc+1 ^Sk
(7.58)
+ l,Sk-hl
for some matrix Y5^4.] and where
(7.59) Combining equation (7.58) and equation (7.55) with j = Sk gives
where U Sfc +
Usfc Ys^+1 l
—
(7.60)
and since this is equation (7.55) with j = sk-\- 1, t h e induction is established for
148
Chapter
7. Look-ahead
methods
For 5fc + l
(7.61)
(as we cannot now assume t h a t Ua = I ) . Since Us^^Sk ^ ^ ^ ^a ^^^ both nonsingular by the induction hypothesis and since TJ Fk is nonsingular from the hypotheses of the theorem it follows from equations (7.59) and (7.61) t h a t b o t h Ug^+i and Ui4.i are nonsingular block upper triangular. This establishes the induction for Sk
7,4. Practical details One of the most difficult problems in numerical analysis is deciding whether or not a small computed quantity should be regarded as zero, and it is effectively this problem t h a t has to be resolved at every iteration of a look-ahead Lanczos
7.4. Practical details
149
or HS method. In this context t h e problem presents itself in Step 3 of t h e lookahead version of t h e basic algorithm (see page 134) as the determination of t h e singularity or otherwise of P f G P ^ , where G is indefinite. There appears t o be no single diagnostic procedure t o cover the general case though specific tests have been devised for particular versions of the algorithm ([180], [36], [37], [150], [151], [117]). Luenberger [180] appears t o have been the first t o propose what is in effect a look-ahead version of a C O method when he suggested the method of hyperbolic pairs. This is applied t o t h e original algorithm of Hestenes and Stiefel when the latter is used to solve A x = b where A is symmetric but indefinite, and is based on t h e proof of Theorem 27. If vector sequences are being generated then, after one step of look-ahead, equation (7.37) may be written, from equation (7.36) with G = A , as D^+l =
pr,4.iAp,^^i pJ^^+iAKAp.^^^i pj;+iAKAp,^+i pJ^^^2Ap,,4.2
If pJ.^_iApg.^_i is exactly zero and K is positive definite then D^+i is nonsingular and the next values of x and p can be determined. However if pj._j_j A p ^ . ^ i is only small, Luenberger in effect set it equal t o zero and computed the new x and p using t h e mutilated value of D^+i. Under certain circumstances, if Di+i is near t o singularity, this can introduce serious inaccuracies and the algorithm in this form was never really satisfactory. Fletcher [110] proposed, in effect, using the true matrix D i + i b u t did not give a specific test for deciding when pj]_^j Ap^..^^ is too small to be used as a simple divisor. He did, though, note an analogy between his algorithm and an algorithm by Bunch and Parlett [54] for factorising a symmetric indefinite matrix which experiences similar problems, and suggested t h a t "it might be possible t o a d a p t the test of Bunch and Parlett t o work here". Fletcher's lack of enthusiasm for a look-ahead version of CG was based on t h e not unreasonable premise t h a t it might be better t o use a method like O D or CR t h a t is incapable of crashing rather t h a n try t o shore up an inherently unstable algorithm. However both O D and C R can stagnate if A is indefinite and t o date nobody else seems t o have taken u p t h e challenge of finding an efficient and effective look-ahead version of C G when applied t o indefinite systems. An algorithm for which this challenge certainly has been taken up is the original version of QMR, but for this algorithm it was found t h a t in order t o achieve reliability it was necessary for three different tests t o be simultaneously satisfiedj^efore taking a normal step. If in equations (7.45) - (7.47) we scale the matrices U^ and V i so t h a t their columns have unit length, denote the scaled versions by Ui and V i and define Y^ and Zi by Yf=(vfu,)
'vf TTr=rr;r
and
Z, = U^ ( v f u ^ ) " '
then since the matrix V^ U j must be nonsingular in order for the algorithm t o function t h e first thing t h a t could be checked is its smallest singular value, ai say. If this is greater t h a n some prescribed tolerance tol (say tol = e^'"^ or tol = e '
150
Chapter
7. Look-ahead
methods
where e is the machine precision - see [204]) then the while-loop in Step 3 can b e safely abandoned. However one problem with this approach is t h a t in some cases (7i is never greater t h a t t h e required tolerance. T h e algorithm remains t r a p p e d in the while-loop which has, in effect, created an artificial incurable breakdown. Even if this ultimate disaster is avoided t h e algorithm as described can show a marked sensitivity t o the value of tol, not a desirable a t t r i b u t e in a "black-box" algorithm. To overcome these problems it was suggested by Freund et al [117] t h a t the above criterion b e relaxed by reducing tol t o s but a t the same time strengthening the procedure by the imposition of two further checks. The improved version of the algorithm escapes from t h e while-loop only if all three of the following conditions are satisfied: (1) ^ i > ^
(2) (3)
j ^ i - l , i ,
and
j=i-l,i,
where n(A) is some function of ||A|| whose detailed calculation is discussed in [117] and where the notation used here is t h a t of equations (7.45) - (7.47). Similar tests have been incorporated into t h e final versions of other algorithms [36], [37], [150], [151], with various degrees of success. Unfortunately there seems t o be very little in the way of comparative testing t o indicate which, if any, method is the best for a particular t y p e of problem and some of the testing t h a t has been done has lacked rigour. Like has occasionally not been tested against like, e.g. a version of algorithm A incorporating look-ahead being compared with a version of algorithm B t h a t lacks this refinement, and this makes evaluation difficult. However t h e lookahead version of Q M R seems t o perform reasonably satisfactorily, and this makes it all the more surprising t h a t we still await an equally satisfactory look-ahead version of CG when the latter is applied t o indefinite systems.
Look-ahead methods As we have already mentioned, under certain circumstances (see Causes of breakdown, page 24) both GODir and GOMin break down because P f GP^, whose inverse is required by the algorithm, is either singular or too badly conditioned for accurate inversion. It is possible in certain circumstances to overcome this problem by taking the appropriate remedial action, the nature of this remedial action being dictated by the precise cause of the singularity or near-singularity. If it is due to the columns of P^ becoming linearly dependent or nearly so then we have one type of instability which is dealt with in the next chapter. If, on the other hand, P^ has full rank and the singularity of P ^ G P ^ is due solely to the indefiniteness of G, then we have serious breakdown, and a different course of action is required. This action is referred to as "look-ahead" and it is this approach that we now describe. We assume in this chapter, therefore, that P^ always has full rank except on termination. It is convenient to use another of our formal algorithms to give an overview of the look-ahead technique but before doing so we first, for purposes of comparison, give a formal definition of the basic algorithm, i.e. without look-ahead. This basic algorithm is just GODir (or GOMin) with Step 2(e) replaced by compute
P^^i =
QfWi^i
where W^+i is either KGP^ or KFj^i as appropriate, and is as follows: The basic algorithm (1) Select G, H, Pi = KFi (2) while PjGPi (a) compute (b) compute (c) compute (d) compute (e) compute
K and Xi. Set Qo = I and i = 1. Compute Fi = G X i - H and is nonsingular D^ = PJGP^ Xi^i = Xi - P^D^^Pf F^ F^+i = F^ - GP^D^-^Pf F^ Q^ = Qi_i - G P ^ D ^ l p f Pi+i = Qf W^^i for the appropriate Wi4.i 133
134
Chapter
7. Look-ahead
methods
(f) set 2 = z -f 1 (3) end while (4) end t h e basic algorithm. This algorithm terminates when PfGPi becomes singular, something which must inevitably occur sooner or later (see Chapter 2). When this happens prematurely it is not possible t o compute t h e new projection matrix Q i , b u t there is no reason why t h e matrices P i + j , j = 1,2,..., should not be computed using t h e old oblique projector Q ^ - i instead of t h e updated one. This at least ensures t h a t the matrices P i + j , j = 0 , 1 , . . . , are conjugate t o t h e matrices P^;, A: = 1,2, ...,z — 1, even though they will not be conjugate t o each other. However this process cannot be allowed t o continue indefinitely since t h e algorithm functions by constructing a sequence of conjugate matrices, t h e linear-algebraic costs of not doing so being prohibitive. T h e obvious way forward is t o keep adjoining t h e matrices Pi-^j t o t h e matrix P i j u s t as in t h e basic algorithm in t h e hope t h a t , sooner or later, Qi as defined by equation (1.46) will become computable. If it does then it may be calculated a n d t h e algorithm can proceed happily on its way. If it does not then we have t h e problem of incurable hreakdovm for which, as t h e name implies, there is no known remedy apart from restarting from a different initial approximate solution. T h e formal representation of t h e look-ahead version is The look-ahead version (1) Select G , H , K and X i . Set Qo = I , z = 1 a n d A: - 1. Compute F i = G X i - H and J » i = K F i (2) set Pk = "Pi (3) while P^GPk is singular (a) if Pfc is rank-deficient, stop, else (b) compute P i ^ i = Q ^ . ^ K C P ^ (c) s e t P A : = \Pk Pi-hi (d) set z = 2-h 1 (4) endwhile (5) c o m p u t e hk
= P^GPA:
(6) compute X ^ + i = XA. - PfcD'^Pl'F^t (7) compute F ^ + i = F ^ - G P ^ D ^ ^ P ^ F ^ (8) c o m p u t e QA^ = QAT-I -
(9) (10) (11) (12)
GPkDj^^Pl
compute P i + i = Q ^ W ^ + i for t h e appropriate W j ^ . ! set z = 2 + 1, A: = A: 4- 1 return t o Step 2 end t h e look-ahead version.
R e m a r k 2 . This algorithm terminates when t h e columns of P ^ become linearly dependent, as they must eventually do. If thej*ank of Pk is zero then t h e solution has been computed. If it is equal t o t h a t of D ^ t h e retrieval of a partial solution
135 should be possible. If it exceeds t h a t of Djt then incurable breakdown has occurred. See Chapter 2, and Chapter 8 for further details of obtaining a partial solution. Note also t h a t in Step 3(b) it is not possible t o replace K G P i by KFfc^-l since Ffc_|_i can only be calculated after the while-loop has terminated. It follows from t h e above description of the look-ahead version t h a t if P f G P ^ is nonsingular for all i t h e while-loop is never entered and the algorithm reduces t o t h e basic algorithm. If, on the other hand, P f G P ^ is singular for some i then the whileloop will be activated and the resulting matrix computed, P ^ , will have dimensions determined by t h e number of times the while-loop is traversed. Suppose t h a t the loop is traversed tk — I times. Then P ^ will consist of tk blocks each consisting of one of t h e matrices P j for appropriate values of i. In fact, if we define ^o = 0 ^^^ Sk, k = 1,2,..., by k
sk = ^tj
(7.1)
it is readily seen t h a t P f c = [ P . . _ . +lP»._>+2... P a J .
(7.2)
It follows from this and equation (1.31) t h a t an alternative expression for P^^ when look-ahead is used is
(7.3)
p . . = P i P2 Now a glance at Step 3 and Step 9 of the algorithm indicates t h a t Pk =
Q^^i^k
for some matrix W ^ while from Step 8, Q, = I - ^ G P , D 7 i p J .
(7.4)
Comparison of these two equations with Steps 2a and 2b of the generalised GramSchmidt algorithm (see page 23) indicates t h a t the look-ahead versions of both GODir and GOMin are simply more sophisticated examples of this algorithm. They thus generate sequences of conjugate matrices \Pi\ and their projectors, by analogy with equations (1.49) and (1.50), satisfy QlPs. = 0
(7.5)
and
T h e matrices Fk^i and X^j+i are computed by Steps 6 and 7 of t h e algorithm only when a new nonsingular D ^ has been computed, and are given by Xfc+i = Xjfc — PjfcD^
^k^k
136
Chapter 7. Look-ahead methods
and F;,^i = Fk- GPkB-'PlFk.
(7.7)
or F^^i = ( l - G P f c D - i p ^ ) F f c . Thus k
F.+i = n ( l - G P , D - ' P j ) F i or, from the conjugacy of the matrices P j , F^^i = Q^Fi
(7.8)
where Qk is defined by equation (7.4) so that PjFk^i
= 0.
l
(7.9)
Equation (7.8) is essentially the same as equation (1.43) except that the rank of Qk is less than that of Qi if any look-ahead steps have been taken. Finally, it is convenient at this point to define T^ by Tfc = P , D - > = [ T . , _ , + i T . , , , + 2 . . . T . J
(7.10)
where the partitioning of T^^ echoes of that of P^; in equation (7.2), and to express Qk in terms of TA:. With TA; thus defined equation (7.7) may be written
F;,+i = Ffc - G I
^
P,Tj I Fk
(7.11)
(note: since D^ is a full matrix rather than a block-diagonal one there is no simple relationship between the matrices Pj and Tj). The projector Qk becomes, since T>k is symmetric,
Qk = I-GY,^,fJ
= l-Gf2^jTj.
(7.12)
7.0.1. Note on notation From equations (7.2) - (7.4) and the fact that that
Q.=i-GP.,(pr^GP.,)"'pr,
PJGPA;
=
O,
j ^ k, it follows
7.1.
The computational
versions
137
and so might be designated Qs^. However in this context Q j only exists for j = Si, i = 1, 2,..., /c, and even for these values of j is not in general given by equation (1.45). For these reasons the designation Q^ is preferred.
7.1. T h e computational versions In just the same way t h a t , if G and K are both symmetric, GODir and GOMin reduce t o t h e block Lanczos and block HS algorithms respectively, similar simplifications may be made t o the look-ahead versions of the basic algorithm in the case when these two matrices are symmetric. Since the matrices P j are mutually conjugate it follows from equation (7.4) t h a t Q^-lGPfc+l = GPfc+1,
1 < i < A: + 1.
(7.13)
Now for 5i_i + 1 < j < 5i the matrices P j are computed by P,=Qr_iW,
(7.14)
SO t h a t premultiplying equation (7.13) by W ^ gives P
GP^+i=W
GP^+i.
This equation holds for b o t h Si^i -\- I < j < Si and i = 1,2, ...,/c + 1, and since SQ = 0 by definition (see text immediately preceding equation (7.1)) it follows t h a t P j G P ; , ^ ! = WjGPfc+i,
1 < J < Sk^i.
(7.15)
Hence, since P f GPjt-t-i = O for 1 < i < k, W j G P ^ + i = O,
1 < J < Sk.
(7.16)
Similarly, by analogy with equation (7.13), equations (7.4) and (7.9) give a _ i F f c + i = Ffc+i,
l
(7.17)
Premultiplying equation (7.17) by W j then yields, from equation (7.14), PjFfc+i = WjFfc+i,
5,_i + 1 < j < 5^.
A parallel argument t o t h a t used to establish equation (7.15) then gives P j F f c + i = WjFk^i,
1 < j < Sk^i
so t h a t , from equations (7.2) and (7.9), WjF;,.^i=.0,
l
(7.18)
Equations (7.16) and (7.18) permeate the whole of conjugate gradient theory. T h e y were used implicitly in Chapter 3 t o derive the short recurrences and will
138
Chapter
7. Look-ahead
methods
be called upon t o perform the same office again in Chapter 8. It is the properties of these equations t h a t , if the matrices Wjt are chosen appropriately, enable us t o obtain short recurrences when computing the matrix sequences { P i } for t h e look-ahead algorithms. J u s t as in the case without look-ahead we have essentially two choices for W/^ leading again t o t h e Lanczos and to the HS versions of t h e algorithms, and as before t h e Lanczos choice leads t o a slightly more complicated algorithm but by a simpler route. T h e analysis for t h e look-ahead methods is inherently more complicated t h a n t h a t for the standard ones as not only do we have to consider the existence of the while-loop (Step 3 of the algorithm) but also the fact t h a t this loop is not executed a fixed number of times. As a result of this the notation required for the lookahead methods is significantly more complicated t h a n t h a t for the standard ones. Moreover, for the HS versions of the algorithms, we are forced t o use two different values of W^^ since the one t h a t essentially defines t h e HS method is not available until after t h e while-loop is terminated. T h e approach here will be similar t o t h a t of Chapter 3. We first obtain look-ahead variants of b o t h main versions, HS and Lanczos, of the block-CG algorithm in terms of G and K . We then substitute in the equations so obtained t h e particular values of these two matrices appropriate t o the specific algorithms. 7.LL
The look-ahead block Lanczos
algorithm
For this version we choose W j by Wj = KGPj-i
(7.19)
for all values of j so equation (7.16) becomes, on replacing j by j -f 1 for clarity, PjGKGPk+i
= O,
l
(7.20)
Suppose now t h a t the matrices P j , j = 1, 2 , . . . , /c, have already been computed, a n d consider the calculation of PAJ+I. Equations (7.14) and (7.19) yield Pi = Q^KGPi_i,
Sk + l
sk^i.
(7.21)
P u t t i n g i = Sk-\- I then gives, from equation (7.12),
or, from equation (7.20) and since P^^ is a block column of P^;,
P , , + i = [ l - f;
T.PjGJKGP,,.
(7.23)
7.1. The computational versions
139
An appeal to equation (7.10) then yields P.,+i= (l-T,,_,Pj;_^G-TfcP^G)KGP,,.
(7.24)
For 5/e. + 2 < 2 < Sk^i equation (7.21) still holds but the matrices Pi_i are now block-columns of Pk-^i. This equation and equation (7.12) give
P^=
JI-X^T,PJG|KGP,_I
and this, from equation (7.20), becomes P^ = (I -
T,,PJ;G)
KGPi_i.
(7.25)
Note that, if the while-loop (Step 3) of the look-ahead version is never invoked, equation (7.24) reverts to the standard block Lanczos algorithm. 7.1.2. Uie Hestenes-Stiefel version In this version, based on GOMin, the matrices W^ are chosen by W^ = KFi
(7.26)
in the absence of breakdown but if, for some i, P f GP^ is singular then it is not possible to compute Fi4.i and we thus cannot choose Wi+i using this equation. We can, though, make the Lanczos choice, i.e. Wi^i = KGP^, so that if we want to apply look-ahead to the HS version we arrive at a hybrid algorithm, owing something to both the Lanczos and the Hestenes-Stiefel ideas. The derivation of the computational form of the look-ahead HS version of the algorithm is, notwithstanding the greater simplicity of the final result, considerably more complicated than that of the Lanczos form. An inspection of the formal lookahead version of the algorithm shows that, just as in the Lanczos form, it generates sequences of matrices \Pi\ and ^ Qi [, defined by equations (7.2) and (7.4), where the matrices P^ are mutually conjugate with respect to G. If the integers Si are as previously defined it follows from the algorithm and equation (7.2) that the block columns of P^ are computed by P j = Q ^ ^ W j where W^- = KFi
(7.27)
for j = Si_i -h 1 and W^ = KGP^-_i.
(7.28)
for 5i_i 4- 2 < j < ^i, where i = 1,2,..., /c. It now follows from these two equations and from equations (7.18) and (7.27) that F f KFfc^i = O
(7.29)
140
Chapter 7. Look-ahead
methods
for 2 = 1,2,..., k, and from equations (7.18) and (7.28) that PjGKF;t+i = O
(7.30)
for all values of j between 1 and Sk — I inclusive except si,S2,^..,Sk-i- We shall now show^ that equation (7.30) holds even for these values of j if certain conditions are satisfied. Substituting i for k in equation (7.11) and transposing gives
Ff+i - Ff = - F f which on post multiplying by KF^^+i gives, from equation (7.29), f2
T , P j I GKF^+i = O,
1 < z < A: - 1.
l
(7.31)
If we now make the assumption (see below) that F f T^. is nonsingular for 1 < i < k — 1 this leads immediately to the conclusion that Pj;GKFfc+i = 0 ,
l
and combining this equation with equation (7.30) gives P j G K F ; t + i = O,
l
(7.32)
This is one of the key results that enable us to obtain a simple computational implementation of the look-ahead version of the block HS algorithm. The other key result is equation (7.20), which is proved for the HS methods by noting that equation (7.16) is simply equation (7.18) with GPA;_^I substituted for F/c+i. Making the same substitution in equations (7.29) - (7.32) and following the same argument used to derive equation (7.32) gives the desired result. With these two equations established we can now obtain the computational form of the HS version. Suppose then that the matrices P^, i = 1,2, .../c, together with Qk and FA;+I have already been computed. The first (and possibly only) block column of P^-f-i is then given by P^.^i^Q^KF;,^!.
(7.33)
or, from equation (7.12),
Sk-^l
I-5^T,PjG
KF.^i.
7.i.
The computational
versions
141
It then follows immediately from equation (7.32) t h a t P,,+i = (l-T,,Pj;G)KFfc+i.
(7.34)
T h e remaining blocks of Pfc+i, if any, are computed by P^ = Q j K G P i _ i ,
sk-^2
Sk^i,
where the matrices P i - i are themselves submatrices of Pfc+i- Equations (7.12) and (7.20) then give Pi = {l-
T , , P J ; G ) KGPi_i,
sk + 2
6'fc+i,
(7.35)
and this equation, together with equation (7.34), defines the final computational form of the algorithm. Note t h a t unlike t h e Lanczos case, the projection matrices of equations (7.34) and (7.35) are identical. We note t h a t in order to obtain these two equations it is necessary t o assume t h a t F ^ T s ^ is nonsingular. For t h e basic version (no look-ahead) a similar assumption is necessary in order t o derive the corresponding formula but in t h a t case we have been able t o relate the nonsingularity of P j F j t t o t h e original Krylov sequence (see Theorem 19 above). Although we conjecture t h a t a similar result holds in t h e look-ahead case we know of no comparable theorem t o Theorem 19. To this extent our analysis of look-ahead methods is incomplete. T h e above analysis completes the derivation of the basic computational forms of the two look-ahead algorithms but before we go on to examine particular examples we prove the following theorem. T h e o r e m 2 7 . Let G be symmetric and nonsingular and K be symmetric positive definite. Assume t h a t P ^ , X ^ + i and F^^+i, /c = L 2 , . . . , i , have been successfully computed by one or other version of t h e look-ahead algorithm and let the first block column of P i + i , Ps^+i say, be computed either by P,,+i = Q f K G F , . (Lanczos) or P^i+i = Qi KFi_^i (HS). Then if P s , + i has full rank and Pj]^_iGPs,4-i = O (serious breakdown), only one look-ahead step will be necessary to restore the algorithm t o normality. Proof. Since Ps^+i is a submatrix of P^+i it follows from conjugacy and equation (7.4) t h a t QiGP,,+i = GP,,+i. Now Psi+2 is computed by Psi+2 = Q f K G P g . + i so t h a t D,+i =
P^
G [ P , , + i P,,.+2]
(7.36)
142
Chapter
7. Look-ahead
methods
or
(7.37)
Di4-1 = and it follows from hypothesis and equation (7.36) t h a t O B B C
where B = P f . _ ^ i G K G P s .^_i. From t h e hypotheses of the theorem, B is nonsingular and hence D i + i is nonsingular so t h a t Xi+2, Fi^2 and Qi+i may be duly computed. • This theorem underlines yet again the crucial importance of the properties of the matrices G and K t o the functioning of conjugate-gradient algorithms and indicates why it is desirable (at least in theory) for at least one of these two matrices to be positive definite. It may be thought t h a t the requirement t h a t Pj]_^2GPs._|_i = O (rather t h a n being merely singular) is a weakness of t h e theorem but in most practical cases, as we shall see, Pj^_^iGP5.+i is either itself a scalar or has the [O a l form r. for some scalar a . It is thus either nonsingular (in which case there is
Ia UJ
no breakdown) or null (when the theorem appUes).
7.2. Particular algorithms We turn now t o the look-ahead versions of some common algorithms, which we obtain from equations (7.24) etc. by making the appropriate substitutions for G and K (see C h a p t e r 3 for the precise values of these two matrices). We begin with the b£isic CG algorithm generating vector sequences and applied t o a symmetric b u t indefinite system A x = b . Thus G = A and is indefinite and p j G p ^ is either zero or nonzero. K is t h e identity matrix and is thus positive definite so t h a t , from Theorem 27, a t most only one step of look-ahead is ever needed. _ _ _ _ If ts^ is the last column (first or second as appropriate) of T ^ where Tjt = ^kJ^k ' equations (7.34) and (7.35) give p,,^i = ( I - t , , p i ; A ) r , + i
(7.38)
and, if ps^_,_2 needs to be computed, P«.+-2 = P . . , , = (I - t . ^ p j . A ) A p . ^ ^ i .
(7.39)
A method similar t o this for overcoming serious breakdown was first proposed by Luenberger [180] who called it the method of hyperbolic pairs. T h e method given above is equivalent to the improved version of Luenberger's algorithm given by Fletcher [110].
7.2. Particular
143
algorithms
Now for all the other look-ahead versions outlined here (BiCG, MRZ, t h e HegediisGalerkin algorithm and Q M R ) , P j has the form of equation (3.26) and G is given either by
G =
O A^ A O
0 I
or
1
O
0 a for some a 0 scalar a so t h a t P ^ G P j is either nonsingular or null, making possible an appeal t o Theorem 27. T h e second is t h a t This impUes two things. The first is t h a t P ^ G P j is equal t o
Us^_i+l 0
0 U, fc-i+2 v,,_j+i 0
0 V,^_^+2
0
V,
Now the matrix P ^ D ^ ^ P ^ is unchanged if P A ^ M is substituted for Pj^, where M is any nonsingular matrix, and since column permutation is equivalent^to postmultipUcation by such a matrix we may assume when evaluating PfcD^ P;, t h a t Ufc 0 ^ 0 V;, where Ufc -= [ U s ^ _ j + i , Us^_^+2, . • - , Us J
with a similar expression for V^. Galerkin algorithms, where G =
Ql
I-TJCZ^VIA O
Thus for the BiCG, MRZ and the HegediisO A T-[ , we have A O
O
(7.40)
I-VfcC-^U^A^
where Ck - V [ A U f c . Let now (7.41) and denote its last column by Zs.. The look-ahead version of BiCG becomes, from equation (7.34), Usfc + l = ( I - ZsfcVj^ A ) Tk+l
(7.42)
while t h e corresponding equation for MRZ is, from equation (7.24), u.,+1 = ( l - ZfcVjA - z,,_,vf^_^ A ) A U , , .
(7.43)
144
Chapter
7. Look-ahead
methods
For both these algorithms, for Sk -\- I < i < Sk+i — 1, Ui+i is given by
(7.44)
Ui+i = ( I - 2 - . v . . A)Au^.
with corresponding equations for v^^-^i and Vi-^i. For the BiCG and MRZ algorithms K is indefinite so t h a t Theorem 27 does not apply, and for these algorithms incurable breakdown is a remote possibility. For the Hegedtis algorithms, on the other hand, K is equal t o the unit matrix so t h a t only one step of look-ahead is necessary [150]. Making the substitutions appropriate to the Hegediis-Galerkin algorithm in equations (7.24) and (7.34) then yields, for the HS version (HG),
and for the Lanczos version (HGL), u , , + i = (I - ZkVlA
- z,,_,vj;_^ A ) A ^ v , ,
while for both versions, if needed, equation (7.35) yields Usf^+2 =
Us
(I-z«fcvJ^A)A^Vs^+i.
We finally derive the equations for t h e look-ahead version of QMR. For this algorithm equation (7.40) becomes
Ql =
i-Ufc(v^Ufc)"'v:
o
o
Vk{uiVky'ui
and if Zfc is now defined by Zfc = U f c ( v r U f c ) * '
(7.45)
with, as before, its last column denoted by Zs^ we obtain u , , + i = ( l - ZfcVf - Zsfc_,vf^_^) A u , ,
(7.46)
and, for Sfe + 1 < j < Sk+i - 1,
Uj+i = (I - z^^vJJ Auj
(7.47)
with corresponding equations for Vs^_^i and Vj_|_i. These equations are somewhat simpler t h a n those given in [117], which take into account certain numerical considerations designed t o make the matrices V ^ U ^ less ill-conditioned. For further details see the original reference. However, since G and K are both indefinite it is still remotely possible t h a t Q M R can suffer from incurable breakdown even if look-ahead is fully implemented.
7.3. More Krylov aspects'''
145
7.3. M o r e Krylov aspects* We now consider the look-ahead algorithms from the point of view of Chapter 4, a n d see how Theorems 16 and 17 have t o be modified for the more sophisticated algorithms of this chapter. We then discuss the implications of the analysis for the look-ahead algorithms. T h e o r e m 2 8 . Let the matrices Sj be defined by equation (4.2) and be nonsingular for j = 5i, 3 2 , . . . , Sk-^i where si < 52 < . . . < s^+i. Then the sequence of matrices { P j } is computable by the look-ahead version of GODir for I < j < Sk-\-i and the matrices P j satisfy
p.- = P . U , -
(7.48)
for some block unit upper triangular, and hence nonsingular, matrix U j . Proof. By induction. We show t h a t if the theorem is true for 1 < j < i, where i is any integer such t h a t s^
Qk^I-GPsA'K
(7.49)
Denote now t h e submatrices of JJj by Urs-, etc. Since equation (7.48) is assumed t o hold for j = i we have, from equations (1.31) and (4.1), i
Pi = ^ ( K G ) ^ - ' P i U j , j=i
so that KGPi = Y, (KG)^ PlUj^•
(7.50)
For the look-ahead version of GODir and for Sk
(7.51)
SO t h a t since, from equations (4.1) and (7.49), Q ^ {KGY
P i = O,
0
(7.52)
it follows t h a t Pi.^i = Q r ^ ( K G ) ^ P i U , , .
(7.53)
146
Chapter
7. Look-ahead
methods
Hence since Ua = I by hypothesis it follows from equation (7.49) t h a t P i + i = Pi,
Yi4.l I
for some Y i + i G W^^^.
But since equation (7.48) is t r u e for j = i this gives
Pi-hi = P i ^ i U i + i where U 14-1
O
(7.54)
I
so t h a t equation (7.48) is t r u e for 2 + 1 , establishing t h e induction. Since t h e theorem is trivially true for j = I ( U i = U n = I) it is true for 1 < j < ^^-^1 — 1, completing t h e proof. • T h e corresponding theorem for GOMin is: T h e o r e m 2 9 . Let Sj be defined by equation (4.2) a n d be nonsingular for j = si^S2,... ,Sk-\-i where ^i < 52 < . . • < ^Ar+i- If t h e matrices Tg. are defined by equation (7.10) a n d t h e matrices Tj.Fi computed by the look-ahead version of GOMin are nonsingular for I < i < k then the sequences of matrices { P j } and {F^} are computable by GOMin for 1 < j < ^^.^i, and t h e matrices P^ satisfy P,- = P , U ,
(7.55)
for some nonsingular block upper triangular matrix U^. Proof. By induction. We show t h a t if the theorem is true for 1 < j < z, where i is any integer such t h a t s^ < i < Sk-^-i — 1, then it is also true for j = i -{-1. Since S^^ is nonsingular from t h e hypotheses of t h e theorem it follows from Lemma 6 and equations (4.2), (7.3) a n d (7.4) t h a t Qk is computable and is given by equation (7.49). Let first i = s^. For this case P.,+1 = QlKFk^i
(7.56)
so t h a t , from equation (7.7), P . , + 1 = Q ^ K ( F ^ - GPfcD^-ip^Ffc) . Now from equations (7.6) a n d (7.56),
QjKFfc = QlQl_,KFk
= QJP.,_,+I = O
since Q ^ P ^ = O and Psfe_i-Hi is a submatrix of P ^ . Hence, from equation (7.10), P..+1 = -QlKGPkTlFk. Now equation (7.55) is assumed to be satisfied for j — Sk so t h a t QjKGP,,
=Q^KGP,,U,,.
(7.57)
7.3. More Krylov
147
aspects*
Hence, from equation (7.52), Q^KGP,, =
O . . . O Q ^ (KG)''^ P i
U Sk
a n d if
u 11
u
u.,= o
u..
l,Sfc
u.
o we have Q^KGP,, =
O...OQnKGrPiU.„.,
where the number of null block columns in t h e right-hand side is Sk — 1. Since Pk is a submatrix of P^^ (see equation (7.3)) this implies t h a t Q ^ K G P i t = [ o . . O Q j ( K G ) ^ ^ P i U Ste,Ste where t h e right-hand side now has Sk — Sk~i — 1 null block columns. Postmultiplying this equation by T ^ then gives, from equation (7.10), QlKGPkTl
= Q ^ (KG)^^ P j U Sk,Sk
-^ Sk
SO t h a t , from equation (7.57),
This equation may be written, from equations (4.1) and (7.49), as Sk+l
Sk + l
Ysfc+1 ^Sk
(7.58)
+ l,Sk-hl
for some matrix Y5^4.] and where
(7.59) Combining equation (7.58) and equation (7.55) with j = Sk gives
where U Sfc +
Usfc Ys^+1 l
—
(7.60)
and since this is equation (7.55) with j = sk-\- 1, t h e induction is established for
148
Chapter
7. Look-ahead
methods
For 5fc + l
(7.61)
(as we cannot now assume t h a t Ua = I ) . Since Us^^Sk ^ ^ ^ ^a ^^^ both nonsingular by the induction hypothesis and since TJ Fk is nonsingular from the hypotheses of the theorem it follows from equations (7.59) and (7.61) t h a t b o t h Ug^+i and Ui4.i are nonsingular block upper triangular. This establishes the induction for Sk
7,4. Practical details One of the most difficult problems in numerical analysis is deciding whether or not a small computed quantity should be regarded as zero, and it is effectively this problem t h a t has to be resolved at every iteration of a look-ahead Lanczos
7.4. Practical details
149
or HS method. In this context t h e problem presents itself in Step 3 of t h e lookahead version of t h e basic algorithm (see page 134) as the determination of t h e singularity or otherwise of P f G P ^ , where G is indefinite. There appears t o be no single diagnostic procedure t o cover the general case though specific tests have been devised for particular versions of the algorithm ([180], [36], [37], [150], [151], [117]). Luenberger [180] appears t o have been the first t o propose what is in effect a look-ahead version of a C O method when he suggested the method of hyperbolic pairs. This is applied t o t h e original algorithm of Hestenes and Stiefel when the latter is used to solve A x = b where A is symmetric but indefinite, and is based on t h e proof of Theorem 27. If vector sequences are being generated then, after one step of look-ahead, equation (7.37) may be written, from equation (7.36) with G = A , as D^+l =
pr,4.iAp,^^i pJ^^+iAKAp.^^^i pj;+iAKAp,^+i pJ^^^2Ap,,4.2
If pJ.^_iApg.^_i is exactly zero and K is positive definite then D^+i is nonsingular and the next values of x and p can be determined. However if pj._j_j A p ^ . ^ i is only small, Luenberger in effect set it equal t o zero and computed the new x and p using t h e mutilated value of D^+i. Under certain circumstances, if Di+i is near t o singularity, this can introduce serious inaccuracies and the algorithm in this form was never really satisfactory. Fletcher [110] proposed, in effect, using the true matrix D i + i b u t did not give a specific test for deciding when pj]_^j Ap^..^^ is too small to be used as a simple divisor. He did, though, note an analogy between his algorithm and an algorithm by Bunch and Parlett [54] for factorising a symmetric indefinite matrix which experiences similar problems, and suggested t h a t "it might be possible t o a d a p t the test of Bunch and Parlett t o work here". Fletcher's lack of enthusiasm for a look-ahead version of CG was based on t h e not unreasonable premise t h a t it might be better t o use a method like O D or CR t h a t is incapable of crashing rather t h a n try t o shore up an inherently unstable algorithm. However both O D and C R can stagnate if A is indefinite and t o date nobody else seems t o have taken u p t h e challenge of finding an efficient and effective look-ahead version of C G when applied t o indefinite systems. An algorithm for which this challenge certainly has been taken up is the original version of QMR, but for this algorithm it was found t h a t in order t o achieve reliability it was necessary for three different tests t o be simultaneously satisfiedj^efore taking a normal step. If in equations (7.45) - (7.47) we scale the matrices U^ and V i so t h a t their columns have unit length, denote the scaled versions by Ui and V i and define Y^ and Zi by Yf=(vfu,)
'vf TTr=rr;r
and
Z, = U^ ( v f u ^ ) " '
then since the matrix V^ U j must be nonsingular in order for the algorithm t o function t h e first thing t h a t could be checked is its smallest singular value, ai say. If this is greater t h a n some prescribed tolerance tol (say tol = e^'"^ or tol = e '
150
Chapter
7. Look-ahead
methods
where e is the machine precision - see [204]) then the while-loop in Step 3 can b e safely abandoned. However one problem with this approach is t h a t in some cases (7i is never greater t h a t t h e required tolerance. T h e algorithm remains t r a p p e d in the while-loop which has, in effect, created an artificial incurable breakdown. Even if this ultimate disaster is avoided t h e algorithm as described can show a marked sensitivity t o the value of tol, not a desirable a t t r i b u t e in a "black-box" algorithm. To overcome these problems it was suggested by Freund et al [117] t h a t the above criterion b e relaxed by reducing tol t o s but a t the same time strengthening the procedure by the imposition of two further checks. The improved version of the algorithm escapes from t h e while-loop only if all three of the following conditions are satisfied: (1) ^ i > ^
(2) (3)
j ^ i - l , i ,
and
j=i-l,i,
where n(A) is some function of ||A|| whose detailed calculation is discussed in [117] and where the notation used here is t h a t of equations (7.45) - (7.47). Similar tests have been incorporated into t h e final versions of other algorithms [36], [37], [150], [151], with various degrees of success. Unfortunately there seems t o be very little in the way of comparative testing t o indicate which, if any, method is the best for a particular t y p e of problem and some of the testing t h a t has been done has lacked rigour. Like has occasionally not been tested against like, e.g. a version of algorithm A incorporating look-ahead being compared with a version of algorithm B t h a t lacks this refinement, and this makes evaluation difficult. However t h e lookahead version of Q M R seems t o perform reasonably satisfactorily, and this makes it all the more surprising t h a t we still await an equally satisfactory look-ahead version of CG when the latter is applied t o indefinite systems.