find_root() is broken when interval borders cannot be evaluated

[sagetrac-mabshoff]

Reported in http://groups.google.com/group/sage-support/browse_thread/thread/40da8039090c3e8a

Hi, I'm trying out SAGE for the first time, so I entered what you 
suggested (see above). 
Now, from the plot, it there seems to be no other roots between 0 and 2 
so I entered 
sage: find_root(x^2*log(x,2)-1,0, 2) 
and got the root = 0.0 
what am I missing here? 
TIA, 
AJG 

But note the following:

sage: find_root(1/(x-1)+1,0, 2) 
0.0 
sage: find_root(1/(x-1)+1,0.00001, 2) 
1.0000000000011564 

Cheers,

Michael

CC: @kcrisman

Component: numerical

Stopgaps: #12709

Author: Eran Assaf

Branch/Commit: e6b7c1e

Reviewer: Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/4942

[sagetrac-mabshoff]

comment:3

This is a critical bug and ought to be fixed in 3.3.

Note that #3870 might be a dupe of this bug.

Cheers,

Michael

[mwhansen]

comment:4

It seems this is a problem with Scipy:

In [16]: def f(x):         
   ....:     return 1.0/(x-1.0)+1.0
   ....: 

In [17]: import scipy.optimize

In [18]: scipy.optimize.brentq(f, 0, 2)
Out[18]: 0.0

In [19]: f(0.001)
Out[19]: -0.0010010010010010895

In [20]: f(2)
Out[20]: 2.0

In [21]: scipy.optimize.brentq(f, 0.001, 2)                                                   
Out[21]: 1.0000000000007283

In [22]: f(1.0000000000007283)
Out[22]: 1373048666882.2488

[sagetrac-cwitty]

comment:5

There are at least a couple of issues here. First, brentq is a variant of a bisection-based solver; if you use any bisection-based solver to find a zero of 1/(x-1) between 0 and 2, it will narrow down and return something very close to 1. So if we don't like that, we should use a different solver (or at least try to check the output; for instance, a simple check that f(x) is "small" would detect this particular problem).

Second, find_root tries to verify that the function evaluates to different signs at the endpoints of the interval (as required by brentq); but it doesn't check the function evaluation results for NaN. In the original test case, fast_float(f)(0) gives NaN.

[sagetrac-mabshoff]

comment:6

Better luck in 3.4.1. Unfortunately this either requires testing of the result of scipy or some deeper surgery in Scipy.

Cheers,

Michael

[williamstein]

comment:7

If we've released for months and months without fixing this, it doesn't make sense to keep it as a blocker.

[jbalakrishnan]

Stopgaps: #12709

[assaferan]

Branch: u/assaferan/find_root___is_broken_when_interval_borders_cannot_be_evaluated

[assaferan]

Author: Eran Assaf

[assaferan]

Commit: 91f14d3

[assaferan]

comment:15

Hi, added two small validity checks:

1. If one of the endpoints is evaluated to NaN we seek a nearby point in the interval which can be evaluated.
2. If the value of the function at the root found is "large", raise an error that we could not find it.

---

New commits:

91f14d3

Added two validity checks - handles the case of NaN values and failure of the algorithm to find a root

[tscrim]

comment:17

I am not sure 1 is necessarily the best solution to this because what if you get a function that always evaluates to NaN as you increase/decrease the endpoints? For instance

sage: f(x) = 0.0 / max(0, x)

will be NaN for infinitely many values. So your current test means this runs forever:

sage: find_root(f, -1, 0)

(before it simply gave a wrong value).

Also, I think for 2 you should raise a NotImplementedError as I think that more accurately reflects the situation.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

8feb5a9

Fixed bug in find_root when full_output=True (size check assumed it was false)

[sagetrac-git]

Changed commit from 91f14d3 to 8feb5a9

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

d67280b

Following the suggestion of tscrim, changed handling of NaN case to finding minimal and maximal values, and raising a NotImplementedError when brentq fails to find a root.

[sagetrac-git]

Changed commit from 8feb5a9 to d67280b

[assaferan]

comment:20

Fixed the bugs and changed behaviour in both cases, as suggested by tscrim

[tscrim]

Reviewer: Travis Scrimshaw

[tscrim]

comment:21

Thanks. Looks better now. A few more little things:

• ticket 4942 -> :trac:`4942` in the documentation.

• Could you add the test from comment:17.

• This change:

          Traceback (most recent call last):
  -           ...
  +       ...
          NotImplementedError: Brent's method failed to find a zero for f on the interval

• Instead of using ... for imprecision, it would be better to use # abs tol (or a # rel tol).

• if statements do not need outer parentheses in Python, so remove them from if (full_output): and the outermost pair from the other if statement 4 lines down.

[sagetrac-git]

Changed commit from d67280b to e6b7c1e

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

e6b7c1e

Some minor repairs - fixed inadequacies in documentation, added example of a function which evaluates to NaN on the entire interval (following review of tscrim)

[tscrim]

comment:23

Thank you. LGTM.

[vbraun]

Changed branch from u/assaferan/find_root___is_broken_when_interval_borders_cannot_be_evaluated to e6b7c1e