A New Algorithm for the Higher-Order $G$-Transformation

TL;DR

This paper introduces the FS/qd-algorithm, a new efficient method combining existing algorithms for higher-order $G$-transformation, improving computational efficiency and applicability to Shanks' transformation.

Contribution

The paper develops the FS/qd-algorithm, merging Ford and Sidi's FS-algorithm with Rutishauser's qd-algorithm, offering a more efficient approach for higher-order $G$-transformation.

Findings

FS/qd-algorithm has fewer operations than rs-algorithm.

FS/qd can implement Shanks' transformation effectively.

FS/qd compares favorably with Wynn's $oldsymbol{ extit{ extbf{ extepsilon}}}$-algorithm.

Abstract

Let the scalars $A^{(j)}_n$ be defined via the linear equations $$A_l=A^{(j)}_n+\sum^n_{k=1}\bar{\alpha}_ku_{k+l-1},\ \ l=j,j+1,\ldots,j+n\ .$$ Here the $A_i$ and $u_i$ are known and the $\bar{\alpha}_k$ are additional unknowns, and the quantities of interest are the $A^{(j)}_n$. This problem arises, for example, when one computes infinite-range integrals by the higher-order $G$-transformation of Gray, Atchison, and McWilliams. One efficient procedure for computing the $A^{(j)}_n$ is the rs-algorithm of Pye and Atchison. In the present work, we develop yet another procedure that combines the FS-algorithm of Ford and Sidi and the qd-algorithm of Rutishauser, and we denote it the FS/qd-algorithm. We show that the FS/qd-algorithm has a smaller operation count than the rs-algorithm. We also show that the FS/qd algorithm can also be used to implement the transformation of Shanks, and compares very favorably with the $\varepsilon$-algorithm of Wynn that is normally used for this purpose.

Equations

$$A_l=A_n^{(j)}+\sum _{k=1}^n{\bar{\alpha }}_k{u}_{k+l-1},l=j,j+1,\dots ,j+n.$$

$$G_{n}(x;h)=\frac{\left|\begin{array}[]{cccc}F(x)&F(x+h)&\cdots&F(x+nh)\\ f(x)&f(x+h)&\cdots&f(x+nh)\\ \vdots&\vdots&&\vdots\\ f(x+(n-1)h)&f(x+nh)&\cdots&f(x+(2n-1)h)\end{array}\right|}{\left|\begin{array}[]{cccc}1&1&\cdots&1\\ f(x)&f(x+h)&\cdots&f(x+nh)\\ \vdots&\vdots&&\vdots\\ f(x+(n-1)h)&f(x+nh)&\cdots&f(x+(2n-1)h)\end{array}\right|},$$

$$F(x)={\int }_a^{x}f(t)dt.$$

$$F(x+ih)=G_n(x;h)+\sum _{k=1}^n{\bar{\alpha }}_{k}f(x+(i+k-1)h),i=0,1,\dots ,n,$$

$$A_l=A_n^{(j)}+\sum _{k=1}^n{\bar{\alpha }}_k{u}_{k+l-1},l=j,j+1,\dots ,j+n,$$

$$r_n^{(j)}=\frac{H_n^{(j)}}{K_n^{(j)}},s_n^{(j)}=\frac{K_{n+1}^{(j)}}{H_n^{(j)}},$$

$$H^{(j)}_{n}=\left|\begin{array}[]{cccc}u_{j}&u_{j+1}&\cdots&u_{j+n-1}\\ u_{j+1}&u_{j+2}&\cdots&u_{j+n}\\ \vdots&\vdots&&\vdots\\ u_{j+n-1}&u_{j+n}&\cdots&u_{j+2n-2}\end{array}\right|,\ \ H^{(j)}_{0}=1,$$

$$K^{(j)}_{n}=\left|\begin{array}[]{cccc}1&1&\cdots&1\\ u_{j}&u_{j+1}&\cdots&u_{j+n-1}\\ u_{j+1}&u_{j+2}&\cdots&u_{j+n}\\ \vdots&\vdots&&\vdots\\ u_{j+n-2}&u_{j+n-1}&\cdots&u_{j+2n-3}\end{array}\right|,\ \ K^{(j)}_{0}=1.$$

$$s_0^{(j)}=1,r_1^{(j)}=u_j,A_0^{(j)}=A_j,j=0,1,\dots .$$

$$s_n^{(j)}=s_{n-1}^{(j+1)}\left(\frac{r_n^{(j+1)}}{r_n^{(j)}}-1\right),r_{n+1}^{(j)}=r_n^{(j+1)}\left(\frac{s_n^{(j+1)}}{s_n^{(j)}}-1\right).$$

