. 12
( 83 .)


paper, we consider the class of methods based on using generalised inverses of vectors, and the formulation speciÿcally
includes the complex case wherever possible. c 2000 Elsevier Science B.V. All rights reserved.

Keywords: Epsilon algorithm; qd algorithm; Padà ; Vector-valued approximant; Wynn; Cross rule; Star identity; Compass
identity; Designant

1. Introduction

A sequence with a limit is as basic a topic in mathematics as it is a useful concept in science
and engineering. In the applications, it is usually the limit of a sequence, or a ÿxed point of its
generator, that is required; the existence of the limit is rarely an issue, and rapidly convergent
sequences are welcomed. However, if one has to work with a sequence that converges too slowly,
the epsilon algorithm is arguably the best all-purpose method for accelerating its convergence. The
algorithm was discovered by Wynn [54] and his review article [59] is highly recommended. The
epsilon algorithm can also be used for weakly diverging sequences, and for these the desired limit
is usually deÿned as being a ÿxed point of the operator that generates the sequence. There are
interesting exceptional cases, such as quantum well oscillators [51], where the epsilon algorithm is
not powerful enough and we refer to the companion paper by Homeier [33] in which the more

Corresponding author.
E-mail address: p.r.graves-morris@bradford.ac.uk (P.R. Graves-Morris).

0377-0427/00/$ - see front matter c 2000 Elsevier Science B.V. All rights reserved.
PII: S 0 3 7 7 - 0 4 2 7 ( 0 0 ) 0 0 3 5 5 - 1
52 P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51“80

powerful Levin-type algorithms, etc., are reviewed. The connections between the epsilon algorithm
and similar algorithms are reviewed by Weniger [50,52].
This paper is basically a review of the application of the epsilon algorithm, with an emphasis on
the case of complex-valued, vector-valued sequences. There are already many reviews and books
which include sections on the scalar epsilon algorithm, for example [1,2,9,17,53]. In the recent past,
there has been progress with the problem of numerical breakdown of the epsilon algorithm. Most
notably, Cordellier™s algorithm deals with both scalar and vector cases [13“16]. This work and its
theoretical basis has been extensively reviewed [26,27]. In this paper, we focus attention on how
the epsilon algorithm is used for sequences (si ) in which si ∈ Cd . The case d = 1 is the scalar case,
and the formulation for si ∈ C is essentially the same as that for si ∈ R. Not so for the vector case,
and we give full details of how the vector epsilon and vector qd algorithms are implemented when
si ∈ Cd , and of the connections with vector Padà approximation. Understanding these connections
is essential for specifying the range of validity of the methods. Frequently, the word “normally”
appears in this paper to indicate that the results may not apply in degenerate cases. The adaptations
for the treatment of degeneracy are almost the same for both real and complex cases, and so we
refer to [25 “27] for details.
In Section 2, we formulate the epsilon algorithm, and we explain its connection with Padà e
approximation and the continued fractions called C-fractions. We give an example of how the
epsilon algorithm works in ideal circumstances, without any signiÿcant loss of numerical precision
(which is an unusual outcome).
In Section 3, we formulate the vector epsilon algorithm, and we review its connection with
vector-valued Padà approximants and with vector-valued C-fractions. There are two major general-
isations of the scalar epsilon algorithm to the vector case. One of them is Brezinski™s topological
epsilon algorithm [5,6,35,48,49]. This algorithm has two principal forms, which might be called the
forward and backward versions; and the backward version has the orthogonality properties associated
with Lanczos methods [8]. The denominator polynomials associated with all forms of the topological
epsilon algorithm have degrees which are the same as those for the scalar case [2,5,8]. By contrast,
the other generalisation of the scalar epsilon algorithm to the vector case can be based on using
generalised inverses of vectors, and it is this generalisation which is the main topic of this paper. We
illustrate how the vector epsilon algorithm works in a two-dimensional real space, and we give a re-
alistic example of how it works in a high-dimensional complex space. The denominator polynomials
used in the scalar case are generalised both to operator polynomials of the same degree and to scalar
polynomials of double the degree in the vector case, and we explain the connections between these
twin generalisations. Most of the topics reviewed in Section 3 have a direct generalisation to the
rational interpolation problem [25]. We also note that the method of GIPAs described in Section 3
generalises directly to deal with sequences of functions in L2 (a; b) rather than vectors Cd ; in this
sense, the vectors are regarded as discretised functions [2].
In Section 4 we review the use of the vector qd algorithm for the construction of vector-valued
C-fractions, and we note the connections between vector orthogonal polynomials and the vector
epsilon algorithm. We prove the cross-rule (4.18), (4.22) using a Cli ord algebra. For real-valued
vectors, we observe that it is really an overlooked identity amongst Hankel designants. Here, the
Cross Rule is proved as an identity amongst complex-valued vectors using Moore“Penrose inverses.
The importance of studying the vector epsilon algorithm lies partly in its potential [20] for
application to the acceleration of convergence of iterative solution of discretised PDEs. For
P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51“80 53

