Differences in the Effects of Rounding Errors in Krylov Solvers for Symmetric Indefinite Linear Systems.
The three-term Lanczos process for a symmetric matrix leads to bases for Krylov subspaces of increasing dimension. The Lanczos basis, together with the recurrence coefficients, can be used for the solution of symmetric indefinite linear systems, by solving a reduced system in one way or another. This leads to well-known methods: MINRES (minimal residual), GMRES (generalized minimal residual), and SYMMLQ (symmetric LQ). We will discuss in what way and to what extent these approaches differ in their sensitivity to rounding errors.In our analysis we will assume that the Lanczos basis is generated in exactly the same way for the different methods, and we will not consider the errors in the Lanczos process itself. We will show that the method of solution may lead, under certain circumstances, to large additional errors, which are not corrected by continuing the iteration process.Our findings are supported and illustrated by numerical examples.
Introduction
We will consider iterative methods for the construction of approximate solutions, starting
with for the linear system A an n by n symmetric matrix, in the
k-dimensional Krylov subspace
with r
With the standard 3-term Lanczos process, we generate an orthonormal basis v 1
process can be recast in matrix formulation
as
in which V j is defined as the n by j matrix with columns v 1
tridiagonal matrix.
Mathematical Institute, Utrecht University, Budapestlaan 6, Utrecht, the Netherlands.
y Institute of Mathematics, Medical University of L-ubeck, Wallstra-e 40, 23560 L-ubeck, Germany.
E-mail: sleijpen@math.uu.nl, vorst@math.uu.nl, modersitzki@informatik.mu-luebeck.de
This assumption does not mean a loss of generality, since the case x0 6= 0 can be reduced to this by a simple
shift
Paige [11] has shown that in finite precision arithmetic, the Lanczos process can be implemented
so that the computed V k+1 and T k satisfy
with, under mild conditions for k,
(u is the machine precision, m 1 denotes the maximum number of nonzeros in any row of A).
we obtain the convenient expression
Popular Krylov subspace methods for symmetric linear systems can be derived with formula
(1) as a starting point: MINRES, GMRES, 2 and SYMMLQ. The matrix T k can be interpreted
as the restriction of A with respect to the Krylov subspace, and the main idea behind these
Krylov solution methods is that the given system replaced by a smaller system
with T k over the Krylov subspace. This reduced system is solved - implicitly or explicitly
- in a convenient way and the solution is transformed with V k to a solution in the original
n-dimensional space. The main differences between the methods are due to a different way of
solution of the reduced system and to differences in the backtransformation to an approximate
solution of the original system. We will describe these differences in relevant detail in coming
sections.
Of course, these methods have been derived assuming exact arithmetic, for instance, the
generating formulas are all based on an exact orthogonal basis for the Krylov subspace. In
reality, however, we have to compute this basis, as well as all other quantities in the methods,
and then it is of importance to know how the generating formulas behave in finite precision
arithmetic. The errors in the underlying Lanczos process have been analysed by Paige [11, 12].
It has been proven by Greenbaum and Strakos [8], that rounding errors in the Lanczos process
may have a delaying effect on the convergence of iterative solvers, but do not prevent eventual
convergence in general. Usually, this type of error analysis is on a worst case scenario, and as
a consequence the error bounds are pessimistic. In particular, the error bounds cannot very
well be used to explain differences between these methods, so as we observe them in practical
situations.
In this paper, we propose a different way of analysing these methods, different in the way
that we do not attempt to derive sharper upper bounds, but that we try to derive upper
bounds for relevant differences between these processes in finite precision arithmetic. This will
not help us to understand why any of these methods converges in finite precision, but it will
give us some insight in answering practical questions such as:
- When and why is MINRES less accurate than SYMMLQ? This question was already posed
in the original publication [14], but the answer in [14, p.625] is largely speculative.
- Is MINRES suspect for ill-conditioned systems, because of the minimal residual approach
(see [14, p.619])? Although hints are given for the reasons of inaccuracies in MINRES, for
MINRES, it is also stated in [14, p. 625] that it is not as accurate as SYMMLQ for the reason
2 GMRES has been designed in combination with Arnoldi's method for unsymmetric systems, but for symmetric
systems Arnoldi's method and Lanczos' method lead, in exact arithmetic, to the same relation (1)
that the minimal residual method is suspect. In [3, p. 43] an explicit relation is suggested
between MINRES and working with A 2 , and it is argued that for that reason sensitivity to
rounding errors of the solution depends on - 2 (A) 2 (it is even stated: 'the squared condition
number of A 2 ', implying - 2 which seems to be a mistake).
- Why and when is SYMMLQ slower than for instance MINRES or GMRES?
- Why does MINRES sometimes lead to rather large residuals, whereas the error in the approximation
is significantly smaller? See, for instance observations on this, made in [14, p.626].
Most important, understanding the differences between these methods will help us in making
a choice.
We will now briefly characterize the different methods in our investigation:
1. MINRES [14]: determine x is minimal. This
minimization leads to a small system with T k , and the tridiagonal structure of T k is
exploited to get a short recurrence relation for x k . The advantage of this is that only
three vectors from the Krylov subspace have to be saved (in fact, MINRES works with
transformed basis vectors; this will be explained in Section 2.3). For the implementation
of MINRES that we have used, see the Appendix.
2. GMRES [16]: This method also minimizes, for y k 2 R k , the residual kb \Gamma Ax k k 2 .
GMRES was designed for unsymmetric matrices, for which the orthogonalisation of the
Krylov basis is done with Arnoldi's method. This leads to a small upper Hessenberg
system that has to be solved. However, when A is symmetric, then, in exact arithmetic,
the Arnoldi method is equivalent to the Lanczos method (see also [7, p.41]). Although
GMRES is commonly presented with an Arnoldi basis, there are various implementations
of it that differ in finite precision, for instance, with Modified Gram-Schmidt, Classical
Gram-Schmidt, Householder, and other variants. We view Lanczos as one way to obtain
an orthogonal basis, and therefore stick to the name GMRES rather than to introduce
a new and possibly confusing acronym. Due to the way of solution in GMRES, all the
basis vectors v j have to be stored, also when A is symmetric.
3. SYMMLQ [14]: determine x such that the error x
Euclidean length. It may come as a surprise that can be minimized without
knowing x, but this can be accomplished by restricting the choice of x k to AK k (A; r 0 ).
Conjugate Gradient approximations can, if they exist, be computed with little effort from
the SYMMLQ information. In the SYMMLQ implementation suggested in [14] this is
used to terminate iterations either at a SYMMLQ iterate or a Conjugate Gradient iterate,
depending on which one is best. For the implementation of SYMMLQ that we have used,
see the Appendix.
Note that these methods can be carried out with exactly the same basis vectors v j and
tridiagonal matrix T j .
Most of our bounds on perturbations in the solutions at the kth iteration step will be
expressed as bounds for corresponding perturbations to the residual in the kth step, relative
to the norm of an initial residual. Since all these iteration methods construct their search
spaces from residual vector information (that is, they all start with kr 0 k 2 ), and since we make
at least errors in the order of u kbk 2 in the computation of the residuals, we may not expect
perturbations of order less than u- 2 (A)kbk 2 in the iteratively computed solutions. So our
bounds can only be expected to show up in the computed residuals, if the errors are larger
than the error induced by the computation of the residuals itself.
Notations: Quantities associated with n dimensional spaces will be represented in bold face,
like A, and v j . Vectors and matrices on low dimensional subspaces are denoted in normal
mode: T , y. Constants will be denoted by Greek symbols, with the exception that we will use
u to denote the relative machine precision.
The absolute value of a matrix refers to elementwise absolute values, that is
2 Differences in round-off error behaviour between MINRES
and GMRES
2.1 The basic formulas for GMRES and MINRES in exact arithmetic
We will first describe the generic formulas for the iterative methods MINRES and GMRES,
and we will assume exact arithmetic in the derivation of these formulas. Without loss of
generality, we may assume that x
The aim is to minimize kb \Gamma Axk 2 over the Krylov subspace, and since
we see that for minimizing must be the linear least squares solution of the
overdetermined system
In GMRES this system is solved with Givens rotations, which leads to an upper triangular
reduction of
in which R k is k by k upper triangular with bandwidth 3, and Q k
is a
orthonormal columns. Using (6), y k can be solved from
and since x
(R
(R
The parentheses have been included in order to indicate the order of computation. In the
original publication [16], GMRES was proposed for unsymmetric A, in combination with
Arnoldi's method for an orthonormal basis for the Krylov subspace. However, when A is
symmetric then Arnoldi's method is equivalent to Lanczos' method, so that (8) describes
GMRES for symmetric A. The well-known disadvantage of this approach is that we have to
store all columns of V k for the computation of x k .
MINRES follows essentially the same approach as GMRES for the minimization of the
residual, but it exploits the banded structure of R k , in order to get short recurrences for x k ,
and in order to save on memory storage.
Indeed, the computations in the generating formula (8) can be reordered as
z k
For the computation of W
k , it is easy to see that the last column of W k is obtained
from the last two columns of W k\Gamma1 and v k . This makes it possible to update x
to x k with a short recurrence, since z k follows from the kth Givens rotation applied to the
vector (z T
This interpretation leads to MINRES.
We see that MINRES and GMRES both use V k , R k , T k , Q k , and z k , for the computation
of x k . Of course, we are not dictated to compute these items in exactly the same way for the
two methods, but there is no reason to compute them differently. Therefore, we will compare
implementations of GMRES and MINRES that are based on exactly the same items in floating
point finite arithmetic. From now on we will study in what way MINRES and GMRES differ
in finite precision arithmetic, given exactly the same set V k , R k , T k , Q k
, and z k (computed
in finite precision too) for the two different methods. Hence, the differences in finite precision
between GMRES and MINRES are only caused by a different order of computation of the
namely
ffl for GMRES: x
ffl for MINRES: x
z k .
Of course, we could have tried to get upper bounds for all errors made in each process, but
this would most likely not reveal the differences between the two methods. If we want to study
the differences between the two methods then we have to concentrate on the two generating
formulas.
2.2 Error analysis for GMRES
In order to understand the difference between GMRES and MINRES, we have to study the
computational errors in V k
. We will indicate actual computation in floating point
finite precision arithmetic by fl, and the result will be denote by a b. Then, according to [5,
p. 89], in floating point arithmetic the computed solution b y
(R
This implies that b
k
so that apart from second order terms in u
Here
is the exact value based on the computed R k and z k . Then we make also
errors in the computation of x k , that is we compute b x
y k ). With the error bounds for
the matrix vector product [10, p.76], we obtain
with Hence, the error \Deltax that can be attributed to
differences between MINRES and GMRES, has two components
This error leads to a contribution \Deltar k to the residual, that is \Deltar k is that part of r k that can
be attributed to differences between MINRES and GMRES (ignoring O(u 2 )
\Deltar
Note that in finite precision we have that AV , and that, because of (3), the
leads to a contribution of O(u 2 ) in \Deltar k . This is also the case in forthcoming situations
where we replace AV k by V k+1 T k in the derivation of upper bounds for error contributions.
Using the bound in (10) and the bound for \Delta 2 , we get (skipping higher order terms in u)
k
3
Here we have used that k jR k
from [21, Th. 4.2]; see Lemma 5.1 for details). The factor - 2 denotes the condition number
with respect to the Euclidean norm. 3
Note that we could bound kV k+1 k 2 by
which is, because of the local orthogonality of the v j , a crude overestimate. According to [15,
p. 267 (bottom)], it may be more realistic to replace this factor
m, where
m denotes the number of times that a Ritz value of T k has converged to an eigenvalue of A.
When solving a linear system, this value of m is usually very modest, 2 or 3 say.
Finally, we note that
R T
It has been shown in [6] that the matrix T k that has been obtained in finite precision arithmetic,
may interpreted as the exact Lanczos matrix obtained from a matrix e
A in which eigenvalues of
A are replaced by multiplets. Each multiplet contains eigenvalues that differ by O(u) 1
4 from
an original eigenvalue of A. 4 With e
k we denote the orthogonal matrix that generates T k , in
exact arithmetic, from e
A. Hence,
e
A T e
A e
3 We also have used that the computed Q k
are orthogonal matrices, with errors in the order of u, i.e.,
O(u). These O(u)-errors lead to O(u) 2 -errors in (13).
4 This order of difference is pessimistic; factors proportional to (u) 1
2 , or even u, are more likely, but have not
been proved [7, Sect.4.4.2].
so that
oe min (R T
A T e
and
oe (R T
A T e
which implies - 2 (R k
(ignoring errors proportional to mild orders of u).
This finally results in the upper bound for the error in the residual due to the difference
between GMRES and MINRES:
Note that, even if there were only rounding errors in the matrix-vector multiplication, then
the perturbation \Deltax to A \Gamma1 b would have been (in norm) in the order of u This
corresponds to an error kA\Deltaxk 2 - u- 2 (A)kbk 2 in the residual. Therefore, the stability of
GMRES cannot essentially be improved.
2.3 Error analysis for MINRES
The differences in finite precision between MINRES and GMRES are reflected by
z k .
We will first analyze the floating point errors introduced by the computation of the columns
of
k . The jth row w j;: of W k satisfies
w
which means that in floating point finite precision arithmetic we obtain the solution b
w j;: of a
perturbed system:
with
Note that the perturbation term \Delta R j depends on j. This gives b
w
when we combine the relations for
c
with
We may replace c
k in (18), because this leads only to O(u 2 ) errors.
Finally, we make errors in the computation of x k because of finite precision errors in the
multiplication of c
with The errors made in c
k and the error term are the only
errors that can be held responsible for the difference between MINRES and GMRES. Added
together, they lead to the \Deltax k related to MINRES:
and this leads to the following contribution to the MINRES residual:
\Deltar
If we use the bound (18) for \Delta W , and use for other quantities bounds similar as for GMRES,
then we obtain
3
3
Here we have also used the fact that
and, with kV k k F -
k, the expression can be further bounded.
This finally results in the following upper bound for the error contribution in the residual
due to the differences in the implementation between MINRES and GMRES:
3k
We see that the different implementation for MINRES leads to a relative error in the residual
that is proportional to the squared condition number of A, whereas for the GMRES implementation
the difference led to a relative error proportional to the condition number only.
This means that if we plot the residuals for MINRES and GMRES then we may expect to
see differences, more specifically, the difference between the computed residuals for the two
methods may be expected to be in the order of the square of the condition number. As soon
as the computed residual of GMRES gets below u difference may be
visible.
2.4 Discussion
In Fig. 1, we have plotted the residuals obtained for GMRES and MINRES. Our analysis
suggests that there may be a difference between both in the order of the square of the condition
number times machine precision relative to kbk 2 . Of course, the computed residuals reflect all
errors made in both processes, and if all these errors together lead to perturbations in the same
order for MINRES and GMRES, then we will not see much difference. However, as we see, all
the errors in GMRES lead to something proportional to the condition number, and now the
effect of the square of the condition number is clearly visible in the error in the residual for
MINRES.
Our analysis implies that one has to be careful with MINRES when solving linear systems
with an ill-conditioned matrix A, specially when eigenvector components in the solution,
corresponding to small eigenvalues, are important.
The residual norm reduction kr k k 2 =kbk 2 for the exact (but unknown) MINRES residual
can be computed efficiently as a product ae k j js 1 of the sines s k of the Givens
rotations. In MINRES (as well as GMRES) this value ae k is used to measure the reduction of
the residual norm: in practical computations, a residual norm is not often computed explicitly
Convergence history MINRES A=Q'*diag(D)*Q, Q Givens
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)
-2Convergence history MINRES , A=Q'*diag(D)*Q, Q Givens
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)
-2Convergence history GMRES A=Q'*diag(D)*Q, Q Givens
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)
-2Convergence history GMRES , A=Q'*diag(D)*Q, Q Givens
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)

Figure

1. MINRES (top) and GMRES (bottom): solid line dotted line
of the estimated residual norm reduction ae k . The pictures show the results for a positive definite
system (the left pictures) and for a non-definite system (the right pictures). For both examples -2
To be more specific: at the left
and G the Givens rotation in the (1; 30)-plane over an angle of at the right diagonal
G the same Givens rotation as for the left example; in both
examples (and others to come) b is the vector with all coordinates equal to 1, and the relative machine
precision
as the kth floating point approximate. Therefore, it is of interest to know
how much the computed ae k may differ from the exact residual norm reduction. The errors
made in the computation of ae k are of order u and can be neglected. Since the computation of
ae k and of b x k are based on the same inexact Lanczos process, (22) implies that
The situation for GMRES is much better: the difference between ae k and the true residual
reduction for GMRES can be bounded by the quantity in the right hand side of (14). In fact,
as observed at the end of x2.2, except for the moderate constant (3
about the most accurate computation that can be expected.
2.5 Diagonal matrices
Numerical Analysts often carry out experiments for (unpreconditioned) iterative solvers with
diagonal matrices, because, at least in exact arithmetic, the convergence behaviour depends
on the distribution of the eigenvalues and the structure of the matrix plays no role in Krylov
solvers. However, the behaviour of these methods for diagonal systems may be quite different
in finite precision, as we will show now, and, in particular for MINRES, experiments with
diagonal matrices may give a too optimistic view on the behaviour of the method.
Rotating the matrix from diagonal to non-diagonal (i.e., diagonal
and Q orthogonal, instead of A = D) has hardly any influence on the errors in the GMRES
residuals (no results shown here). This is not the case for MINRES: experimental results (cf.
Fig. 2) indicate that the errors in the MINRES residuals for diagonal matrices are of order
(A), as for GMRES. This can be understood as follows.
If we neglect O(u 2 ) terms, then, according to (15), the error, due to the inversion of R k ,
in the jth coordinate of the MINRES-x k is given by
k
When A is diagonal with (j; j)-entry - j , the error in the jth coordinate of the MINRES
residual is equal to (use (1) and (6))
k
k
Therefore, in view of (16), and including the error term for the multiplication with c
(cf. (19)), we have for MINRES applied to a diagonal matrix:
which is the same upper bound as for the errors in the GMRES residuals in (14).
The perturbation matrix \Delta R j
depends on the row index j. Since, in general, \Delta R j
will
be different for each coordinate j, (23) cannot be expected to be correct for non-diagonal
matrices. In fact, if orthogonal matrix, then errors of order
in the jth coordinate of x k can be transferred by Q to an mth coordinate
and may not be damped by a small value j- m j. More precisely, if \Gamma is the maximum size of
the off-diagonal elements of A that "couple" small diagonal elements of A to large ones, then
the error in the MINRES residual will be of order \Gamma
(R
we recover the bound (22).
2.6 The errors in the approximations
In exact arithmetic we have that
. Assuming that, in finite
precision, this also gives about the right order of magnitude, then the errors related to differences
between MINRES and GMRES, for the approximate solutions in (11) and (20) can be
bounded by essentially the same upper bound:
. (3
Convergence history MINRES with A=diag(D)
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)
-2Convergence history MINRES with A=diag(D)
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)

Figure

2. MINRES: solid line dotted line (\Delta \Delta \Delta) log 10 of the estimated
residual norm reduction ae k . The pictures show the results for a positive definite diagonal system (the left picture)
and for a non-definite diagonal system (the right picture). Except for the Givens rotation, the matrices in these
examples are equal to the matrices of the examples in Fig. 1: here
This may come as a surprise since the bound for the error contribution to the residual for
MINRES is proportional to
Based upon our observations for numerical experiments, we think that this can be explained
as follows. The error in the GMRES approximation has mainly large components in the
direction of the small singular vectors of A. These components are relatively reduced by
multiplication with A, and then have less effect to the norm of the residual. On the other
hand the errors in the MINRES approximation are more or less of the same magnitude over
the spectrum of singular values of A and multiplication with A will make error components
associated with larger singular values more dominating in the residual.
We will support our viewpoint by a numerical example. The results in Fig. 3 are obtained
with a positive definite matrix with two tiny eigenvalues. For b we took a random perturbation
of Ay in the order of 0:01: This example mimics the situation
where the right-hand side vector is affected by errors from measurements. The solution x of
the equation has huge components in the direction of the two singular vectors with
smallest singular value. In the other directions x is equal to y plus a perturbation of less than
one percent. The coordinates of the vector y in our example form a parabola, which makes
the effects easier visible.
The convergence history of GMRES and of MINRES (not shown here) for this example
with is comparable to the ones in the left pictures of Fig. 1, but, because of a higher
condition number, the final stagnation of the residual norm in the present example takes place
on a higher level (- 3
Fig. 3 shows the solution x k as computed at the 80th step of GMRES (top pictures) and
of MINRES (bottom pictures); the right pictures show the component of x k orthogonal to the
two singular vectors with smallest singular value, while the left pictures show the complete
x k . Note that kx k k . The curve of the projected GMRES solution (top-right picture)
is a slightly perturbed parabola indeed (the irregularities are due to the perturbation p). The
computational errors from the GMRES process are not visible in this picture: these errors are
-0.050.050.150.25x_{GMRES} proj on span(V(3:n))
sing. vectors V with increasing sing. values
-0.50.51.5x_{MINRES} proj on span(V(3:n))
sing. vectors V with increasing sing. values

Figure

3. The pictures show the solution x of computed with 80 steps of GMRES (top pictures)
and of MINRES (bottom pictures). The ith coordinate of xk (along the vertical axis) is plotted against i
(along the horizontal axis).
sin
0:01. The right
pictures show the component of xk orthogonal to the two singular vectors with smallest singular value, while the
left pictures show the complete xk .
mainly in the direction of the two small singular vectors. In contrast, the irregularities in the
MINRES curve (bottom-right) are almost purely the effect of rounding errors in the MINRES
process.
In SYMMLQ we minimize the norm of x which means that y k is
the solution of the normal equations
This system can be further simplified by exploiting the Lanczos relations (1):
A stable way of solving this set of normal equations is based on an L e
Q decomposition of T T
and this is equivalent to the transpose of the Q k
R k decomposition of T k (see (6)), which is
constructed for GMRES and MINRES:
This leads to
from which the basic generating formula for SYMMLQ is obtained:
with
k . We will further assume that x
The actual implementation of SYMMLQ [14] is based on an update procedure for V k+1 Q k ,
and on a three term recurrence relation for kr
Note that SYMMLQ can be carried out with exactly the same computed values for V k+1 , Q k ,
R k , and r 0 , as for GMRES and MINRES. In fact, there is no good reason for using different
values for each of the algorithms. Therefore, differences because of round-off, between the
three methods, must be attributed to the additional rounding errors made in the evaluation
of the right-hand side of (25).
The largest factor in the upper bound for these additional rounding errors in the construction
of the SYMMLQ approximation x k is caused by the inversion of L k . The multiplication
with
and the assembly of x k , leads to a factor k
k in the upper bound (similar as for MINRES
and GMRES). In order to simplify the much more complicated analysis for SYMMLQ,
we have chosen to study only the effect of the errors introduced by the inversion of L k . The
resulting error \Deltax k is written as
where g k represents the exact solution and b g k is the value obtained in finite precision arithmetic.
We likewise the coordinates of bg k =kr 0 k 2 are denoted by
These coordinates can be written as
manipulation leads to
where
From (25) it follows that
Hence, the error in the SYMMLQ residual r ME
k can be written as
The first term can be treated as in GMRES:
We define
By combining (29), (27), and the definition for t k , we conclude that
and because of the orthogonality of v k and v k+1 , we have that
The computed residual reduction k b t k k 2 is usually used for monitoring the convergence, in a
stopping criterion. In actual computations with SYMMLQ, no residual vectors are computed.
Expression (30) can now be bounded realistically byp
3
Here we have used that k jL k
Hence
3
A straight-forward estimate is
3
which is much larger than the first term in (33). Experiments indicate that k b t
towards 0 (even below the value u- 2 (A)). Below, we will explain why this is to be expected
(cf. (49)). Fig. 4 illustrates that the upper bound in (33), with k b t
Accuracy. From (33) it follows that
3
is the SYMMLQ residual with respect to the computed SYMMLQ approximate and
r k is the SYMMLQ residual for the exact SYMMLQ approximate (for the finite precision
Lanczos). Apparently, assuming that kr k increases, SYMMLQ is rather accurate
since, for any method, errors in the order u should be expected anyway.
Convergence history SYMMLQ A=Q'*diag(D)*Q, Q Givens
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)
-55Convergence history SYMMLQ , A=Q'*diag(D)*Q, Q Givens
log10(|r|)
(solid
line),
log10(|rho|)
(dotted
line)

Figure

4. SYMMLQ: solid line dotted line (\Delta \Delta \Delta) log 10 of the estimated
residual norm reduction k b tk k2 . The pictures show the results for the positive definite system (the left picture)
and for the non-definite system (the right picture) of Fig. 1. Both systems have condition number
Convergence. It is not clear yet whether the convergence of SYMMLQ is insensitive to rounding
errors. This would follow from (33) if both t k and b t k would approach 0. It is unlikely that
will be (much) larger than k b t k k 2 , that is, it is unlikely that the inexact process converges
faster than the process in exact arithmetic. Therefore, when it is observed that k b t k k 2 is small
(of order u- 2 (A)), it may be concluded that the speed of convergence has not been affected
seriously by rounding errors. In experiments, we see that b t k approaches zero if k increases.
For practical applications, assuming that kt k k 2 . k b t k k 2 , it is useful to know that the
computable value k b t k k 2 informs us on the accuracy of the computed approximate and on a
possible loss of speed of convergence. However, it is of interest to know in advance whether
the computed residual reduction will decrease to 0. Moreover, we would like to know whether
course, it is impossible to prove that SYMMLQ will converge for any
symmetric problem: one can easily construct examples for which kr k k 2 will be of order 1 for
any k ! n. But, as we will analyse in the next subsection, the interesting quantities can be
bounded in terms of the MINRES residual. That result will be used in order to show that
the will be relatively unimportant as soon as MINRES has converged to some
degree.
3.1 A relation between SYMMLQ and MINRES residual norms
In this section we will assume exact arithmetic, in particular the Lanczos process is assumed
to be exact. The residuals r MR
k and r ME
k denote the residuals of MINRES and SYMMLQ,
respectively.
The norm of the residual b \Gamma Ax b , with x b the best approximate of x in K k
can be bounded in terms of the norm of the
residual r MR
kr MR
This follows from the observation that r MR
k where x MR
k is from the same subspace
from which the best approximate x b has been selected, and furthermore that kb \Gamma Ax b k 2 -
Unfortunately, SYMMLQ selects its
approximation x k from a different subspace, namely AK k (A; r 0 ). This makes a comparison
less straight forward.
The following lemma will be used for bounding the SYMMLQ error in terms of the MINRES
error. Its proof uses the fact that r MR
connects K k+1
spanned by r MR
k and AK k (A; r 0 ).
Lemma 3.1 For each z 2 K k+1
kr MR
2: (36)
Proof. For simplicity we will assume that x
By construction x ME
z in the space AK k (A; r 0 ). Hence
implies that
By construction we have that x ME
as a consequence:
From Pythogoras' theorem, with (37), we conclude that
and (36) follows by combining this result with (38).
Unfortunately, a combination of (36) with
k and the obvious estimate jff k j kr MR
, from (37) does not lead to a useful result. An interesting result follows from an
upper bound for jff k j that can be obtained from a relation between two consecutive MINRES
residuals and a Lanczos basis vector. This result is formulated in the next theorem.
Theorem 3.2
Proof. We use the relation
r MR
where r CG
k is the kth Conjugate Gradient residual. The scalars s and c represent the Givens
transformation used in the kth step of MINRES. This relation is a special case of the slightly
more general relation between GMRES and FOM residuals, formulated in [2, 22]. For symmetric
A, GMRES is equivalent with MINRES, and FOM is equivalent with CG. Since
r CG
r MR
kr MR
kr MR
and
. Moreover, since
r MR
k .
Therefore, with e ME
kr MR
r MR
kr MR
kr MR
r MR
kr MR
kr MR
and hence
kr MR
kr MR
A combination of (42) and (36) with
k+1 leads to
kr MR
kr MR
With
and using the minimal residual property kr MR
we obtain the following recursive
upper bound from (43):
kr MR
A simple induction argument shows that fi k - k+1 , and the definition of fi k implies
kr MR
which completes the proof.
For our analysis of the additional errors in SYMMLQ, we also need a slightly more general
result, formulated in the next theorem.
Theorem 3.3 Let y.
For the best approximation y ME
k of y in AK k (A; r 0 ), and for y MR
is the best approximation of c in AK k (A; r 0 ), with - k as in (39), we have
kr MR
i-k
kr MR
Proof. The proof comes along the same lines as the proof of Theorem 3.2.
Replace the quantities x and x MR
k by y and y MR
k . Since the y quantities fulfill the same
orthogonality relations, (36) is valid also in the y quantities. This is also the case for the
upper bound for jff k j kr MR
Hence, with e ME
j , we have
kr MR
kr MR
If we define b
we find that
kr MR
which implies (45).
For the relations between SYMMLQ and MINRES we have assumed exact arithmetic, that is
we have assumed an exact Lanczos process as well as an exact solve of the systems with L k .
However, we can exclude the influence of the Lanczos process by applying Theorem 3.2 right
away to a system with a Lanczos matrix Tm and initial residual kr 0 k 2 e 1 . In this setting, we
have, for k ! m, that ([2, 22])
kr MR
with s j the sine in the jth Givens rotation for the QR decomposition of T is the estimated
reduction of the norms of the MINRES residuals. From relation (44) in combination with (31)
we conclude that
Note that inequality (47) is correct for any symmetric tri-diagonal extension e
Tm of T
(47) holds with e
Tm instead of Tm . It has been shown in [6] that there is an extension e
Tm
of which any eigenvalue is in a O(u) 1
4 -neighborhood of some eigenvalue of A, and therefore
in fairly good precision. This leads to our upper bound
In x3.1.1, we will show that
The upper bound in (49) contains a square of the condition number. However, in the interesting
situation where ae k decreases towards 0, the effect of the condition number squared will be
annihilated eventually.
Remark 3.4 Except for the constants 'k
the estimates (48)
and (49), respectively, appear to be sharp (see Fig. 5).
Although the maximal values of the ratio of kt Fig. 5 exhibit slowly growing
behavior, the growth is not of order k 3 . In the proof of (49) (cf. x3.1.1), upper bounds as
in (48) are used in a consecutive number of steps. In view of the irregular convergence of
SYMMLQ, the upper bound (48) will be sharp for at most a few steps. By exploiting this
observation, one can show that a growth of order k 2 , or even less, will be more likely.
versus MINRES
log10
of
quotient of the estimates of the residual norms: SYMMLQ / MINRES
-22(e_k^t(L+Delta)\e1)./rho_k, |Delta|<eps*|L|, eps=2.958e-13
log10
of
perturbations in SYMMLQ

Figure

5. Results for the non-definite matrix with condition number (as in the right pictures) of
Fig. 1 and Fig. 4. The left picture shows log 10 of the ratio k b tk k2 =ae k of the estimated residual norm reduction
of SYMMLQ with the one of MINRES, the right picture models k b tk \Gamma tk k2 =ae k : it shows the log 10 of e T
3.1.1 SYMMLQ recurrences
In this section we derive the upper bound (49).
Suppose that the jth recurrence for the fl i 's is perturbed by a relatively small ffi and all
other recurrence relation are exact:
The resulting perturbed quantities are labeled as e.
Then
For is a multiple of the SYMMLQ residual for the Tm -system (m ?
as in the proof of inequality (48), Theorem 3.2 could be applied for estimating k e t . For
the situation where j 6= 1, Theorem 3.3 can be used.
To be more precise, with we have (in the notation of
Theorem 3.3), for
y
and
ae k
with c j the cosine in the jth Givens rotation. Therefore, by Theorem 3.3,
ae k
For this specific situation, the estimate for fi k in the last paragraph of the proof of Theorem
3.2 can be improved. It can be shown that fi j - 1 if fi k - k\Gammaj . Therefore, the - k+1 in
(54) can be replaced by - k\Gammaj .
A combination of (51) with (54) gives
Using the definition of M j and the recurrence relations for the fl j , we can express
\Gamma' jj
Therefore, from (48), we have that
Hence (cf. (50))
and, with (55), this gives gives
Because the recurrences are linear, the effect of a number of perturbations is the cumulation
of the effects of single perturbations. If each recurrence relation is perturbed as in (50) then the
estimate (49) appears as a cumulation of bounds as in (57). The vector b t k in (49) represents
the result of these successive perturbations due to finite precision arithmetic.
Finally, we will explain that the effect of rounding errors in solving L can be described
as the result of successively perturbed recurrence relations (50), with
First we note that the efl k 's resulting from the perturbation
are the same as those resulting from the perturbation
which means that a perturbation to the second term in the jth recurrence relation can also be
interpreted as a similar perturbation to the first term in the (j \Gamma 1)st recurrence relation.
Now we consider perturbations that are introduced in each recurrence relation due to finite
precision arithmetic errors. Let b actually computed
and this can be rewritten, with different - and - 0 , as
Since the perturbation to the second term in this jth recurrence relation can be interpreted as
a similar perturbation to the first term in the (j \Gamma 1)st recurrence relation (which was already
perturbed with a factor (1 + 3-)), we have that the computed b fl j can be interpreted as the
result of perturbing each leading term with a factor (1
4 Discussion and Conclusions
In Krylov subspace methods there are two main effects of floating point finite precision arithmetic
errors. One effect is that the generated basis for the Krylov subspace deviates from the
exact one. This may lead to a loss of orthogonality of the Lanczos basis vectors, but the main
effect on the iterative solution process is a delay in convergence rather than mis-convergence.
In fact, what happens is that we try to find an approximated solution in a subspace that is
not as optimal, with respect to its dimension, as it could have been.
The other effect is that the determination of the approximation itself is perturbed with rounding
errors, and this is, in our view a serious point of concern; it has been the main theme of
this study. In our study we have restricted ourselves to symmetric indefinite linear systems
b. Before we review our main results, it should be noted that we should expect upper
bounds for relative errors in approximations for x that contain at least the condition number
of A, simply because we can in general not compute Ax k exactly. We have studied the effects
of perturbations to the computed solution through their effect on the residual, because the
residual (or its norm) is often the only information that we get from the process. This residual
information is often obtained in a cheap way from some update procedure, and it is not
uncommon that the updated residual may take values far beyond machine precision (relative
to the initial residual). Our analysis shows that there are limits on the reduction of the true
residual because of errors in the approximated solution.
In view of the fact that we may expect at least a linear factor - 2 (A), when working with
Euclidean norms, GMRES (x2.2) and SYMMLQ (x3) lead to acceptable approximate solutions.
When these methods converge then the relative error in the approximate solution is, apart from
modest factors, bounded by u - 2 (A). SYMMLQ is attractive since it minimizes the norm of
the error, but it does so with respect to A times the Krylov subspace, which may lead to a
delay in convergence with respect to GMRES (or MINRES), by a number of iterations that is
necessary to gain a reduction by in the residual, see Theorem 3.2. For ill-conditioned
systems this may be considerable.
As has been pointed out in [14], the Conjugate Gradient iterates can be constructed with
little effort from SYMMLQ information if the they exist. For indefinite systems the Conjugate
Gradient iterates are well-defined for at least every other iteration step, and they can be used
to terminate the iteration if this is advantageous. However, the Conjugate Gradient process
has no minimization property (as for the positive definite case) when the matrix is indefinite
and so there is no guarantee that any of these iterates will be sufficiently close to the desired
solution before SYMMLQ converges.
For indefinite symmetric systems we see that MINRES may lead to large perturbation
errors: for MINRES the upper bound contains a factor This means that if the
condition number is large, then the methods of choice are GMRES or SYMMLQ. Note that
for the symmetric case, GMRES can be based on the three-term recurrence relation, which
means that the only drawback is the necessity to store all the Lanczos vectors. If storage is at
premium then SYMMLQ is the method of choice.
If the given system is well-conditioned, and if we are not interested in very accurate solu-
tions, then MINRES may be an attractive choice.
Of course, one may combine any of the discussed methods with a variation on iterative
refinement: after stopping the iteration at some approximation x k , we compute the residual
possible in higher precision, and we continue to solve
solution z j of this system is used to correct x . The procedure could be
repeated and eventually this leads to approximations for x so that the relative error in the
residual is in the order of machine precision (for more details on this, see [20]). However, if we
would use MINRES then, after restart, we have to carry out at least a number of iterations
for the reduction by a factor equal to the condition number, in order to arrive at something of
the same quality as GMRES, which may make the method much less effective than GMRES.
For situations where
u, MINRES may be even incapable of getting at a sufficient
reduction for the iterative refinement procedure to converge.
It is common practice, among numerical analysts, to test the convergence behavior of
Krylov subspace solvers for symmetric systems with well-chosen diagonal matrices. This gives
often a quite good impression of what to expect for non-diagonal matrices with the same
spectrum. However, as we have shown in our x2.5, for MINRES this may lead to a too
optimistic picture, since floating point error perturbations with MINRES lead to errors in the
residual (and the approximated solution) that are a factor smaller as for non-diagonal
matrices.



--R

Templates for the solution of linear sys- tems:building blocks for iterative methods
A theoretical comparison of the Arnoldi and GMRES algorithms
A survey of preconditioned iterative methods
Polynomial Based Iteration Methods for Symmetric Linear Systems
Matrix Computations
Behavior of slightly perturbed Lanczos and conjugate-gradient recurrences
Iterative Methods for Solving Linear Systems

Methods of conjugate gradients for solving linear systems
Accuracy and Stability of Numerical Algorithms
analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix
Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem
Approximate solutions and eigenvalue bounds from Krylov subspaces
Solutions of sparse indefinite systems of linear equations
The symmetric eigenvalue problem
GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems

Reliable updated residuals in hybrid Bi-CG methods
Relaxiationsmethoden bester Strategie zur L-osung linearer Gleichungssysteme
Efficient High Accuracy Solutions with GMRES(m)

The superlinear convergence behaviour of GMRES
--TR
GMRES;MINRES;linear systems;iterative methods;SYMMLQ;stability