Superlinear Convergence of Conjugate Gradients.
We give a theoretical explanation for superlinear convergence behavior observed while solving large symmetric systems of equations using the conjugate gradient method or other Krylov subspace methods. We present a new bound on the relative error after $n$ iterations. This bound is valid in an asymptotic sense when the size $N$ of the system grows together with the number of iterations. The bound depends on the asymptotic eigenvalue distribution and on the ratio $n/N$. Under appropriate conditions we show that the bound is asymptotically sharp.Our findings are related to some recent results concerning asymptotics of discrete orthogonal polynomials. An important tool in our investigations is a constrained energy problem in logarithmic potential theory.The new asymptotic bounds for the rate of convergence are illustrated by discussing Toeplitz systems as well as a model problem stemming from the discretization of the Poisson equation.
Introduction
. The Conjugate Gradient (CG) method is widely used for solving
systems of linear equations with a positive denite symmetric matrix A.
The CG method is popular as an iterative method for large systems, stemming e.g.
from the discretisation of boundary value problems for elliptic PDEs. The rate of
convergence of CG depends on the distribution of the eigenvalues of A. A well-known
upper bound for the error e n in the A-norm after n steps is
where e 0 is the initial error and the condition number is the ratio
of the two extreme eigenvalues of A. In practical situations, this bound is often too
pessimistic, and one observes an increase in the convergence rate as n increases. This
phenomenon is known as superlinear convergence of the CG method. It is the purpose
of this paper to give an explanation for this behavior in an asymptotic sense.
The error bounds are derived from the following polynomial minimization prob-
lem. For any compact set S  R, we dene
p2Pn
where Pn is the set of polynomials p of degree at most n with 1. The standard
convergence analysis of the CG method leads to
Laboratoire d'Analyse Numerique et d'Optimisation, UFR IEEA { M3, UST Lille, F-59655
Villeneuve d'Ascq CEDEX, France, e-mail: bbecker@ano.univ-lille1.fr
y Department of Mathematics, Katholieke Universiteit Leuven, Celestijnenlaan 200 B, B-3001
Leuven, Belgium, e-mail: arno@wis.kuleuven.ac.be
2 B. BECKERMANN AND A. B. J. KUIJLAARS
where (A) is the spectrum of A. The usual way to analyze (1.3) is to include the
spectrum into a 'continuous' compact set S, so that
The quantity En (S) can be estimated using notions from potential theory, since
lim
log En
where gS (z) is the Green function for the complement of S with pole at 1. Thus one
arrives at
as an upper estimate for the error. For example, if one chooses
the Green function evaluated at 0 is known to be
min
min
which leads to the asymptotic estimate
in terms of the condition number which is in agreement with (1.1).
We refer the reader to the survey paper [DTT98] of Driscoll, Toh, and Trefethen for an
excellent account on the interaction between iterative methods in Numerical Linear
Algebra and logarithmic potential theory.
The estimate (1.7) is typically accurate at early stages of the iteration. The
reason for this is that for small n, a polynomial p 2 Pn that is small on (A) has
to be uniformly small on the full interval [ min ; max ] as well. When n gets larger,
however, a better strategy for p is to have some of its zeros very close to some of the
eigenvalues of A, thereby annihilating the value of p at those eigenvalues, while being
uniformly small on a subcontinuum of S only. Then the right-hand side of (1.7) may
become a great overestimation of the error. This eect is the reason for the superlinear
convergence behavior of the CG iteration, observed in practical situations.
As an illustration we look at the case of a matrix A with 100 equally spaced
100g. The error curve computed for this example is the
solid line in Figure 1. See also [DTT98, page 560]. The classical error bound given
by (1.1) with is the straight line in Figure 1. For smaller values of n, the
classical error bound gives an excellent approximation to the actual error. The other
curve (the one with the dots) is the new asymptotic bound for the error that we nd
in Corollary 3.2 below. This curve follows the actual error especially well in the region
of superlinear convergence (for n  40).
The phenomenon of superlinear convergence has been understood for compact
operators, see [Win80, Mor97, Nev93]. Also, the above heuristic for the convergence
behavior of CG for large matrices has been discussed and further analyzed by several
authors [VSVV86, Gre89, SlVS96, DTT98]. To our knowledge, a formula for
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 3
-5log(rel.
iterations
CG error energy norm
Classical bound
Our asymptotic bound
Fig. 1. The CG error curve versus the two upper bounds for the system
solution x, and initial residual r Our new asymptotic
bound is given in formula (3.11).
the relative error improving (1.7) and explaining the superlinear convergence is still
lacking.
Our goal in this paper is to provide a better understanding of the superlinear
convergence of CG iteration, and in particular to explain the form of the error curve
as seen in Figure 1. We will argue that for a large NN matrix A, the error En
in the polynomial minimization problem (1.2) is approximatelyn
log En ((A)) < t
decreasing family of sets, depending
on the distribution of the eigenvalues of A. The sets S() have the following inter-
pretation: S() is the subcontinuum of [ where the optimal polynomial of
degree uniformly small.
From (1.3) and (1.8) we nd the improved estimate
with
Note that  t depends on n, since . As the sets S() are decreasing as
increases, their Green functions g S() (0), evaluated at 0, increase with  . Hence the
numbers  t decrease with increasing n (see also Remark 2.3 below), and this explains
the eect of superlinear convergence.
4 B. BECKERMANN AND A. B. J. KUIJLAARS
The phenomenon of superlinear convergence may also occur for other Krylov
subspace methods applied to a system A is no longer symmetric
positive denite. For instance, for symmetric but not positive denite A one usually
employs iterative methods like MINRES. Also, the method GMRES may be applied
in case of a general matrix A. Supposing that A is diagonalizable, i.e.,
with D a diagonal matrix containing the (possibly complex) eigenvalues of A, the nth
relative residual may be bounded for these methods by
(see, e.g., [Saa96, Proposition 6.15]). In particular, for symmetric or more generally
normal matrices, V is unitary, and thus again we may give bounds for the relative
residual by describing the (asymptotic) behavior of En ((A)). Indeed, for the ease of
presentation our results are stated for real spectra, but they remain equally valid for
complex spectra (see also Remark 2.4 below).
The paper is organized as follows: In x2, we describe the (sequence of) matrices
AN under considerations. We explain the potential-theoretic origin of our sets S(t),
and establish in Theorem 2.1 the estimate (1.8). Under some stronger assumption
concerning the clustering of eigenvalues, we prove in Theorem 2.2 that estimate (1.8)
is sharp. Section 3 contains a description of eigenvalue distributions where our sets
are explicit intervals. Subsequently, we give an analysis of the plot of Figure 1.
In x4 it is shown that our assumptions are valid for a large class of symmetric positive
denite Toeplitz matrices. Our ndings are illustrated by considering some Toeplitz
matrix occurring in time series analysis. The discretized two-dimensional Poisson
equation on a uniform grid is analyzed in x5. Finally, a lemma used in the proof of
Theorem 2.1 is proved in the appendix.
We should mention that our results concerning the two applications above are
more of theoretical nature since in the present paper neither preconditioning nor
nite precision arithmetic are considered. The main aim of this paper is to illustrate
that some recent results in logarithmic potential theory may help to understand better
a classical phenomenon in Numerical Linear Algebra (see also [BeSa98, Kuij99]).
2. The main result. Properly speaking, the concept of superlinear convergence
for the CG method applied to a single linear system does not make sense. Indeed, in
the absence of roundo errors, the iteration will terminate after N steps if N is the
size of the system. Also the notion that the eigenvalues are distributed according to
some continuous distribution is problematic when considering a single matrix.
Therefore we are not going to consider a single matrix A, but instead a sequence
(AN ) N of symmetric positive denite (or more generally invertible symmetric) matri-
ces. The matrix AN has size NN , and we are interested in asymptotics for large N .
These matrices need to have an asymptotic eigenvalue distribution. By this we mean
that there exists a positive Borel measure  with compact support supp() such that
the following condition is satised.
Condition (i) The spectra (AN ) are all contained in a xed compact set S  R,
and for every function f continuous on S we have
lim
Z
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 5
This condition is equivalent to the weak  convergence of the normalized eigenvalue
counting measures N dened by
where   is the unit point mass at , to . As all the AN have spectra contained in
S, the measure  is supported on S. The total mass kk is at most one, and it can be
strictly less than one if the matrices AN have many coinciding eigenvalues. Note that
in the sum in (2.1) each  in (AN ) is taken only once, regardless of its multiplicity
(see also Remark 2.5 below).
For the use of the potential theory in what follows, we need to impose a condition
on . The logarithmic potential of a Borel measure  with compact support is the
function
U
Z
logj
This is a superharmonic function on C taking values in (1;1]. In particular it is
lower semi-continuous. We refer the reader to [Ran95, SaTo97] for detailed accounts
of logarithmic potential theory. Our assumption is the following.
Condition (ii) The logarithmic potential U  of the measure  from (i) is a continuous
real-valued function on C .
The condition (ii) is not very restrictive. For example, if  is absolutely continuous
with respect to Lebesgue measure with a bounded density then (ii) is satised. It
is also satised if the density has only logarithmic-type or power-type singularities
at a nite number of points. On the other hand, condition (ii) is not satised if
has point masses. A consequence of (ii) is that for any measure  satisfying   ,
the potential U  is also continuous. Indeed, U  is lower semi-continuous, and since
continuous and U   lower semi-continuous, U  is also
upper semi-continuous; hence U  is continuous.
There is a third condition we impose on the sequence (AN ) N .
Condition (iii) The limit (2.1) also holds for
Notice that (iii) follows from (i) if 0 62 S, or even if the (in modulus) small eigenvalues
of AN do not approach zero too fast. If (iii) would not hold, then the matrices AN
are too ill-conditioned and the estimate (2.9) given below may very well fail.
In many practical applications, the family (AN ) N of matrices appears as discretizations
of a continuous operator, and then (i){(iii) are natural conditions, see for
instance the discussion in x4 and x5 below.
The sets S(t) that were announced in (1.8) depend only on the asymptotic eigenvalue
distribution . They are determined by the solution of an energy minimization
problem which we describe now.
The logarithmic energy of a Borel measure  with compact real support is the
double integral
Z
U  ()
ZZ
logj
For every t 2 (0; kk), we dene the class
is a Borel probability measure on R :
6 B. BECKERMANN AND A. B. J. KUIJLAARS
and we let  t be the unique measure minimizing the logarithmic energy (2.2) in the
class M(t; ) (compare [Rak96], [DrSa97, Theorem 2.1]). Thus
The minimizer  t depends on t and . The minimization problem (2.3) is a constrained
problem, since measures in M(t; ) are dominated by the constraint =t. It
is known that the minimizer  t is characterized by the following variational conditions
associated with (2.3). There exists a constant F t such that
see [Rak96, Theorem 3] and [DrSa97, Theorem 2.1]. From these variational conditions
one obtains
Finally, the sets S(t) which are crucial in our ndings are dened by
The extremal problem (2.3) has been studied before in connection with the
asymptotic behavior of discrete orthogonal polynomials, see, e.g., [Rak96], [DrSa97],
[KuVA98], [Beck98], [BeSa98], and [Joh99]. In particular, the monic analogue of (1.2)
is covered by these results, that is, the study of
where P
n denotes the class of monic polynomials of degree n. Notice that if S  [0; 1)
(as for instance for symmetric positive denite matrices), then E
n (S) and En (S) are
realized (up to scaling) by the same polynomial, namely the generalized Chebyshev
polynomial.
Our main result is the following.
Theorem 2.1. Let (AN ) N be a sequence of symmetric invertible matrices, AN
of size N N , satisfying the conditions (i), (ii) and (iii) for some measure . Let the
measures  t , the constants F t , and the sets S(t) be dened by (2.3), (2.4){(2.5), and
(2.7), respectively. Then for t 2 (0; kk), we have
lim sup
n=N!tn
log En ((AN
Proof. Let t 2 (0; kk), and let depend on N in such a way that
Our goal is to construct for every large N a polynomial pN in
Pn which is suciently small on (AN ), so as to obtain the estimate (2.9).
We x  > 0, and dene
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 7
Since U  t is a continuous function (cf. the discussion following condition (ii)), the set
K  is closed. It is disjoint from S(t) because of (2.4) and (2.7). Thus (K
By choosing a smaller  if necessary, we may assume that (@K
possible to nd for every large N , a set ZN such that
(3) for all continuous functions f ,
lim
Z
For the proof that this is indeed possible, we refer to Lemma A.1 in the Appendix.
We write for N large,
Y
Then pN 2 Pn by property (1). We want to estimate max 2(AN ) jp N ()j. Let
be such that
Since pN vanishes on (AN ) \ K  by property (2) and the denition (2.13), we have
and the latter is a bounded set. Passing to a subsequence, if necessary, we may assume
that the sequence (N ) converges as N !1 with limit
We have by (2.13),n
log jp N (N )j =n
log
log jj: (2.16)
From property (3), we have that the normalized counting measures of ZN , i.e.,
converge in weak  sense to t t . The principle of descent, see [SaTo97, Theorem I.6.8],
and (2.15) then imply that
U N (N
t, this gives
lim sup
log
8 B. BECKERMANN AND A. B. J. KUIJLAARS
and thus by (2.15)
lim sup
log
The principle of descent also implies
By property (2) we have  N  N , where N is the normalized counting measure of
(AN ). Since N !  by condition (i), we nd that ( N  N ) N is a sequence of
positive measures that converges to  t t in weak  sense. Applying the principle of
descent once more, we obtain
U N N (0): (2.19)
Also the condition (iii) gives
U
The relations (2.18){(2.20) easily imply that
U N (0)
and this is equivalent to
lim
log
Combining (2.16) with (2.17) and (2.21), we obtain
lim sup
log jp N (N )j  F t
By (2.14) and the denition (1.2) of En , we then see
lim sup
En ((AN
The number  > 0 can be chosen arbitrarily close to 0. Hence (2.9) follows.
To obtain (2.10) we need to show that
To establish this and the inclusion property
claimed in the introduction, we recall the connection of the constrained minimization
problem (2.3) with the energy problem in the presence of an external eld. For a
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 9
continuous function Q sucient growth at 1, the logarithmic energy
of a measure  in the presence of the external eld Q is
I
ZZ
log 1
Z
The minimizer  Q
s for the extremal problem
Z
exists and is unique if
log jj (2.24)
and it is characterized by the conditions
U  Q
s
U  Q
for some constant G s , see, e.g., [SaTo97, Theorem I.1.3]. Buyarov and Rakhmanov
[BuRa99, Theorem 2] proved the following formula for  Q
s
where !S is the equilibrium measure for the set S. In fact, the authors consider
external elds where the limit on the right-hand side of (2.24) equals +1, and s 2
(0; +1). However, from their proof (see the last paragraph of [BuRa99, Section 2])
it becomes clear that (2.27) remains valid as long as (2.24) holds.
Now, in our situation with the constraint , we take
By comparing the conditions (2.4){(2.5) with (2.25){(2.26) we can easily check that
for s; t > 0 with s
In particular, ( Q
converges in weak  sense to  for s ! kk. Then the Buyarov{
Z t! S() d: (2.29)
From (2.29) we obtain the inequality    t t for  < t, and thus (2.23) holds. In
order to show (2.22), notice that the Green function is connected with the potential
of the equilibrium measure by the formula
A. B. J. KUIJLAARS
where cap denotes the logarithmic capacity. Combining this with (2.29), we obtain
for
Z tlog cap (S()) d: (2.30)
For  2 S(t), the left-hand side of (2.30) vanishes according to (2.4). Also, by (2.23),
each  2 S(t) belongs to S() for all  < t, so that the integral in (2.30) involving the
Green functions vanishes for  2 S(t). Consequently,
Z tlog cap (S()) d;
and the equation (2.22) follows from (2.30). This completes the proof of Theorem 2.1.
Under additional conditions the inequality (2.9) can be improved to give equality
lim
n=N!tn
log En ((AN
These additional conditions are related to the separation of the eigenvalues. If many
eigenvalues are very close to each other then the inequality (2.9) may be strict. For the
extremal problem (2.8), various separation conditions were considered by Rakhmanov
[Rak96], Dragnev{Sa [DrSa97], Kuijlaars{Van Assche [KuVA98], and Beckermann
[Beck98], see also [KuRa98] for a survey.
If one of these conditions holds in the present situation, the limit (2.31) can be
proved. Indeed, according to Theorem 2.1, we only require a sharp lower bound for
En ((AN )). For sets S with positive capacity, lower bounds for En (S) are usually
obtained by applying the Bernstein-Walsh inequality. In our discrete setting, some
analogue of the Bernstein-Walsh inequality in terms of the extremal measure  t exists,
see [KuVA98, Lemma 8.1 and Corollary 8.2] and [Beck98, Theorem 1.4(c)], implying
(2.31).
We will give here a proof using the separation condition of Beckermann [Beck98],
which was rst conjectured by Rakhmanov [KuRa98]. For a nite subset Z  C , we
introduce
I  (Z) :=(#Z) 2
logj
which may be thought of as the discrete energy of a system of #Z particles each
having a charge 1=#Z. Beckermann's condition is:
Condition (iv) With I() as in (2.2) we have
lim
I  ((AN
It can be shown that lim inf I  ((AN
already follows from condition (i).
Notice also that the separation conditions of [Rak96, DrSa97] imply (iv).
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 11
It is not dicult to prove using for instance [Beck98, Lemma 2.2(b)] that condition
(iv) is equivalent to the fact that
I  (XN
whenever (XN ) is a sequence of sets satisfying XN  (AN ) for each N , lim(#XN )=N
In this form the condition (iv) will be used.
Theorem 2.2. Suppose that the assumptions of Theorem 2.1 are satised and
that in addition the condition (iv) holds. Then for every t 2 (0; kk), the limit (2.31)
holds.
Proof. Let t 2 (0; kk). As in the proof of Theorem 2.1 we assume that
depends on N in such a way that n=N ! t.
For every N 2 N, let N be a set of n points in (AN ). That is, N
has denoted by  maximizes the product
Y
among all n 1-point subsets of (AN ). Equivalently, N minimizes the discrete
energy I  (N ). Since N  (AN ), it is clear that
En ((AN
Our rst goal is to show that the normalized counting measures of the Fekete
points tend to  t , that isn
as N ! 1. Since the sets N are all contained in the compact S, Helly's theorem
asserts that from any subsequence of the sequence of normalized counting measures
of the Fekete points, we may extract a further subsequence having a weak  limit
(which clearly is an element of M(t; )). Our claim (2.35) follows by showing that
. According to (2.33), we nd that along an appropriate subsequence we then
have
I  (N
Let (ZN ) N be a sequence of sets satisfying
It follows from Lemma A.1 in the Appendix that such a sequence exists. Again by
(2.33), we nd that
I  (ZN
12 B. BECKERMANN AND A. B. J. KUIJLAARS
by (2.3), and I  (ZN )  I  (N ) by the denition of Fekete
points, we may conclude that I(  by the uniqueness of
the minimizer in (2.3). This proves the claim (2.35).
Next, we dene for N 2 N and
Y
Then P k;N has degree n, and any polynomial p 2 Pn can be written in the form
a k
with coecients a k satisfying
k=0 a
and
a k
P k;N (0)
P k;N (0)
ng be such that it maximizes
P k;N (0)
among all k 2 f0; ng. Then it follows from (2.38) that
Since this holds for every p 2 Pn , we nd
We write shorter
~
with normalized counting measure
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 13
Because of (2.35) we see that (N ) has the weak  limit  t for N ! 1. From the
principle of descent [SaTo97, Theorem I.6.8] it follows that
lim sup
log jP kN ;N
Also, we will show below that
lim inf
log jP kN ;N ( kN ;N )j  F t (2.41)
(compare with [Beck98, Lemma 2.6]). Combining (2.40) and (2.41) with (2.34) and
(2.39), we may conclude that
lim inf
log En ((AN
log
which in view of Theorem 2.1 is the inequality required for the proof of Theorem 2.2.
Finally, in order to establish (2.41), we note that by the denition of Fekete points
we have for every  2 (AN ),
Taking logarithms, and adding the inequalities for  2 (AN
N , we obtain
log jP kN ;N ()j
and therefore
log jP kN ;N ( kN ;N )j  1
log 1
One easily veries that the right-hand of (2.42) side equals2
which according to (2.33) converges to2t 2
Z
where for the last equality we have used the variational condition (2.4). Since
(#(AN assertion (2.41) follows from (2.42), and Theorem 2.2
is proved.
Remark 2.3. We have shown in Theorem 2.1 that, for n; N !1, the quantity
log En ((AN )) is asymptotically bounded above by
log  n
14 B. BECKERMANN AND A. B. J. KUIJLAARS
and this bound is sharp (under some additional assumptions) according to Theorem
2.2. This conrms our claims (1.9) and (1.10) of the introduction. The graph
of Nt log  t for xed N and varying is drawn in the plots of Figures 1, 2,
and 4. From (2.43) one sees that N t log  t is dierentiable, up to at most a countable
number of points, with derivative
d
Thus it follows that (2.43) is decreasing. Also because of (2.23) one sees that g S(n=N) (0)
is increasing with n, and therefore the graph of (2.43) is concave.
If S is a compact set containing all the spectra (AN ), then S(t)  S, for every
one easily checks that
In other words, the bound (1.9) is sharper than (1.6). The equality
holds if and only if which again is true if and only if the
equilibrium distribution !S of S is less then or equal to =t. This may be translated
by saying that, roughly, about tN out of the eigenvalues of AN are asymptotically
distributed like the equilibrium distribution of S.
Remark 2.4. Theorems 2.1 and 2.2 are equally valid for complex discrete
sets (AN ), here supp() may be a subset of the complex plane. Indeed, the energy
problems with constraint have been studied in a complex setting (see, e.g., [DrSa97]),
and it is possible to show that the representation of F t U  t in terms of Green
functions remains equally true. Furthermore, all other arguments used in the proofs
of Theorems 2.1 and 2.2 still apply for complex sets (AN ). As a consequence, our
Theorems can also be used for bounding the relative residual while solving systems of
linear equations with normal matrices AN via methods like MINRES or GMRES.
Remark 2.5. In many applications (as for instance for symmetric Toeplitz ma-
trices) it is dicult to know in advance the multiplicities of the eigenvalues, and one
only obtains a measure ~
dened by a modication of condition (i) where multiple
eigenvalues are counted according their multiplicities. We will refer to this modi-
cation as condition (i'). Condition (ii) with this (possibly) new measure ~
will be
called (ii'), and accordingly (iii) becomes (iii'), where again we count multiplicities.
Notice that Theorem 2.1 remains valid if assumptions (i),(ii),(iii) are replaced by
(i'),(ii'),(iii') (and  is replaced by the new measure ~ ).
In case of, e.g., real S, conditions (i') and (iii') have an interesting interpretation
in terms of asymptotics of determinants: Indeed,N
log jdet(I N AN
log 1
(in this formula we count multiplicities), and from logarithmic potential theory we
know that relation (2.1) holds for every function f continuous on S if and only if
lim
log jdet(I N AN
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 15
for all  2 C n S. Furthermore, it is sucient that (2.44) holds for  2  where
C has a nite accumulation point outside of S. Notice that condition (iii') may
be rewritten as (2.44) with Finally, condition (i') is known to hold i
lim
trace
for all  2 ,   C as above. Using (2.45), one can for instance easily show that
condition (i') remains valid (with the same measure) if AN is perturbed by some
matrix BN , with sup N kBN k < 1, and rank(BN )=N ! 0.
3. Equidistant eigenvalues. In order to apply Theorems 2.1 and 2.2 we have
to calculate the sets S(t) from the eigenvalue distribution . This is a problem in
itself. In general, the sets S(t) can have a complicated form. They may consist of
several intervals, or even have a Cantor-like structure. The easiest case would be
if all S(t) are single intervals. This would also be the most convenient case for the
computation of the Green function at 0, since for an interval [a; b], we have
a
a
Lemma 3.1. Suppose that  is supported on the interval [a; b] and has a density
w() with respect to Lebesgue measure. We write
~
w() :=
(a) Suppose ~
w is strictly increasing on (a; b). Then S(t) is an interval containing
b for every t 2 (0; kk). More precisely, we have
and
is the unique solution in (a; b) of the equation
Z r
a
r
r
w() d: (3.2)
(b) Suppose ~
w is strictly decreasing on (a; b). Then S(t) is an interval containing
a for every t 2 (0; kk). More precisely, we have
w(b
and
w(b
is the unique solution in (a; b) of the equation
r
r
a
r
w() d:
B. BECKERMANN AND A. B. J. KUIJLAARS
(c) Suppose ~
w is symmetric with respect to the midpoint m := (a + b)=2 and
strictly decreasing on (m; b). Let t 2 (0; kk). Then
w(b
and
w(b
is the unique solution in (0; (b a)=2) of the equation
a
w()d:
Proof. (a) We consider as in the proof of Theorem 2.1 the external eld
a
log
Let  Q
s be the extremal measure with external eld Q and normalization s 2 (0; kk)
(cf. the paragraph preceding formula (2.26)). In [KuDr99, Theorem 2] it was proved
that the support of  Q
s is an interval of the form [r; b] if Q and w are related as in (3.3)
and if ~
w() increases on [a; b]. In [KuDr99] this is stated under the assumption that
Q is dierentiable with a Holder continuous derivative. An inspection of the proof,
however, shows that this assumption is not necessary. It was also assumed that s = 1.
This is also not essential. Since
s by (2.28), it thus
follows that S(t) is an interval containing b for every t.
We show that
w(a+). For t  ~
w(a+), we have from
the fact that ~
w is strictly increasing,
(b )( a)
Thus the equilibrium measure ! [a;b] of [a; b], i.e.,
(b )( a)
d
belongs to the class M(; t). Since ! [a;b] minimizes the energy (2.2) among all probability
measures on [a; b], it is then also the minimizer over M(; t). Thus
and it follows from (3.4) that
Conversely, if is a probability measure on [a; b] whose potential
is constant on [a; b] by (2.4). This implies that  Hence t! [a;b]   and (3.4)
holds. Thus t  ~
w(a+).
For the rest of the proof we assume that t 2 ( ~
w(a+); kk). Then
with a < r(t) < b. From [DrSa97, Corollary 2.15] we know that t
is the balayage (see [SaTo97, Section II.4]) of  onto the interval [r; b].
Consequently, according to [SaTo97, Eqn. (II.4.47)], t t has the density
w() if  2 (a; r),
Z r
a
s
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 17
For  2 (r; b), we rewrite the density as
Z r
a
Since a < r < b and 0  v()  w() < 1 for  2 (a; b), we must have
lim
In view of (3.5), the relation (3.2) follows.
To show that there is only one r satisfying (3.2), we rewrite the right-hand side
of (3.2) as
Z r
a
r
r
Z r
a
~
w()
d
Z =2~
where for the second equality, we used the change of variables
Since ~
w is strictly increasing, it is then clear that (3.6) strictly increases for r 2 (a; b).
This completes the proof of part (a).
(b) The proof of part (b) is similar.
(c) Part (c) follows using a quadratic transformation, compare with [Kuij99, Proof
of Theorem 5.1].
Lemma 3.1 allows us to determine the sets S(t) in a number of situations. We
consider here the case of equidistant eigenvalues.
Suppose AN has N equidistant eigenvalues . Multiplying the matrix
by a positive constant does not change the numbers En ((AN )). We multiply AN by
1=N and so we consider instead matrices with spectrum
These matrices have an asymptotic eigenvalue distribution
and the conditions (i){(iv) are satised.
The explicit solution of the energy minimization problem (2.3) with  given by
(3.7) is due to Rakhmanov [Rak96]. We show how the sets S(t) can be determined
from Lemma 3.1. The assumptions of Lemma 3.1(c) are clearly satised with a = 0,
Therefore we have for 0 < t < 1,
d
=1=2 r
B. BECKERMANN AND A. B. J. KUIJLAARS
Thus
and
11p
Using (3.1) and (3.8) we nd after a little calculation
log
Hencet
Z tlog
d
Theorem 2.2 and (3.10) then give the following result.
Corollary 3.2. For every t 2 (0; 1) we have
lim
n=N!tn
log En (f1;
Corollary 3.2 gives the theoretical justication for our CG bound in the case of
equidistant eigenvalues as given in Figure 1. Notice that, already for 100, the
approximation for log En (f1; obtained by multiplying the right-hand side
of (3.11) by n is quite accurate.
To conclude this section, we note that Lemma 3.1(c) also applies to the case of
ultraspherical eigenvalue distributions. The corresponding sets S(t) were determined
in [Kuij99].
4. Applications to Toeplitz matrices. Toeplitz matrices provide interesting
examples for our results. Toeplitz systems arise in a variety of applications, such
as signal processing and time series analysis, see [ChNg96] and the references cited
therein.
be an integrable function with Fourier coecients
Z
()e ik d;
We assume  is bounded and not equal to a constant. The Nth Toeplitz matrix with
symbol  is given by
TN
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 19
Then TN () is a Hermitian matrix and it is well known that ( inf ;  sup ) is the smallest
interval containing the spectrum of TN () for every N , where  inf and  sup denote
the essential inmum and essential supremum of , respectively. Thus, since  is
non-negative, all eigenvalues  1;N   2;N      N;N of TN () are strictly positive,
and the matrix TN () is positive denite.
A classical result of Szeg}o, see e.g. [GrSz84, pp. 63-65], [BoSi99, Theorem 5.10
and Corollary 5.11], says that
lim
Z
for every continuous function f on [ min ; max ]. It follows that the sequence (TN ())
satises the condition (i') (see Remark 2.5) with the measure  given by
Z
Z
Then  is a probability measure and its support is equal to the essential range of .
Another result of Szeg}o (see [Sz67, Eqn. (12.3.3)] or [GrSz84, p. 44 and p. 66]) is
that
lim
Z
log () d
provided that  satises the Szeg}o condition
Z
log
Notice that this condition can be rewritten as U  (0) < +1. It follows from (4.2),
(4.3) that
lim
log det TN ()
Z
log
Z
log  d() 2 R;
and the condition (iii') is satised.
Consequently, for Toeplitz matrices TN () with non-negative, integrable, and
bounded symbol  and continuous real-valued potential U  , the conditions (i')-(iii')
are satised, and we may apply Theorem 2.1.
We will discuss an example of Kac, Murdock and Szeg}o [KaMuSz53, p. 783]
with
1). Toeplitz matrices with this symbol (or with a multiple of this
arise as covariance matrices of rst-order autoregressive processes [ChNg96,
Section 4.6.1]. The corresponding Fourier coecients are given by
Suppose without loss of generality that
> 0. Then the measure  from (4.2) has
support [a; b] where
20 B. BECKERMANN AND A. B. J. KUIJLAARS
Since  is even we have2
Z
Making the substitution  = (), we obtain after some calculations2
Z
a
d
Thus the measure  has density
a <  < b: (4.4)
with respect to Lebesgue measure. From (4.4) it is easy to show that U  is continuous,
so that Theorem 2.1 applies in this case.
Now we apply Lemma 3.1(b) in order to compute S(t). Notice that for r 2 [a; b),
r
r
a
r
r
d
r a
r
Consequently, by Lemma 3.1(b), we have
r a
a
if a < t < 1:
In particular, we get from (1.10) and (3.1) the convergence rate
log
whereas for a < t < 1, we have
t log
a
log
a
a
d
a log(
a
log
d
a log(
a
log
It is quite interesting that, in the superlinear range, we obtain (up to some linear
transformation) the same function as for equidistant nodes.
Numerical experiments for the symmetric positive denite Toeplitz matrix T 200 of
order 200 of Kac, Murdock and Szeg}o are given in Figure 2. The four dierent plots
correspond to the choices
19=20g of the parameter. Notice that the
CG error curve (solid line) of the last two plots is clearly aected by rounding errors
leading to loss of orthogonality, whereas the GMRES relative residual curves (dotted
line) behave essentially like predicted by our theory. In particular, the classical bound
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 21
log(
rel.
residual
log(
rel.
residual
iterations
iterations
Fig. 2. The error curve of CG (solid line) and GMRES (dotted line) versus the classical upper
bound (crosses) and our asymptotic upper bound (circles) for the system T 200
solution x, and initial residual r Here TN is the Kac, Murdock and Szeg}o matrix of
x4, with parameter
19=20g.
(1.1), (1.11) (crosses) does no longer describe correctly the size of the relative residual
of GMRES for n  20 and
19=20g. Experimentally we observe that the
range of superlinear convergence starts in the dierent examples approximately at
the iteration indices  50, 30, 20, and 10, respectively. This has to compared with
the predicted quantity N  a which for the dierent choices of
approximately takes
the values 66, 40, 29, and 5, respectively. Though theses numbers dier slightly,
we observe that the new bound (1.9), (1.10) re
ects quite precisely the shape of the
relative residual curve, and in particular allows to detect the ranges of linear and of
superlinear convergence.
Let us nally mention the Toeplitz matrices occuring in the context of the rst-
order moving average process [ChNg96, Section 4.6.1], where the symbol is given by
Here the eigenvalues are asymptotically distributed like the equilibrium distribution
on
and therefore there will be no superlinear convergence in
22 B. BECKERMANN AND A. B. J. KUIJLAARS
this case.
5. The Model Problem. Consider the two dimensional Poisson equation
for (x; y) in the unit square 0 < x; y < 1, with Dirichlet boundary conditions on the
boundary of the square. The usual ve-point nite dierence approximation on the
uniform grid
leads to a linear system of size N N where After rescaling, the coecient
matrix of the system may be written as a sum of Kronecker products
Bm
where
. 1
mm
and I m is the identity matrix of order m. It is well known and easy to verify that the
eigenvalues of Bm are
and that the eigenvalues  j;k of AN are connected with the eigenvalues of Bm via
most of the eigenvalues have multiplicity at least 2. Also,  j;m+1 j
and the eigenvalue 4 has multiplicity m. We suspect that
which is conrmed by our numerical
experiments presented below.
To calculate the asymptotic distribution of the eigenvalues  j;k as
rst note that the eigenvalues  k of Bm are in [0; 4] and have the asymptotic density
(4
1. The asymptotic density of the  is then given by the
convolution of v with itself, i.e.,
SUPERLINEAR CONVERGENCE OF CONJUGATE GRADIENTS 23
Fig. 3. Bar chart of the eigenvalue distribution of the matrix A 1600 (without counting multi-
plicities) resulting from discretizing the 2D Poisson equation on a uniform grid with
points. The solid line corresponds to the asymptotic density function.
where the factor 1=2 is added in accordance with the multiplicities of the eigenvalues
of AN . The density w is symmetric around 4. For  2 (0; 4), we have from (5.4) and
In (5.6) we put make the change of variables
to obtain
w(4 4x) =8 2
Z 11
By the Euler integral representation for hypergeometric functions, (5.7) is
w(4 4x) =8 F (1=2; 1=2;
is a Pochhammer symbol. It is interesting to observe
that 4 2 w(4 4x) equals the complete elliptic integral of the rst kind, evaluated at
Eqn. 7.3.2.(75)]. Since w is symmetric around 4 and the
right-hand side of (5.8) is even in x, the formula (5.8) holds for 1 < x < 0 as well.
B. BECKERMANN AND A. B. J. KUIJLAARS
In the series in the right-hand side of (5.8) each term is clearly decreasing as
increases. Thus w() is increasing for  2 (0; 4), which is also clear from

Figure

3. Then also
increases for  2 (0; 4), and therefore the assumptions
of Lemma 3.1(c) are satised. From Lemma 3.1(c) we then conclude that,
for every t 2 (0; 1=2), the set S(t) associated with w()d is
with r 2 (0;
w() d: (5.9)
Putting in (5.9) we have
x
w(4 4x) dx: (5.10)
Inserting the series (5.8) for w(4 4x) and interchanging integration and summation,
we nd
x
For each k, the integral is easily transformed to a beta-integral, and it follows that
x
see also [Kuij99]. Inserting this in (5.11), we obtain
This is a known series expansion for the arccos function
see [PrBrMa90, Eqn. 7.3.2.(76)]. Inverting this we obtain the remarkably simple
and so
-5log(rel.
iterations
2D Poisson discretized, 150 2 inner gridpoints
CG error energy norm
Classical bound
Our asymptotic bound
Fig. 4. The CG error curve versus the two upper bounds for the system AN resulting
from discretizing the 2D Poisson equation on a uniform grid with points. We have
chosen a random solution x, and initial residual As predicted by (5.15), we obtain
superlinear convergence from the beginning. Notice that the classical upper bound for CG is far too
pessimistic for larger iteration indices. Similar plots are obtained for other mesh sizes.
Finally, after a small calculation using (1.10) and (3.1) we obtain the convergence
factor
log  t =t
Z tlog(tan(
Notice that, for small t, the set S(t) of (5.14) approximately equals the set obtained for
equidistant eigenvalues on [0; 8], compare with x3. This observation is in accordance
with the behavior of the eigenvalues of AN at the endpoints of
One might be curious about what CG error curve is obtained if other boundary
conditions are imposed. In this case, we need to modify O(m) rows of AN , and
such \small rank" perturbations have been covered in Remark 2.5. However, since
multiplicities are in general not preserved by such modications, we need to have a
closer look in order to obtain sharp error bounds.
In our case we can be more precise since again the eigenvalues can be computed
explicitly for a number of congurations (see, e.g., [ChEl89]). For instance, in case of
periodic boundary conditions, most of the eigenvalues are of multiplicity 8. Thus, in
accordance with [ChEl89], the convergence behavior for Dirichlet boundary conditions
with mesh size h is similar to the one obtained for periodic boundary conditions with
mesh size h=2. In case of \no-
ow" Neumann boundary conditions on the vertical
boundaries discretized by a rst order scheme, the corresponding eigen-values
are given by (3.4) plus the m eigenvalues of Bm . Here we may expect the same
convergence behavior as for Dirichlet boundary conditions.


Appendix

A. In the appendix we state and prove a lemma that is used in the
proof of Theorem 2.1.
26 B. BECKERMANN AND A. B. J. KUIJLAARS
Lemma A.1. Let  be a nite Borel measure on R with compact support. Suppose
(N ) N is a sequence of sets, all contained in a xed compact set, such that
lim
Z
for all continuous functions f on R.
let  be a Borel probability measure such that t  . Let
#N such that n=N ! t. Then there exists a sequence of sets (ZN ) N such
that
(a)
(b) ZN  N , and
(c) for all continuous functions f ,
lim
Z
Furthermore, if K is a closed set such that then the
sets ZN can be chosen such that in addition to (a), (b) and (c), we also have for N
large enough,
(d) N \ K  ZN .
Proof. We have to prove that for some sets ZN satisfying (a) and (b) the normalized
counting measures
converge in weak  sense to t. To show this, we proceed in three steps.
Step 1 Suppose we have a nite partition of R consisting of measurable sets
j. Since the normalized counting
measures of the sets N tend to , we then have for every j,
lim
then possible to choose, for every j and N , a subset
ZN;j  N \ U j such that
lim
The sets ZN;j are disjoint and for their union
Z
we have Z
Hence
lim
#(Z
Then also
lim
lim
#(Z
so that #Z
1. The set Z
N may not have exactly n elements.
By adding or deleting o(N) elements, we obtain from Z
N a set ZN with exactly n
elements. If we add elements, we choose them from N . Then ZN  N and the
limits
lim
hold.
Now assume we have a nite collection U j , of measurable
sets such that (@U j j. The sets U j are not necessarily disjoint. For each
I
I =@ \
U jA \@ \
(R
The sets V I with I ranging over all subsets of a partition of R. By
Step 1, see (A.1), there exist sets ZN such that
lim
Since every U j is a nite disjoint union of some of the V I , it also follows that
lim
be a basis for the topology of R, chosen such
that (@U j j. From Step 2 we get for each k, a sequence of sets (Z (k)
such that #Z (k)
N  N and
lim
see (A.2). Then by a diagonal argument, it is possible to nd a sequence (k N ) tending
to innity, such that the sets ZN dened by
lim
We also have
28 B. BECKERMANN AND A. B. J. KUIJLAARS
so that (a) and (b) hold.
if  N is the normalized counting measure of ZN , then by (A.3) we have for
every j,
lim
Since the U j form a basis for the open sets, it follows that the measures N tend in
sense to t. Thus (c) holds.
Next, assume that K is a closed set such that
Then
lim
#(N \ K)
and since also we have because of (c)
lim
#(ZN \ K)
we then have
Then we modify ZN by adding the elements of (N n ZN ) \ K to ZN and removing
o(N) arbitrary elements from ZN n K. This is always possible for N large enough.
Then clearly (d) is satised, while (a), (b) and (c) continue to hold.
This completes the proof of the lemma.



--R

On a conjecture of
The sensitivity of least squares polynomial approxima- tion

Families of equilibrium measures with external
Fourier analysis of iterative methods for elliptic problems
Conjugate Gradient Methods for Toeplitz systems
Constrained energy problems with applications to orthogonal polynomials of a discrete variable
From potential theory to matrix iteration in six steps
Comparisions of splittings used with the conjugate gradient algorithm



Which eigenvalues are found by the Lanczos method?
Equilibrium problems associated with fast decreasing polynomials
Zero distributions for discrete orthogonal poly- nomials
Extremal polynomials on discrete sets
A note on the superlinear convergence of GMRES
Convergence of Iterations for Linear Equations
Integrals and Series
Potential theory in the complex plane
Equilibrium measure and the distribution of zeros of the extremal polynomials of a discrete variable
Iterative methods for sparse linear systems
Logarithmic potentials with external
Further results on the convergence behavior of conjugate-gradients and Ritz values

The rate of convergence of conjugate gradients
Some superlinear convergence results for the conjugate gradient method
--TR

--CTR
A. L. Levin , D. S. Lubinsky, Green equilibrium measures and representations of an external field, Journal of Approximation Theory, v.113 n.2, p.298-323, December 2001
S. Helsen , M. Van Barel, A numerical solution of the constrained energy problem, Journal of Computational and Applied Mathematics, v.189 n.1, p.442-452, 1 May 2006
superlinear convergence;conjugate gradients;toeplitz systems;krylov subspace methods;logarithmic potential theory