Reverse function numerically

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Reverse function numerically

stn021
Hello,

I have 2 functions:

y1 = f1( x1,x2 )
y2 = f2( x1,x2 )

also there is the reverse function

x1 = g1( y1,y2 )
x2 = g2( y1,y2 )

in octave syntax:
[ y1,y2 ] = f( [ x1,x2 ] )
[ x1,x2 ] = g( [ y1,y2 ] )

Both functions lead to exactly one distinct pair of results for each
pair of input variables. That means that always

f ( g ( [x1,x2] ) ) == [ x1,x2 ]


My problem is this: I can calculate f(x1,x2) but I cannot calculate g(y1,y2).
Meaning that f( [x1,x2] ) cannot be algebraically reversed.

I am looking for a way to calculate g( [y1,y2] ).
The obvious solution would be some kind of approximation.
(fft looks like a good choice)

So far I could not piece together how to do that.
Could you please give me a hint ?

THX,stn


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Re: Reverse function numerically

stn021
Hello,

this is the same question I asked in my previous mail about an hour ago.
I have added a code-example to illustrate what I mean.

See at the bottom of this email.
This is only meant as an example. It can be solved algebraically, so
please simply assume that it cannot be solved,


My question is this:
I have 2 functions:

y1 = f1( x1,x2 )
y2 = f2( x1,x2 )

also there is the reverse function

x1 = g1( y1,y2 )
x2 = g2( y1,y2 )

in octave syntax:
[ y1,y2 ] = f( [ x1,x2 ] )
[ x1,x2 ] = g( [ y1,y2 ] )


Both functions lead to exactly one distinct pair of results for each
pair of input variables.
That means that within predefined limits ( im my example all variables
are >=0 and <=1 )

this is true :   f ( g ( [x1,x2] ) ) == [ x1,x2 ]
and this also:   ( f(a,b) == f(c,d) )  <=>  ( a==c and b==d )

I am not a mathmatician so I hope I got this one right :-)


My problem is this: I can calculate f(x1,x2) but I cannot calculate g(y1,y2).
Meaning that f( [x1,x2] ) cannot be algebraically reversed.

I am looking for a way to calculate g( [y1,y2] ).
The obvious solution would be some kind of approximation.
(fft looks like a good choice)

So far I could not piece together how to do that.
Could you please give me a hint ?

THX,stn


code-example:

# function [x1 x2] = g(y1,y2)
#    ...unclear...
# end

function [y1 y2] = f( x1,x2 )
  y1 = x1.^2 .* (2-x2) / 2 ;
  y2 = (2-x1) .* x2.^2 / 2 ;
end


x = 0:.025:1 ;
[ x1 x2 ] = meshgrid( x,x ) ;
[ y1 y2 ] = f( x1 , x2 ) ;

plot3( x1,x2,y1,".g" ) ; hold on ;
plot3( x1,x2,y2,".b" ) ;
xlabel( "x1" ) ; ylabel( "x2" ) , zlabel( "y1 y2" ) ;


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Re: Reverse function numerically

sdorsher
A classic example of this problem is x=e^x. This can be generalized to vectors. This is called a transcendental equation, because it has no algebraic solution. These can be solved by picking an initial guess and iteratively improving the guess by feeding it back into the right hand side of the function. In this case it works better if the equation is written in a different form. It may take some manipulation. The technique is called the relaxation method, and can be found in many introductory numerical methods textbooks. I haven’t looked this up before summarizing, so you might want to. Best of luck!

There are other methods for solving transcendental equations, but I can’t remember their names or details without looking them up. Please ask again if you need my help with this.

S. Dorsher

On Jan 28, 2018, at 5:35 PM, stn021 <[hidden email]> wrote:

Hello,

this is the same question I asked in my previous mail about an hour ago.
I have added a code-example to illustrate what I mean.

See at the bottom of this email.
This is only meant as an example. It can be solved algebraically, so
please simply assume that it cannot be solved,


My question is this:
I have 2 functions:

y1 = f1( x1,x2 )
y2 = f2( x1,x2 )

also there is the reverse function

x1 = g1( y1,y2 )
x2 = g2( y1,y2 )

in octave syntax:
[ y1,y2 ] = f( [ x1,x2 ] )
[ x1,x2 ] = g( [ y1,y2 ] )


Both functions lead to exactly one distinct pair of results for each
pair of input variables.
That means that within predefined limits ( im my example all variables
are >=0 and <=1 )

this is true :   f ( g ( [x1,x2] ) ) == [ x1,x2 ]
and this also:   ( f(a,b) == f(c,d) )  <=>  ( a==c and b==d )

