on epsilon extrapolation

Overview

Wynn’s $ϵ$-algorithm turns a scalar sequence into a triangular table of nonlinear extrapolants (Wynn, 1956). This note builds the method from Aitken’s ${\Delta }^2$ process, explains the Shanks transform as exact elimination of finitely many geometric error modes, and then gives the main structural fact: when the input sequence consists of Taylor partial sums, the even $ϵ$-entries are Padé approximants (Shanks, 1955; Wynn, 1966; Gragg, 1972). In this sense epsilon extrapolation is a sequence-side algorithm for rational approximation.

🏷️ The Extrapolation Problem

Suppose we know a sequence

$$s_0,s_1,s_2,\dots$$

which is expected to converge to an unknown limit $s$. A convergence acceleration method constructs a new sequence

$$t_n=T(s_n,s_{n+1},\dots ,s_{n+m})$$

which should have the same limit but smaller error. In numerical analysis this is useful because the original sequence often comes from partial sums, mesh refinements, fixed-point iterations, or truncations of an asymptotic expansion.

The guiding principle is model-based. If the error has a simple form, then finitely many sequence values may determine the missing limit exactly. For example, Richardson extrapolation assumes an error expansion in powers of a small step size $h$. The $ϵ$-algorithm assumes, locally, that the tail is well approximated by a short sum of geometric modes:

$$s_n=s+c_1{\lambda }_1^n+\cdots +c_k{\lambda }_k^n.$$

This model is natural for linearly convergent fixed-point iterations and for partial sums whose remainders are controlled by nearby singularities. It is also the bridge to rational approximation: a rational function is exactly what one obtains after summing finitely many geometric tails.

We use the forward difference notation

$$\Delta s_n=s_{n+1}-s_n,{\Delta }^2{s}_n=s_{n+2}-2s_{n+1}+s_n.$$

🏷️ Why Rational Functions Appear

The geometric model is already a rational model in disguise. If the error is

$$r_n=s_n-s=\sum _{j=1}^k{c}_j{\lambda }_j^n,$$

then its generating function is

$$\sum _{n\ge 0}{r}_n{w}^n=\sum _{j=1}^k\frac{c_j}{1-{\lambda }_{j}w},$$

which is rational with denominator ${\prod }_j(1-{\lambda }_{j}w)$. Thus removing a few geometric modes from a sequence is the same algebraic operation as fitting a denominator to a formal series. This is the conceptual reason the $ϵ$-algorithm meets Padé approximation rather than ordinary polynomial approximation.

🏷️ Aitken’s Process as the First Case

Assume first that the error has one geometric mode:

$$s_n=s+c{\lambda }^n,|\lambda |<1.$$

Then

$$\Delta s_n=c(\lambda -1){\lambda }^n,{\Delta }^2{s}_n=c(\lambda -1{)}^2{\lambda }^n.$$

Therefore

$$\frac{(\Delta s_n{)}^2}{{\Delta }^2{s}_n}=c{\lambda }^n=s_n-s.$$

Solving for $s$ gives Aitken’s ${\Delta }^2$ transform:

$$\mathcal{A}(s_n)=s_n-\frac{(\Delta s_n{)}^2}{{\Delta }^2{s}_n}=\frac{s_n{s}_{n+2}-s_{n+1}^2}{s_{n+2}-2s_{n+1}+s_n}.$$

Interpretation

Aitken’s transform is exact on sequences of the form $s+c{\lambda }^n$. It estimates the dominant geometric error and subtracts it. This is why it can be extremely effective for linearly convergent iterations, but also why it may fail when ${\Delta }^2{s}_n$ is small or when the leading error is not geometric.

The rightmost formula is the first Shanks transform. The higher Shanks transforms remove several geometric modes at once.

🏷️ The Shanks Transform

For a fixed integer $k$, suppose that near index $n$ the sequence is modeled by

$$s_m=s+\sum _{r=1}^k{c}_r{\lambda }_r^m,m=n,n+1,\dots ,n+2k,$$