example, Gauss“Seidel iteration generates sequences of vectors which often converge too slowly
to be useful. SOR, multigrid and Lanczos methods are alternative approaches to the problem which
are currently popular, but the success of the techniques like CGS and LTPMs (see [31] for an expla-
nation of the techniques and the acronyms) indicates the need for continuing research into numerical
methods for the acceleration of convergence of vector-valued sequences.
To conclude this introductory section, we recall that all algorithms have their domains of validity.
The epsilon algorithm fails for logarithmically convergent sequences (which converge too slowly) and
it fails to ÿnd the ÿxed point of the generator of sequences which diverge too fast. For example, if
+ O(n’2 ); C = 0;
sn ’ s =
the sequence (sn ) is logarithmically convergent to s. More precisely, a sequence is deÿned to con-
verge logarithmically to s if it converges to s at a rate governed by
sn+1 ’ s
lim = 1:
n’∞ sn ’ s

Not only does the epsilon algorithm usually fail for such sequences, but Delahaye and Germain-Bonne
[18,19] have proved that there is no universal accelerator for logarithmically convergent
Reviews of series transformations, such as those of the energy levels of the quantum-mechanical
harmonic oscillator [21,50,51], and of the Riemann zeta function [34], instructively show the in-
adequacy of the epsilon algorithm when the series coe cients diverge too fast. Information about
the asymptotic form of the coe cients and scaling properties of the solution is exploited to create
purpose-built acceleration methods. Exotic applications of the -algorithm appear in [55].

2. The epsilon algorithm

The epsilon algorithm was discovered by Wynn [54] as an e cient implementation of Shanks™
method [47]. It is an algorithm for acceleration of convergence of a sequence
S = (s0 ; s1 ; s2 ; : : : ; si ∈ C) (2.1)
and it comprises the following initialisation and iterative phases:
Initialisation: For j = 0; 1; 2; : : :
( j)
= 0 (artiÿcially); (2.2)

( j)
= sj : (2.3)

Iteration: For j; k = 0; 1; 2; : : :
( j) ( j+1) ( j+1) ( j) ’1
= +[ ’ k] : (2.4)
k+1 k’1 k

The entries k j) are displayed in the epsilon table on the left-hand side of Fig. 1, and the initialisation

has been built in.
54 P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51“80

Fig. 1. The epsilon table, and a numerical example of it.

Example 2.1. Gregory™s series for tan’1 z is
z3 z5 z7
tan’1 z = z ’ + ’ + ··· : (2.5)
3 5 7
This series can be used to determine the value of by evaluating its MacLaurin sections at z = 1:
sj := [4 tan’1 (z)]0 ; j = 0; 1; 2; : : : : (2.6)

Nuttall™s notation is used here and later on. For a function whose MacLaurin series is
(z) = + 1z + 2z + ···;

its sections are deÿned by
(z)]k i
[ = iz for 06j6k: (2.7)

In fact, sj ’ as j ’ ∞ [2] but sequence (2.6) converges slowly, as is evidenced in the column
k = 0 of entries sj = 0j) in Fig. 1. The columns of odd index have little signiÿcance, whereas the

columns of even index can be seen to converge to , which is the correct limit [2], increasingly
( j)
fast, as far as the table goes. Some values of 2k are also shown on the bar chart (Fig. 2). Notice
(2) (0)
that 2 = 3:145 and 4 = 3:142 cannot be distinguished visually on this scale.

