Best way to invert numerically a function and to evaluate it on a large array

Hello, I’m a new user of Julia.

I want to invert numerically a function and to compute the inverse on a ‘large’ array.
For the moment I studied only two options: using a Newton-Rhapson algorithm (I know the first derivative analytically) and using interpolation.
The function I have to invert is monotonically increasing and it is called luminosity_distance_at_z in the example below.

using BenchmarkTools
using Interpolations
using QuadGK

function NewtonRhapson(f, df, y0; x_guess=1.0, tol=1e-6, maxIter = 1000)
    # found x0 : f(x0) = y0
    ff(x) = f(x) - y0                 
    x = x_guess
    ffx = ff(x)
    iter = 0

    while abs(ffx) > tol && iter < maxIter
        x -= ffx/df(x)        # Iteration
        ffx = ff(x)           # Precompute ff(x)
        iter += 1
    return x

function inverse_interpolation(f, y0s)
    xs  = vcat( range(start = 0, stop = 9.99e-9, length = 125),
                10 .^(range(start = -8, stop = log10(7.99), length = 1000)),
                10 .^(range(start = log10(8), stop = 8, length = 125)) )    

    ys = f.(xs)
    inv = interpolate((ys,), xs, Gridded(Linear()))
    return inv.(y0s)
function E_at_z(z::Real)
    sqrt(0.3*(1+z)^3 + 0.69991 + 9e-5*(1+z)^4)

function integral_of_E_inv(z::Real)
    QuadGK.quadgk(z -> 1/E_at_z(z), 0., z)[1]  

const c_light = 299792.458 # km/s

function comoving_radial_distance_at_z(z::Real)
    dH = c_light/67
    return dH * integral_of_E_inv(z)  

function differential_luminosity_distance_at_z(z::Real)
    dH = c_light/67
    comoving_radial_distance_at_z(z) + (dH/E_at_z(z))*(1+z)

function luminosity_distance_at_z(z::Real)

When I have to evaluate the inverse function on a single point, the NR algorithm wins.
When I have to evaluate the inverse on a ‘large’ set of points, the interpolation algorithm seems to perform better

I think that the difference is due to the fact that the interpolation method computes f(x) 1250 times (in this example), while the NR algorithm computes f(x) two times each step needed to evaluate each point of the array (and a few steps per points are usually needed).

Since the NR is more precise, are there any ways to make it as fast as the interpolation one when evaluated on a large set of points?
Since I’m new I thought that maybe I’m using wrongly the vectorization of the NR algorithm via broadcast and there could be more efficient ways to vectorize it.

I would go the opposite direction and make the interpolation more accurate.

In particular, if you have a smooth function and evaluate at Chebyshev points (not equally spaced), then polynomial interpolation is exponentially accurate—for a simple monotonic function (not too wiggly) you will probably achieve machine precision with only a handful of function evaluations. Once you have a polynomial interpolation accurate to machine precision, you can invert it easily. (e.g. you could then find the Chebyshev interpolant of the inverse function.)

There are various packages to help you with this. For example, ApproxFun or FastChebInterp.


Orthogonal to @stevengj answer, you may want to use Cosmology.jl instead of implementing luminosity_distance from scratch.
And Roots.jl for function inversion instead of coding your own newton methods.


Sorry but I’m not sure I have understood correctly.

To use Chebyshev interpolation as an inversion algorithm, I should first approximate the function I want to invert with a polynomial. Then, what do I have to do to invert such polynomial ‘easily’? Evaluate it on the Chebyshev points of the function co-domain and then?

Calculate the inverse on the Chebyshev points (appropriately transformed) of the codomain of the original problem. In other words, just treat the inverse like any other function you would approximate. You can calculate it with a Newton method or any other rootfinding algorithm for these points.

1 Like

So, something like that?

using FastChebInterp

f(x) = GW_distance_at_z(c,x) # the function to invert
g(x) = f(x) - dist # I want to solve g(x) = 0

yi = f(0)
yf = f(4.5)

y = chebpoints(10, yi, yf) # Compute Chebyshev points on the original co-domain
inv_inter = chebinterp(g.(y), yi, yf) # Chebyshev interpolation of g(x)

# then I have to find the roots of g(x) = 0?

If so, what is the advantage over a traditional NR algorithm and the the naive interpolation I implemented in the first post?

The assumption here is that your original function f(x) is expensive, so that a polynomial interpolant is vastly cheaper to evaluate. So, you:

  1. Evaluate f(x) at a set of Chebyshev points to obtain an interpolating polynomial p(x). This is the only time you evaluate your expensive function f, and you only need a relatively small number of evaluations because of the exponential convergence.
  2. In the codomain y, for a set of Chebyshev points y_k, you solve p(x) = y_k to obtain x = p^{-1}(y_k) by whatever method, e.g. Newton. Then you interpolate these p^{-1}(y_k) values with a Chebyshev polynomial q(y). This is fast because it only involves polynomial evaluations.
  3. Your q(y) is now an approximation for f^{-1}(y). Evaluate q(y) as desired. This is fast because it only involves polynomial evaluations.

Compared to piecewise-polynomial interpolation ala Interpolations.jl, this should require many fewer function evaluations in step (1) to obtain a high-accuracy interpolant; you can easily approach the limits of roundoff errors. (This is only for smooth f that you can evaluate at arbitrary points, as is the case here. Piecewise interpolation is useful for non-smooth functions or data that is provided at a fixed set of points that you can’t change to Chebyshev points.)


Thank you so much! Now it’s clear and worked perfectly.