Richardson package for Julia

The Richardson package provides a function extrapolate that extrapolates any function f(x) to f(x0), evaluating f only at a geometric sequence of points > x0 (or optionally < x0) or at a given sequence of points. f(x) can return scalars, vectors, or any type implementing a normed vector space.

The key algorithm is Richardson extrapolation using a Neville–Aitken tableau, which adaptively increases the degree of an extrapolation polynomial until convergence is achieved to a desired tolerance (or convergence stalls due to e.g. floating-point errors). This allows one to obtain f(x0) to high-order accuracy, assuming that f(x0+h) has a Taylor series or some other power series in h. (See e.g. these course notes by Prof. Flaherty at RPI.)

Usage

extrapolate(f, h; contract=0.125, x0=zero(h), power=1,
                  atol=0, rtol=atol>0 ? 0 : sqrt(ε), maxeval=typemax(Int), breaktol=2)

Extrapolate f(x) to f₀ ≈ f(x0), evaluating f only at x > x0 points (or x < x0 if h < 0) using Richardson extrapolation starting at x=x₀+h. It returns a tuple (f₀, err) of the estimated f(x0) and an error estimate.

The return value of f can be any type supporting ± and norm operations (i.e. a normed vector space). Similarly, h and x0 can be in any normed vector space, in which case extrapolate performs Richardson extrapolation of f(x0+s*h) to s=0⁺ (i.e. it takes the limit as x goes to x0 along the h direction).

On each step of Richardson extrapolation, it shrinks x-x0 by a factor of contract, stopping when the estimated error is < max(rtol*norm(f₀), atol), when the estimated error increases by more than breaktol (e.g. due to numerical errors in the computation of f), when f returns a non-finite value (NaN or Inf), or when f has been evaluated maxeval times. Note that if the function may converge to zero, you may want specify a nonzero atol (which cannot be set by default because it depends on the scale/units of f); alternatively, in such cases extrapolate will halt when it becomes limited by the floating-point precision. (Passing breaktol=Inf can be useful to force extrapolate to continue shrinking h even if polynomial extrapolation is initially failing to converge, possibly at the cost of extraneous function evaluations.)

If x0 = ±∞ (±Inf), then extrapolate computes the limit of f(x) as x ⟶ ±∞ using geometrically increasing values of h (by factors of 1/contract).

In general, the starting h should be large enough that f(x0+h) can be computed accurately and efficiently (e.g. without severe cancellation errors), but small enough that f does not oscillate much between x0 and x0+h. i.e. h should be a typical scale over which the function f varies significantly.

Technically, Richardson extrapolation assumes that f(x0+h) can be expanded in a power series in h^power, where the default power=1 corresponds to an ordinary Taylor series (i.e. assuming f is analytic at x0). If this is not true, you may obtain slow convergence from extrapolate, but you can pass a different value of power (e.g. power=0.5) if your f has some different (Puiseux) power-series expansion. Conversely, if f is an even function around x0, i.e. f(x0+h) == f(x0-h), so that its Taylor series contains only even powers of h, you can accelerate convergence by passing power=2.

extrapolate(fh_itr; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf)

Similar to extrapolate(f, h), performs Richardson extrapolation of a sequence of values f(h) to h → 0, but takes an iterable collection fh_itr of a sequence of (f(h), h) tuples (in order of decreasing |h|).

There is no contract keyword argument since the contraction factors are determined by the sequence of h values (which need not contract by the same amount). The tolerances atol and rtol both default to 0 so that by default it examines all of the values in the fh_itr collection. Otherwise, the keyword arguments have the same meanings as in extrapolate(f, h).

extrapolate!(fh::AbstractVector; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf)

Similar to extrapolate(fh_itr), performs Richardson extrapolation on an array fh of (f(h), h) tuples (in order of decreasing |h|), but overwrites the array fh in-place with intermediate calculations.

(Thus, the array fh must be a vector of Tuple{T,H} values, where H<:Number is the type of h and T is the type of the extrapolated f(0) result. This T should be a floating-point type, i.e. fh should contain float(f(h)) if the function you are extrapolating is not already floating-point-valued.)

Examples

For example, let's extrapolate sin(x)/x to x=0 (where the correct answer is 1) starting at x=1, printing out the x value at each step so that we can see what the algorithm is doing.

