Richardson.jl: Richardson extrapolation in Julia

Package Announcements
1

Richardson extrapolation is a technique for extrapolating a given function f(x) to f(x_0), taking a one-sided x \to x_0 limit without evaluating f(x_0) explicitly, to high-order accuracy using adaptive polynomial extrapolation.

I couldn’t find any library in Julia for this, so I put together a little package:

I plan to register it soon (perhaps under JuliaMath). (Please let me know if this duplicates functionality already implemented elsewhere!)

24 Likes
2

There’s ODE solver versions of Richardson extrapolation in OrdinaryDiffEq: https://docs.sciml.ai/latest/solvers/ode_solve/#Parallelized-Explicit-Extrapolation-Methods-1 but while it’s the same general method, that seems to be a fairly different use case from what you’re doing here.

3

I coded simpler (and MUCH wronger) versions of this general code countless times…

Thank you Steven for the marvellous work! register it ASAP, someone much wiser than me said: “its usefulness for practical computations can hardly be overestimated.” :smiley:

EDIT 1:
It can be useful to protect tricky computations in AD:

using Richardson
using ForwardDiff: Dual
f(x) = cos(x + sqrt(x) * x)

A,_=extrapolate(Dual(0.1,0.0),x0=Dual(0.0,1.0)) do x
    @show x
    f(x)
end
#=
x = Dual{Nothing}(0.1,2.0)
x = Dual{Nothing}(0.0125,1.125)
x = Dual{Nothing}(0.0015625,1.015625)
x = Dual{Nothing}(0.0001953125,1.001953125)
x = Dual{Nothing}(2.44140625e-5,1.000244140625)
x = Dual{Nothing}(3.0517578125e-6,1.000030517578125)
=#
B=f( Dual(0.0,1.0) )

@show A # Dual{Nothing}(0.9999999999995401,1.740011327412078e-8)
@show B # Dual{Nothing}(1.0,NaN))

The following is now fixed/supported! :smiley:

#EDIT 2: via super ugly, surely broken approach it seems to work for vector valued functions, is there a way to truly support vector functions? 
#```julia
using LinearAlgebra: norm
import Base.abs
import Base.isless
Base.abs(x::AbstractVector) = norm(x)
Base.isless(y::Number,x::AbstractVector) = y<norm(x)
extrapolate(-0.1,x0=0) do x
   @show x
   [x/sqrt(-x),sin(x)/x]
end
# ([-4.76091483644921e-9, 1.0], 8.704985825570316e-9)
#```

#EDIT 3: limits in R^N can be easily supported i think

# ```julia
using Richardson

@inline function Richardson.extrapolate(f,h::AbstractArray{T};x0::AbstractArray{T},args...) where T
   @boundscheck size(h)==size(x0)
   @inline g(t)=@inbounds x0.+t.*h
   extrapolate(T(1),x0=T(0);args...) do t
       f(g(t))
   end
end

extrapolate([0.1,-0.2],x0=[0.0,1.0]) do x
   @show x
   sin(x[1])/x[1]*log(x[2])/(x[2]-1)
end
#=
x = [0.1, 0.8]
x = [0.0125, 0.975]
x = [0.0015625, 0.996875]
x = [0.0001953125, 0.999609375]
x = [2.44140625e-5, 0.999951171875]
=#

#(1.000000000000058, 1.4038002982275088e-9)
#```

#EDIT 4: for some functions it can go into endless loops (valid up to commit: https://github.com/JuliaMath/Richardson.jl/commit/c828d05a0b07fde094092496eecdea0e1c4588c9)

#```julia
using Richardson
extrapolate(0.2,x0=1.0) do x
   @show x # this is there to slow it down just enough to be able to abort the computation easily
   log(x)/sqrt(x-1)
end
#```
4

Thanks for the comments, but please file bug reports and feature requests as github issues rather than posting here.

5

There are two issues here. First, if the limit could be zero then you need to specify an absolute tolerance (atol); we can’t have a default atol because it is dimensionful.

Second, Richardson extrapolation with polynomials assumes an analytic function, i.e. a Taylor series, which does not hold for \log(x) / \sqrt{x-1}, so you will get rather slow convergence (and eventually NaN values). If you change variables to x = 1 + u^2, and hence use \log(1+u^2) / u or log1p(u^2) / u, then it is analytic in u and you recover rapid convergence.

Technically, we could do Richardson extrapolation for an arbitrary Puiseux series if you specify the power law, but I haven’t implemented this yet.

3 Likes
6

ofcourse, the example was chosen on purpose \log(1+\epsilon)/\sqrt\epsilon \propto \sqrt\epsilon, slow convergence was expected, a parametrized cutoff on the number of iterations is enough i think. fixed! :smiley:

i thank you again for filling this gap in the ecosystem, i already modified some of my code to use Richardson.jl, (i was using a fixed order extrapolation) and with your package the accuracy is much better :slight_smile:

sorry if i posted all that shebang here, the first post started of as a simple thank you.

1 Like
7

That’s a reasonable feature, I’ll add something like that.

However, I’ve fixed it to terminate on your example (basically, it needs to terminate if it hits a NaN). I also fixed support for vector-valued functions (or any type that supports a norm).

4 Likes
8

Now implemented. Just pass power=0.5 and you get decent convergence for your example function (which has a Puiseux series in \sqrt{h}):

julia> extrapolate(0.2, x0=1.0, power=0.5) do x
           @show x
           log(x)/sqrt(x-1)
       end
x = 2.2
x = 1.15
x = 1.01875
x = 1.00234375
x = 1.00029296875
x = 1.00003662109375
x = 1.0000045776367188
x = 1.00000057220459
x = 1.0000000715255737
x = 1.0000000089406966
(-2.1650306436709847e-13, 2.2837517952902017e-13)

As an added bonus, if you have an even/symmetric function f(x_0 + h) = f(x_0 - h), which has only even powers of h in its Taylor series, then you can pass power=2 to accelerate convergence.

7 Likes