Residual Replacement Strategies for Krylov Subspace Iterative Methods for the Convergence of True Residuals.
In this paper, a strategy is proposed for alternative computations of the residual vectors in Krylov subspace methods, which improves the agreement of the computed residuals and the true residuals to the level of O(u) ||A|| ||x||. Building on earlier ideas on residual replacement and on insights in the finite precision behavior of the Krylov subspace methods, computable error bounds are derived for iterations that involve occasionally replacing the computed residuals by the true residuals,  and they are used to monitor the deviation of the two residuals and hence to select residual replacement steps, so that the recurrence relations for the computed residuals, which control the convergence of the method, are perturbed within safe bounds. Numerical examples are presented to demonstrate the effectiveness of this new residual replacement scheme.
Introduction
Krylov subspace iterative methods for solving a large linear system typically consist of
iterations that recursively update approximate solutions x n and the corresponding residual vectors
They can be written in a general form as follows.
Algorithm 1. Template for Krylov subspace Method:
Input: an initial approximation x
For convergence
Generate a correction vector q n by the method;
(the vector x n does not occur in other statements)
End for
Department of Mathematics, Utrecht University, P.O. Box 80010, NL-3508 Utrecht, The Netherlands E-mail:
vorst@math.uu.nl
y Department of Mathematics, University of Manitoba, Winnipeg, Manitoba, Canada R3T 2N2. E-mail:
ye@gauss.amath.umanitoba.ca Research supported by grants from University of Manitoba Research Development
Fund and from Natural Sciences and Engineering Research Council of Canada
Most Krylov subspace iterative methods, including the conjugate gradient method (CG) [12], the
bi-conjugate gradient method (Bi-CG) [4, 13], CGS [19], and BiCGSTAB [22], fit in this framework
(see [2, 11, 16] for other methods).
In exact arithmetic, the recursively defined r n in Algorithm 1 is exactly the residual for the
approximate solution x In a floating
point arithmetic, however, the round-off patterns for x n and r n will be different. It is important
to note that any error made in the computation of x n is not reflected by a corresponding error
in r n , or in other words, computational errors to x n do not force the method to correct, since x n
has no influence on the iteration process. This leads to the well known situation that b \Gamma Ax n and
r n may differ significantly. This phenomenon has been extensively discussed in the literature, see
[10, 11, 18] and the references cited there. Indeed, if we denote the computed results of x
respectively (but we still use q n to denote the computed update vector of the algorithm),
then we have
where f l(z) denotes the computed result of z in finite arithmetic, the absolute value and inequalities
on vectors are componentwise, and u is the machine roundoff unit. The vectors / n and
rounding error terms, and they can be bounded by a straightforward error analysis (see Section
3 for details). In particular, the relations (1) and (2) show that / n and j n depend only on the
iteration vectors -
r n , and q n .
We will call b \Gamma A-x n the true residual for the approximation -
x n and call -
r n , as obtained by
recurrence formula (2), the computed residual (or the updated residual). Then the difference between
the two satisfies (using the finite precision recurrences (1) and (2))
where we assume for now that b \Gamma Ax Hence, the difference between the true and the
updated residuals is a result of accumulated rounding errors. In particular, a significant deviation
of
r n may be expected, if there is a -
r i with large norm during the iteration (a
not uncommon situation for Bi-CG and CGS). On the other hand, even when all / i or j i are small
(as is common for CG), but if it takes a relatively large number of iterations for convergence, the
sheer accumulation of / i and j i could also lead to a nontrivial deviation.
What makes all this so important is that, in a finite precision implementation, the sequence - r n
satisfies, almost to machine precision u, its defining recurrence relation, and as was observed for
many Krylov subspace methods, this is the driving force behind convergence of - r n [10, 15, 18, 20].
Indeed, residual bounds have been obtained in [20] for CG and Bi-CG, which show that even
a significantly perturbed recurrence relation (with perturbations much larger than the machine
precision) usually still leads to eventual convergence of the computed residuals. This theoretical
insight has been our motivation and justification for the residual replacement scheme to be presented
in Section 2.1. On the other hand, the true residual b \Gamma A-x n itself has no self-correcting mechanism
for convergence, mainly because any perturbation made to - x n does not have an effect on the
iteration parameters, whereas errors in -
immediately lead to other iteration parameters.
Thus, in a typical convergent iteration process, - r n converges to a level much smaller than u
eventually, but the true residual b \Gamma A-x n can only converge to the level dictated by \Sigma n
since
Usually, when -
r n is still bigger than the accumulated error \Sigma n
agrees well
with -
r n in magnitude, but when - r n has converged to a level that is smaller than the accumulated
error, then
just the accumulated error and has no agreement at all
with -
r n . In summary, a straightforward implementation would reduce the true residuals at best to
bound for this has been obtained in [10] and it is called the attainable accuracy.
We note that this term could be significant even if only one of / i or j i is large, or if n is large.
The above problems become most serious in methods such as CGS and Bi-CG where intermediate
x n and - r n can have very large norm, and this may result in a large / n or j n . Several popular
methods, such as BiCGSTAB [22], BiCGSTAB(') [17], QMR [7], TFQMR [5], and composite step
BiCG [1], have been developed to reduce the norm of - r n (see [6] for details). We note that controlling
the size of k-r n k only does not solve the deviation problem in all situations, as, for instance,
the accumulation of tiny errors over a long iteration may still result in a nontrivial deviation.
A simple approach for solving the deviation problem is to replace the computed residuals by the
true residuals at some iteration step to restore the agreement. Then the deviation at subsequent
steps will be the error accumulation after that iteration only. This includes a complete replacement
strategy that simply computes r n by b \Gamma Ax n at every iteration, and a periodic replacement strategy
that updates r n by b \Gamma Ax n only at intervals of the iteration count. While such a strategy maintains
agreement of the two kinds of residuals, it turns out that the convergence of the r n may deteriorate
(as we will see, it may result in unacceptably large perturbations to the lanczos recurrence
relation for the residual vectors that steers the convergence, see Section 2.3). Recently, Sleijpen
and van der Vorst [18], motivated by suggestions made by Neumaier (see [11, 18]), introduced a
very sophisticated replacement scheme that includes a so-called flying-restart procedure. It was
demonstrated that this new residual replacement strategy can be very effective in the sense that
it can improve the convergence of the true residuals by several orders of magnitude. For practical
implementations, such a strategy is very useful because it leads to meaningful residuals and this is
important for stopping the iteration process at the right point. Of course, one could, after termination
of the classical process, simply test the true residual, but the risk is that the true residual
stagnated already long before termination, so that much work has been done in vain.
The present paper will follow the very same idea of replacing the computed residual by the true
residual at selected steps, in order to maintain close agreement between the two residuals, but we
propose a simpler strategy so that the replacement is done only when it is necessary and at phases
in the iteration where it is harmless, that is that convergence mechanism for -
r n is not destroyed.
Specifically, we shall present a rigorous error analysis for iterations with residual replacement and
we will propose computable bounds for the deviation between the computed and true residuals.
This will be used to select the replacement phases in the iteration in such a way that the Lanczos
recurrence among -
r n is sufficiently well maintained. For the resulting strategy, it will be shown
that, provided that the computed residuals converge, the true residual will converge to the level
O(u)kAkkxk, the smallest level that one can expect for an approximation.
The paper has been organized as follows. In Section 2, we develop a refined residual replacement
strategy and we discuss some strategies that have been reported by others. We give an error analysis
in Section 3, and we derive some bounds for the deviation to be used in the replacement condition.
We present a complete implementation in Section 4. It turns out that the residual replacement
strategy can easily be incorporated in existing codes. Some numerical examples are reported in
Section 5, and we finish with remarks in Section 6.
The vector norm used in this paper is one of the 1, 2, or 1-norm.
Residual Replacement Strategy
In this section, we develop a replacement strategy that maintains the convergence of the true
residuals. A formal analysis is postponed to the next section. The specific iterative method can
be any of those that fit in the general framework of Algorithm 1. Throughout this paper, we shall
consider only iteration processes for which the computed residual - r n converges to a sufficiently
small level.
As mentioned in Section 1, we follow the basic idea to replace the computed residual -
r m by
the true residual f selected steps We will refer to such
an iteration step as one where residual replacement occurs. Hence, the residual generated at an
arbitrary step n could be either the usual updated residual - r or the true
residual depending on whether replacement has taken place or not at step n. In order
to distinguish the two possible formulations, we denote by r n the residual obtained at step n of the
process with the replacement strategy, that is
With the residual replacement at step m residual deviation is immediately
reduced to
and it can be shown (see Lemma 1 of Section 2.2) that For the
subsequent iterations n ? m, but before the next replacement step, we clearly have that
Therefore, the accumulated deviation before step m has no effect to the deviation after updating
(n ? m). However, in order for such a strategy to succeed, two conditions must be met, namely,
ffl the computed residual r n should preserve the convergence mechanism of the original process
that has been steered by the -
ffl from the last updating step m to the termination step K, the accumulated error \Sigma K
should be small relative to u(jr which is the upperbound for j- m j.
We discuss in the next two subsections how to satisfy these two objectives.
2.1 Maintaining convergence of computed residuals
In order that r n maintains the convergence mechanism of the original updated residuals, it should
preserve the property that gives rise to the convergence of the original - r n . We therefore need to
identify the properties that lead to convergence of the iterative method in finite precision arithmetic.
While this may be different for each individual method, it has been observed for several Krylov
subspace methods (including CG [10, 20], Bi-CG [20], CGS, BiCGSTAB, and BiCGSTAB(') [18]),
that the recurrence r and a similar one for q n is satisfied almost to machine
precision and this small local error is one of the properties behind the convergence of the computed
residuals. Furthermore, the analysis of [20] suggests that convergence is well maintained even when
the recurrence equations are perturbed with perturbations that are significantly greater than the
machine precision. This latter property is the basis for our residual replacement strategy. Therefore,
we briefly discuss this perturbation phenomenon for Bi-CG (or CG), as presented in [20].
Consider the Bi-CG iteration which contains r In
finite
which denote the computed results of r n and p n , respectively, satisfy
the perturbed recurrence
are rounding error terms that can be bounded in terms of u. Combining these
two equations, we obtain the following perturbed matrix equation in a normalized form
r n+1
where T n is an invertible tridiagonal matrix 1 , ff 0
with
ff
We note that (4) is just an equation satisfied by an exact Bi-CG iteration under a perturbation F n .
In particular, detailed bounds on - n and j n will, under some mild assumptions, lead to F n - O(u).
The main result of [20] states that if a sequence -
r n satisfies (4) and Z n+1 has full rank, then we
have
where
n . The case F reduces to the
known theoretical bound for the exact BiCG residuals [1]. Therefore, even when - r n and its exact
counterpart are completely different, their norms are bounded by similar quantities and are usually
comparable. Of course, in both cases, the bounds depend on the quality of the constructed basis.
More importantly, a closer examination of the bound reveals that even if the perturbation F n
is in magnitude much larger than u, the quantities in the bound, and thus k-r n+1 k, may not be
significantly affected. Indeed, in [20] numerical experiments were presented, where relatively large
artificial random perturbations had been injected to the recurrence for r n ; yet it did not significantly
affect the convergence mechanism.
An implication of this analysis is that, regardless of how -
r n is generated but as long as it
satisfies (4), its norm can be bounded by (6). Hence, we can replace -
r n by r
We assume that no breakdowns of the iteration process have occurred
when are not too large relative to kr n k and kr (see (5)), and we may
still expect it to converge in a similar fashion. Indeed, this criterion explains why the residual
replacement strategies like r but do not work always (see Section
2.3). Here, it will be used to determine when it is safe to replace -
r n by r
note that the above discussion is for Bi-CG, but the phenomenon it reveals seems to be valid for
many other methods, especially for those methods that are based on Bi-CG (CGS, BiCGSTAB,
and others).
Now we consider the case that residual replacement is carried out at step m, that is r
It follows from the definition of ffi m and -
r m that
. So, the updated residual r m satisfies
Thus depending on the magnitude of kj 0
k relative to kr m k and kr m\Gamma1 k, the use of r
may result in large perturbations to the recurrence relation. Therefore, a residual replacement
strategy should ensure that the replacement is only done when kj 0
kg is not
too large.
In a typical iteration, as the iteration proceeds, kffi n k, and hence kj 0
increases while k-r n k
decreases. Replacement will reduce ffi n but, in order to maintain the recurrence relation, it should
be carried out before kj 0
becomes too large relative to k-r n k. For this reason, we propose to set a
threshold ffl and carry out a replacement when kj 0
reaches the threshold. To be precise, we
replace the residual at step n by r
We note that, in principle, residual replacement can be carried out for all steps up to where
reaches certain point. However, from the stability point of view, it is preferred to generate the
residual by the recurrence as much as possible, since kj 0
n k is generally bigger than the recurrence
rounding error kj n k (of order u).
2.2 Groupwise solution updating to reduce error accumulations
From the discussions of Section 2.1, we learn that residual replacement should only be carried out
up to certain point. In this subsection, we will discuss how to maintain, after the last replacement,
the deviation at the order of ujAjjx n j, in which case x n is a backward stable solution. Note that,
for any x n , ukAkkx n k is the lowest value one can expect for its residual. This is simply because
even with the exact solution x, both
is the last updating step, which menas that we are in the final phase of the iteration
process, then, because of (3), the deviation at step n ? m is
From our updating condition, we have that kr n k - kj 0
is chosen not too close
to u, kr n k is small and -
m. We now discuss the three different parts of ffi n . The
discussion here is only to motivate the groupwise updating strategy; a more rigorous analysis will
be given in the next section.
we have that \Sigma n
ffl For the / i part, j/
m)ukAkkxk. If large, the accumulation of errors over steps can be significant.
We note that this is the same type of error accumulation in evaluating a sum
of small numbers by direct recursive additions, which can fortunately be corrected through
appropriately grouping the arithmetic operations as
with terms of similar order of magnitude in the same group S i
\Delta. In this way, the rounding errors associated with a large number of
additions inside a group S i is of the magnitude of uS i , which can be much smaller than uS.
The same technique can be adopted for computing x n as
Specifically, the recurrence for x n can be carried out in the following equivalent form
Groupwise Solution Update:
For convergence
End for
Such a groupwise update scheme has been suggested by Neumaier, and it has been worked
out by Sleijpen and van der Vorst (see [18] for both references). By doing so, the error in the
local recurrence is reduced. Indeed, for
(instead of ujx i j). Hence, \Sigma n
In summary, with groupwise updating of the approximated solution, all three parts of ffi n can
be maintained at the level of ujAjjxj. We mention that groupwise updating can also be used to
obtain better performance of a code for modern architectures, because it allows for level-3 BLAS
operations. This has been suggested in [21, page 52, note 5].
2.3 Some other residual replacement strategies
We briefly comment on some other residual replacement strategies.
For the naive strategy of "replacing always" (the residuals are computed always as b\GammaAx n ) or for
"periodic replacement" (update periodically at every ' steps), replacement is carried out throughout
the iteration, even when kr n k is very small. This, as we know, may result in large perturbations to
the recurrence equations relative to kr n k, since jj 0
j is at least j- n j - ujAjjx n j, see (7). In that case,
as kr n k decreases, the recurrence relation may be perturbed too much and hence the convergence
property deteriorates. This is the typical behaviour observed in such implementations.
We note that if - n can be made to decrease as kr n k does, then replacement can be carried out at
later stages of the iterations. This leads to the strategy of "flying-restart" of Sleijpen and van der
Vorst [18], which significantly reduces - n , and hence j 0
n , at a replacement step. In the flying-restart
strategy b is replaced by f at some but not all of the residual replacement steps (say
addition to the residual replacement r The advantage of
this is that, at the flying-restart step n i+1 , the residual is updated by r n i+1
(noting that b / r n i
. Then
which decreases as r n i
and -
decrease. This is the term that determines the perturbation to
the recurrence and can be kept small relative to r n . However, the deviation satisfies
(assuming x n i+1
). Namely, the deviation at each flying-restart step carries forward
to the later steps. Therefore flying-restart should only be used at carefully selected steps where
However, it is not easy to identify a condition to monitor this. It
is also necessary to have two different conditions for the residual replacement and flying-restart.
Fortunately, our discussion in the last two subsections shows that carrying out replacement carefully
at some selected steps, in combination with groupwise update, is usually sufficient. We shall not
pursue the flying-restart idea further in this paper.
Analysis of the Residual Replacement Scheme
In this section, we formally analyze the residual replacement strategy as developed in Section 2.1
(and presented in Algorithm 2 below). In particular, we develop a computable bound for kffi n k and
n k, that can be used for the implementation of the residual replacement condition.
We first summarize residual replacement strategy in the following algorithm, written in a form
that identifies relevant rounding errors for later theoretical analysis.
Algorithm 2: Iterative Method with Residual Replacement:
Given an initial approximation floating point vector);
set -
For convergence
Generate a correction vector q n by the method;
if residual replacement condition (8) holds
else
(denote but not compute x
End for
Note that x n and ffi n are theoretical quantities as defined by the formulas and are not to be
computed. The vectors / due to finite precision arithmetic

At step n of the iterative method, q n is computed in finite precision arithmetic by the algorithm.
However, the rounding errors involved in the computation of q n are irrelevant for the deviation of
the two residuals, which solely depends on the different treatment of q n in the recurrences for r n
and x n .
Throughout this paper, we assume that A is a floating point matrix. Our error analysis is
based on the following standard model for roundoff errors in basic matrix computations [8, p.66]
(all inequalities are componentwise).
where are floating point vectors, N is a constant associated with the matrix-vector
multiplication (for instance, the maximal number of nonzero entries per row of A).
It is easy to show that
Using this, the following lemma, which includes (1) and (2), is obtained.
Lemma 1 The error terms in the computed recurrence of Algorithm 2 are bounded as follows:
For a step at which a residual replacement is carried out:
Proof From (9), we have that j/ n j - uj-x This leads to the bound for
j/ n j: For a residual replacement step, the updated z is x n by definition, that is x
Therefore, The bounds for j n and - n follow similarly.
be the number of step at which a residual replacement is carried out and let
later step, but still before the next replacement step. Then, we have that
and
Proof The first bound follows directly from Lemma 1. For we have that q
Noting that -
Similarly,
We now consider the deviation of the two residuals.
be the number of an iteration step at which residual replacement is carried out
and let n ? m denote a later iteration step, still before the next replacement step. Then, we have
that
Proof At step m, by the definition of xm in Algorithm 2,
z with z being the
updated z-vector and -
Therefore . Hence, for the range of n ? m, and before the next residual replacement
step:
Lemma 2, we obtain the following computable bound on ffi n .
Lemma be the number of an iteration step at which residual replacement is carried out
and let n ? m denote a later iteration step, still before the next replacement step. Then, we have
Proof The bound for kffi m k follows from that for - m , see (14). From Lemma 2 and Lemma 3, it
follows that
which leads to the bound for ffi n in terms of norms.
We note that it is possible to obtain a sharper bound by accumulating the vectors in the bound
for jffi n j. Our experiments do not show any significant advantage of such an approach. We next
consider the perturbation to the recurrence.
Theorem 1 Consider step n of the iteration and let m ! n be the last step before n, at which a
residual replacement is carried out. If replacement is also done at step n, then let x 0
be the computed approximate solution and r 0
the residual. Then the residual r 0
n satisfies the following approximate recurrence
Proof First, in the notation of Alg. 2, x 0
where we have used that b \Gamma Ax . Furthermore, by
Lemma 3,
Also, kAi
Combining these three, and using that r 0
O(u), the
bound on kj 0
n k is obtained as in Lemma 4.
Note that bound (16) is computable at each iteration step. Therefore, we can implement the
residual replacement criterion (8) with this bound instead of kj 0
k. We note that the factor 2 in
the bound comes from the bound for q i in Lemma 2, which is pessimistic since q i -
x i . Therefore,
we can use the following d n as an estimate for kj 0
Hence, we shall use the following residual replacement criterion, that is residual replacement is
done if
With this strategy, the replaced residual vector r n satisfies the recurrence equation (15) with
k. With this property, we consider situations where r n converges. We now discuss
convergence of the true residual.
Theorem 2 Consider Algorithm 2 with the residual replacement criterion (18), and assume that
the algorithm terminates at step be the number of the last
residual replacement iteration step before termination. If
then
Proof From (17), we have dK ? k. Furthermore, at the
termination step, we have kr is the
last updating step, we have for n - m, d n ? fflkr n k as otherwise there would be another residual
replacement after m. That implies kr
~
which is an upper bound for kffi n k (Lemma 4) and ~
where
~
which implies
~
Thus the bound follows from
We add two remarks with respect to this theorem.
Remark 1: If the main condition (19) is satisfied, then the deviation, and hence the true residual,
will remain at the level of uNkAkkxK k at termination. Such an approximate solution is backward
stable and it is best one can expect. The condition suggests that ffl should not be chosen too small.
Otherwise, the replacement strategy will be terminated too early so that the accumulation after the
last replacement might become significant. As can be expected, however, the theoretical condition
is more restrictive than practically necessary and our numerical experience suggests that ffl can be
much smaller than what (19) dictates, without destroying the conclusion of the theorem.
Remark 2: On the other hand, in Section 2.1 we have seen that ffl controls perturbations to
the recurrence of r n , and for this reason it is desirable to choose it as small as possible. In our
experience, there is a large range of ffl around p u that balances the two needs.
Reliable Implementation of Iterative Methods
In this section, we summarize the main results of the previous sections into a complete implemen-
tation. We also address some implementation issues.
It is easy to see from the definition of d n (see (17)) that it increases except at the residual
replacement steps when it is reset to u(NkAkkxm k Our residual replacement strategy
is to reduce d n whenever necessary (as determined by the replacement criterion) so as to keep it
at the level of uNkAkkxK k at termination. With the use of criterion (18), however, there are
situations where the residual replacement is carried out in consecutive steps while d n remains
virtually unchanged, namely when kr n k stays around d n =ffl - uNkAkkx n k=ffl. ?From the stability
point of view, it is preferred not to replace the residuals in such situations. To avoid unnecessary
replacement in such cases, we impose as an additional condition that residual replacement is carried
out only when d n has a nontrivial increase from the dm of the previous replacement step m.
Therefore, we propose d n ? 1:1d m as a condition in addition to (18) for the residual replacement.
The following scheme sketches a complete implementation.
Algorithm 3: Reliable Implementation of Algorithm 1.
Input an initial approximation residual replacement threshold ffl; an estimate of NkAk;
For convergence
Generate a correction vector q n by the Iterative Method;
if d
End for
Remark: In this reliable implementation, we need estimates for N (the maximal number of
nonzero entries per row of A) and kAk. In our experience with sparse matrices, the simple choice
still leads to a practical estimate d n for kffi n k. In any case, we note that precise estimates
are not essential, because the replacement threshold ffl can be adjusted. We also need to choose
this ffl. Our extensive numerical testing (see section 5) suggests that ffl -
p u is a practical criterion.
However, there are examples where this choice leads to stagnating residuals at some unacceptable
level. In such cases, choosing a smaller ffl will regain the convergence to O(u).
The presented implementation requires one extra matrix-vector multiplication when an replacement
is carried out. Since only a few steps with replacement are required, this extra cost is marginal
relative to the other costs. However, some savings can be made by selecting a slightly smaller ffl and
carrying out residual replacement at the step next to the one for which the residual replacement
criterion is satisfied (cf [18]). It also requires one extra vector storage for the groupwise solution up-date
(for z) and computation of a vector norm k-x n k for the update of d n (kr n k is usually computed
in the algorithm for stopping criteria).
5 Numerical Examples
In this section, we present some numerical examples to show how Algorithm 3 works and to demonstrate
its effectiveness. We present our testing results for CG, Bi-CG and CGS. All tests are carried
out in MATLAB on a SUN Sparc-20 workstation, with
In all examples, unless otherwise specified, the replacement threshold ffl is chosen to be 10 \Gamma8 .
kAk1 is explicitly computed and N is set to 1. In Examples 1 and 2, we also compare d n and the
deviation kffi n k.
Example 1: The matrix is a finite-difference discretization on a 64 \Theta 64 grid for
with a homogeneous Dirichlet boundary condition. a(x; y. We apply
CG and Reliable CG (i.e. Alg. 3) to solve this linear system and the convergence results are given
in

Figure

1.
In

Figure

(and similarly in Figures 2 and 3 for the next example), we give in (a) the convergence
history of the (normalized) computed residual for CG (solid line), the (normalized) true residuals
for CG (dashed line) and for reliable CG (dotted line). In (b), we also give the (normalized)
deviations of the two residuals kffi (dash-dotted line) and for reliable
CG (dotted line) and the bound d n for reliable CG (in x-mark).
Example 2: The matrix is a finite-difference discretization on a 64 \Theta 64 grid for the following
convection diffusion equation
with a homogeneous Dirichlet boundary condition. The function f is a constant. We consider Bi-
CG and CGS for solving the linear systems with
The results are shown in Figure 2 for Bi-CG, and in Figure 3 for CGS.
In the above examples, we have observed the following typical convergence behaviour. For
the original implementations, the deviation increases and finally stagnates at some level, which
is exactly where the true residual stagnates, while the computed residual continues to converge.
With the reliable implementations, when the deviation increases to a certain level relative to r n , a
residual replacement is carried out and this reduces the error level. Eventually, the deviation and
hence the true residual arrive at the level of ukAkkxk. We also note that the bound d n captures
the behaviour of kffi n k very closely, although it may be an overestimate for ffi n by a few orders of
magnitude. In all three cases, the final residual norms for the reliable implementation are smaller
than the ones as obtained by the MATLAB function Anb.
Example 3: In this case, we have tested the algorithm for Bi-CG (or CG if symmetric definite)
and CGS, on the Harwell-Boeing collection of sparse matrices [3]. We compare the original imple-
mentations, the reliable implementations and the implementations of Sleijpen and van der Vorst
[18] (based on their replacement criteria (16) and (18)). In Table 1, we give the results for those
matrices for which the computed residuals converge to a level smaller than ukAkkxk so that there
is a deviation of the two residuals. For those cases where b is not given, we choose it such that a

Figure

1: Example 1 (CG) (a): solid - computed residual of CG; dashed - true residual of CG; dotted
true residual of reliable CG; (b): dash-dotted - of CG, dotted - nk of reliable
CG; x - dn of reliable CG
(a) Convergence History
iteration number
normalized
residual
norm
(b) Residual Deviation and Bound
iteration number
deviation
given random vector is the solution. We note that for some matrices, it may take 10n iterations
to achieve that, which is not practical. However, we have included these results in order to show
that even with excessive numbers of iterations, we still arrive at small true residuals eventually. We
list the normalized residuals res attained by the three implementations
and by Gaussian elimination with partial pivoting (MATLAB Anb). We also list the number of
residual replacements (n r ) for our reliable implementations and the number of flying-restart (n f )
and the number of residual replacements (n r ) for the implementations of Sleijpen and van der Vorst
(SvdV). There are two cases for which the computed residuals do not converge to O(u)kbk with the
choice of 8. For those cases, a slightly smaller ffl will recover the stability and the results
are listed in the last row of the table.
We see that in all cases, the reliable implementation reduces the normalized residual to O(u)
and res2 is the smallest among the three implementations, even smaller than MATLAB Anb. The
improvement on the true residual is more apparent in CGS than in Bi-CG (or CG). Except in a
few cases, both the reliable implementation presented here and the implementation of Sleijpen and
van der Vorst work well and are comparable. So the main advantage of the new approach is its
simplicity and an occasional improvement in accuracy.

Figure

2: Example 2 (Bi-CG) (a): solid - computed residual of Bi-CG; dashed - true residual of Bi-
CG; dotted - true residual of reliable Bi-CG; (b): dashed - of Bi-CG, dotted -
of reliable Bi-CG; x - dn of reliable Bi-CG
(a) Convergence History
iteration number
normalized
residual
norm
(b) Residual Deviation and Bound
iteration number
deviation
6 Concluding Remarks
We have presented a new residual replacement scheme for improving the convergence of the true
residuals in finite precision implementations of Krylog subspace iterative methods. By carefully
monitoring the deviation of the computed residual and the true residual and incorporating the
earlier ideas on residual replacement, we obtain a reliable implementation that preserves the convergence
mechanism of the computed residuals, as well as sufficiently small deviations. An error
analysis shows that this approach works under certain conditions, and numerical tests demonstrate
its effectiveness. Comparison with an earlier approach shows that the new scheme is simpler and
easier to implement as an add-on to existing implementations for iterative methods.
We point out that the basis for the present work is the understanding that the convergence
behaviour (of computed residuals) in finite precision arithmetic is preserved under small perturbations
to the recurrence relations. Such a supporting analysis is available for Bi-CG (and CG)
but it is still an empirical observation for most other Krylov subspace methods. It would be
interesting to derive a theoretical analysis confirming this phenomenon for those methods as well.

Acknowledgements

We would like to thank Ms. Lorrita McKnight for assistance in carrying
out the tests on Harwell-Boeing matrices.

Figure

3: Example 2 (CGS) (a): solid - computed residual of CGS; dashed - true residual of CGS;
dotted - true residual of reliable CGS; (b): dashed - of CGS, dotted - nk of
reliable CGS; x - dn of reliable CGS
(a) Convergence History
iteration number
normalized
residual
norm
(b) Residual Deviation and Bound
iteration number
deviation



--R

An Analysis of the Composite Step Biconjugate Gradient Algorithm for Solving nonsymmetric Systems
Templates for the solution of linear systems: Building blocks for iterative methods
Sparse Matrix Test Problems
Conjugate Gradient Methods for Indefinite Systems
A transpose-free quasi-minimal residual algorithm for non-Hermitian linear systems SIAM J
Iterative solutions of linear systems Acta Numerica

Matrix Computations
Behavior of Slightly Perturbed Lanczos and Conjugate-Gradient Recurrences
Estimating the attainable accuracy of recursively computed residual methods

Methods of Conjugate Gradients for solving linear systems
Solution of Systems of Linear Equations by Minimized Iterations
On the convergence rate of the conjugate gradients in presence of rounding errors
Accuracy and effectiveness of the Lanczos algorithm for the Symmetric eigenproblem
Iterative Methods for Sparse Linear Systems PWS Publishing
BICGSTAB(L) for linear equations involving unsymmetric matrices with complex spectrum Electronic Trans.
Reliable updated residuals in hybrid Bi-CG methods Computing <Volume>56</Volume>:<Pages>144-163</Pages> (<Year>1996</Year>)

Analysis of the Finite Precision Bi-Conjugate Gradient algorithm for Nonsymmetric Linear Systems
The performance of FORTRAN implementations for preconditioned conjugate gradients on vector computers

--TR

--CTR
Stefan Rllin , Martin H. Gutknecht, Variations of Zhang's Lanczos-type product method, Applied Numerical Mathematics, v.41 n.1, p.119-133, April 2002
residuals;finite precision;residual replacement;krylov subspace methods