Aitken's delta-squared process

Properties

Aitken's delta-squared process is an acceleration of convergence method and a particular case of a nonlinear sequence transformation.

A sequence $X=(x_n)$ that converges to a limiting value $\ell$ is said to converge linearly, or more technically Q-linearly, if there is some number $\mu \in (0,1)$ for which

$${\displaystyle \underset{n\to \mathrm{\infty }}{lim}\frac{|x_{n+1}-\ell |}{|x_n-\ell |}=\mu .}$$

This means that asymptotically, the distance between the sequence and its limit shrinks by nearly the same proportion, ${\displaystyle \mu ,}$ on every step and the ratio of reduction becomes closer and closer to that proportion. This is also sometimes called "geometric convergence," since it is a characteristic property for geometric series, or "exponential convergence," since it is convergence like ${\displaystyle {\mu }^n=\exp (n\ln \mu ).}$

Aitken's method will accelerate the convergence of a sequence ${\displaystyle X}$ if ${\displaystyle A[X]=(a_n),}$ with terms defined above, satisfies $\underset{n\to \mathrm{\infty }}{lim}\frac{a_n-\ell }{x_n-\ell }=0.$

${\displaystyle A}$ is not a linear operator on sequences, but it is linear with respect to addition of constant sequences: $A[X-C]=A[X]-C,$ if ${\displaystyle C}$ is any constant sequence ${\displaystyle C=(c)}$, constant for all ${\displaystyle n.}$ This is clear from the expression of ${\displaystyle A[X]}$ in terms of the finite difference operator ${\displaystyle \mathrm{\Delta }.}$

The new process does not in general converge quadratically, but for an iterated function sequence satisfying ${\displaystyle x_{n+1}=f(x_n)}$ for some function ${\displaystyle f}$ converging to a fixed point, the accelerated sequence's convergence is quadratic. In this case, the technique is known as Steffensen's method.

Empirically, the A-operation eliminates the "most important error term". One can check this by considering a sequence of the form ${\displaystyle x_n=\ell +a^n+b^n}$, where ${\displaystyle 0<b<a<1}$: The sequence ${\displaystyle A[X]}$ will then go to the limit $\ell$ like ${\displaystyle b^n}$ goes to zero.

Geometrically, the graph of an exponential function ${\displaystyle f(t)}$ that satisfies ${\displaystyle f(n)=x_n}$, ${\displaystyle f(n+1)=x_{n+1}}$ and ${\displaystyle f(n+2)=x_{n+2}}$ has an horizontal asymptote at ${\displaystyle \frac{x_n{x}_{n+2}-x_{n+1}^2}{x_n-2x_{n+1}+x_{n+2}}}$ (if ${\displaystyle x_n-2x_{n+1}+x_{n+2}\ne 0}$).

One can also show that if a sequence ${\displaystyle X}$ converges to its limit ${\displaystyle \ell }$ at a rate strictly greater than 1, ${\displaystyle A[X]}$ does not have a better rate of convergence. (In practice, one rarely has e.g. quadratic convergence which would mean over 30 (respectively 100) correct decimal places after 5 (respectively 7) iterations (starting with 1 correct digit); usually no acceleration is needed in that case.)

In practice, ${\displaystyle A[X]}$ often converges much faster to the limit than ${\displaystyle X}$ does, as demonstrated by the example calculations below. Usually, it is much cheaper to calculate ${\displaystyle A[X]}$ (involving only calculation of differences, one multiplication and one division) than to calculate many more terms of the sequence ${\displaystyle X}$. Care must be taken, however, to avoid introducing errors due to insufficient precision when calculating the differences in the numerator and denominator of the expression.

Example calculations

Example 1: The value of ${\displaystyle \sqrt{2}\approx 1.4142136}$ can be approximated by assuming an initial value for ${\displaystyle x_0}$ and iterating the following sequence, called Heron's method: $${\displaystyle x_{n+1}=\frac{x_n+\frac{2}{x_n}}{2}.}$$ Starting with ${\displaystyle x_0=1:}$

It is worth noting here that Aitken's method does not save the cost of calculating two iterations here; computation of the first three $A[X]$ values required the first five $X$ values. Also, the second $A[X]$ value is less accurate than the 4th $X$ value, which is not surprising due to the fact that Aitken's process is best suited for sequences that converge linearly, rather than quadratically, and Heron's method for calculating square roots converges quadratically.

Example 2: The value of ${\displaystyle \frac{\pi }{4}}$ may be calculated as an infinite sum via the Leibniz formula for π:

$${\displaystyle \frac{\pi }{4}=\sum _{n=0}^{\mathrm{\infty }}\frac{(-1{)}^n}{2n+1}\approx 0.785398}$$

