Accurate conjugate gradient methods for families of shifted systems.
We consider the solution of the linear system real values of . This family of shifted systems arises, for example, in Tikhonov regularization and computations in lattice quantum chromodynamics. For each single shift  this system can be solved using the conjugate gradient method for least squares problems (CGLS). In literature various implementations of the, so-called, multishift CGLS methods have been proposed. These methods are mathematically equivalent to applying the CGLS method to each shifted system separately but they solve all systems simultaneously and require only two matrix-vector products (one by A and one by AT) and two inner products per iteration step. Unfortunately, numerical experiments show that, due to roundoff errors, in some cases these implementations of the multishift CGLS method can only attain an accuracy that depends on the square of condition number of the matrix A. In this paper we will argue that, in the multishift CGLS method, the impact on the attainable accuracy of rounding errors in the Lanczos part of the method is independent of the effect of roundoff errors made in the construction of the iterates. By making suitable design choices for both parts, we derive a new (and efficient) implementation that tries to remove the limitation of previous proposals. A partial roundoff error analysis and various numerical experiments show promising results.
Introduction
In various scientific computations the problem arises to compute solutions to
for various values of #. The matrix I denotes the identity matrix. Krylov subspace methods
are iterative solution methods for solving linear systems. These methods, with zero initial
guess, are characterized by the fact that they construct their approximations in step j from
the so-called j dimensional Krylov subspace defined as K j (A, b) # span{b, Ab . , A j-1 b}.
An important property of Krylov subspaces is that they are shift invariant, that is K
By exploiting this property, Equation (1.1) can be solved for various values of
Mathematical Institute, Utrecht University, P.O. Box 80.010, NL-3508 Utrecht, The Netherlands.
the shift # by constructing a basis for the Krylov subspace only once. This observation has
led to many e#cient implementations of known Krylov subspace methods that can handle
multiple shifts simultaneously. We refer the interested reader for further information and
applications to [4, 6, 11, 9, 17, 10, 22, 8, 21]. The multishift variants, in general, require the
number of matrix-vector products and inner products of the original method applied to a
single system and for the solution of each additional shifted system only a few extra vector
updates are needed.
In this paper we focus on the numerical solution of the system
for several real values of # 0. The solution of this family of systems plays an important role
in Tikhonov regularization [10] and in the computation of the overlap operator in simulations
in quantum chromodynamics [18]. Despite the fact that this system is sensitive to the e#ects
of using computer arithmetic, since A T A+ #I can be ill conditioned, there is overwhelming
numerical evidence that accurate solutions to this shifted system can be obtained by applying
a su#cient number of iterations of the CGLS method, a variant of the CG method designed
for solving least squares problems. Unfortunately, previously proposed implementations of
multishift-type for solving (1.2), i.e., multishift CGLS methods, can in some cases only achieve
a final precision for the shifted systems that depends on the square of the condition number
of the matrix due to roundo# errors. In the present paper, we will demonstrate this
with a numerical example (cf. Section 5.2). The main goal of this paper is to derive a new
implementation that tries to overcome this limitation.
The use of computer arithmetic a#ects the standard CGLS method in two ways: it alters
the convergence properties of the method, and restricts the accuracy of the method that
can eventually be achieved. This is also the case for the multishift versions of the CGLS
method. In this paper we mainly focus on the second aspect. Our implementation tries
to improve previous proposals in this area. The ultimate attainable accuracy of CG-type
methods is often investigated, e.g., [14, 3], by considering the residual gap which is defined as
the di#erence between the recursively computed approximation to residual and the unknown
true residual. For the CGLS method this is fairly straightforward. For the multishift variants
the situation is more complicated and in order to understand the influence of finite precision
arithmetic on the attainable accuracy of the multishift CGLS method, we will discriminate
between rounding errors made in the construction of the basis for the Krylov subspace (i.e.,
the Lanczos part) and the computation of the approximate solution vector from the Krylov
subspace (the inversion part). The subdivision into a Lanczos and inversion part is not
immediately visible in the standard implementations of the CGLS method and is also not
reflected by the analysis of the attainable accuracy through an inspection of the residual gap.
We stress that in the multishift context these two parts are necessarily independent and it
is a priori not clear if it is even possible to develop multishift versions of the CGLS method
that are able to obtain approximations with similar precisions as a direct application of the
CGLS method to each system separately. In this paper we propose an implementation of
the multishift CGLS method by making suitable design choices for the implementation of the
Lanczos part and for the computation of the iterates for the shifted systems. Confidence in
the success of our method will be given by a partial rounding error analysis (for the important
situations that partially, by numerical experiments.
This paper has the following structure. In Section 2 we review the CG method and
its variant for least squares problems (CGLS) and we discuss some well-known results on
their attainable accuracy. An abstract formulation of the multishift CGLS method is given in
Section 3. The implementation of the 'Lanczos part' is the subject of Section 4. We review an
important result by Paige on the finite precision behavior of this method. As an alternative
for computing an orthogonal basis for the Krylov subspace, we propose to use the CGLS
recurrences. Section 5 deals with the influence of rounding errors made in the 'alternative'
Lanczos method on the attainable accuracy of the multishift CGLS method. The topic of
Section 6 is the accurate computation of the iterates for the shifted systems. Finally, we show
by several numerical experiments that, if both main ingredients are chosen properly, we can
achieve high accuracy for the shifted systems.
Conjugate gradient methods
In this section we review two variants of the conjugate gradient method from the paper of
Hestenes and Stiefel [16] and some of their properties. The conjugate gradient method is an
iterative solution method for solving linear systems when the matrix A is symmetric positive
definite. In this method the iterates, x j and their corresponding residuals, r are
computed for using the recurrence relations
with the coe#cients given by
and, initially, (Norms in this paper are Euclidean.) The approximate solution
then follows using the recurrence
In practice nonzero starting vectors are sometimes used but, for future convenience, we
will assume that the initial guess, x 0 , is zero here. A key characterization of the CG method is
that its iterates, x j , minimize the error in the energy norm (that is an A-weighted norm) over
all approximations from the j-th Krylov subspace K j (A, b), see e.g., [16]. As a consequence
of this, the residuals r i for an orthogonal basis for K j (A, b). Another
useful property is the fact that # i r T
which follows from Equation (5.2) in
[16] and the positivity of the # i . This shows, in combination with (2.2), that
In the following proposition we summarize an important result, due to Greenbaum, that helps
to provide insight into the attainable accuracy of the standard conjugate gradient method.
Proposition 2.1 (Greenbaum [14]). The di#erence between the true residual b-Ax k and
the computed residual r k satisfies
with
#x#
# is the unit roundo#, which for double precision computations is in the order of 10 -16 and c
is a constant that depends on the error in the matrix-vector product.
The quantity # k is di#cult to bound in computer arithmetic. However, in exact arithmetic
it immediately follows from (2.3) that # k # 1. The argument to bound this quantity used in
[14] depends on the fact that, also in exact arithmetic, the errors for the CG method are in
Euclidean norm monotonically decreasing [16, Theorem 6.3] and, therefore,
2. (2.4)
Unfortunately, the argument for proving the reduction of the error in Euclidean norm in [16]
uses the orthogonality of the residuals which is in general lost in finite precision computations.
However, using a relation of the CG method with an exact CG method applied to a larger
matrix and specific right-hand side, the author argues in [14, Section 3.1] that (2.4) might
also hold approximately in the finite precision context.
The relevance of Proposition 2.1 for explaining the attainable precision of the conjugate
gradient method depends on the experimental observation that, in practical computations,
the computed vector r k converges to zero or, at least, stagnates when its norm is many orders
of magnitude smaller than the machine precision #. With these assumptions it, therefore,
follows from Proposition 2.1 that for the iterate, x k , we essentially have for su#ciently large
k that
which implies the following bound on the relative error:
For least squares problems the CG method can be directly applied to the normal equations.
Nevertheless, it is common practice to use an alternative, but mathematically equivalent,
variation of CG, known as CGLS, in this case [16, Section 10]. For method
is defined by the following recurrence relations:
with
c j-1
and z computed as in (2.2).
The advantage of this method, compared to applying CG directly to the normal equations,
is that, here the least squares residuals, z are directly available. Furthermore,
it was shown in [14, Section 3.3] and [3], with similar arguments as used in the proof of
Proposition 2.1 for CG, that recurring the residuals for the least squares problem, z j , improves
the attainable accuracy of the method. For future convenience, we state here a result from [3].
Algorithm 1: CGLS implementation for solving (1.2).
Proposition 2.2 (Bj-orck, Elfving and Strako-s [3]). The di#erence between the true least
squares residual, b -Ax k , and the computed least squares residual, z k , satisfies
with
#x#
c is a constant that depends on the error in the matrix-vector product with A and z = b-Ax.
Again, to apply this result, the authors have to bound # k , which can be accomplished
with similar heuristics as used for the CG method. Assuming, furthermore, that the CGLS
method is terminated at a point such that #z k - z# d#A#x#, the authors conclude that
the CGLS method is eventually expected to achieve an approximate solution that is as good
as any forward stable method for solving the least squares problem. We return to this in
more detail in Section 5.2.
The damped least squares problem (1.2) is equivalent to the least squares problem
#I
Hence, the shifted system can be solved, for one shift, by applying the CGLS method to
this augmented matrix and right-hand side. However, due to the special structure of this
matrix and right-hand side some computational work can be saved in the CGLS method, see,
for example, Algorithm 6 in [10]. For convenience of the reader the resulting algorithm is
summarized in Algorithm 1.
3 An abstract formulation of the multishift CGLS method
In the previous section, we reviewed the celebrated conjugate gradient least squares method
from [16]. In exact arithmetic, a mathematical equivalent method can be obtained based on
the Lanczos method, e.g., [12, Chapter 9], applied to the matrix A T A with starting vector
A T b. This method constructs an orthonormal basis for the Krylov subspace K k
After k iterations, the process can summarized by the Lanczos relation:
where the columns v 0 , . , v k of V k+1 , form an orthonormal basis for K k+1 and the
symmetric k-k tridiagonal matrix T k collects the coe#cients computed during the execution
of the Lanczos algorithm. It is often convenient to include also # k-1 into one
tridiagonal matrix T k by adding an additional row, e T
. The vector e k denotes the
k-th standard basis vector, i.e., is the
vector with all components one and, similarly, # 0 is the vector with all components zero. The
dimension of these vectors should be apparent from the context. The approximation in step
k for the corresponding CG process is equal to
and the vector r is the initial residual. The Equations (3.1) and (3.2) together give
us an abstract formulation of the CGLS method.
It follows from (3.1) that the Lanczos relation for the shifted system reads
k . (3.3)
This shows that, if the Lanczos relation is known for one system, the iterates of the CGLS
method for the damped least squares problem (2.8) can be directly computed from this relation
and are, in fact, equal to
A multishift CGLS method constructs the orthonormal basis once and at the same time
computes (3.4) for the required values of #, of course without storing all Lanczos vectors.
Since there are many, mathematically equivalent, ways to compute the Lanczos vectors and
tridiagonal matrix T k and just as many ways to compute the vectors in (3.4), the number
of possible implementations is countless. Two specific implementations are presented in [17]
and [10]. However, as discussed in the introduction, the obtainable precision of the computed
iterates for these implementations can be limited. We propose a new (and e#cient)
implementation by choosing a suitable algorithm for the Lanczos part and for computing the
iterates (3.4) and we will discuss the relationship of previously proposed multishift CG-type
methods to our implementation in the course of this paper.
4 The implementation for the Lanczos part
As discussed in the previous section, the first key ingredient of a multishift CGLS method
is the construction of an orthonormal, or orthogonal, basis for the Krylov subspaces and
the tridiagonal matrix containing the orthogonalization coe#cients. The obvious choice is
to apply the standard Lanczos method, e.g., [12, Algorithm 9.2.1] to the matrix A T A with
starting vector A T b. It is not di#cult to see that this is not an optimal choice in the context
of solving least squares problems. We will give detailed arguments in Section 5 when we
discuss the impact of rounding errors in the Lanczos part on the attainable accuracy of the
multishift CGLS method.
It is well known that the coupled two-term recurrences of the CG method can be used
as an alternative to the application of the Lanczos method for constructing an orthonormal
basis for the Krylov subspace, see [7, 1]. Numerical experiments in these papers showed
that alternative recurrences could in some cases considerably improve the robustness of the
methods. In a similar spirit, the recurrences in (2.6) can be used to build an orthonormal
basis for the Krylov subspace K j b). In this section we consider the use of the CGLS
recurrences as alternative Lanczos-type method which we will refer to as the CGLS-Lanczos
method.
First, a little remark about notational conventions: with R k we denote the n - k matrix
with columns r 0 , . , r k-1 . Similarly, other capitals will be used to group together the corresponding
vectors. Now, the relations in (2.6) can be summarized by the following matrix
formulations
where
. -# k-2#
and J k
lower bidiagonal matrix with 1 and -1 on, respectively, the diagonal and
sub-diagonal. Substitution yields the residual relation:
Apart from some di#erence in scaling, this is precisely the Lanczos relation given in (3.1).
We can obtain the quantities of the Lanczos method in its standard form by using a simple
diagonal scaling # k with diagonal elements
, which shows that
holds with
Since the # j are available in the CGLS method, the recurrences in (2.6) can be used as an
alternative to applying the Lanczos method to A T A with starting vector A T b at virtually the
same cost. 1
We have argued that the CGLS recurrences can be used as an alternative to the standard
Lanczos method and, therefore, the iterates for the shifted systems can be computed using
(3.4) with the quantities given in (4.2). An, equivalent, alternative that avoids the scaling in
(4.2), is to directly compute
with S k defined as the upper k - k block of S k
given in (4.1).
1 Notice that the vectors v j (i.e., the columns of V k ) are plus or minus the Lanczos vectors of the standard
Lanczos method.
4.1 Perturbed Lanczos relations in computer arithmetic
In finite precision computations, the computed Lanczos vectors V k and tridiagonal matrix
do not satisfy the Lanczos relation exactly. For the standard Lanczos method, Paige [19]
proved that, instead, the computed quantities now satisfy a perturbed Lanczos relation. In
this section we give an analogous result for the CGLS-Lanczos method which we need in the
remainder of this paper.
We assume the standard rules for floating point arithmetic with machine precision #, see,
e.g., [12, Section 2.4.2],
Here, a and b are floating point numbers and # stands for any basic operation like addi-
tion, subtraction, multiplication and division. Furthermore, we assume that the errors in the
matrix-vector product of A and A T with some vector y (of appropriate length) are, respec-
tively, bounded by
#y#.
With this notation, the result of Paige, for a certain implementation of the Lanczos
method, applied to the matrix A T A and starting vector A T b is summarized by the following
proposition.
Proposition 4.1 (Paige [19]). The Lanczos method in computer arithmetic results in a
perturbed Lanczos relation
where, ignoring higher order terms,
We now turn our attention to CGLS-Lanczos, the alternative Lanczos type method discussed
in the previous section. Let f c
j denote the perturbation in the computation of c j caused
by the use of computer arithmetic and assume similar notation for the other perturbations.
Using the standard model for floating point arithmetic we get (cf. (2.6) and the initialization
computations assuming that all quantities in step minus one have length zero)
#f z
#f r
Combining this with the notation from the previous section, we find that the quantities computed
by the CGLS recurrences in finite precision arithmetic satisfy the perturbed relations
given by
k , R
k .
Finally, a straightforward sequence of substitutions and multiplications with suitable matrices
yields a perturbed relation which is summarized by the following lemma.
Lemma 4.1. The CGLS recurrences given in (2.6) in finite precision computations lead to a
perturbed relation of the form
where S k is given by (4.1) and the perturbation is of the form
Upper bounds on the norms of the columns of the perturbations on the right are given by
(4.5)-(4.8).
This lemma does not explicitly provide a bound on the norm of the perturbation as in
Proposition 4.1 for the standard Lanczos method. We note that the norm of this perturbation
can be much larger than the norm of the perturbation term for the standard Lanczos method
(taking into account di#erence in scaling). This can be explained from the fact that #p j #r j #
can become very large which typically occurs when the matrix A has a very small isolated
eigenvalue. However, the perturbation term in Lemma 4.1 has an interesting structure. In
the next section we will exploit this special structure to investigate the attainable accuracy
of the multishift CGLS method based on CGLS-Lanczos.
Finally, we conclude this section by giving, for future use, an explicit expression for the
iterate in step k of the CGLS method. Assuming that no roundo# errors are made in (2.2),
this approximation is given by
x CGLS
(R k U
The first expression follows from (2.2) and the fact that J -1
5 The impact of errors in the Lanczos part on multishift CGLS
We discuss the e#ects that rounding errors in the Lanczos part have on the attainable accuracy
of the multishift CGLS method. The first choice that we mentioned for the Lanczos part was
an application of the standard Lanczos method applied to the matrix A T A and starting vector
A T b. This is not a very good choice as will be argued using the same line of arguments as
given in [3, Section 4.2] to explain the failure of the LSCG method. The main observation is
that we have to startup the Lanczos method by computing A T b to obtain the first Lanczos
vector. Since the vector b does not appear in the Lanczos part, it follows that the multishift
CGLS method based on the standard Lanczos method at best computes a solution to
This shows that, in the nearly consistent case and # = 0, the forward error is, at best, much
less than optimal. For more details on this argument consult [3, Section 4.2].
We now focus on the alternative Lanczos procedure, CGLS-Lanczos. We assume that the
multishift iterates are computed as
defined as in (4.1) and, moreover, for the moment no rounding errors are considered
in the computation of this vector, i.e., x #
is assumed to be computed exactly. In the
multishift CGLS method the iterates are computed independently of the Lanczos part. As a
consequence this will generate additional sources of errors that are not present in the CGLS
method while having impact on the attainable accuracy.
At this point, we warn the reader for some unconventional notation. If we apply a matrix
with k columns to an #-vector with # k, then we assume the vector to be expanded with
zeros if necessary (we do the same with other operations and equalities). In the case # > k,
the matrix is assumed to be applied to the first k elements of the vector and the remaining
elements are assumed to be unchanged by the operation.
For future convenience we define
and we note that
Using this relation and the relations given in Section 4.1, the true residual, corresponding to
the normal equations, can be written as
We plug in the expression for the perturbation term given by Lemma 4.1 and we arrive at
our main relation:
with
This expression plays a similar role in this section as Proposition 2.2 plays for the analysis
of the attainable accuracy of the CGLS method. It shows that, if the vector z k eventually
approaches z, and, therefore, the first term becomes very small, then the true residual corresponding
to (1.2) stagnates at a level determined by the three w-vectors in (5.1). To be more
precise, we assume that we terminate, at an iteration step k, such that
#A#x# d#, (5.2)
for some constant d. Although there is no rigorous proof that this condition, eventually,
always can be achieved, there is according to [3] "overwhelming experimental evidence" to
justify this assumption, see also the numerical experiments in [3].
The attained forward error is given by multiplying the expression for the true residual
(5.1) from the left with This gives
#w (1)
(5.
We can get a similar expression for the least squares residual by multiplying from the left by
We notice that the special structure of the perturbation in Lemma 4.1 has
played an important role in the derivation of this expression. In the next section we bound
the size of these w-vectors for the important case that (notice that already for this
simple case, multishift CGLS based on the standard Lanczos method fails). In Section 5.2 we
combine the obtained bounds with (5.3) and compare the attainable accuracy of the CGLS
method and the multishift CGLS method based on the CGLS-Lanczos method. We will,
furthermore, discuss the situation of more general #.
5.1 Bounding the size of the w-vectors for
We bound the size of the vectors w (1)
and w (3)
for the case of # equal to zero. Since
the expression for the first w-vector can be rewritten to
Therefore, the estimates in (4.5)-(4.8) lead to the following bounds on the w-vectors for
w (2)
To further bound these quantities, it turns out to be convenient to use a relation with an
exact conjugate gradient process applied to the 1)-matrix A and right-hand
side b that generates the same matrix S k (and therefore the same coe#cients as in the CGLS
process). We define
with
. (5.4)
The elements of the diagonal scaling # k are chosen such that the k - k upper block of T k is
symmetric. The symmetric matrix A is now defined by taking its upper left block equal to
1)-element is arbitrary but we will assume that the matrix A is positive
definite which requires that the computed #-coe#cients are all positive. The vector b equals
It is interesting to notice that for the matrix A more advanced completions could have
been chosen is in this case of higher dimension). For example, Greenbaum [13] shows that
the matrix T k can be completed such that the eigenvalues of A are in tiny clusters around the
eigenvalues of A. This gives important insight into the finite precision convergence behavior of
the conjugate gradient method. For our purposes the precise completion is not of importance
and we work with a simple k 1-dimensional matrix.
It is not di#cult to see that the exact conjugate gradient method applied to the matrix
A with right-hand side b, indeed, yields our computed matrix S k . The vectors of this hypothetical
conjugate gradient process will be denoted with non-bold characters, in contrast to
the bold characters denoting the computed vectors in the CGLS-Lanczos method. So, for
example,
appended with k +1- j additional zeros, where
is the upper j - j block of A. The following lemma gives some useful relations.
Lemma 5.1. We have that
and
Proof. Let # i and # denote variables such that |# i |, |#
varies). We have, using the standard model of floating point arithmetic (4.3),
From this it follows that
For the second inequality we use that
Recursive application of this estimate leads to
#1/2 +O(#).
Plugging the expression (5.5) into this expression we find (5.6). #
Now we continue with our original goal, which is to bound #w (1)
k #. We have that
where in the last inequality we have used (2.3). For the second quantity that appears in w (1)
Proposition 2.2 shows that
Furthermore, by combining Lemma 5.1 and (4.9), it is not di#cult to see that
#x CGLS
All together we get the following upper bound:
#w (1)
In order to estimate the size of the vector w (3)
, the problem is reduced to bounding:
For this purpose, the following lemma is of use.
Lemma 5.2. Let j < k. We have that
Proof. First observe that #r j #|e T
| and the first inequality follows.
From [16, Theorem 5.3] we know that #p j #r j #r
The value # j is essentially the norm of the minimal residual approximation corresponding
to the approximation x MR
j with only components in the first j elements of its vector, thus
and the second inequality in (5.10) follows. Here, we used that r j # x MR
j and the fact that
the errors for the CG method are in 2-norm monotonically decreasing [16, Theorem 6.3]. #
Notice that the obvious estimates |e T
could have
been used in the proof of the lemma. This would have given the crude bound
k #) 1/2
which is very similar to the estimate (9.13) in [23] for the CG method in finite precision
computations. Since it contains the term #T -1
would not have been su#cient for our
purposes. In the proof we used the fact that a large value of #r j # is canceled against a smaller
1-th element in the vector x k .
Combining these expressions for all j < k, we finally find
#w (3)
5.2 The attainable accuracy of the multishift CGLS method
Using our bounds for the w-vectors and the expression for the true residual (5.1), we are now
in a position to give some discussion on the ultimate attainable accuracy of the multishift
CGLS method based on the CGLS-Lanczos method. Recall that we neglect, for the moment,
rounding errors made in the computation of the x #
The condition number for the least squares problem, e.g., [2, Section 1.4.3], can be written
as
where #(A) is defined as #A#A denotes the pseudo-inverse of the matrix A given
by A . This means that if
# x is computed using a (normwise) backward stable
method, the forward error is essentially bounded by
. (5.12)
This provides us with some idea what we can expect.
In case of # = 0, we have for the first term in (5.3) that
To bound the contribution of w (1)
and w (3)
to the forward error we define
Although in exact arithmetic this expression does not exceed one, its finite precision value
is unknown. In the following we assume that this quantity is modest. We stress, however,
that the quantity #x j #x#, as in Proposition 2.1, is argued to be bounded in finite precision
computations using similar assumptions.
Using (5.8), (5.11) and the fact that #b#z#A#x#, we finally arrive at the following
result.
Theorem 5.1. Assume that (5.2) holds and let For the accuracy of the multishift
CGLS approximation, we have
and for the least squares residual:
Notice that the forward error appears to be smaller than the upper bound on the error of
a normwise backward stable method given in (5.12). However, as remarked in [3], the number
of iterations k to reach this state might also depend on the conditioning of the matrix A.
With a direct application of the CGLS method, it follows from (4.9) that x CGLS
k .
Therefore, the di#erence in attainable accuracy between the multishift CGLS method and a
direct application of the CGLS method is determined by the vector w (3)
. Our analysis clearly
shows that this term is not expected to form an essential di#erence. In addition, it shows
that the residual gap for the CGLS method is not influenced by the third w-vector.
In the preceding, we restricted our attention to the zero shift case. The generalization for
# > 0 is not straightforward. We want to provide some, preliminary, insight into using the
CGLS-Lanczos method for the Lanczos part of the multishift CGLS method. If #, then
we have:
which yields
Notice that only a perturbation of the right-hand side b of size #b# leads to a forward error
of this size and, therefore, the error of the multishift CGLS method is as small as might be
expected.
In