$$A_n^{(j)}=\frac{r_n^{(j)}{A}_{n-1}^{(j+1)}-r_n^{(j+1)}{A}_{n-1}^{(j)}}{r_n^{(j)}-r_n^{(j+1)}}.$$

$$A_l=A_n^{(j)}+\sum _{k=1}^n{\bar{\alpha }}_k{g}_k(l),l=j,j+1,\dots ,j+n,$$

$$A_n^{(j)}=\frac{R_n^{(j)}{A}_{n-1}^{(j+1)}-R_n^{(j+1)}{A}_{n-1}^{(j)}}{R_n^{(j)}-R_n^{(j+1)}},$$

$$|u_{1}(j)\,u_{2}(j)\,\cdots\,u_{n}(j)|=\left|\begin{array}[]{cccc}u_{1}(j)&u_{2}(j)&\cdots&u_{n}(j)\\ u_{1}(j+1)&u_{2}(j+1)&\cdots&u_{n}(j+1)\\ \vdots&\vdots&&\vdots\\ u_{1}(j+n-1)&u_{2}(j+n-1)&\cdots&u_{n}(j+n-1)\end{array}\right|,$$

$$G_n^{(j)}=|g_1(j)g_2(j)\cdots g_n(j)|,G_0^{(j)}=1.$$

$$f_n^{(j)}(b)=|g_1(j)g_2(j)\cdots g_n(j)b(j)|,f_0^{(j)}(b)=b(j).$$

$${\psi }_n^{(j)}(b)=\frac{f_n^{(j)}(b)}{G_{n+1}^{(j)}}.$$

$$A_n^{(j)}=\frac{{\psi }_n^{(j)}(a)}{{\psi }_n^{(j)}(I)},$$

$${\psi }_n^{(j)}(b)=\frac{{\psi }_{n-1}^{(j+1)}(b)-{\psi }_{n-1}^{(j)}(b)}{D_n^{(j)}},D_n^{(j)}=\frac{G_{n+1}^{(j)}{G}_{n-1}^{(j+1)}}{G_n^{(j)}{G}_n^{(j+1)}}.$$

$$D_n^{(j)}=\frac{H_{n+1}^{(j)}{H}_{n-1}^{(j+1)}}{H_n^{(j)}{H}_n^{(j+1)}}=e_n^{(j)}.$$

$$q_n^{(j)}=\frac{H_{n-1}^{(j)}{H}_n^{(j+1)}}{H_n^{(j)}{H}_{n-1}^{(j+1)}}.$$

$$e_n^{(j)}=q_n^{(j+1)}-q_n^{(j)}+e_{n-1}^{(j+1)},q_{n+1}^{(j)}=\frac{e_n^{(j+1)}}{e_n^{(j)}}{q}_n^{(j+1)},j\ge 0,n\ge 1,$$

$$\begin{array}[]{cccccc}&q^{(0)}_{1}&&&&\\ e^{(1)}_{0}&&e^{(0)}_{1}\\ &q^{(1)}_{1}&&q^{(0)}_{2}\\ e^{(2)}_{0}&&e^{(1)}_{1}&&e^{(0)}_{2}\\ &q^{(2)}_{1}&&q^{(1)}_{2}&&\ddots\\ e^{(3)}_{0}&&e^{(2)}_{1}&&e^{(1)}_{2}\\ &q^{(3)}_{1}&&q^{(2)}_{2}&&\ddots\\ e^{(4)}_{0}&\vdots&e^{(3)}_{1}&&e^{(2)}_{2}\\ \vdots&&\vdots&q^{(3)}_{2}&\vdots&\ddots\\ &&&\vdots&\end{array}$$

$$e_0^{(j)}=0,q_1^{(j)}=\frac{u_{j+1}}{u_j},M_0^{(j)}=\frac{A_j}{u_j},N_0^{(j)}=\frac{1}{u_j}.$$

$$e_n^{(j)}=q_n^{(j+1)}-q_n^{(j)}+e_{n-1}^{(j+1)},q_{n+1}^{(j)}=\frac{e_n^{(j+1)}}{e_n^{(j)}}{q}_n^{(j+1)},$$

$$M_n^{(j)}=\frac{M_{n-1}^{(j+1)}-M_{n-1}^{(j)}}{e_n^{(j)}},N_n^{(j)}=\frac{N_{n-1}^{(j+1)}-N_{n-1}^{(j)}}{e_n^{(j)}}.$$

$$A_n^{(j)}=\frac{M_n^{(j)}}{N_n^{(j)}}.$$

