Question about best practices with numerical functions

Question about best practices with numerical functions

[mockturtle]

Dear.all, 
I have a question about the best way to manage a potential loss of precision in a numerical function.  This is a doubt that came to my mind while writing a piece of software; now I solved the specific problem, but the curiosity remains.

Let me explain.

Recently I needed to write an implementation of the Lambert W function (is the function that given y finds x such that x*exp(x)=y).  This function cannot be expressed with elementary functions and the algorithm I found basically solves the equation in an iterative way.  Of course, if you fix the maximum number of iterations, it can happen that the convergence is not fast enough and you obtain a result that is potentially less precise than what you would expect.  

I was wondering how to manage such a non convergence case. Please note that I am supposing that I am writing a "general" function that could be used in many different programs.  If the function was specific for a single program, then I would choose the line of action (e.g., ignore, log a warning or raise an exception) depending on the needs of the specific program.

(Incidentally, it turned out that the implementation converges nicely for any value of interest; nevertheless, the curiosity remains...)

I can see few line of actions that would make sense

   [1] Raise an exception.  

Maybe this is a bit too drastic since there are cases where a moderate loss of precision does not matter (this was my case, i just needed one or two decimal digits)     

   [2] Let the function have an optional "precision" parameter and raise an exception if the precision goes below that

   [3] Let the function return a record with a field Value with the actual result and a field Error with the estimated precision.

This would make the code a bit heavier since instead of calling

X := Lambert(Y);

you would say

X := Lambert(Y).Value;

Not really a huge deal, however...

   [4] Print a warning message to standard error or some logging system and go on.

This sounds like the worst option to me.  The message could be overlooked and, moreover, it supposes there is some logging facilities or that  the standard error is available for logging... Remember that the function should be general, to be used in any program.

   [5] ??? 

Any suggestions?

Thank you in advance

Riccardo

Re: Question about best practices with numerical functions

[Dmitry A. Kazakov]

On 04/07/2020 07:30, mockturtle wrote:

> I have a question about the best way to manage a potential loss of precision in a numerical function.  This is a doubt that came to my mind while writing a piece of software; now I solved the specific problem, but the curiosity remains.

[...]

>     [5] ???
> 
> Any suggestions?

[5] Interval computations is the best way to handle rounding errors:

    https://en.wikipedia.org/wiki/Interval_arithmetic

An Ada implementation is here:

    http://www.dmitry-kazakov.de/ada/intervals.htm

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Re: Question about best practices with numerical functions

[Nasser M. Abbasi]

On 7/4/2020 12:30 AM, mockturtle wrote:
> Dear.all,
> I have a question about the best way to manage a potential loss of precision in a numerical function.  This is a doubt that came to my mind while writing a piece of software; now I solved the specific problem, but the curiosity remains.
> 
> Let me explain.
> 
> Recently I needed to write an implementation of the Lambert W function (is the function that given y finds x such that x*exp(x)=y).  This function cannot be expressed with elementary functions and the algorithm I found basically solves the equation in an iterative way.  Of course, if you fix the maximum number of iterations, it can happen that the convergence is not fast enough and you obtain a result that is potentially less precise than what you would expect.
> 
> I was wondering how to manage such a non convergence case. Please note that I am supposing that I am writing a "general" function that could be used in many different programs.  If the function was specific for a single program, then I would choose the line of action (e.g., ignore, log a warning or raise an exception) depending on the needs of the specific program.
> 
> (Incidentally, it turned out that the implementation converges nicely for any value of interest; nevertheless, the curiosity remains...)
> 
> I can see few line of actions that would make sense
> 
>     [1] Raise an exception.
> 
> Maybe this is a bit too drastic since there are cases where a moderate loss of precision does not matter (this was my case, i just needed one or two decimal digits)
> 
>     [2] Let the function have an optional "precision" parameter and raise an exception if the precision goes below that
> 
>     [3] Let the function return a record with a field Value with the actual result and a field Error with the estimated precision.
> 
> This would make the code a bit heavier since instead of calling
> 
> X := Lambert(Y);
> 
> you would say
> 
> X := Lambert(Y).Value;
> 
> Not really a huge deal, however...
> 
>     [4] Print a warning message to standard error or some logging system and go on.
> 
> This sounds like the worst option to me.  The message could be overlooked and, moreover, it supposes there is some logging facilities or that  the standard error is available for logging... Remember that the function should be general, to be used in any program.
> 
>     [5] ???
> 
> Any suggestions?
> 
> Thank you in advance
> 
> Riccardo
> 


