<< Computer Algorithms Use  Ctrl +  or  –  to enlarge or reduce text size.

Fast Function Inversion
Introduction

Problem:  You are given an invertable function F of a single variable over a finite interval. This function might take a long time to compute. Write a program that computes a good approximation to the function’s inverse at any point, and that is always fast. You are allowed one arbitrarily long setup computation.

In other words, write a program that takes a given y in the function’s range and quickly computes x so that F(x) = y to a good approximation. And this even though F might take a long time to compute.

The solution is to use function interpolation on the inverse, using some algorithm to compute the inversion table. Any of the usual root finding algorithms (bisection, secant, Brent’s, or – if the derivative is known – Newton’s) could be used for the latter. Since this computation needn’t be very fast we use the bisection method, which is always reliable. (If you use one of the other root finding algorithms, note that F is strictly monotonic because it is invertable, which might simplify the computation.)

Using the bisection algorithm (or any other) by itself will not solve the original problem because it takes even longer to compute than the function. The method shown below, after its long initial computation, computes the inverse at once for any given y.

Before we get to the method first a description of the two algorithms it depends on.

Bisection algorithm for solving F(x) = 0

The accuracy is determined by the constant tol (tolerance). No upper bound has been placed on the number of iterations. In order to avoid an infinite loop tol must not be smaller than the machine precision used for x.

Two initial numbers must be given that bracket the solution
a  ≤  x  ≤  b

We assume the function is invertable over the given domain which means there is only one root.

The program below and the ones given later are written for PBCC and can be modified for Visual Basic, FORTRAN and any other BASIC-like compiler.

```Function F(x As Single) As Single
Function = x*x - 2   'simple example, the solution in the domain  1 < x < 2  is  Sqr(2)
End Function
'----------------------------

Function PBMain
Local i As Long
Local x, x1, x2 As Single
Local tol As Single

x1 = 1               'bracket solution
x2 = 2
tol = .000001        'fixed
i = 0                'count loops, i not used otherwise

Do
Incr i
x = (x1 + x2) * .5
If F(x) = 0 Or Abs(x1 - x2) < tol Then Exit Do
If Sgn(F(x)) = Sgn(F(x2)) Then
x2 = x
Else
x1 = x
End If
Loop

Print
Print tol; "tolerance"
Print i; "iterations"
Print " Solution: "; x
Print " Actual solution: "; Sqr(2)
WaitKey\$
End Function
```

Function Interpolation

Given evaluations of G(y) at evenly spaced y over an interval, find a good approximation of G(y) for any y in the interval by linear interpolation.

This is useful when G takes a lot of time to compute and great accuracy is unnecessary. You compute G at a finite number of values once and for all, then use this table to approximate G at any value any number of times.

The greater the finite number of these values (n) or the shorter the interval over which they are spread (interval), the more accurate the interpolation.

In order to minimize truncation error we put up with a trivial time penalty by using
... * n / interval

Usually y is considered a function of x, but here the relation is the reverse, because we have in mind an application to the inverse of a function.

In this example the function is G(y) = Sqr(y).  You can change this.

The interpolation will be exact at the y() and less accurate in between. For example, with the fixed parameters in the program as is, where
n = 50,  y0 = 0,  yn = 5,  G(y) = Sqr(y)
consider  y = 4:
Ginterpolate(4)  =  2
G(4)  =  2
With 4 it’s exact, but at  y = 3.95
Ginterpolate(3.95)  =  1.987340
G(3.95)  =  1.987461