$${\epsilon }_{k+1}^{(j)}={\epsilon }_{k-1}^{(j+1)}+\frac{1}{{\epsilon }_k^{(j+1)}-{\epsilon }_k^{(j)}},j,k=0,1,\dots ,$$

Full text

[TABLE]

Here the $A_{i}$ and $u_{i}$ are known and the $\bar{\alpha}_{k}$ are additional unknowns, and the quantities of interest are the $A^{(j)}_{n}$. This problem arises, for example, when one computes infinite-range integrals by the higher-order $G$-transformation of Gray, Atchison, and McWilliams. One efficient procedure for computing the $A^{(j)}_{n}$ is the rs-algorithm of Pye and Atchison. In the present work, we develop yet another procedure that combines the FS-algorithm of Ford and Sidi and the qd-algorithm of Rutishauser, and we denote it the FS/qd-algorithm. We show that the FS/qd-algorithm has a smaller operation count than the rs-algorithm. We also show that the FS/qd algorithm can also be used to implement the transformation of Shanks, and compares very favorably with the $\varepsilon$-algorithm of Wynn that is normally used for this purpose.

Mathematics Subject Classification 2000 : 65B05, 65B10, 65D30, 41A21.

1 The Higher Order $G$-Transformation

The $G$-transformation was designed by Gray and Atchison [5] as an extrapolation method for evaluating infinite integrals of the form $\int^{\infty}_{a}f(t)\,dt\equiv I[f]$. It was later generalized in different ways in Atchison and Gray [1] and Gray and Atchison [6], the ultimate generalization being given in Gray, Atchison, and McWilliams [7]. This generalization was denoted the higher-order $G$-transformation. The way it is defined in [7], this transformation produces approximations to $I[f]$ that are of the form

[TABLE]

where

[TABLE]

The approximations produced by the $G$-transformation of [5] are simply the $G_{1}(x;h)$.

It follows from (1) that $G_{n}(x;h)$ is also the solution of the linear system

[TABLE]

where $\bar{\alpha}_{k}$ are additional unknowns.

It has been shown in [7] that the kernel of the higher-order $G$-transformation is the set of functions $f(x)$ that are integrable at infinity in the sense of Abel and that satisfy linear homogeneous ordinary differential equations of order $n$ with constant coefficients. Thus $f(x)$ is in this kernel if it is of the form $f(x)=\sum^{r}_{k=1}p_{k}(x)e^{c_{k}x}$, where the $c_{k}\neq 0$ are distinct and $\Re c_{k}\leq 0$, $p_{k}(x)$ are polynomials. If $\mu_{k}$ is the degree of $p_{k}(x)$ for each $k$, and if $\sum^{r}_{k=1}(\mu_{k}+1)=n$, then $G_{n}(x;h)=I[f]$ for all $x$ and $h$. On the basis of this result it was concluded in Levin and Sidi [11] that the higher-order $G$-transformation will be effective on functions of the form $f(x)=\sum^{s}_{k=1}e^{c_{k}x}v_{k}(x)$, where $v_{k}(x)\sim\sum^{\infty}_{i=0}\alpha_{ki}x^{\gamma_{k}-i}$ as $x\to\infty$, with arbitrary $\gamma_{k}$, provided $h$ is of a suitable size.

In the present work, we are concerned with the actual computation of the $G_{n}(x;h)$. Of course, it is not desirable to compute $G_{n}(x;h)$ via the determinantal representation in (1). Direct solution of the linear systems in (3) is expensive too. A very efficient and elegant procedure for computing the $G_{n}(x;h)$ was given by Pye and Atchison in [12], and it has been denoted the rs-algorithm in Brezinski and Redivo Zaglia [3]. The derivation of the rs-algorithm makes use of the representation in (1). In the present work, we develop yet another procedure for computing the $G_{n}(x;h)$, and we call this new procedure the FS/qd-algorithm. We show that the FS/qd-algorithm is more efficient than the rs-algorithm. With proper substitutions to be discussed in Section 3, the FS/qd-algorithm can also be used to implement the transformation of Shanks [15], which is normally implemented via the well-known $\varepsilon$-algorithm of Wynn [17]. We show that, when used for implementing the transformation of Shanks, the FS/qd-algorithm compares very favorably with the $\varepsilon$-algorithm.

2 Algorithms for the Higher Order $G$-Transformation

We start with a review of the algorithm of Pye and Atchison [12]. Actually these authors consider the more general problem in which one would like to compute the quantities $A^{(j)}_{n}$ defined via the linear equations

[TABLE]

