Hi,
with my students, this afternoon, I encountered the following annoying behaviour of Scipy: >>> from scipy import optimize >>> optimize.newton(lambda x: x**3-x+1, 0) 0.999999989980082 This is very weird, because the equation is rather simple with simple coefficients, and the initial point is simple also. Any teacher or student could fall into this case when trying to use optimize.newton. Of course, the answer is absolutely wrong. First, some facts: * using the function with the derivative gives the correct answer (we are not in some tricky case); * my students wrote pieces of code with Newton-Raphson algorithm returning the expected answer even without the derivative; they guessed the value of the derivative either by some (f(x+eps)-f(x-eps))/(2*eps) or (f(x+eps)-f(x))/eps, etc. All these methods worked well. I had a glance in the source code and found that: * the function is actually thinking the answer is correct (the answer has nothing to do with iterating too many times) because two consecutive values of x are close "enough". * the answer is returned after only 3 iterations (correct answer requires about 20 iterations). I think it has to do with the way the derivative is computed by Scipy (in order to call the f function only once in each iteration, the secant is computed by using f(x) for consecutive values of x (rather than by using some epsilon); furthermore, in order to start this method the function has to provide an arbitrary "second" value of x). Is it a bug? Regards. -- Thomas Baruchel _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
On 24 February 2016 at 19:43, Thomas Baruchel <[hidden email]> wrote:
> with my students, this afternoon, I encountered the following annoying > behaviour of Scipy: > >>>> from scipy import optimize >>>> optimize.newton(lambda x: x**3-x+1, 0) > > 0.999999989980082 > > This is very weird, because the equation is rather simple with simple > coefficients, > and the initial point is simple also. Somehow the initial guess is problematic: In [1]: from scipy import optimize In [2]: optimize.newton(lambda x: x**3 - x + 1, 0) Out[2]: 0.999999989980082 In [3]: optimize.newton(lambda x: x**3 - x + 1, 0.1) Out[3]: -1.324717957244745 In [4]: optimize.newton(lambda x: x**3 - x + 1, -0.1) Out[4]: -1.3247179572447458 These last two are the correct answer: In [6]: from numpy import roots In [7]: roots([1, 0, -1, 1]) Out[7]: array([-1.32471796+0.j , 0.66235898+0.56227951j, 0.66235898-0.56227951j]) The only thing that jumps out to me is that you've picked an initial guess at which the function has a zero second derivative. I'm not sure why that would cause the method to fail but Newton solvers are sensitive to various singularities in the input function. -- Oscar _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
First of all it sounds like your students are receiving a very good educational experience! typing help(optimize.newton) returns the following: fprime : function, optional The derivative of the function when available and convenient. If it is None (default), then the secant method is used. A look in https://en.wikipedia.org/wiki/Secant_method#Convergence returns the following: "If the initial values are not close enough to the root, then there is no guarantee that the secant method converges." So it looks like the behavior of scipy.optimize.newton in your example is neither unexpected nor a bug. On Thu, Feb 25, 2016 at 7:36 AM, Oscar Benjamin <[hidden email]> wrote: On 24 February 2016 at 19:43, Thomas Baruchel <[hidden email]> wrote: _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
Sorry - you are talking about the fact that it returns a wrong answer without indication that ti's wrong. Don't know if that's a "bug" but it's certainly incorrect behavior inherent to the secant method. On Thu, Feb 25, 2016 at 8:21 AM, David Mikolas <[hidden email]> wrote:
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by David Mikolas
Brent's method (scipy.optimize.brentq) is probably the best general
scalar (as opposed to multidimensional) root-finding algorithm, it's the scipy equivalent of MATLAB's fzero. Worth a read of Numerical Recipes to understand the different root-finding algorithms. The previous (2nd) edition is readable for free online: http://apps.nrbook.com/c/index.html and root-finding starts on p. 347. Chandrupatla's algorithm appears to be "better" than Brent's in the sense that both are robust and Chandrupatla's usually converges faster by a factor of 2-3x (and if not then follows Brent's within an extra iteration or two) -- but it's not well-known. I have a writeup here: http://www.embeddedrelated.com/showarticle/855.php On Wed, Feb 24, 2016 at 5:21 PM, David Mikolas <[hidden email]> wrote: > First of all it sounds like your students are receiving a very good > educational experience! > > typing help(optimize.newton) returns the following: > > fprime : function, optional > The derivative of the function when available and convenient. If it > is None (default), then the secant method is used. > > A look in https://en.wikipedia.org/wiki/Secant_method#Convergence returns > the following: > > "If the initial values are not close enough to the root, then there is no > guarantee that the secant method converges." > > So it looks like the behavior of scipy.optimize.newton in your example is > neither unexpected nor a bug. > > > > On Thu, Feb 25, 2016 at 7:36 AM, Oscar Benjamin <[hidden email]> > wrote: >> >> On 24 February 2016 at 19:43, Thomas Baruchel <[hidden email]> wrote: >> > with my students, this afternoon, I encountered the following annoying >> > behaviour of Scipy: >> > >> >>>> from scipy import optimize >> >>>> optimize.newton(lambda x: x**3-x+1, 0) >> > >> > 0.999999989980082 >> > >> > This is very weird, because the equation is rather simple with simple >> > coefficients, >> > and the initial point is simple also. >> >> Somehow the initial guess is problematic: >> >> In [1]: from scipy import optimize >> >> In [2]: optimize.newton(lambda x: x**3 - x + 1, 0) >> Out[2]: 0.999999989980082 >> >> In [3]: optimize.newton(lambda x: x**3 - x + 1, 0.1) >> Out[3]: -1.324717957244745 >> >> In [4]: optimize.newton(lambda x: x**3 - x + 1, -0.1) >> Out[4]: -1.3247179572447458 >> >> These last two are the correct answer: >> >> In [6]: from numpy import roots >> >> In [7]: roots([1, 0, -1, 1]) >> Out[7]: >> array([-1.32471796+0.j , 0.66235898+0.56227951j, >> 0.66235898-0.56227951j]) >> >> The only thing that jumps out to me is that you've picked an initial >> guess at which the function has a zero second derivative. I'm not sure >> why that would cause the method to fail but Newton solvers are >> sensitive to various singularities in the input function. >> >> -- >> Oscar >> _______________________________________________ >> SciPy-User mailing list >> [hidden email] >> https://mail.scipy.org/mailman/listinfo/scipy-user > > > > _______________________________________________ > SciPy-User mailing list > [hidden email] > https://mail.scipy.org/mailman/listinfo/scipy-user > SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
I use scipy.optimize/brentq frequently and so far have had no problems... except: if one of the limits is np.nan (which can happen without you necessarily knowing about it) then there is trouble. I'm having difficulty finding my report. I wrote this in Stackoverfow but for some reason it seems to be missing. Here is a residual copy on a "mirror" Briefly - it sometimes returns an incorrect/invalid answer but indicates convergence=True. It can also sometimes take time running to the iteration limit before giving an incorrect/inavlid answer. With the first limit set to np.nan, it returns a wrong answer (0.0), but indicates convergence: With the second limit set to np.nan, and disp = False, it runs to the iteration limit, then returns nan as the answer and indicates convergence: With the second limit set to np.nan, and disp = True, it runs to the iteration limit, then raises: On Thu, Feb 25, 2016 at 8:29 AM, Jason Sachs <[hidden email]> wrote: Brent's method (scipy.optimize.brentq) is probably the best general _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
Ah - I see what happened - 365 days with zero votes. Can be accessed and reopened if signed in http://stackoverflow.com/questions/27792180/brentq-gives-wrong-results-and-reports-converged-if-limits-contain-a-nan-how On Thu, Feb 25, 2016 at 9:05 AM, David Mikolas <[hidden email]> wrote:
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Oscar Benjamin
On Wed, 24 Feb 2016, Oscar Benjamin wrote:
> The only thing that jumps out to me is that you've picked an initial > guess at which the function has a zero second derivative. I'm not sure > why that would cause the method to fail but Newton solvers are > sensitive to various singularities in the input function. First, I thank you all for your answers. I am not absolutely sure it has to do with initial guess; as I said, the method works correctly when the derivative is provided or when it is computed by various methods. Furthermore (I forgot to tell it in my initial message), the function in Scipy would return the correct answer if some more terms were computed! By replacing if abs(p - p1) < tol: with if abs(p - p1) < tol and abs(p - p0) < tol: in Scipy's code, the correct answer is returned: of course this would imply computing one more term in most of the cases, but it wouldn't return too early in some weird cases where convergence is "fictive". Regards, -- Thomas Baruchel _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
On 25 February 2016 at 06:32, Thomas Baruchel <[hidden email]> wrote:
> On Wed, 24 Feb 2016, Oscar Benjamin wrote: >> >> The only thing that jumps out to me is that you've picked an initial >> guess at which the function has a zero second derivative. I'm not sure >> why that would cause the method to fail but Newton solvers are >> sensitive to various singularities in the input function. > > > First, I thank you all for your answers. I am not absolutely sure it has to > do > with initial guess; as I said, the method works correctly when the > derivative is > provided or when it is computed by various methods. > > Furthermore (I forgot to tell it in my initial message), the function in > Scipy > would return the correct answer if some more terms were computed! > > By replacing > > if abs(p - p1) < tol: > with > if abs(p - p1) < tol and abs(p - p0) < tol: > > in Scipy's code, the correct answer is returned: of course this would imply > computing one more term in most of the cases, but it wouldn't return too > early > in some weird cases where convergence is "fictive". Maybe that is worth doing then. However as David says these methods are only really expected to converge if the algorithm begins close to the root. So it might be that the proposed improvement helps in some cases. It's also possible that there are other cases where it could cause problems. Generally these kinds of solvers are not robust at globally searching for a root. The Newton/secant methods are about rapid improvement of an estimate of the root so you should locate the root approximately first. (I would use this as a lesson for my students in understanding that the need to provide an initial guess is not just a formality). In your particular example the function has a single global root but it also has two stationary points. The initial guess that you provide is between the two stationary points and the root is outside of that interval. There aren't any theorems I know of that give any guarantee that either method will converge in this sort of situation. If you want guaranteed convergence bracket the root and use bisect: this is guaranteed to converge in a fixed number of iterations (barring floating point weirdness). -- Oscar _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
Free forum by Nabble | Edit this page |