(Since f is passed as the first argument to extrapolate, we can use Julia's do syntax to conveniently define a multi-line anonymous function to pass.)

extrapolate(1.0, rtol=1e-10) do x
    @show x
    sin(x)/x
end

giving the output:

x = 1.0
x = 0.125
x = 0.015625
x = 0.001953125
x = 0.000244140625
x = 3.0517578125e-5
(1.0000000000000002, 2.0838886172214188e-13)

That is, it evaluates our function sin(x)/x for 6 different values of x and returns 1.0000000000000002, which is accurate to machine precision (the error is ≈ 2.2e-16). The returned error estimate of 2e-13 is conservative, which is typical for extrapolating well-behaved functions.

Since sin(x)/x is an even (symmetric) function around x=0, its Taylor series contains only even powers of x. We can exploit this fact to accelerate convergence for even functions by passing power=2 to extrapolate:

extrapolate(1.0, rtol=1e-10, power=2) do x
    @show x
    sin(x)/x
end

gives

x = 1.0
x = 0.125
x = 0.015625
x = 0.001953125
x = 0.000244140625
(1.0, 0.0)

which converged to machine precision (in fact, the exact result) in only 5 function evaluations (1 fewer than above).

Infinite limits

Using the x0 keyword argument, we can compute the limit of f(x) as x ⟶ x0. In fact, you can pass x0 = Inf to compute a limit as x ⟶ ∞ (which is accomplished internally by a change of variables x = 1/u and performing Richardson extrapolation to u=0). For example:

extrapolate(1.0, x0=Inf) do x
    @show x
    (x^2 + 3x - 2) / (x^2 + 5)
end

gives

x = 1.0
x = 8.0
x = 64.0
x = 512.0
x = 4096.0
x = 32768.0
x = 262144.0
(1.0000000000000002, 1.2938539128981574e-12)

which is the correct result (1.0) to machine precision.

Extrapolating series

One nice application of infinite limits is extrapolating infinite series. If we start with an integer x, then the default contract=0.125 will increase x by a factor of 8.0 on each iteration, so x will always be an exact integer and we can use it as the number of terms in a series.

For example, suppose we are computing the infinite series 1/1² + 1/2² + 1/3² + ⋯. This is the famous Basel problem, and it converges to π²/6 ≈ 1.644934066848226…. If we compute it by brute force, however, we need quite a few terms to get high accuracy:

julia> sum(n -> 1/n^2, 1:100) - π^2/6
-0.009950166663333482

julia> sum(n -> 1/n^2, 1:10^4) - π^2/6
-9.99950001654426e-5

julia> sum(n -> 1/n^2, 1:10^9) - π^2/6
-9.999985284281365e-10

Even with 10⁹ terms we get only about 9 digits. Instead, we can use extrapolate (starting at 1 term):

julia> val, err = extrapolate(1, x0=Inf) do N
           @show N
           sum(n -> 1/n^2, 1:Int(N))
       end
N = 1.0
N = 8.0
N = 64.0
N = 512.0
N = 4096.0
N = 32768.0
(1.6449340668482288, 4.384936858059518e-12)

julia> (val - π^2/6)/(π^2/6)
1.4848562646983628e-15

By 32768 terms, the extrapolated value is accurate to about 15 digits.

Numerical derivatives

A classic application of Richardson extrapolation is the accurate evaluation of derivatives via finite-difference approximations (although analytical derivatives, e.g. by automatic differentiation, are typically vastly more efficient when they are available). In this example, we use Richardson extrapolation on the forward-difference approximation f'(x) ≈ (f(x+h)-f(x))/h, for which the error decreases as O(h) but a naive application to a very small h will yield a huge cancellation error from floating-point roundoff effects. We differentiate f(x)=sin(x) at x=1, for which the correct answer is cos(1) ≈ 0.5403023058681397174009366..., starting with h=0.1

extrapolate(0.1, rtol=0) do h
    @show h
    (sin(1+h) - sin(1)) / h
end

Although we gave an rtol of 0, the extrapolate function will terminate after a finite number of steps when it detects that improvements are limited by floating-point error:

h = 0.1
h = 0.0125
h = 0.0015625
h = 0.0001953125
h = 2.44140625e-5
h = 3.0517578125e-6
(0.5403023058683176, 1.7075230118734908e-12)

The output 0.5403023058683176 differs from cos(1) by ≈ 1.779e-13, so in this case the returned error estimate is only a little conservative. Unlike the sin(x)/x example, extrapolate is not able to attain machine precision (the floating-point cancellation error in this function is quite severe for small h!), but it is able to get surprisingly close.

Another possibility for a finite-difference/Richardson combination was suggested by Ridders (1982), who computed both f'(x) and f''(x) (the first and second derivatives) simultaneously using a center-difference approximation, which requires two new f(x) evaluations for each h. In particular, the center-difference approximations are f'(x) ≈ (f(x+h)-f(x-h))/2h and f''(x) ≈ (f(x+h)-2f(x)+f(x-h))/h², both of which have errors that go as O(h²). We can plug both of these functions simultaneously into extrapolate (so that they share f(x±h) evaluations) by using a vector-valued function returning [f', f'']. Moreover, since these center-difference approximations are even functions of h (identical for ±h), we can pass power=2 to extrapolate in order to exploit the even-power Taylor expansion. Here is a function implementing both of these ideas:

# returns (f'(x), f''(x))
function riddersderiv2(f, x, h; atol=0, rtol=atol>0 ? 0 : sqrt(eps(typeof(float(real(x+h))))), contract=0.5)
    f₀ = f(x)
    val, err = extrapolate(h, atol=atol, rtol=rtol, contract=contract, power=2) do h
        f₊, f₋ = f(x+h), f(x-h)
        [(f₊-f₋)/2h, (f₊-2f₀+f₋)/h^2]
    end
    return val[1], val[2]
end

(This code could be made even more efficient by using StaticArrays.jl for the [f', f''] vector.) The original paper by Ridders accomplishes something similar in < 20 lines of TI-59 calculator code, by the way; so much for high-level languages!

For example,

julia> riddersderiv2(1, 0.1, rtol=0) do x
           @show x
           sin(x)
       end
x = 1
x = 1.1
x = 0.9
x = 1.05
x = 0.95
x = 1.025
x = 0.975
x = 1.0125
x = 0.9875
x = 1.00625
x = 0.99375
(0.5403023058681394, -0.841470984807975)

evaluates the first and second derivatives of sin(x) at x=1 and obtains the correct answer (cos(1), -sin(1)) to about 15 and 13 decimal digits, respectively, using 11 function evaluations.

Handling problematic convergence

It is useful to consider a finite-difference approximation for the derivative of the function 1/x at some x ≠ 0: i.e. computing the limit of f(h) = (1/(x+h) - 1/x) / h as h goes to zero similar to above.

This function f(h) has a pole at h=-x, i.e. f(-x) blows up. This means that the Taylor series of f(h) only converges for h values small enough to avoid this pole, and in fact for |h| < |x|. Since Richardson extrapolation is essentially approximating the Taylor series, this means that the extrapolation process doesn't converge if the starting h is too large, and extrapolation will give up and halt with the wrong answer.

This lack of convergence is easily observed: set x=0.01 (where the correct derivative of 1/x is -10000) and consider what happens for a starting h that is too large compared to |x|:

julia> extrapolate(1.0) do h
           @show h
           x = 0.01
           (1/(x+h) - 1/x) / h
       end
h = 1.0
h = 0.125
h = 0.015625
(-832.4165749908325, 733.4066740007335)

Before reaching an |h| < 0.01 where the power series could begin to converge, extrapolate gave up and returned a wrong answer (with a large error estimate to let you know that the result is garbage)! In contrast, if we start with a small enough h then it converges just fine and returns the correct answer (-10000) to nearly machine precision:

julia> extrapolate(0.01) do h
           @show h
           x = 0.01
           (1/(x+h) - 1/x) / h
       end
h = 0.01
h = 0.00125
h = 0.00015625
h = 1.953125e-5
h = 2.44140625e-6
h = 3.0517578125e-7
(-10000.000000000211, 4.066770998178981e-6)

Of course, if you know that your function blows up like this, it is easy to choose good initial h, but how can we persuade extrapolate to do a better job automatically?

The trick is to use the breaktol keyword argument. breaktol defaults to 2, which means that extrapolate gives up if the best error estimate increases by more than a factor of 2 from one iteration to the next. Ordinarily, this kind of breakdown in convergence arises because you hit the limits of floating-point precision, and halting extrapolation is the right thing to do. Here, however, it would converge if we just continued shrinking h. So, we simply set breaktol=Inf to force extrapolation to continue, which works even for a large initial h=1000.0:

julia> extrapolate(1000.0, breaktol=Inf) do h
           @show h
           x = 0.01
           (1/(x+h) - 1/x) / h
       end
h = 1000.0
h = 125.0
h = 15.625
h = 1.953125
h = 0.244140625
h = 0.030517578125
h = 0.003814697265625
h = 0.000476837158203125
h = 5.9604644775390625e-5
h = 7.450580596923828e-6
h = 9.313225746154785e-7
h = 1.1641532182693481e-7
(-10000.000000029328, 5.8616933529265225e-8)

Not that it continues extrapolating until it reaches small h values where the power-series converges, and in the end it again returns the correct answer to nearly machine precision (and would reach machine precision if we set a smaller rtol). (The extrapolate function automatically discards the initial points where the polynomial extrapolation fails.)