Quantcast

[SciPy-User] Is this weird answer a bug in optimize.newton?

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[SciPy-User] Is this weird answer a bug in optimize.newton?

Thomas Baruchel
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

Oscar Benjamin
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

David Mikolas
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.


"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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

David Mikolas
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:
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.


"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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

Jason Sachs
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

David Mikolas
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
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


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

David Mikolas
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:
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
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



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

Thomas Baruchel
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is this weird answer a bug in optimize.newton?

Oscar Benjamin
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
Loading...