```Global y0, yn As Single            'first and last points
Global interval As Single          'length of interval from y0 to yn
Global n As Long                   'number of divisions plus one
Global x(), y() As Single          'x = G(y),  y() equally spaced from y0 to yn
'----------------------------

' The larger n is (see Sub SetupG) the more accurate the interpolation.

Function G(y As Single) As Single  'defined from y0 to yn
Function = Sqr(y)                 'simple example
End Function
'----------------------------

Sub SetupG
n = 50             'fixed

Dim x(n), y(n) As Single
Local i As Long

' chose the first and last points so the line through them on the graph
' does not intersect a zero of the function
y0 = 0             'first point, fixed
yn = 10            'last point, fixed
interval = yn - y0

For i = 0 To n - 1
y(i) = y0 + interval * i / n      ' y0 + s*i  where  s = (yn - y0) / n
x(i) = G(y(i))
Next
y(n) = yn          'possibly more accurate than putting in loop
x(n) = G(yn)
End Sub
'----------------------------

Function Ginterpolate(y As Single) As Single
Local k As Long
k = Int((y - y0)*n / interval)     ' Int((y - y0) / s)  where  s = (yn - y0) / n
Function = x(k) + (y - y(k)) * ( x(k+1) - x(k) ) * n / interval    ' ...*n/interval = .../s
End Function
'----------------------------

Function PBMain
Local y As Single
Local a As String

SetupG

Do
Line Input "Enter a number:  "; a  : If Trim\$(a) = "" Then End
y = Val(a)
If y < y0 Or y > yn Then    'check range
Print "Number must be between" + Str\$(y0) + " and" + Str\$(yn)
Iterate
End If
Print "Interpolated: "; Ginterpolate(y); "     From function: "; G(y)
Print
Loop
End Function
```

Function Inversion

For the purpose of demonstrating the idea described in the introduction the particular function to invert below is fast, simply x squared, and has a known simple solution. You can change the function in Function F, then either change or delete the “Actual inverse” in the output.

```Global y0, yn As Single            'first and last points
Global interval As Single          'length of interval from y0 to yn
Global n As Long                   'number of interpolation points plus one
Global x(), y() As Single          'x = G(y),  y() equally spaced from y0 to yn
Global tol As Single
'----------------------------

Function F(x As Single) As Single
Function = x*x    'simple example for the purpose of demonstration
End Function
'----------------------------

Function G(y As Single) As Single  'the inverse of F
Local x, x1, x2 As Single
Local y1, y2 As Single
x1 = y0
x2 = yn
Do
x = (x1 + x2) * .5
If F(x) = y Or Abs(x1 - x2) < tol Then Exit Do
If Sgn(F(x) - y) = Sgn(F(x2) - y) Then
x2 = x
Else
x1 = x
End If
Loop
Function = x
End Function
'----------------------------

Sub SetupG               'called only once
Dim x(n), y(n) As Global Single
Local i As Long
interval = yn - y0
For i = 0 To n - 1
y(i) = y0 + interval * i / n
x(i) = G(y(i))
Next
y(n) = yn
x(n) = G(yn)
End Sub
'----------------------------

Function Ginterpolate(y As Single) As Single
Local k As Long
k = Int((y - y0)*n / interval)
Function = x(k) + (y - y(k)) * ( x(k+1) - x(k) ) * n / interval
End Function
'----------------------------

Function PBMain
Local y As Single
Local a As String
y0 = 1                                      '<--- change to suit
yn = 10                                     '<--- change to suit
tol = .000001                               'must be no smaller than precision of x
SetupG
Print " Range:  "; Format\$(y0);" to ";Format\$(yn)
Print " Number of interpolation points:  "; Format\$(n)
Print " Tolerance:  ";  Format\$(tol,"0.######")
Print
Do
Line Input " Enter a number:  "; a  : If Trim\$(a) = "" Then End
y = Val(a)
If y < y0 Or y > yn Then                   'check range
Print " Number must be between" + Str\$(y0) + " and" + Str\$(yn)
Iterate
End If
Print " Interpolate:"; Ginterpolate(y)
Print " Our inverse:"; G(y)                '<--- not used in practice
Print " Actual inverse:"; Sqr(y)           '<--- not used in practice
Print
Loop
End Function
```