Trying to build inverse function

Question

* If you are doing intro2cs in huji be careful... I know why you're here *

I'm trying to build a function that get a function (assume continuity of monotone) as argument and return its inverse function. From mathematics I know that I need to reflect the function to y=x. but this isn't going very well.

I wrote something that will give me the x0 that makes f(x0)=0.

def solve(f, x0=-10000, x1=10000, epsilon=EPSILON):
"""return the solution to f in the range between x0 and x1"""

while x1-x0>EPSILON:

    if (f(x0)*f((x1+x0)/2))<0:
        x1=((x1+x0)/2)
    elif (f(x0)*f((x1+x0)/2))>0:
        x0=((x1+x0)/2)
    elif (f(x0)*f(x1))==0:
        if f(x0)==0:
            return x0
        else:
            return x1
if f(x0)*f(x1)<0:
    return (x1+x0)/2
else:
    return None

The problem is that I don't know what will be the range's input on the function. I think I need to work on small range at first, and if I won't find a solution i'll expand the range exponentially.

Any ideas?

UPDATE: ok so I tried to write it as @Spektre suggested and I succeed:

def inverse(g, epsilon=EPSILON):
    """return f s.t. f(g(x)) = x"""

    def helper(y):
        x0=-2
        x1=2
        while x1<10000:
            if solve((lambda x:g(x)-y),x0,x1,epsilon):
                return solve((lambda x:g(x)-y),x0,x1,epsilon)
            else:
                x0=math.pow(x0,3)
                x1=math.pow(x1,3)
    return helper

Answer 1

I just published a python package that does this precisely, and have to go through many of those caveats. You may want to borrow some ideas from it: https://pypi.python.org/pypi/pynverse

It essentially follows this strategy:

1. Figure out if the function is increasing or decreasing. For this two reference points ref1 and ref2 are needed:
   • In case of a finite interval, the points ref points are 1/4 and 3/4 through the interval.
   • In an infinite interval any two values work really.
   • If f(ref1) < f(ref2), the function is increasing, otherwise is decreasing.
2. Figure out the image of the function in the interval.
   • If values are provided, then those are used.
   • In a closed interval just calculate f(a) and f(b), where a and b are the ends of the interval.
   • In an open interval try to calculate f(a) and f(b), if this works those are used, otherwise it will be assume to be (-Inf, Inf).
3. Built a bounded function with the following conditions:
   • bounded_f(x):
     • return -Inf if x below interval, and f is increasing.
     • return +Inf if x below interval, and f is decreasing.
     • return +Inf if x above interval, and f is increasing.
     • return -Inf if x above interval, and f is decreasing.
     • return f(x) otherwise
4. If the required number y0 for the inverse is outside the image, raise an exception.
5. Find roots for bounded_f(x)-y0, by minimizing (bounded_f(x)-y0)**2, using the Brent method, making sure that the algorithm for minimising starts in a point inside the original interval by setting ref1, ref2 as brackets. As soon as if goes outside the allowed intervals, bounded_f returns infinite, forcing the algorithm to go back to search inside the interval.
6. Check that the solutions are accurate and they meet f(x0)=y0 to some desired precision, raising a warning otherwise.

Answer 2

without the algebraic expression of y=f(x)

There is no way to know the actual x or y interval.

As you mentioned you need continuous and monotonous interval only

So if you do not know it You can try to find extremal points (where first derivation cross zero) with some search with dynamic step dependent on the derivation value (linearly or logarithmic) but there will always be a risk of missing some small bumps. If you create the list of these points then the intervals between them should be monotonous. After this if the x=g(y) is function you can use binary search. If small bumps are present then approximation search is better.

For unbound monotone function y=f(x)

you can change the approximation search from above link to something like this:

1. initial guess

   x0=0
   

   there is not much we can do without more info on g,f. if x0=0 is not an option use any other safe value.

2. detect search step dx

   dx=0.0001 // or dx=function(|f(x0),f(x0+0.0001)|)
   if (|f(x0-dx)-y|<|f(x0+dx)-y|) dx=-dx;
   

   so find which way is closer to the solution. The 0.0001 constant should be selected to some meaningful manner. Too small will slow the process too big will miss bumps. For dynamic step you can use the magnitude of first derivation or the distance to solution itself like dx=sign(dx)*abs(f(x)-y) and if dx==0 stop ...

3. search the closest match

   while (|f(x0)-y|>|f(x0+dx)-y|) x0+=dx;
   

   just stop on the closest match. If you want the dynamic step then add also the dx=function(|f(x0),f(x0+0.0001)|) inside loop

4. now you have bound search <x0-dx,x0+dx>

   so you can use binary approximation search from above link as is.

5. if not close enough solution found

   then you got some missed bumps or too big step or initial guess point is extremal point with symetric function shape causing the initial direction estimate to fail. So you can change the initial x0 constant.

Answer 3

I would suggest working on small range first, than expanding it exponentially. If you use ranges like [-1, 1], [-2, 2], [-4, 4] and so on, the overhead will be constant even in worst case.