Figure

1 we have plotted upper bounds on the quantities in (5.3) for two very ill-conditioned
test problems from [15] for general values of #. For this purpose, we ran the
CGLS-Lanczos method for 100 iterations and computed upper bounds on the perturbation
terms in Lemma 4.1 using the estimates given in (4.5)-(4.8) evaluated using the computed
quantities in the CGLS-Lanczos method. Furthermore,
as a reference we have included two upper bounds on the error that are smaller than the upper
s
s

Figure

1: Upper bounds, as a function of #, on the quantities
#w (1)
# and
#b# (both dotted). Left picture:
HEAT(100). Right picture: URSELL(100).
bound predicted by the condition number for the least squares problem for the damped least
squares problem (2.8). This picture shows the behavior that we observed for all our test
problems: for small # the error is dominated by the contribution of w (1)
and for very large
values of #, the contribution of w (2)
becomes dominant. The impact of the vector w (3)
is
always relatively small. The upper bound on the contribution of w (2)
increases for increasing #
until in reaches it maximum value after which it appears to decay with a speed proportional to
coincides with the the bound #(A T A+#I)
#b#. In our numerical experience,
this quantity never becomes much larger than the upper bound on the contribution of w (1)
Unfortunately, we do not have precise expressions that bound these quantities from above for
general #.
6 The solution of the shifted systems
The second main component of the multishift CGLS method is the computation of the approximate
solutions x #
k as
or using the equivalent formulation given by (3.4).
We start this section by summarizing an approach that is used at several places in literature
in multishift versions of the (Bi-)CG method based on coupled two-term recurrences, which
also requires the computation (6.1), see [17, 8]. For more details consult these references.
An important observation is that the residuals for shifted systems are colinear for CG-
type methods, that is, there exist constants #
j such that r #
. Writing out the
three-term recurrence of the residuals r #
k (similar to (4.1)) and comparing terms reveals, with
and a recurrences for the iterates and search directions
with initially p #= r 0 and x #= 0.
Scaling
by #
j leads to a version with the same stability properties
j-1 . (6.3)
If a stable method is applied for computing the inverse in (6.1), then this leads to an
approximate solution that is usually su#ciently accurate when dealing with ordinary linear
systems. However, this might not be the case for normal equations, since then a dependence of
the attainable precision on the square of the condition number of A may have been introduced.
Therefore, a point of concern of the approach given by (6.2)-(6.3) is that it implicitly forms
the ill-conditioned tridiagonal matrix S k (cf., (4.1)) in the computation of the #
in (6.2). An
example of the failure of this algorithm in the multishift CGLS context is given at the end of
this section (-6.2).
6.1 An alternative implementation
The goal is to give an alternative algorithm which prevents the formation of the ill-conditioned
In case the CGLS-Lanczos method is used for the Lanczos part, then it is clear
from (4.2), that also the L T DL factorization of T k is directly available. The quotient-di#erence
algorithms introduced by Rutishauser [20], provide a means to construct an L T DL factorization
of the shifted matrix T k +#I directly from the factors of T k . These algorithms construct
the factors, D #
and L #
, in a step-by-step fashion such that
where D #
is diagonal with diagonal elements
and L #
is upper bidiagonal with
diagonal elements one and upper diagonal elements l #
k-2 . If we take the di#erential
form of the stationary qd transformation (dstqds) as presented in [5, Algorithm 4.2], then
we have, with t #
#, the following recurrence relations for computing the elements of the
factors D #
and L #
.
Just as for the conjugate gradient method, the construction of the vector x #
can be accomplished
e#ciently by introducing the auxiliary vectors p #
defined by the relation P #
Starting with #
.
A simple scaling of the p #
denoting #
leads to a
slightly more e#cient, but equally stable, variant
j-1 . (6.4)
The Lanczos part
The inversion part
Algorithm 2: Multishift CGLS implementation for solving families of the form (1.2).
j-1 . (6.5)
We have discussed the two components of our implementation of the multishift CGLS method
which consists of combining (6.4) and (6.5) with the CGLS recurrences in (2.6). The resulting
algorithm is summarized in Alg. 2. Comparing (6.5) and (6.3) reveals that the search
directions
j and the scalars #
coincide (in exact arithmetic). So, the recursions (6.4) can
be seen as a reformulation of the equivalent recursion given in (6.2), which is hopefully more
stable. As a consequence we, furthermore, have that the norm of the residual of the shifted
system (for the normal equation) is equal to #r j #
In Section 5.2, we addressed the influence of rounding errors in the Lanczos part on the
attainable accuracy of the multishift CGLS method based on the CGLS recurrences. We
now discuss the second main source of errors in the multishift method which is the di#erence
between x #
, computed using (6.4) and (6.5), and the exact expression R k In
[5, Section 4.3], Dhillon and Parlett present a roundo# error analysis of the dstqds algorithm
which shows that the outcome of this algorithm is relatively close to the outcome of an exact
transformation applied to factors relatively close to the original input. This is also expected
to hold for our scaled recursions (6.4). Alternatively, this stability can also be understood
by a close inspection of the recursions in (6.4): only elements are added with the same sign
and therefore the computed value of #
is relatively close to the exact value (since there
is no cancellation). For this reason, we expect that the influence of rounding errors in our
recursion for the #
j is small. On the other hand, the three-term recurrence in (6.2) requires
the computation of #
which may be sensitive to roundo# errors due to cancellation
in case the two quantities are relatively close.
Iteration
relative
error
-2Iteration
relative
error