where the ${\lambda }_r$ are distinct and different from $1$. The idea is to find coefficients $q_0,\dots ,q_k$ such that the forward differences are annihilated:

$$\sum _{j=0}^k{q}_j\Delta s_{n+i+j}=0,i=0,\dots ,k-1.$$

After normalizing by

$$\sum _{j=0}^k{q}_j=1,$$

define

$$e_k(s_n)=\sum _{j=0}^k{q}_j{s}_{n+j}.$$

If the model is exact, then the same coefficients annihilate the error terms and hence $e_k(s_n)=s$. For $k=1$ this reduces to Aitken’s formula.

Equivalently, $e_k(s_n)$ can be written as a ratio of determinants:

$$e_k(s_n)=\frac{\left|\begin{array}{cccc}{s}_n& s_{n+1}& \cdots & s_{n+k}\\ \Delta s_n& \Delta s_{n+1}& \cdots & \Delta s_{n+k}\\ \Delta s_{n+1}& \Delta s_{n+2}& \cdots & \Delta s_{n+k+1}\\ ⋮& ⋮& & ⋮\\ \Delta s_{n+k-1}& \Delta s_{n+k}& \cdots & \Delta s_{n+2k-1}\end{array}\right|}{\left|\begin{array}{cccc}1& 1& \cdots & 1\\ \Delta s_n& \Delta s_{n+1}& \cdots & \Delta s_{n+k}\\ \Delta s_{n+1}& \Delta s_{n+2}& \cdots & \Delta s_{n+k+1}\\ ⋮& ⋮& & ⋮\\ \Delta s_{n+k-1}& \Delta s_{n+k}& \cdots & \Delta s_{n+2k-1}\end{array}\right|}.$$

Exactness of the Shanks Transform

If $s_m=s+{\sum }_{r=1}^k{c}_r{\lambda }_r^m$ for $m=n,\dots ,n+2k$, with distinct ${\lambda }_r\ne 1$, and if the denominator determinant is nonzero, then $e_k(s_n)=s$.

Proof

Let $p(x)={\prod }_{r=1}^k(x-{\lambda }_r)=p_0+p_{1}x+\cdots +p_k{x}^k$. The error sequence $r_m=s_m-s$ satisfies the finite recurrence

$$\sum _{j=0}^k{p}_j{r}_{m+j}=0.$$

Taking a forward difference shows that the differences $\Delta s_m=\Delta r_m$ satisfy the same recurrence. Thus the coefficient vector solving the annihilation equations is proportional to $(p_0,\dots ,p_k)$. Since $p(1)\ne 0$, it can be normalized so that ${\sum }_j{q}_j=1$. Then

$$\sum _{j=0}^k{q}_j{s}_{n+j}=s\sum _{j=0}^k{q}_j+\sum _{j=0}^k{q}_j{r}_{n+j}=s.$$

The determinant expression is Cramer’s rule applied to the annihilation equations and the normalization. ▮

This theorem explains both the power and the limitation of the method. It is exact for a finite geometric error model. When the real error is only approximately of this form, the transform may still be a good local extrapolation; when the model is poor or the determinant is nearly singular, it can amplify noise.

🏷️ Wynn’s Epsilon Algorithm

The determinant formula is conceptually useful but computationally expensive. Wynn’s observation was that all Shanks transforms can be generated by a short nonlinear recurrence. Define the $ϵ$-table by

$${ϵ}_{-1}^{(n)}=0,{ϵ}_0^{(n)}=s_n,$$

and

$${ϵ}_{k+1}^{(n)}={ϵ}_{k-1}^{(n+1)}+\frac{1}{{ϵ}_k^{(n+1)}-{ϵ}_k^{(n)}}.$$

The even entries are the approximants:

$${ϵ}_{2k}^{(n)}=e_k(s_n).$$

The odd entries are auxiliary quantities needed by the recurrence. For example,

$${ϵ}_2^{(n)}=s_{n+1}+\frac{1}{\frac{1}{s_{n+2}-s_{n+1}}-\frac{1}{s_{n+1}-s_n}}=s_n-\frac{(\Delta s_n{)}^2}{{\Delta }^2{s}_n},$$