where the $A_{i}$ and $u_{i}$ are known scalars, while the $\bar{\alpha}_{k}$ are not necessarily known. Comparing (3) with (4) we can draw the analogy $A_{i}\leftrightarrow F(x+ih)$, $u_{i}\leftrightarrow f(x+ih)$, and $A^{(j)}_{n}\leftrightarrow G_{n}(x+jh;h)$.

The algorithm of [12] computes the $A^{(j)}_{n}$ with the help of two sets of auxiliary quantities, $r^{(j)}_{n}$ and $s^{(j)}_{n}$. These quantities are defined by

[TABLE]

where $H^{(j)}_{n}$, the Hankel determinant associated with $\{u_{s}\}$, and $K^{(j)}_{n}$ are given as in

[TABLE]

and

[TABLE]

The rs-algorithm computes the $r^{(j)}_{n}$, and $s^{(j)}_{n}$ simultaneously by efficient recursions. Once these have been computed the $A^{(j)}_{n}$ can be computed via a separate recursion.

The rs-Algorithm

Set

[TABLE] 2. 2.

For $j=0,1,\ldots,$ and $n=1,2,\ldots,$ compute recursively

[TABLE] 3. 3.

For $j=0,1,\ldots,$ and $n=1,2,\ldots,$ set

[TABLE]

Before proceeding further, we note that the equations in (4) are the same as

[TABLE]

with $g_{k}(l)=u_{k+l-1}$, $l=0,1,\ldots\ .$ Linear systems as in (8) arise also in the definition of a generalized Richardson extrapolation method. This suggests that the E-algorithm of Schneider [14] and the FS-algorithm of Ford and Sidi [4] can be used for computing the $A^{(j)}_{n}$. Different derivations of the E-algorithm were given by Håvie [8] and Brezinski [2]. Of course, direct application of these algorithms without taking into account the special nature of the $g_{k}(l)$ is very uneconomical. By taking the nature of the $g_{k}(l)$ into consideration, it becomes possible to derive fast algorithms for the $A^{(j)}_{n}$.

Now the E-algorithm produces the $A^{(j)}_{n}$ by a recursion relation of the form

[TABLE]

the most expensive part of the algorithm being the determination of the $R^{(j)}_{n}$. Comparing the known expression for $R^{(j)}_{n}$ when $g_{k}(l)=u_{k+l-1}$ with that of $r^{(j)}_{n}$, we realize that $R^{(j)}_{n}=(-1)^{n-1}r^{(j)}_{n}$. We thus conclude that the rs-algorithm is simply the E-algorithm in which the $R^{(j)}_{n}$, whose determination forms the most expensive part of the E-algorithm, are computed by a fast recursion. For this point and others, see [3, Section 2.4].

In view of the close connection between the rs- and E-algorithms, it is natural to investigate the possibility of designing another algorithm that is related to the FS-algorithm. This is worth the effort as the FS-algorithm is more economical than the E-algorithm to begin with. To this end we start with a brief description of the FS-algorithm and refer the reader to [4] for details. For a comprehensive summary, see also Sidi [16].

Let us first define the short-hand notation

[TABLE]

and set

[TABLE]

Next, let us agree to denote the sequence $\{b(l)\}^{\infty}_{l=0}$ by $b$ for short, and define

[TABLE]

Finally, let us define

[TABLE]

Then we have

[TABLE]

where $a$ and $I$ denote the sequences $\{a(l)=A_{l}\}^{\infty}_{l=0}$ and $\{I(l)=1\}^{\infty}_{l=0}$ respectively. The FS-algorithm computes the quantities $\psi^{(j)}_{n}(b)$ by a recursion of the form

[TABLE]

By this recursion, first the $\psi^{(j)}_{n}(a)$ and $\psi^{(j)}_{n}(I)$ are computed and then $A^{(j)}_{n}$ is determined via (13).

We recall that the most expensive part of the FS-algorithm is the (recursive) determination of the quantities $D^{(j)}_{n}$, and we would like to reduce the cost of this part. Fortunately, this can be achieved once we realize that, with $G^{(j)}_{n}$ and $H^{(j)}_{n}$ defined as in (10) and (6) respectively, and with $g_{k}(l)=u_{k+l-1}$, we have $G^{(j)}_{n}=H^{(j)}_{n}$ in the present case. From this and (14), we obtain the surprising result that

[TABLE]

Here $e^{(j)}_{n}$ is a quantity computed by the famous qd-algorithm of Rutishauser [13], a clear exposition of which can also be found in Henrici [9]. Actually, the qd-algorithm computes along with the $e^{(j)}_{n}$ also the quantities $q^{(j)}_{n}$ that are given as in

[TABLE]