In this example, Aitken's method is applied to a sublinearly converging series and accelerates convergence considerably. The convergence is still sublinear, but much faster than the original convergence: the first $A[X]$ value, whose computation required the first three $X$ values, is closer to the limit than the eighth $X$ value.

Example pseudocode for Aitken extrapolation

The following is an example of using the Aitken extrapolation to help find the limit of the sequence ${\displaystyle x_{n+1}=f(x_n)}$ when given some initial ${\displaystyle x_0,}$ where the limit of this sequence is assumed to be a fixed point ${\displaystyle f}$ (say ${\displaystyle \alpha =f(\alpha )}$). For instance, if the sequence is given by $x_{n+1}=\frac{1}{2}\left(x_n+\frac{2}{x_n}\right)$ with starting point ${\displaystyle x_0=1,}$ then the function will be $f(x):=\frac{1}{2}\left(x+\frac{2}{x}\right),$ which has ${\displaystyle \alpha :=\sqrt{2}}$ as a fixed point (see Methods of computing square roots); it is this fixed point whose value will be approximated.

This pseudo code also computes the Aitken approximation to ${\displaystyle f^{\mathrm{\prime }}(\alpha )}$. The Aitken extrapolates will be denoted by aitkenX. During the computation of the extrapolate, it is important to check if the denominator becomes too small, which could happen if we already have a large amount of accuracy; without this check, a large amount of error could be introduced by the division. This small number will be denoted by epsilon. Because the binary representation of the fixed point could be infinite (or at least too large to fit in the available memory), the calculation will stop once the approximation is within tolerance of the true value.