which recovers Aitken’s process.

Singular and Nearly Singular Updates

The recurrence requires division by ${ϵ}_k^{(n+1)}-{ϵ}_k^{(n)}$. If this denominator vanishes, the table entry is undefined. If it is merely small, the algorithm may produce a large spurious value. In the Padé interpretation below, this is the same numerical phenomenon as encountering a nearby pole or a nearly canceled pole-zero pair.

🏷️ Padé Approximation Preliminaries

Let

$$f(z)=\sum _{m=0}^{\infty }{a}_m{z}^m$$

be a formal power series. The Taylor polynomial $S_N(z)={\sum }_{m=0}^N{a}_m{z}^m$ is a polynomial approximation. A Padé approximant replaces this polynomial by a rational function

$$[L/M{]}_f(z)=\frac{P_L(z)}{Q_M(z)},$$

where

$$\mathrm{deg}{P}_L\le L,\mathrm{deg}{Q}_M\le M,Q_M(0)=1,$$

and the Taylor expansion matches as far as possible:

$$Q_M(z)f(z)-P_L(z)=O(z^{L+M+1}).$$

Thus $[L/M{]}_f$ uses the same local coefficient data as a Taylor polynomial of degree $L+M$, but stores it in a rational form. This matters because rational functions can represent poles, geometric tails, and continued fractions much more efficiently than polynomials (Baker & Graves-Morris, 1996).

The denominator is determined by the homogeneous linear equations obtained from the coefficients of $z^{L+1},\dots ,z^{L+M}$:

$$\sum _{r=0}^M{q}_r{a}_{\ell -r}=0,\ell =L+1,\dots ,L+M,$$

where $Q_M(z)={\sum }_{r=0}^M{q}_r{z}^r$ and $q_0$ is then used for normalization. Once $Q_M$ is known, $P_L$ is the truncation of $Q_{M}f$ through degree $L$.

🏷️ The Padé Meaning of the Epsilon Table

Now take the input sequence to be the Taylor partial sums

$$s_n(z)=S_n(z)=\sum _{m=0}^n{a}_m{z}^m.$$

Then

$$\Delta s_n(z)=a_{n+1}{z}^{n+1}.$$

The Shanks annihilation equations become

$$\sum _{j=0}^k{q}_j\Delta S_{n+i+j}(z)=0,i=0,\dots ,k-1.$$

After substituting $\Delta S_m(z)=a_{m+1}{z}^{m+1}$, the equations take the form

$$\sum _{j=0}^k{q}_j(z)a_{n+i+j+1}{z}^j=0,i=0,\dots ,k-1.$$

Reversing the coefficients converts these into the usual Hankel equations for the denominator of the Padé approximant with numerator degree at most $n+k$ and denominator degree at most $k$. When the relevant denominator is nonzero, the resulting identity is

$${ϵ}_{2k}^{(n)}(z)=e_k(S_n(z))=[n+k/k{]}_f(z).$$

Main Principle

The $ϵ$-algorithm is a sequence transformation for arbitrary scalar data, but on Taylor partial sums it becomes a Padé-table algorithm. The row index $n$ controls how many numerator degrees are allowed, and the transform order $k$ controls the denominator degree. Degenerate cases correspond to the same defects one sees in the Padé table: zero denominator determinants, common numerator-denominator factors, or nearby poles.

This explains why the method can sum beyond what a Taylor polynomial sees. A Taylor polynomial stores finitely many local coefficients as a polynomial. The $ϵ$-algorithm reorganizes the same coefficients into a rational function, so it can represent the geometric patterns created by nearby singularities.

🏷️ Example: The First Padé Approximant of $e^z$

For

$$e^z=1+z+\frac{z^2}{2}+\cdots ,$$

the first three Taylor partial sums are

$$S_0=1,S_1=1+z,S_2=1+z+\frac{z^2}{2}.$$

Applying the first nontrivial $ϵ$ step gives