I am not a mathmatician so I hope I got this one right :-)


My problem is this: I can calculate f(x1,x2) but I cannot calculate g(y1,y2).
Meaning that f( [x1,x2] ) cannot be algebraically reversed.

I am looking for a way to calculate g( [y1,y2] ).
The obvious solution would be some kind of approximation.
(fft looks like a good choice)

So far I could not piece together how to do that.
Could you please give me a hint ?

THX,stn


code-example:

# function [x1 x2] = g(y1,y2)
#    ...unclear...
# end

function [y1 y2] = f( x1,x2 )
 y1 = x1.^2 .* (2-x2) / 2 ;
 y2 = (2-x1) .* x2.^2 / 2 ;
end


x = 0:.025:1 ;
[ x1 x2 ] = meshgrid( x,x ) ;
[ y1 y2 ] = f( x1 , x2 ) ;

plot3( x1,x2,y1,".g" ) ; hold on ;
plot3( x1,x2,y2,".b" ) ;
xlabel( "x1" ) ; ylabel( "x2" ) , zlabel( "y1 y2" ) ;


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------



-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Re: Reverse function numerically

stn021
Hello Steven,

 thank you for the idea. That ist exactly what I have been doing so far.

First I precalculated a few thousand pairs of values and then used Interpolation to get the starting value. 

Then I used leasqr() to do the rest.

That works fine but it is horribly slow.
150ms per value pair.
Now I am hoping to find a much faster way.

THX, stn




Am 29.01.2018 01:09 schrieb "Steven Dorsher" <[hidden email]>:
A classic example of this problem is x=e^x. This can be generalized to vectors. This is called a transcendental equation, because it has no algebraic solution. These can be solved by picking an initial guess and iteratively improving the guess by feeding it back into the right hand side of the function. In this case it works better if the equation is written in a different form. It may take some manipulation. The technique is called the relaxation method, and can be found in many introductory numerical methods textbooks. I haven’t looked this up before summarizing, so you might want to. Best of luck!

There are other methods for solving transcendental equations, but I can’t remember their names or details without looking them up. Please ask again if you need my help with this.

S. Dorsher

On Jan 28, 2018, at 5:35 PM, stn021 <[hidden email]> wrote:

Hello,

this is the same question I asked in my previous mail about an hour ago.
I have added a code-example to illustrate what I mean.

See at the bottom of this email.
This is only meant as an example. It can be solved algebraically, so
please simply assume that it cannot be solved,


My question is this:
I have 2 functions:

y1 = f1( x1,x2 )
y2 = f2( x1,x2 )

also there is the reverse function

x1 = g1( y1,y2 )
x2 = g2( y1,y2 )

in octave syntax:
[ y1,y2 ] = f( [ x1,x2 ] )
[ x1,x2 ] = g( [ y1,y2 ] )


Both functions lead to exactly one distinct pair of results for each
pair of input variables.
That means that within predefined limits ( im my example all variables
are >=0 and <=1 )

this is true :   f ( g ( [x1,x2] ) ) == [ x1,x2 ]
and this also:   ( f(a,b) == f(c,d) )  <=>  ( a==c and b==d )

I am not a mathmatician so I hope I got this one right :-)


My problem is this: I can calculate f(x1,x2) but I cannot calculate g(y1,y2).
Meaning that f( [x1,x2] ) cannot be algebraically reversed.

I am looking for a way to calculate g( [y1,y2] ).
The obvious solution would be some kind of approximation.
(fft looks like a good choice)

So far I could not piece together how to do that.
Could you please give me a hint ?

THX,stn


code-example:

# function [x1 x2] = g(y1,y2)
#    ...unclear...
# end

function [y1 y2] = f( x1,x2 )
 y1 = x1.^2 .* (2-x2) / 2 ;
 y2 = (2-x1) .* x2.^2 / 2 ;
end


x = 0:.025:1 ;
[ x1 x2 ] = meshgrid( x,x ) ;
[ y1 y2 ] = f( x1 , x2 ) ;

plot3( x1,x2,y1,".g" ) ; hold on ;
plot3( x1,x2,y2,".b" ) ;
xlabel( "x1" ) ; ylabel( "x2" ) , zlabel( "y1 y2" ) ;


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------




-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Re: Reverse function numerically

Montgomery-Smith, Stephen
Could you use fsolve?
From: Help-octave <help-octave-bounces+stephen=[hidden email]> on behalf of stn021 <[hidden email]>
Sent: Sunday, January 28, 2018 6:24:15 PM
To: Steven Dorsher
Cc: stn021; Help GNU Octave
Subject: Re: Reverse function numerically
 
