bugGNU Scientific Library - Bugs: bug #39713, roots/secant.c "derivative...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #39713: roots/secant.c "derivative value is not finite" for a good guess

Submitter:  Max <nikulin>
Submitted:  Wed 07 Aug 2013 03:52:26 AM UTC
   
 
Category:  Runtime error Severity:  3 - Normal
Operating System:  Linux Status:  Fixed
Assigned to:  gladman Open/Closed:  Open
Release:  bzr

Jump to the original submission

Sat 05 Oct 2013 01:14:00 PM UTC, comment #17: 

Thanks for testing the fix on Linux, which allows me to mark it as fixed.

The form df_new = df * (1 - f_new / f) is better than the earlier expression because it avoids the severe loss of precision when (x_new - x) is used in the demominator of the expression for df_new (this was the primary cause of the problems you experienced).

In fact I have just made a further change to:

   df_new = df * ((f - f_new) / f);

which avoids the possible loss of precision in using (1 - f_new / f).

Brian Gladman <gladman>
Group Member
Sat 05 Oct 2013 11:45:00 AM UTC, comment #16: 

I can confirm that my examples work with 4848
revision (2013-10-03 08:58:48) of roots/secant.c.

On the other hand I am in doubt if the form
  df_new = df * (1 - f_new / f);
may lead to a problem with some functions
(slowly converging third order root or something
else).

Max <nikulin>
Sat 05 Oct 2013 10:05:56 AM UTC, comment #15: 

With a small correction to the value of r (= M_PI / M_E), your modified test program (gsl-secant64.c) gives a correct result on Windows x64 using the version of secant.c now in the GSL repository:

# i   x_i-r      x_i - x_p status
  0 -7.509e-003  4.822e-002 -2
  1  0.000e+000  7.509e-003 -2
  2  0.000e+000  0.000e+000 0
# Converged
# f(x_i) =  2.2204e-019

Have you tested this revised version of secant.c?

Brian Gladman <gladman>
Group Member
Sat 05 Oct 2013 08:19:47 AM UTC, comment #14: 

Concerning API changes, in my opinion, the implementation
of the secant method should be rewritten to avoid
requirement of the function derivative, that is actually
computed on the first step only. It is enough to
supply two starting points in the beginning. It may be
hard to determine the derivative or to implement
its calculation.

Finally, a note on Python. Some people do not trust it
for the following reason. Supporting of a long-lasting
project e.g. in scientific research may be a real pain
due to API changes.


Max <nikulin>
Sat 05 Oct 2013 08:08:52 AM UTC, comment #13: 

Sorry for not participating in the discussion
in the middle of the week. Let us return to
the original issue. The key to reproduce it is
the operation system (linux and it was specified
in the original report). Concerning the version,
the last commit to secant.c was in 2007.

The difficulty with reproducing is related to difference
of the default FPU precision on Linux and on Windows
(80 and 64 bit accordingly).

My problem is not connected to the case when
the function evaluates exactly to zero.
It may happen that for some x_0
f(x_0) < 0 and f(x_0 + DBL_EPSILON) > 0
even for a linear function. So any test
for zero would not help.

If the function may be modified as (attached)
double test_f (double x, void *params)
{
   (void) params; // unused
   return -M_E*x + M_PI + 0.001*GSL_DBL_EPSILON;
}
then the problem may be reproduced with -mpc64
compile option, so I expect the same should happen
on Windows as well

$ cc -mpc64 -g -Wall -W   gsl-secant64.c -lgsl -lgslcblas -lm -o gsl-secant64
# i   x_i-r      x_i - x_p status
  0  2.830e-01  4.822e-02 -2
  1  2.905e-01  7.509e-03 -2
gsl: secant.c:104: ERROR: derivative value is not finite
Default GSL error handler invoked.