Most Fortran Lapack use INFO code.

"
All documented routines have a diagnostic argument INFO that
indicates the success or failure of the computation, as follows:

INFO = 0: successful termination
INFO < 0: illegal value of one or more arguments -- no computation performed
INFO > 0: failure in the course of computation
"

https://www.netlib.org/lapack/lug/node138.html

So you could follow that.

Another option is to throw an exception. So your program will look like

     try
       X= Lambert(Y,....);
       ... if you get here, then no error
    catch:
       ..... error happend. handle it
    end try

--Nasser
    

Re: Question about best practices with numerical functions

[Shark8]

On Friday, July 3, 2020 at 11:30:53 PM UTC-6, mockturtle wrote:
> 
>    [5] ??? 
> 
> Any suggestions?
One way to do this would be to use fixed-point types:
(1) Convert your input to the fixed-point that has a "good enough" delta for the precision you want.
(2) Run the algorithm.
(3) Convert back to your normal value-type.
This assumes you're using floating-point or integers, but one nice thing about fixed-point is that it has a bounded error when dealing with operations, unlike floating-point. -- I remember some years ago seeing a bug report dealing with floating-point, where the particular error simply couldn't have happened with fixed-point.

Re: Question about best practices with numerical functions

[Dmitry A. Kazakov]

On 06/07/2020 19:02, Shark8 wrote:

> One way to do this would be to use fixed-point types:
> (1) Convert your input to the fixed-point that has a "good enough" delta for the precision you want.
> (2) Run the algorithm.
> (3) Convert back to your normal value-type.
> This assumes you're using floating-point or integers, but one nice thing about fixed-point is that it has a bounded error when dealing with operations, unlike floating-point. -- I remember some years ago seeing a bug report dealing with floating-point, where the particular error simply couldn't have happened with fixed-point.

Rounding error is bounded in both cases. Fixed-point has same error 
regardless the values involved in the operations. Floating-point has 
error depending on the values.

I would say that floating-point error would be roughly same for 
addition, subtraction and multiplication, provided fixed-point does not 
overflow. It will be hugely better for division.

Using fixed-point arithmetic has only sense for a few marginal cases of 
rounding.

Furthermore converting many algorithms to fixed-point might turn quite 
non-trivial as you will have to ensure absence of overflows and 
underflows. Where floating-point computation just would lose some 
precision, fixed-point will catastrophically fail.

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Re: Question about best practices with numerical functions

[mockturtle]

Also in my case the loss of precision is not due to approximation introduced by  operations, but to a loss of convergence because (maybe) for that input value the convergence is too slow and the bound on the number of iterations is not large enough.

On Tuesday, July 7, 2020 at 12:23:09 AM UTC+2, Dmitry A. Kazakov wrote:
> On 06/07/2020 19:02, Shark8 wrote:
> 
> > One way to do this would be to use fixed-point types:
> > (1) Convert your input to the fixed-point that has a "good enough" delta for the precision you want.
> > (2) Run the algorithm.
> > (3) Convert back to your normal value-type.
> > This assumes you're using floating-point or integers, but one nice thing about fixed-point is that it has a bounded error when dealing with operations, unlike floating-point. -- I remember some years ago seeing a bug report dealing with floating-point, where the particular error simply couldn't have happened with fixed-point.
> 
> Rounding error is bounded in both cases. Fixed-point has same error 
> regardless the values involved in the operations. Floating-point has 
> error depending on the values.
> 
> I would say that floating-point error would be roughly same for 
> addition, subtraction and multiplication, provided fixed-point does not 
> overflow. It will be hugely better for division.
> 
> Using fixed-point arithmetic has only sense for a few marginal cases of 
> rounding.
> 
> Furthermore converting many algorithms to fixed-point might turn quite 
> non-trivial as you will have to ensure absence of overflows and 
> underflows. Where floating-point computation just would lose some 
> precision, fixed-point will catastrophically fail.
> 
> -- 
> Regards,
> Dmitry A. Kazakov
> http://www.dmitry-kazakov.de