Reverse complex 2D lookup table

I have some function which maps some input to the output . The output is a complex number. What I'm actually interested in is the inverse function . But since this inversion can't be done in a analytical way I need to do it with a numerical approximation.

Since is computationally expensive my idea was to use a lookup table approach. I can generate a 2D lookup table with the dimensions (forward lookup table), but what I actually need is the inverse of this lookup table—yielding based on a given .

For the inversion of the lookup table the simplest approach I can think of is using the entries of the forward lookup table as vertices and interpolate between them in a regular grid yielding the inverse lookup table. In case the inverse lookup table would be to large for the required precision, I would generate coarse table and use the values as starting values for an optimization algorithm. Are there any easier approaches I'm missing?

Where , are constant and , are .

for f(x,y) you can use f(x,y)->(a,b) 2D LUT (Look Up Table)

But the stored grid points must be selected so dense so that there is max one bump per grid rectangle, otherwise interpolation would not work correctly.

If you want to use linear interpolation then the local min/max must be points inside LUT, for polynomial interpolation of higher order is this not always needed. I would use 4 point cubic interpolation
how to compute g(a,b)->(x,y)
- first is reverse mapping possible?
- so how many (x,y) points return the same (a,b)=f(x,y) ?
- it is the same as question: Is f function or not?
if f is not function then you got a problem and you can not solve this, unless somehow subdivide the range to sub ranges where f is function and then you will have to select the proper range according to some rules dependent on what are you trying to do. So let assume that f is function

So how to compute (x,y)=g(a,b) that f(x,y)=(a,b)?
1. I would start with approximation of the result. So try enough (x,y) values along the whole range and store the closest point to the desired output so that |f(x,y)-(a,b)| is minimal.
2. Then start this again but not on full range but around this point instead
  - recursively increasing accuracy up to a point you need.
  - look here Increasing accuracy of solution of transcendental equation there I am computing similar problem have 1D LUT of 2D points (a(t),y(t)) and need the inverse 3D point (a0,y0,z0) You can use my approximation class there

Nesting of approximations is done like this:

int n=5;  // recursions
double e; // Error Of Solution Variable
approx ax,ay;
//            min    max   step   
for (ax.init(-100.0,+100.0,10.0,n,&e);!ax.done;ax.step())
for (ay.init(-100.0,+100.0,10.0,n,&e);!ay.done;ay.step())
    {
    e=|f(ax.a,ay.a)-(a,b)|;
    }
// here (ax.a,ay.a) should hold your solution for input point `(a,b)`

initial step should be as small so there is no bump missed
if your g(a,b) shape is too complex then this will probably not work properly

From this you can compute the inverse LUT table ...

Each recursion divides the step by 10 so choose n wisely.

For 2D and singular point is the performance of this not that bad O((log(N))^2). I am doing this on 3D O((log(N))^3) with 100 points per e computation and that is painfully slow (about 35 seconds)