If the patch from the file file #28779 is applied
$ cc -mpc64 -g -Wall -W   gsl-secant64.c ../gsl-roots/libgsl-roots.a  -lgsl -lgslcblas -lm -o gsl-secant64
$ ./gsl-secant64
# i   x_i-r      x_i - x_p status
  0  2.830e-01  4.822e-02 -2
  1  2.905e-01  7.509e-03 -2
  2  2.905e-01  0.000e+00 0
# Converged
# f(x_i) =  2.2204e-19

Without -mpc the results are essentially the same
(# f(x_i) =  6.7273e-18 last case)

$ gcc --version
gcc (Ubuntu/Linaro 4.7.3-1ubuntu1) 4.7.3
Linux, i686 kernel,
Intel Core 2 Duo processor.

(file #29299)

Max <nikulin>
Thu 03 Oct 2013 09:10:07 AM UTC, comment #12: 

I have fixed this by returning GSL_SUCCESS if the solver input is already at a root value (as suggested below in order to maintain the existing API). 

Brian Gladman <gladman>
Group Member
Wed 02 Oct 2013 05:49:17 PM UTC, comment #11: 

I'm not objecting to a quick fix but I don't wan't to put a lot of effort into making a contribution to GSL if we can't also make improvements.  And I certainly consider that an API that in common situations forces twice as many function evaluations as are necessary to find a root is in need of improvement (in fact I really thought that I must be wrong about this and I still hope that I am).

I am sorry if this is not the right place to discuss such wider issues but I am new here and I wanted to understand the design as in my experience I have found it quite dangerous to 'fix' things that I don't understand.

My experience with having a compatibility define in MPIR (a fork of GMP) is quite the opposite to your own.  But, as you suggest, that discussion is for another place.

Brian Gladman <gladman>
Group Member
Wed 02 Oct 2013 05:03:17 PM UTC, comment #10: 

Agreed on GSL 2.0 permitting API changes.

I'm not aware of an API-breaking #define sitting anywhere.  Having seen these employed extensively along with goofy backward-compatible macros within the HDF5 APIs, I hate them.  An example is http://www.hdfgroup.org/HDF5/doc/RM/RM_H5T.html#Datatype-ArrayCreate.  This is a disgusting, disgusting thing that screws up package management (e.g. is the hdf5 API the new or old one) which percolates into every dependent package.  It is consequently wildly unpleasant to attempt to use HDF5 packages in Debian-ish distros.

AFAIK, GSL 1.x is to be backwards compatible forever, including ABI as much as possible.  Much of this sort of conversation should be had on the developer list so others can correct me where I'm wrong.

What's the objection to having the secant solver shortcircuiting on f == 0.0 and then returning GSL_SUCCESS?  Would that solve the reported issue?  I'm agreeable to discussing fixing or improving the solver API behavior in a feature ticket or on the mailing list, but we're too far afield here for a bug tracker.

Rhys Ulerich <rhysu>
Group administrator
Wed 02 Oct 2013 03:58:09 PM UTC, comment #9: 

Not if we have to maintain the existing API's since we would need s new item in the state for the functiopn value.  And if you are saying that we have to maintain the existing API's for ever more, we cannot expect to improve GSL in areas where experience shows that API improvments are desirable.  I like the Python approach where they are willing to break legacy code every few years in order to accomodate improvements built on real experience (i.e GSL v2.0 would allow API changes).

One way round the legacy issue is to have two returns from the solvers SOLVER_SUCCESS and SOLVER_CONTINUE and to map these both onto GSL_SUCCESS for legacy code and onto GSL_SUCCESS and GSL_CONTINUE otherwise (do we have a define for detecting a user desire to use a legacy API?).

Brian Gladman <gladman>
Group Member
Wed 02 Oct 2013 03:43:14 PM UTC, comment #8: 


> As for duplicated evalution, gsl_root_fsolver_root returns the root without requiring an additional function evaluation. The root doesn't need to be exposed directly to the user through the iteration API as an accessor is available to get it at near-zero cost. Again, one could imagine a different API but this is what's in-place now. You can see this routine exercised in the examples and the source code in fsolver.c.


I apologize and retract all of this.  I mixed up root and function evaluation.  Sorry.  Could one build something similar atop the solver-specific state?

Rhys Ulerich <rhysu>
Group administrator
Wed 02 Oct 2013 03:35:46 PM UTC, comment #7: 

I agree that one could design a more crisp API (like you suggest) that'll better handle conditions like this, but you can't make the solver iteration routine return something other than GSL_SUCCESS without breaking legacy code.  Returning GSL_CONTINUE is nonzero and you'd break existing logic looking for failure from the iteration routine.

To stay within the current API/ABI, shortcircuiting on f == 0.0 and returning GSL_SUCCESS may be the only feasible change.  Anything else that's backwards compatible come to mind?

As for duplicated evalution, gsl_root_fsolver_root returns the root without requiring an additional function evaluation.   The root doesn't need to be exposed directly to the user through the iteration API as an accessor is available to get it at near-zero cost.  Again, one could imagine a different API but this is what's in-place now.  You can see this routine exercised in the examples and the source code in fsolver.c.

Rhys Ulerich <rhysu>
Group administrator
Wed 02 Oct 2013 03:20:10 PM UTC, comment #6: 

In fact the problem arises because there is actually one condition where the solver can detect and defend against a root that has already been found and this is when the previous function value is exactly zero.  And it is the failure of the root solver to detect exactly this condition that is the cause of the problem.

I am coming to the conclusion that solvers should return one of two 'normal' values - GSL_SUCCESS when entered with a root that has already produced a funtion value of zero, and GSL_CONTINUE if they have not detected an end to the iteration in this way.

This would allow a user to react to success immediately or to use a convergence test if continuation is suggested. 

There is also another issue I am thinking about - to use gsl_root_test_residual() for termination the user has to have a function value but function values are not returned to the user.  So they are forced to evaluate the function on each loop and this doubles the number of function evaluations involved in finding the root.

If I have this right, this doesn't seem sensible to me.

Brian Gladman <gladman>
Group Member
Wed 02 Oct 2013 02:50:48 PM UTC, comment #5: 


> So, in my view, the real question raised by this issue is what should the solver do if it is entered at a position where the previous iteration has already found the root and the user has not terminated the iteration?


This feels like a usage error on the user's part.  How one can have secant_iterate defend against having already found the root when "the solver found the root" is a problem- and tolerance-specific statement.  See below, however, ...

> To do nothing risks an infinite loop. But to do what is done now indicates failure for what is really success! Maybe a warning return value such as "root already found'?


Maybe, when the root is obtained identically (i.e. state->f == 0 on entry to secant_iterate in secant.c), have the method immediately return GSL_SUCCESS without attempting to re-evaluate the function.  I say identically because only in that case should all sane tolerancing schemes agree that the root has been found.

Presumably the user not checking for success will either (a) notice an infinite loop in their code or (b) notice that their solve procedure hits the maximum number of iterations if they control for that.

Thanks for digging into this,
Rhys

Rhys Ulerich <rhysu>
Group administrator
Wed 02 Oct 2013 01:36:19 PM UTC, comment #4: 

When the secant method is used to find a root with a linear function, its root prediction is essentially perfect and it will give a function value that is zero (or very close). In this situation, when gsl_root_test_delta() is used to test for termination without also using gsl_root_test_residual(), they will get a return which indicates the need to continue the iteration even though this is not necessary.  And this WILL lead to the "derivative value is not finite" error.

So, in my view, the real question raised by this issue is what should the solver do if it is entered at a position where the previous iteration has already found the root and the user has not terminated the iteration?  To do nothing risks an infinite loop. But to do what is done now indicates failure for what is really success!  Maybe a warning return value such as "root already found'?


 


Brian Gladman <gladman>
Group Member
Wed 02 Oct 2013 10:17:00 AM UTC, comment #3: 

The section within secant.c in which the root position, the function value and the derivative value are updated is:

  x_new = x - (f / df);

  f_new = GSL_FN_FDF_EVAL_F(fdf, x_new) ;
  df_new = (f_new - f) / (x_new - x) ;

  *root = x_new ;

Although it involves an additional multiplication, I think this may be safer:

  x_new = x - (f / df);

  f_new = GSL_FN_FDF_EVAL_F(fdf, x_new) ;
  df_new = df * (1.0 - f_new / f) ;

  *root = x_new;

since it avoids the potential problem of x_new and x being the same because (f / df) is too small to change x when it is subtracted.  Also, since the algorithm shouldn't presumably be entered with f = 0, the division f_new / f should not cause a divide by zero problem.

Brian Gladman <gladman>
Group Member
Tue 01 Oct 2013 07:29:54 PM UTC, comment #2: 

Running the program provided on Windows x64 gives:

# i   x_i-r      x_i - x_p status
  0 -1.110e-016 -2.347e-001 -2
  1  0.000e+000  1.110e-016 0
# Converged
# f(x_i) =  0.0000e+000

so the problem appears to be system dependent. 

On which systems does this error show up? 

Can anyone else confirm this as a problem?



All bug reports should include:

       The version number of GSL
       The hardware and operating system
       The compiler used, including version number and
compilation options
       A description of the bug behaviour
       A short program which reproducibly exercises the bug

It is also useful if you can report whether the same problem occurs
when the library is compiled without optimization.

Thank you.

Brian Gladman <gladman>
Group Member
Wed 07 Aug 2013 05:30:38 AM UTC, comment #1: 

My bad. Of course, the relative change of root approximation
matters, not the absolute one. The updated patch is attached.

(file #28783)

Max <nikulin>
Wed 07 Aug 2013 03:52:26 AM UTC, original submission:  

Implementation of the secant method for root finding fails
for an (almost) linear function. If at a certain iteration
the root approximation is within DBL_EPSILON, then
the next iteration fails with "derivative value is not finite"
error. Actually it means that a good approximation has been found.

It may be tested with the following function
(full program source is attached as well)

double test_f (double x, void *params)
{
   return -M_PI*x + M_E;
}
double test_df (double x, void *params)
{
   return -M_PI;
}

# i   x_i-r      x_i - x_p status
  0  0.000e+00 -2.347e-01 -2
gsl: secant.c:104: ERROR: derivative value is not finite
Default GSL error handler invoked.

I think that an addition test should be added and
derivative should be updated only if the step is
significantly greater then DBL_EPSILON. Otherwise
rounding error or even zero division may be the issue.

The proposed patch is attached. When it is applied
the programs works correctly.
# i   x_i-r      x_i - x_p status
  0  0.000e+00 -2.347e-01 -2
  1  0.000e+00  0.000e+00 0
# Converged
# f(x_i) =  2.9707e-17

It seems that all 82 tests supplied in roots directory
are passed correctly with this patch.


Max <nikulin>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #29299:  gsl-secant64.c added by nikulin (2KiB - text/x-csrc - test case for 64 bit FPU mode (Windows))
file #28783:  gsl-roots-secant-v2.patch added by nikulin (803B - text/x-patch - corrected patch)
file #28779:  gsl-secant.c added by nikulin (2KiB - text/x-csrc - A program to demonstrate the issue and the patch)
file #28780:  gsl-roots-secant.patch added by nikulin (795B - text/x-patch - A program to demonstrate the issue and the patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rhysu (Posted a comment)
  • -email is unavailable- added by gladman (Posted a comment)
  • -email is unavailable- added by nikulin (Submitted the item)
  •  

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-10-05 gladman StatusReady For Test Fixed
    2013-10-05 nikulin Attached File- Added gsl-secant64.c, #29299
    2013-10-03 gladman Assigned toNone gladman
    2013-10-03 gladman StatusNone Ready For Test
    2013-08-07 nikulin Attached File- Added gsl-roots-secant-v2.patch, #28783
    2013-08-07 nikulin Attached File- Added gsl-secant.c, #28779
        Attached File- Added gsl-roots-secant.patch, #28780

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code