Figure

2: Relative error as function of k for Alg. 1 (solid), Algorithm 6 from [10] (dashed),
multishift CGLS based on CGLS-Lanczos with (6.5) (dash-dot) and (6.3) (dotted) for two
di#erent shifts:
Another issue is the impact of roundo# errors in the vector updates (6.5). The impact of
these errors can be analyzed in detail using the techniques in Section 5.1. If we define x #
k as
the iterate computed in an exact conjugate gradient method applied to the matrix A
with starting vector b, as defined in Section 5.1, then we can show that
Here, we have used, essentially, the same technique as used for bounding the vector w (3)
in
the case that Hence, we do not expect that rounding errors in the computation of
the vector x #
have a significant influence on the attainable accuracy of the multishift CGLS
method.
6.2 A numerical comparison
Numerical experiments suggest that a multishift CGLS method (based on the CGLS-Lanczos
method) combined with (6.2)-(6.3) for the inversion part, is often remarkably accurate and,
in most situations, as accurate as with our alternative given in Alg. 2. Nevertheless, there are
examples where di#erences are clear. We show this for two simple systems. The matrix A has
eigenvalues {1/250, 240, 241, . , 250} (n = 12), the orthonormal eigenvector basis is random
and the right-hand side has equal components in all eigenvector directions. The results for
solving the system (1.2) are presented in Figure 2 for 1. In this picture
we have also presented the results for Alg. 1 as a reference to the other methods. In this case
coupled recurrences (6.4) clearly gives more accurate results than the three-term recurrence
. For larger values of #, and very small values, the di#erences become small.
The implementation of the multishift CGLS method presented in [10, Alg. 6] uses a variant
of the standard Lanczos method for the Lanczos part. This implementation is based on
a three-term recurrence for the least squares residuals and results in a tridiagonal matrix in
standard form. So, even if these changes in the Lanczos part are improving the attainable
accuracy of the method, the accuracy is expected to be limited by the inversion of the tridiagonal
since the ill-conditioned tridiagonal is not given in factorized form. The dashed lines
in