The $q^{(j)}_{n}$ and $e^{(j)}_{n}$ serve the construction of the regular $C$-fractions, hence the Padé approximants, associated with the formal power series $\sum^{\infty}_{k=0}u_{k}z^{k}$. Regular $C$-fractions are continued fractions of a special type. (See Jones and Thron [10].) The qd-algorithm computes the $q^{(j)}_{n}$ and $e^{(j)}_{n}$ via the recursions

[TABLE]

with the initial conditions $e^{(j)}_{0}=0$ and $q^{(j)}_{1}=u_{j+1}/u_{j}$ for all $j\geq 0$. The quantities $q^{(k)}_{n}$ and $e^{(k)}_{n}$ can be arranged in a two-dimensional array as in Figure 1.

This observation enables us to combine the FS- and qd-algorithms to obtain the following economical implementation, the FS/qd-algorithm, for the higher-order $G$-transformation. For simplicity of notation, we will let $\psi^{(j)}_{n}(a)=M^{(j)}_{n}$ and $\psi^{(j)}_{n}(I)=N^{(j)}_{n}$ in the FS-algorithm.

The FS/qd-Algorithm

Set

[TABLE] 2. 2.

For $j=0,1,\ldots,$ and $n=1,2,\ldots,$ compute recursively

[TABLE]

[TABLE] 3. 3.

For $j=0,1,\ldots,$ and $n=1,2,\ldots,$ set

[TABLE]

It seems that, when a certain number of the $u_{i}$ are given, it is best to compute the associated qd-table columnwise. (Note that the quantities in each of the recursions for $e^{(j)}_{n}$ and $q^{(j)}_{n+1}$ in step 2 of the FS/qd-algorithm above form the four corners of a lozenge in Figure 1.) Following that we can compute the $M^{(j)}_{n}$, $N^{(j)}_{n}$, and $A^{(j)}_{n}$ columnwise as well.

3 Comparison of the rs- and FS/qd-Algorithms

Let us now compare the operation counts of the two algorithms. First, we note that the $r^{(j)}_{n}$ and $s^{(j)}_{n}$ in the rs-algorithm can be arranged in a table similar to the qd-table of the $e^{(j)}_{n}$ and $q^{(j)}_{n}$. Thus, given $A_{0},A_{1},\ldots,A_{L}$, and $u_{0},u_{1},\ldots,u_{2L}$, we can compute $A^{(j)}_{n}$ for $0\leq j+n\leq L$. Now the number of the $e^{(j)}_{n}$ in the relevant qd-table is $L^{2}+O(L)$ and so is that of the $q^{(j)}_{n}$. A similar statement can be made concerning the $r^{(j)}_{n}$ and $s^{(j)}_{n}$. The number of the $A^{(j)}_{n}$ is $L^{2}/2+O(L)$, and so are the numbers of the $M^{(j)}_{n}$ and the $N^{(j)}_{n}$. Consequently, we have the following operation counts.

[TABLE]

In case only the $A^{(0)}_{n}$ are needed (as they have the best convergence properties), the number of divisions in the FS/qd-algorithm can be reduced from $5L^{2}/2+O(L)$ to $2L^{2}+O(L)$. In any case, we see that the operation count of the rs-algorithm is about 30% more than that of the FS/qd-algorithm.

Finally, we observe that when $u_{k}=\Delta A_{k}=A_{k+1}-A_{k}$ the higher-order $G$-transformation reduces to the Shanks [15] transformation, and therefore the rs- and FS/qd-algorithms can be used for computing the $A^{(j)}_{n}$ in this case too. Of course, the most famous and efficient implementation of the Shanks transformation is via the $\varepsilon$-algorithm of Wynn [17], which reads

[TABLE]

with the initial values $\varepsilon^{(j)}_{-1}=0$ and $\varepsilon^{(j)}_{0}=A_{j}$ for all $j\geq 0$. Then $\varepsilon^{(j)}_{2n}=A^{(j)}_{n}$ for all $j$ and $n$. Thus, given $A_{0},A_{1},\ldots,A_{2L}$, we can compute $\varepsilon^{(j)}_{2n}$ for $0\leq j+2n\leq 2L$. (In particular, we can compute the diagonal approximations $\varepsilon^{(j)}_{2n}$, $n=0,1,\ldots,L,$ that have the best convergence properties.) Since there are $2L^{2}+O(L)$ of the $\varepsilon^{(j)}_{k}$ to compute, the operation count of this computation is $4L^{2}+O(L)$ additions and $2L^{2}+O(L)$ divisions and no multiplications. It is seen from the table above that the FS/qd-algorithm compares very favorably with the $\varepsilon$-algorithm as well.