<span id="mw8g"></span><span class="c" id="mw8w">%These choices depend on the problem being solved</span> <span class="n" id="mw9A">x0</span><span class="w" id="mw9Q"> </span><span class="p" id="mw9g">=</span><span class="w" id="mw9w"> </span><span class="mi" id="mw-A">1</span><span class="w" id="mw-Q"> </span><span class="c" id="mw-g">%The initial value</span> <span class="n" id="mw-w">f</span><span class="p" id="mw_A">(</span><span class="n" id="mw_Q">x</span><span class="p" id="mw_g">)</span><span class="w" id="mw_w"> </span><span class="p" id="mwAQA">=</span><span class="w" id="mwAQE"> </span><span class="p" id="mwAQI">(</span><span class="mi" id="mwAQM">1</span><span class="o" id="mwAQQ">/</span><span class="mi" id="mwAQU">2</span><span class="p" id="mwAQY">)</span><span class="o" id="mwAQc">*</span><span class="p" id="mwAQg">(</span><span class="n" id="mwAQk">x</span><span class="w" id="mwAQo"> </span><span class="o" id="mwAQs">+</span><span class="w" id="mwAQw"> </span><span class="mi" id="mwAQ0">2</span><span class="o" id="mwAQ4">/</span><span class="n" id="mwAQ8">x</span><span class="p" id="mwARA">)</span><span class="w" id="mwARE"> </span><span class="c" id="mwARI">%The function that finds the next element in the sequence</span> <span class="n" id="mwARM">tolerance</span><span class="w" id="mwARQ"> </span><span class="p" id="mwARU">=</span><span class="w" id="mwARY"> </span><span class="mi" id="mwARc">10</span>^<span class="o" id="mwARg">-</span><span class="mi" id="mwARk">10</span><span class="w" id="mwARo"> </span><span class="c" id="mwARs">%10 digit accuracy is desired</span> <span class="n" id="mwARw">epsilon</span><span class="w" id="mwAR0"> </span><span class="p" id="mwAR4">=</span><span class="w" id="mwAR8"> </span><span class="mi" id="mwASA">10</span>^<span class="o" id="mwASE">-</span><span class="mi" id="mwASI">16</span><span class="w" id="mwASM"> </span><span class="c" id="mwASQ">%Do not divide by a number smaller than this</span> <span class="n" id="mwASU">maxIterations</span><span class="w" id="mwASY"> </span><span class="p" id="mwASc">=</span><span class="w" id="mwASg"> </span><span class="mi" id="mwASk">20</span><span class="w" id="mwASo"> </span><span class="c" id="mwASs">%Do not allow the iterations to continue indefinitely</span> <span class="n" id="mwASw">haveWeFoundSolution</span><span class="w" id="mwAS0"> </span><span class="p" id="mwAS4">=</span><span class="w" id="mwAS8"> </span><span class="nb" id="mwATA">false</span><span class="w" id="mwATE"> </span><span class="c" id="mwATI">%Were we able to find the solution to within the desired tolerance? not yet</span> <span class="k" id="mwATM">for</span><span class="w" id="mwATQ"> </span><span class="nb" id="mwATU">i</span><span class="w" id="mwATY"> </span><span class="p" id="mwATc">=</span><span class="w" id="mwATg"> </span><span class="mi" id="mwATk">1</span><span class="w" id="mwATo"> </span><span class="p" id="mwATs">:</span><span class="w" id="mwATw"> </span><span class="n" id="mwAT0">maxIterations</span><span class="w" id="mwAT4"> </span> <span class="w" id="mwAT8"> </span><span class="n" id="mwAUA">x1</span><span class="w" id="mwAUE"> </span><span class="p" id="mwAUI">=</span><span class="w" id="mwAUM"> </span><span class="n" id="mwAUQ">f</span><span class="p" id="mwAUU">(</span><span class="n" id="mwAUY">x0</span><span class="p" id="mwAUc">)</span> <span class="w" id="mwAUg"> </span><span class="n" id="mwAUk">x2</span><span class="w" id="mwAUo"> </span><span class="p" id="mwAUs">=</span><span class="w" id="mwAUw"> </span><span class="n" id="mwAU0">f</span><span class="p" id="mwAU4">(</span><span class="n" id="mwAU8">x1</span><span class="p" id="mwAVA">)</span> <span class="w" id="mwAVE"> </span><span class="k" id="mwAVI">if</span><span class="w" id="mwAVM"> </span><span class="p" id="mwAVQ">(</span><span class="n" id="mwAVU">x1</span><span class="w" id="mwAVY"> </span><span class="o" id="mwAVc">~=</span><span class="w" id="mwAVg"> </span><span class="n" id="mwAVk">x0</span><span class="p" id="mwAVo">)</span> <span class="w" id="mwAVs"> </span><span class="n" id="mwAVw">lambda</span><span class="w" id="mwAV0"> </span><span class="p" id="mwAV4">=</span><span class="w" id="mwAV8"> </span><span class="n" id="mwAWA">absoluteValue</span><span class="p" id="mwAWE">((</span><span class="n" id="mwAWI">x2</span><span class="w" id="mwAWM"> </span><span class="o" id="mwAWQ">-</span><span class="w" id="mwAWU"> </span><span class="n" id="mwAWY">x1</span><span class="p" id="mwAWc">)</span><span class="o" id="mwAWg">/</span><span class="p" id="mwAWk">(</span><span class="n" id="mwAWo">x1</span><span class="w" id="mwAWs"> </span><span class="o" id="mwAWw">-</span><span class="w" id="mwAW0"> </span><span class="n" id="mwAW4">x0</span><span class="p" id="mwAW8">))</span><span class="w" id="mwAXA"> </span><span class="c" id="mwAXE">%OPTIONAL: Computes an approximation of |f'(fixedPoint)|, which is denoted by lambda</span> <span class="w" id="mwAXI"> </span><span class="k" id="mwAXM">end</span> <span class="w" id="mwAXQ"> </span><span class="n" id="mwAXU">denominator</span><span class="w" id="mwAXY"> </span><span class="p" id="mwAXc">=</span><span class="w" id="mwAXg"> </span><span class="p" id="mwAXk">(</span><span class="n" id="mwAXo">x2</span><span class="w" id="mwAXs"> </span><span class="o" id="mwAXw">-</span><span class="w" id="mwAX0"> </span><span class="n" id="mwAX4">x1</span><span class="p" id="mwAX8">)</span><span class="w" id="mwAYA"> </span><span class="o" id="mwAYE">-</span><span class="w" id="mwAYI"> </span><span class="p" id="mwAYM">(</span><span class="n" id="mwAYQ">x1</span><span class="w" id="mwAYU"> </span><span class="o" id="mwAYY">-</span><span class="w" id="mwAYc"> </span><span class="n" id="mwAYg">x0</span><span class="p" id="mwAYk">);</span> <span class="w" id="mwAYo"> </span><span class="k" id="mwAYs">if</span><span class="w" id="mwAYw"> </span><span class="p" id="mwAY0">(</span><span class="n" id="mwAY4">absoluteValue</span><span class="p" id="mwAY8">(</span><span class="n" id="mwAZA">denominator</span><span class="p" id="mwAZE">)</span><span class="w" id="mwAZI"> </span><span class="o" id="mwAZM"><</span><span class="w" id="mwAZQ"> </span><span class="n" id="mwAZU">epsilon</span><span class="p" id="mwAZY">)</span><span class="w" id="mwAZc"> </span><span class="c" id="mwAZg">%To avoid greatly increasing error, do not divide by too small of a number</span> <span class="w" id="mwAZk"> </span><span class="nb" id="mwAZo">print</span><span class="p" id="mwAZs">(</span><span class="s" id="mwAZw">'WARNING: denominator is too small'</span><span class="p" id="mwAZ0">)</span> <span class="w" id="mwAZ4"> </span><span class="k" id="mwAZ8">break</span><span class="w" id="mwAaA"> </span><span class="c" id="mwAaE">%Leave the loop</span> <span class="w" id="mwAaI"> </span><span class="k" id="mwAaM">end</span> <span class="w" id="mwAaQ"> </span><span class="n" id="mwAaU">aitkenX</span><span class="w" id="mwAaY"> </span><span class="p" id="mwAac">=</span><span class="w" id="mwAag"> </span><span class="n" id="mwAak">x2</span><span class="w" id="mwAao"> </span><span class="o" id="mwAas">-</span><span class="w" id="mwAaw"> </span><span class="p" id="mwAa0">(</span><span class="w" id="mwAa4"> </span><span class="p" id="mwAa8">(</span><span class="n" id="mwAbA">x2</span><span class="w" id="mwAbE"> </span><span class="o" id="mwAbI">-</span><span class="w" id="mwAbM"> </span><span class="n" id="mwAbQ">x1</span><span class="p" id="mwAbU">)</span>^<span class="mi" id="mwAbY">2</span><span class="w" id="mwAbc"> </span><span class="p" id="mwAbg">)</span><span class="o" id="mwAbk">/</span><span class="n" id="mwAbo">denominator</span> <span class="w" id="mwAbs"> </span> <span class="w" id="mwAbw"> </span><span class="k" id="mwAb0">if</span><span class="w" id="mwAb4"> </span><span class="p" id="mwAb8">(</span><span class="n" id="mwAcA">absoluteValue</span><span class="p" id="mwAcE">(</span><span class="n" id="mwAcI">aitkenX</span><span class="w" id="mwAcM"> </span><span class="o" id="mwAcQ">-</span><span class="w" id="mwAcU"> </span><span class="n" id="mwAcY">x2</span><span class="p" id="mwAcc">)</span><span class="w" id="mwAcg"> </span><span class="o" id="mwAck"><</span><span class="w" id="mwAco"> </span><span class="n" id="mwAcs">tolerance</span><span class="p" id="mwAcw">)</span><span class="w" id="mwAc0"> </span><span class="c" id="mwAc4">%If the value is within tolerance</span> <span class="w" id="mwAc8"> </span><span class="nb" id="mwAdA">print</span><span class="p" id="mwAdE">(</span><span class="s" id="mwAdI">"The fixed point is "</span><span class="p" id="mwAdM">,</span><span class="w" id="mwAdQ"> </span><span class="n" id="mwAdU">aitkenX</span><span class="p" id="mwAdY">))</span><span class="w" id="mwAdc"> </span><span class="c" id="mwAdg">%Display the result of the Aitken extrapolation</span> <span class="w" id="mwAdk"> </span><span class="n" id="mwAdo">haveWeFoundSolution</span><span class="w" id="mwAds"> </span><span class="p" id="mwAdw">=</span><span class="w" id="mwAd0"> </span><span class="nb" id="mwAd4">true</span> <span class="w" id="mwAd8"> </span><span class="k" id="mwAeA">break</span><span class="w" id="mwAeE"> </span><span class="c" id="mwAeI">%Done, so leave the loop</span> <span class="w" id="mwAeM"> </span><span class="k" id="mwAeQ">end</span> <span class="w" id="mwAeU"> </span><span class="n" id="mwAeY">x0</span><span class="w" id="mwAec"> </span><span class="p" id="mwAeg">=</span><span class="w" id="mwAek"> </span><span class="n" id="mwAeo">aitkenX</span><span class="w" id="mwAes"> </span><span class="c" id="mwAew">%Update x0 to start again </span> <span class="k" id="mwAe0">end</span> <span class="k" id="mwAe4">if</span><span class="w" id="mwAe8"> </span><span class="p" id="mwAfA">(</span><span class="n" id="mwAfE">haveWeFoundSolution</span><span class="w" id="mwAfI"> </span><span class="o" id="mwAfM">==</span><span class="w" id="mwAfQ"> </span><span class="nb" id="mwAfU">false</span><span class="p" id="mwAfY">)</span><span class="w" id="mwAfc"> </span><span class="c" id="mwAfg">%If we were not able to find a solution to within the desired tolerance</span> <span class="w" id="mwAfk"> </span><span class="nb" id="mwAfo">print</span><span class="p" id="mwAfs">(</span><span class="s" id="mwAfw">"Warning: Not able to find solution to within the desired tolerance of "</span><span class="p" id="mwAf0">,</span><span class="w" id="mwAf4"> </span><span class="n" id="mwAf8">tolerance</span><span class="p" id="mwAgA">)</span> <span class="w" id="mwAgE"> </span><span class="nb" id="mwAgI">print</span><span class="p" id="mwAgM">(</span><span class="s" id="mwAgQ">"The last computed extrapolate was "</span><span class="p" id="mwAgU">,</span><span class="w" id="mwAgY"> </span><span class="n" id="mwAgc">aitkenX</span><span class="p" id="mwAgg">)</span> <span class="k" id="mwAgk">end</span>