Figure

confirm this.
7 Numerical experiments
In this section, we compare the attainable accuracy of the CGLS method for damped least
squares problems, Alg. 1, to the attainable accuracy of our version of the multishift CGLS
method, Alg. 2. The 'exact' solution, x # , was computed using a singular value decomposition
of the matrix A and we report the relative error given by
#,
is the computed approximation with either method. The number of iterations, k,
was chosen such that the error for the particular method was minimal. The results for various
test problems from [15] are given in Table 1.
Alg. 1 4.9(-13) 5.1(-15) 6.6(-16) 6.5(-16)
Alg. 2 4.8(-13) 5.2(-15) 6.1(-16) 7.0(-16)
Alg. 1 6.7(-14) 2.9(-15) 2.9(-16) 2.6(-16)
Alg. 2 8.7(-14) 3.3(-15) 2.5(-16) 2.7(-16)
Alg. 1 2.2(-13) 3.0(-15) 3.7(-16) 6.7(-16)
Alg. 2 2.7(-13) 3.0(-15) 3.7(-16) 7.3(-16)
Alg. 1 9.2(-13) 1.8(-14) 1.2(-15) 6.7(-16)
Alg.

Table

1: Attained relative errors for various problems and various choices for #.
The results in this table confirm that the proposed implementation of the multishift CGLS
method achieves a comparable accuracy to applying the CGLS method directly to the regularized
system. However, there are a few interesting di#erences between both methods that
occur now and then and which are not apparent from this table. One property of Alg. 1 is
that, for large shifts, the method appears to have the tendency to diverge after reaching its
maximal precision. An interesting observation is that the multishift version of CGLS does not
have this behavior. This is illustrated in the left figure in Figure 3. The computational costs
per step are much lower for the multishift version of CGLS (no matrix-vector multiplication,
no inner products, less vector updates for solving the shifted problems) than the direct application
of CGLS. However, it is remarkable that, in addition, the multishift version sometimes
needs less iteration steps. An example of this is given in the right picture in Figure 3.
Iteration
relative
error
-1Iteration
relative
error