In Example 2.1, convergence can be proved and the rate of convergence is also known [2]. From
the theoretical viewpoint, Example 2.1 is ideal for showing the epsilon algorithm at its best. It
is noticeable that the entries in the columns of odd index are large, and this e ect warns us to
beware of possible loss of numerical accuracy. Like all algorithms of its kind (which use reciprocal
di erences of convergent sequences) the epsilon algorithm uses (and usually uses up) numerical
precision of the data to do its extrapolation. In this case, there is little loss of numerical precision
using 16 decimal place (MATLAB) arithmetic, and 22 = almost to machine precision. In this
case, the epsilon algorithm converges with great numerical accuracy because series (2.5) is a totally
oscillating series [4,7,17,59].
To understand in general how and why the epsilon algorithm converges, whether we are referring
( j) ( j)
to its even columns ( 2k ; j = 0; 1; 2; : : : ; k ÿxed) or its diagonals ( 2k ; k = 0; 1; 2; : : : ; j ÿxed) or any
P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51“80 55

( j)
Fig. 2. Values of 2k for Example 2.1, showing the convergence rate of the epsilon algorithm using n + 1 = 1; 2; 3; 4; 5
terms of the given sequence.

other sequence, the connection with Padà approximation is essential [1,2,56]. Given a (possibly
formal) power series
f(z) = c0 + c1 z + c2 z 2 + · · · ; (2.8)
the rational function
A(z)B(z)’1 ≡ [˜=m](z) (2.9)
is deÿned as a Padà approximant for f(z) of type [˜=m] if
(i) deg{A(z)}6˜; deg{B(z)}6m; (2.10)

(ii) f(z)B(z) ’ A(z) = O(z ˜+m+1 ); (2.11)

(iii) B(0) = 0: (2.12)
The Baker condition
B(0) = 1 (2.13)
is often imposed for reliability in the sense of (2.14) below and for a deÿnite speciÿcation of A(z)
and B(z). The deÿnition above contrasts with the classical (Frobenius) deÿnition in which axiom
(iii) is waived, and in this case the existence of A(z) and B(z) is guaranteed, even though (2.14)
below is not. Using speciÿcation (2.10) “ (2.13), we ÿnd that
f(z) ’ A(z)B(z)’1 = O(z ˜+m+1 ); (2.14)
56 P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51“80

Fig. 3. Relative location of Padà approximants.

provided that a solution of (2.15) below can be found. To ÿnd B(z), the linear equations corre-
sponding to accuracy-through-orders z ˜+1 ; z ˜+2 ; : : : ; z ˜+m in (2.11) must be solved. They are
® ®  ® 
c˜’m+1 ::: c˜ bm c˜+1
. .  .  .
 . . »° . » = ’° . »: (2.15)
° . . . .
c˜ ::: c˜+m’1 b1 c˜+m
The coe cients of B(z) = m bi z i are found using an accurate numerical solver of (2.15). By
contrast, for purely theoretical purposes, Cramer™s rule is applied to (2.15). We are led to deÿne
c˜’m+1 c˜’m+2 ::: c˜+1
c˜’m+2 c˜’m+3 ::: c˜+2
. . .
. . .
q[˜=m] (z) = (2.16)
. . .
c˜ c˜+1 ::: c˜+m
zm z m’1 ::: 1
and then we ÿnd that
B[˜=m] (z) = q[˜=m] (z)=q[˜=m] (0) (2.17)
is the denominator polynomial for the Padà approximation problem (2.9) “ (2.15) provided that
q (0) = 0.
The collection of Padà approximants is called the Padà table, and in Fig. 3 we show ÿve neigh-
e e
bouring approximants in the table.
These approximants satisfy a ÿve-point star identity,
[N (z) ’ C(z)]’1 + [S(z) ’ C(z)]’1 = [E(z) ’ C(z)]’1 + [W (z) ’ C(z)]’1 ; (2.18)
called Wynn™s identity or the compass identity. The proof of (2.18) is given in [1,2], and it is
also a corollary (in the case d = 1) of the more general result (3.59) that we prove in the next
section. Assuming (2.18) for the moment, the connection between Padà approximation and the epsilon
P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51“80 57

Fig. 4. Some artiÿcial entries in the Padà table are shown circled.

algorithm is given by connecting the coe cients of f(z) with those of S with
c0 = s0 ; ci = si ’ si’1 ; i = 1; 2; 3; : : : ;
and by


. 12
( 83 .)