$${ϵ}_2^{(0)}=S_0-\frac{(S_1-S_0{)}^2}{S_2-2S_1+S_0}=1-\frac{z^2}{\frac{z^2}{2}-z}=\frac{1+\frac{z}{2}}{1-\frac{z}{2}}.$$

This is exactly the Padé approximant

$$[1/1{]}_{e^z}(z)=\frac{1+\frac{z}{2}}{1-\frac{z}{2}}.$$

The calculation is small, but it captures the entire philosophy: three polynomial partial sums have been reorganized into a rational approximant with a denominator. Higher even entries of the $ϵ$-table build the near-diagonal Padé approximants $[n+k/k]$ in the same way.

🏷️ When It Works

The $ϵ$-algorithm is most reliable when the remainder is well described by a few geometric components:

$$s_n-s\approx c_1{\lambda }_1^n+\cdots +c_k{\lambda }_k^n.$$

This includes many linearly convergent fixed-point iterations, many alternating or nearly geometric series, and Taylor series whose coefficients are dominated by a few nearby singularities. In the Padé setting, the method is especially natural for meromorphic functions, because poles are exactly the singularities rational functions are designed to represent.

It is less automatic for algebraic or logarithmic tails, such as

$$s_n=s+\frac{c}{n^{\alpha }}+o(n^{-\alpha }),$$

unless the local finite window happens to be well approximated by geometric modes. In such cases Richardson-type, Levin-type, or problem-specific transformations may be more appropriate (Brezinski & Redivo Zaglia, 1991).

Practical Stability

The transformed sequence may be worse than the original one if the input data contain noise, if the denominator determinant is small, or if the Padé approximant introduces a pole near the evaluation point. The algorithm should be monitored by comparing neighboring even entries, not by trusting a single high-order transform.

🏷️ Notes

• The first transform ${ϵ}_2^{(n)}$ is Aitken’s process; the higher transforms are recursive Shanks transforms.
• The even entries ${ϵ}_{2k}^{(n)}$ are the approximants. The odd entries are computational scaffolding.
• The determinant formula clarifies the exactness theorem, while Wynn’s recurrence is the efficient way to compute the same quantities.
• The Padé connection is not a metaphor: for Taylor partial sums, ${ϵ}_{2k}^{(n)}$ is exactly $[n+k/k]$ whenever the relevant denominators are nonzero.
• Epsilon extrapolation should be understood as rational extrapolation. It competes with polynomial extrapolation methods when the error has a rational or geometric structure.

🔗 See Also

• on Sylvester’s determinantal identity and Schweinsian expansion — The determinant identities behind condensation explain why ratios of neighboring determinants can be updated recursively instead of recomputed from scratch.
• on Q-D recursive Runge-Kutta methods — Extrapolation-based high-order constructions use related rational and recursive table ideas, especially in the comparison between order growth and stage count.
• on best uniform approximation from subspaces — Padé approximation is local moment matching, while best uniform approximation is global minimax approximation; the contrast is useful for judging what the $ϵ$-algorithm is and is not optimizing.

📚 References

🐻  Baker, G.A. & Graves-Morris, P. 1996. Padé Approximants 2nd ed., Cambridge University Press,p.

🐻  Brezinski, C. & Redivo Zaglia, M. 1991. Extrapolation Methods: Theory and Practice, North-Holland,p.

🐻  Gragg, W.B. 1972. The Padé Table and its Relation to Certain Algorithms of Numerical Analysis. SIAM Review 14(1), 1–62.

🐻  Shanks, D. 1955. Non-linear Transformations of Divergent and Slowly Convergent Sequences. Journal of Mathematics and Physics 34, 1–42.

🐻  Wynn, P. 1956. On a Device for Computing the e_m(Sₙ₎ Transformation. Mathematical Tables and Other Aids to Computation 10(54), 91–96.

🐻  Wynn, P. 1966. Upon Systems of Recursions which Obtain among the Quotients of the Padé Table. Numerische Mathematik 8(3), 264–269.