Figure

3: The relative error (as a function of k) of Alg. 1 (solid) and Alg. 2 (dash-dot). Left:
problem FOXGOOD(100) with
8 Summary and outlook
In this paper we have proposed a new implementation of the multishift CGLS method for
solving families of damped least squares problems. This work was motivated by the observation
that, in some cases, previous proposals can only attain a precision that depends on the
square of the condition number. The first key ingredient of our implementation is the use
of coupled recurrences for the construction of an orthogonal basis for the Krylov subspace.
We showed that this leads to a perturbed Lanczos-type relation with a perturbation that has
a desirable structure. The size of this perturbation can be relatively large compared to the
perturbation term of the standard Lanczos process. Nevertheless, using an extensive rounding
error analysis for which exploits the special structure of the perturbation term,
we showed that this alternative Lanczos method is desirable in multishift CGLS methods. A
second advantage of the use of CGLS recurrences for the Lanczos part, is the fact that the
tridiagonal matrix is available in a factorized form. This is important information for the
second key ingredient: the construction of the iterates for the shifted systems. We showed,
by exploiting the factorized form of the tridiagonal using the stationary qd transformation in
[5, Algorithm 4.2], that the iterates can be accurately and e#ciently constructed.
In future work, we plan to extend our analysis for the Lanczos part to the situation of
more general #. This should also for more general # theoretically confirm the suitability
of the CGLS-Lanczos method for the Lanczos part. Another interesting extension of this
paper is to consider the implementation of the multishift CGLS method based on Lanczos
bidiagonalization e.g., [12, Section 9.3.3].
Furthermore, we believe that our analysis provides the key ingredients for analyzing the
advantages of coupled two-term recurrence Lanczos in the QMR method as experimentally
shown in [7]. In the QMR method there is also a clear separation of the Lanczos and inversion
part, as for the multishift CGLS method. It could also help to analyze the advantages of
alternative Lanczos methods in other applications, e.g., [1].