Hello Steven,

 thank you for the idea. That ist exactly what I have been doing so far.

First I precalculated a few thousand pairs of values and then used Interpolation to get the starting value. 

Then I used leasqr() to do the rest.

That works fine but it is horribly slow.
150ms per value pair.
Now I am hoping to find a much faster way.

THX, stn




Am 29.01.2018 01:09 schrieb "Steven Dorsher" <[hidden email]>:
A classic example of this problem is x=e^x. This can be generalized to vectors. This is called a transcendental equation, because it has no algebraic solution. These can be solved by picking an initial guess and iteratively improving the guess by feeding it back into the right hand side of the function. In this case it works better if the equation is written in a different form. It may take some manipulation. The technique is called the relaxation method, and can be found in many introductory numerical methods textbooks. I haven’t looked this up before summarizing, so you might want to. Best of luck!

There are other methods for solving transcendental equations, but I can’t remember their names or details without looking them up. Please ask again if you need my help with this.

S. Dorsher

On Jan 28, 2018, at 5:35 PM, stn021 <[hidden email]> wrote:

Hello,

this is the same question I asked in my previous mail about an hour ago.
I have added a code-example to illustrate what I mean.

See at the bottom of this email.
This is only meant as an example. It can be solved algebraically, so
please simply assume that it cannot be solved,


My question is this:
I have 2 functions:

y1 = f1( x1,x2 )
y2 = f2( x1,x2 )

also there is the reverse function

x1 = g1( y1,y2 )
x2 = g2( y1,y2 )

in octave syntax:
[ y1,y2 ] = f( [ x1,x2 ] )
[ x1,x2 ] = g( [ y1,y2 ] )


Both functions lead to exactly one distinct pair of results for each
pair of input variables.
That means that within predefined limits ( im my example all variables
are >=0 and <=1 )

this is true :   f ( g ( [x1,x2] ) ) == [ x1,x2 ]
and this also:   ( f(a,b) == f(c,d) )  <=>  ( a==c and b==d )

I am not a mathmatician so I hope I got this one right :-)


My problem is this: I can calculate f(x1,x2) but I cannot calculate g(y1,y2).
Meaning that f( [x1,x2] ) cannot be algebraically reversed.

I am looking for a way to calculate g( [y1,y2] ).
The obvious solution would be some kind of approximation.
(fft looks like a good choice)

So far I could not piece together how to do that.
Could you please give me a hint ?

THX,stn


code-example:

# function [x1 x2] = g(y1,y2)
#    ...unclear...
# end

function [y1 y2] = f( x1,x2 )
 y1 = x1.^2 .* (2-x2) / 2 ;
 y2 = (2-x1) .* x2.^2 / 2 ;
end


x = 0:.025:1 ;
[ x1 x2 ] = meshgrid( x,x ) ;
[ y1 y2 ] = f( x1 , x2 ) ;

plot3( x1,x2,y1,".g" ) ; hold on ;
plot3( x1,x2,y2,".b" ) ;
xlabel( "x1" ) ; ylabel( "x2" ) , zlabel( "y1 y2" ) ;


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------




-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Fwd: Reverse function numerically

Juan Pablo Carbajal-2
As said before, you are essentially searching for zeros of the
function g(x1,x2) - f(y1,y2), i.e.

g(x1,x2) - f(y1,y2) = 0

You can solve this with with any zero finding routine: Gauss-Seidel,
Picard, Newton-Rapson, etc... or you can approximate it with a
minimizing problem, e.g.

[x1,x2] = argmin ||g(x1,x2) - f(y1,y2)||

where ||.|| is some distance (or squared norm). The iterative methods
can be fast or not depending on the structure of your functions, the
minimization method will in general be slow.
To sped up you can resource to interpolation, as you did already,
solve the problem for a bunch of relevant points, then interpolate.
The solving and the creation of the interpolant is done off-line and
only once, To use your inverse, you just evaluate your interpolant
(usually very cheap). The interpolation can be done with Polynomials
(since it is a 2D domain), with Special functions or generalized
Fourier series (truncated, of course), or with Gaussian processes
(kriging or co-kriging).

There are packages in octave to do any f these thigns.

The particular choice depends on your function, is the example you
send the actual function you want to invert?


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Re: Reverse function numerically

stn021
Hello all,

thanks for the ideas.

I have been busy so it took me a while to get back to you.

So far I have checked out the hint of Stephen Montgomery-Smith, which
was fsolve().

That essentially does the same as I have previously done with
leasqr(), only a bit more efficient.

