GNU 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: |
Brian Gladman <gladman>![]() ![]() |
Sat 05 Oct 2013 11:45:00 AM UTC, comment #16: I can confirm that my examples work with 4848
On the other hand I am in doubt if the form
|
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
Have you tested this revised version of secant.c? |
Brian Gladman <gladman>![]() ![]() |
Sat 05 Oct 2013 08:19:47 AM UTC, comment #14: Concerning API changes, in my opinion, the implementation
Finally, a note on Python. Some people do not trust it
|
Max <nikulin> |
Sat 05 Oct 2013 08:08:52 AM UTC, comment #13: Sorry for not participating in the discussion
The difficulty with reproducing is related to difference
My problem is not connected to the case when
If the function may be modified as (attached)
$ cc -mpc64 -g -Wall -W gsl-secant64.c -lgsl -lgslcblas -lm -o gsl-secant64
If the patch from the file file #28779 is applied
Without -mpc the results are essentially the same
$ gcc --version
|
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>![]() ![]() |
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>![]() ![]() |
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>![]() |
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>![]() ![]() |
Wed 02 Oct 2013 03:43:14 PM UTC, comment #8:
|
Rhys Ulerich <rhysu>![]() |
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>![]() |
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>![]() ![]() |
Wed 02 Oct 2013 02:50:48 PM UTC, comment #5:
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 Ulerich <rhysu>![]() |
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>![]() ![]() |
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) ;
*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) ;
*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>![]() ![]() |
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
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
It is also useful if you can report whether the same problem occurs
Thank you.
|
Brian Gladman <gladman>![]() ![]() |
Wed 07 Aug 2013 05:30:38 AM UTC, comment #1: My bad. Of course, the relative change of root approximation
|
Max <nikulin> |
Wed 07 Aug 2013 03:52:26 AM UTC, original submission:
Implementation of the secant method for root finding fails
It may be tested with the following function
double test_f (double x, void *params)
# i x_i-r x_i - x_p status
I think that an addition test should be added and
The proposed patch is attached. When it is applied
It seems that all 82 tests supplied in roots directory
|
Max <nikulin> |
Depends on the following items: None found
Items that depend on this one: None found
Follow 7 latest changes.
Date | Changed by | Updated Field | Previous Value | => | Replaced by |
---|---|---|---|---|---|
2013-10-05 | gladman | Status | Ready For Test | ![]() |
Fixed |
2013-10-05 | nikulin | Attached File | - | ![]() |
Added gsl-secant64.c, #29299 |
2013-10-03 | gladman | Assigned to | None | ![]() |
gladman |
2013-10-03 | gladman | Status | None | ![]() |
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 |
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).