Acknowledgment

. The research of J. van den Eshof was financially supported by the
Dutch scientific organization (NWO) through project 613.002.035. We are thankful to the
two referees. Their remarks helped us improve the presentation of this paper.



--R

A symmetric band Lanczos process based on coupled recurrences and some applications


Arnoldi methods for large Sylvester-like observer matrix equa- tions

Solution of shifted linear systems by quasi-minimal residual iterations
An implementation of the QMR method based on coupled two-term recurrences
BICGSTAB(l) for families of shifted linear systems
Restarted GMRES for shifted linear systems
Fast CG-based methods for Tikhonov-Phillips regularization
Many masses on one stroke: Economic computation of quark
Matrix computations
Behavior of slightly perturbed Lanczos and conjugate-gradient recur- rences

a Matlab package for analysis and solution of discrete ill-posed problems
Methods of conjugate gradients for solving linear systems
Krylov space solvers for shifted linear systems
Numerical challenges in Lattice Quantum Chro- modynamics (Berlin
analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix
Der Quotienten-Di#erenzen-Algorithmus
Restarted full orthogonalization method for shifted linear systems
On the numerical solution of
On error estimation in the conjugate gradient method and why it works in finite precision computations
--TR
An implementation of the QMR method based on coupled two-term recurrences
Matrix computations (3rd ed.)
Estimating the Attainable Accuracy of Recursively Computed Residual Methods
Restarted GMRES for Shifted Linear Systems
Stability of Conjugate Gradient and Lanczos Methods for Linear Least Squares Problems
Fast CG-Based Methods for Tikhonov--Phillips Regularization
A Symmetric Band Lanczos Process Based on Coupled Recurrences and Some Applications
On the Numerical Solution of $(\lambda^2 A Application to Structural Dynamics
finite precision arithmetic;tikhonov regularization;iterative methods;shifted systems;accuracy