With leasqr() it takes about 150ms to calculate one single value. With
fsolve() it gets done in 90ms.
That is an improvement of course but I had hoped to find something
more significant.
For comparison: the "forward"-function takes about 5 microseconds per value.

@Juan Pablo Carbajal
The function in my example is not the actual function but it has a
very similar shape.
Specifically the two intersecting (almost-)planes are similar and the
limit of all variables in the range of 0..1 is also given.
So I believe that any solutions for the example would work on the
original functions as well
The actual modeling-function involves lots of sums of nonlinear
functions and would IMO not help at this point.

I continue work on the matter ...

THX for now, Stefan


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Re: Reverse function numerically

Juan Pablo Carbajal-2
The modeling function is what you want to find.

If f(x1,x2) is continuously differentiable (as your example) then
split the domain where det(J f) == 0 and obtain local inverses using
root finding (this exploits the existence of the inverse function
based on the implicit function theorem).

For each domain, make a unstructured triangular mesh and use 2D
polynomial interpolation. However since your inverse looks like a
rational function, maybe splines are a better option.


-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------
Reply | Threaded
Open this post in threaded view
|

Re: Reverse function numerically

Montgomery-Smith, Stephen
In reply to this post by stn021
On 01/28/18 17:35, stn021 wrote:

> Hello,
>
> this is the same question I asked in my previous mail about an hour ago.
> I have added a code-example to illustrate what I mean.
>
> See at the bottom of this email.
> This is only meant as an example. It can be solved algebraically, so
> please simply assume that it cannot be solved,
>
>
> My question is this:
> I have 2 functions:
>
> y1 = f1( x1,x2 )
> y2 = f2( x1,x2 )
>
> also there is the reverse function
>
> x1 = g1( y1,y2 )
> x2 = g2( y1,y2 )
>
> in octave syntax:
> [ y1,y2 ] = f( [ x1,x2 ] )
> [ x1,x2 ] = g( [ y1,y2 ] )
>
>
> Both functions lead to exactly one distinct pair of results for each
> pair of input variables.
> That means that within predefined limits ( im my example all variables
> are >=0 and <=1 )
>
> this is true :   f ( g ( [x1,x2] ) ) == [ x1,x2 ]
> and this also:   ( f(a,b) == f(c,d) )  <=>  ( a==c and b==d )
>
> I am not a mathmatician so I hope I got this one right :-)
>
>
> My problem is this: I can calculate f(x1,x2) but I cannot calculate g(y1,y2).
> Meaning that f( [x1,x2] ) cannot be algebraically reversed.
>
> I am looking for a way to calculate g( [y1,y2] ).
> The obvious solution would be some kind of approximation.
> (fft looks like a good choice)
>
> So far I could not piece together how to do that.
> Could you please give me a hint ?
>
> THX,stn
>
>
> code-example:
>
> # function [x1 x2] = g(y1,y2)
> #    ...unclear...
> # end
>
> function [y1 y2] = f( x1,x2 )
>    y1 = x1.^2 .* (2-x2) / 2 ;
>    y2 = (2-x1) .* x2.^2 / 2 ;
> end
>
>
> x = 0:.025:1 ;
> [ x1 x2 ] = meshgrid( x,x ) ;
> [ y1 y2 ] = f( x1 , x2 ) ;
>
> plot3( x1,x2,y1,".g" ) ; hold on ;
> plot3( x1,x2,y2,".b" ) ;
> xlabel( "x1" ) ; ylabel( "x2" ) , zlabel( "y1 y2" ) ;


You can put the problem into Wolfram alpha:

https://www.wolframalpha.com/input/?i=Solve%5B%7By1%3D%3D1%2F2x1%5E2(2-x2),y2%3D%3D1%2F2(2-x1)x2%5E2%7D,%7Bx1,x2%7D%5D

then you will see it involves solving a quintic.  As far as I know, the
best way to solve a quintic is using a numerical method, like fsolve.
So I think you are stuck using fsolve, or some other numerical method.

Another way I have seen is to use Neural Nets to find the coefficients
of a polynomial or rational function approximation to the inverse.  That
would end up being rather fast once you have computed the coefficients.
I personally have no idea how to use Neural Nets, so I couldn't help you
any more than that.  You could use some other optimization method to
compute the coefficients.

You would want to graph the inverse function to get some idea of what
the form of the polynomial or rational function should be, before you
start fitting coefficients.



-----------------------------------------
Join us March 12-15 at CERN near Geneva
Switzerland for OctConf 2018.  More info:
https://wiki.octave.org/OctConf_2018
-----------------------------------------