Repository: scipy/scipy
Pull request: BUG: ENH: Check if guess for newton is already zero before checking derivative
* fixes #8904 
* adds a test to see if `RuntimeWarning` is raised for case where `fval(p0)` is already zero where `p0` is the initial guess
* checks in `newton` if `fval(p0)` is already zero, before iterating for any method or evaluating derivatives `fprime` (or `fprime2`) if given, and returns `_ECONVERGED` instead of raising a `RuntimeWarning`
* fails the current test `scipy/optimize/tests/test_zeros.py::TestBasic::test_deriv_zero_warning` since these two are in opposition.
* see [post on SciPy Dev mailing list](https://mail.python.org/pipermail/scipy-dev/2018-June/022861.html) for more info

Problem statement:
BUG: if zero derivative at root, then Newton fails with RuntimeWarning
As indicated in [this discussion](../pull/8357#discussion_r192612946) by @pvanmulbregt

>... if the derivative is 0 at the root? E.g. `f(x)=x**3 - x**2`, starting at `x0=0` (immediately at the root), ... It looks like this might be considered a failure


### Reproducing code example:
```
from scipy import optimize

# use secant
r = optimize.newton(lambda x: x**3 - x**2, x0=0)
# -1.3552527156068808e-20  close enough to zero to be the correct answer

# now try Newton
r = optimize.newton(lambda x: x**3 - x**2, x0=0, fprime=lambda x: 3*x**2 - 2*x)
# raises RuntimeWarning but solution is zero

# or Halley's
r = optimize.newton(lambda x: x**3 - x**2, x0=0, fprime=lambda x: 3*x**2 - 2*x, fprime2=lambda x: 6*x - 2)
# also raises RuntimeWarning even though solution should be zero
```

### Error message:
```
~/miniconda/lib/python3.6/site-packages/scipy/optimize/zeros.py:169: RuntimeWarning: derivative was zero.
  warnings.warn(msg, RuntimeWarning)
```

### Scipy/Numpy/Python version information:
```
1.1.0 1.14.3 sys.version_info(major=3, minor=6, micro=5, releaselevel='final', serial=0)
```

Hints from issue discussion:
So I added a test for this:

```python
@pytest.mark.filterwarnings("error")
def test_gh8904_zeroder_at_root_fails():
    """Test that Newton or Halley don't warn if zero derivative at root"""

    # a function that has a zero derivative at it's root
    def f_zeroder_root(x):
        return x ** 3 - x ** 2

    # should work with secant
    r = zeros.newton(f_zeroder_root, x0=0)
    assert_allclose(r, 0, atol=1e-16)

    r = zeros.newton(f_zeroder_root, x0=[0]*10)
    assert_allclose(r, 0, atol=1e-16)

    # 1st derivative
    def fder(x):
        return 3 * x ** 2 - 2 * x

    # should work with newton and halley
    r = zeros.newton(f_zeroder_root, x0=0, fprime=fder)
    assert_allclose(r, 0, atol=1e-16)
    r = zeros.newton(f_zeroder_root, x0=0, fprime=fder,
                     fprime2=lambda x: 6 * x - 2)
    assert_allclose(r, 0, atol=1e-16)

    r = zeros.newton(f_zeroder_root, x0=[0]*10, fprime=fder)
    assert_allclose(r, 0, atol=1e-16)
    r = zeros.newton(f_zeroder_root, x0=[0]*10, fprime=fder,
                     fprime2=lambda x: 6 * x - 2)
    assert_allclose(r, 0, atol=1e-16)
```

and, if you remove the warning filter that raises warnings as errors, then it passes in both [master](../tree/67672b0f8) and #8357 because it just raises a `RuntimeWarning`, but since the initial guess is zero, it passes the assertions.

I think the issue that the only way to test for convergence is to calculate the Newton step and compare it to the tolerance. However if the Newton step is infinite, then there is no way to tell if the solution has been found. So I'm inclined to believe this is a feature and not a bug, because even if the result is correct, there's no way to know, and the outcome is that the solver failed to converge.

* In the scalar case as of #8803 if an infinite Newton step is discovered and the `full_output` flag is `True`, then it will yield an error code of -2 indicating a "convergence error"
* in the vector case, if the `failure_idx_flag` is true, then the indices of elements with an infinite Newton step will be `True` in `zero_der`, but they are also marked as failures. Therefore if both are true this would be equivalent to an scalar-case error code of -2, but if only failures is true, then one can conclude that it reached the maximum number of iterations.

One proposal to consider in #8357, is that we could change `failures` to `converged`, _ie_: take the opposite `converged = ~failures`? In my mind that implies a certain ambiguity which may be desired here, _ie_ it hasn't necessarily failed, it just didn't converge. Or we could change `failures` to `convergence_errors` which would match the root results object in the scalar case, although I don't like that option as much, because if the maximum iterations are reached it doesn't necessarily mean it wouldn't have still converged, so it's not an error. 
Hmm, maybe I'll walk back and say since we're searching for zeros, we can always compare the initial guess versus zero first? How about that? 
I think that is right, moving the evaluation of `fval = func(p0, *args)` prior to `fder = fprime(p0, *args)` and only evaluating `fprime ` if `fval != 0`,  would allow avoiding failure at roots with multiplicity.  If locally `f(x) ~ Cx**n`, then `f'(x) ~ nCx**(n-1) ~ (n/x)f(x)`, so one would expect `|f'(x)|`to be bigger than `f(x), except at the root itself when both are 0. 
So I can fix this issue by testing `fval` first:

```diff
diff --git a/scipy/optimize/zeros.py b/scipy/optimize/zeros.py
index 7272ff0fb..eb7c9fdea 100644
--- a/scipy/optimize/zeros.py
+++ b/scipy/optimize/zeros.py
@@ -187,6 +187,9 @@ def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50,
     # Multiply by 1.0 to convert to floating point.  We don't use float(x0)
     # so it still works if x0 is complex.
     p0 = 1.0 * x0
+    # check if initial guess is the solution
+    if func(p0, *args) == 0:
+        return _results_select(full_output, (p0, 0, 0, _ECONVERGED))
     funcalls = 0
     if fprime is not None:
         # Newton-Raphson method
```

... but then it causes [`scipy/optimize/tests/test_zeros.py::TestBasic::test_deriv_zero_warning`](../tree/2ce7116d2/scipy/optimize/tests/test_zeros.py#L99-L102) to fail. _NB_ `test_deriv_zero_warning` was written 8 years ago by @WarrenWeckesser  and explicitly tests the condition when `fval` and `fder` are both zero, which we are now considering to change to not raising `RuntimeWarning`.  So now I'm inclined to wait on #8904 and not block #8357 in order to allow the SciPy community to weigh in on whether they want to change when `RuntimeWarning` is raised. [I've raised the question on the SciPy Dev mailing list to get more feedback.](https://mail.python.org/pipermail/scipy-dev/2018-June/022861.html)

Here is the traceback:
<details>

```python
self = <scipy.optimize.tests.test_zeros.TestBasic object at 0x7f2062993898>

    def test_deriv_zero_warning(self):
        func = lambda x: x**2
        dfunc = lambda x: 2*x
>       assert_warns(RuntimeWarning, cc.newton, func, 0.0, dfunc)

dfunc      = <function TestBasic.test_deriv_zero_warning.<locals>.<lambda> at 0x7f20629956a8>
func       = <function TestBasic.test_deriv_zero_warning.<locals>.<lambda> at 0x7f20629958c8>
self       = <scipy.optimize.tests.test_zeros.TestBasic object at 0x7f2062993898>

scipy/optimize/tests/test_zeros.py:102:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
../../../../../venv/lib/python3.5/site-packages/numpy/testing/nose_tools/utils.py:1731: in assert_warns
    return func(*args, **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <contextlib._GeneratorContextManager object at 0x7f2062993dd8>, type = None, value = None, traceback = None

    def __exit__(self, type, value, traceback):
        if type is None:
            try:
>               next(self.gen)
E               AssertionError: No warning raised when calling newton

self       = <contextlib._GeneratorContextManager object at 0x7f2062993dd8>
traceback  = None
type       = None
value      = None

/usr/lib/python3.5/contextlib.py:66: AssertionError
```
</details>

In my discussions there were two examples associated with the function `f(x)=x**3 - x**2`. 
The first was `x0=0` where the initial value is the root. 
The second was `x0=0.5`, and this is the more important of the two.  One step of N-R takes `x0=0.5` to `x1=0`, a root with multiplicity > 1.

If the starting point is already a root, the user can determine that before starting. 
On the other hand, an intermediate iterate landing on a multiple root is something internal to the algorithm/implementation.
So is there _anyone_ who is **not** okay with changing the existing tests?
> NB test_deriv_zero_warning was written 8 years ago by @WarrenWeckesser and explicitly tests the condition when fval and fder are both zero

To be honest, I don't recall the full context of those changes, but it looks like I was doing some clean up (e.g. change `import *` to explicit imports), and added a test that reflected the current behavior of the implementation.  It does seem more reasonable that if you are at a root, then the derivative being zero shouldn't matter because there is no need for another iteration, and therefore no need for a warning.  (Clearly, though, if a root hasn't been found yet and the derivative is zero, it is impossible perform another "pure" Newton iteration.  But in that case, you need more than a warning.  If the algorithm is supposed to be a pure Newton method than a `RuntimeError` should be raised.)  
> If the algorithm is supposed to be a pure Newton method than a RuntimeError should be raised.

optimizers don't raise on convergence failure. A RuntimeError makes it difficult to see where the algorithm stopped.
That's true.  On the other hand, a warning that is unnoticed or ignored ("It's just a warning...") allows the code to continue on with a return value that is garbage.  Different use cases have different API requirements.
Don't they all have return values to indicate successful completion?

I only use brentq and fsolve in statsmodels, and I try them one after the other if one fails.
e.g. after fsolve
```
            if ier == 1 and np.abs(fval) < 1e-4 :
                success = 1
            else:
               <try something else>
```

Commit message: TST: remove last test for halley since fder isn't zero, not right test

Patch to review:
diff --git a/scipy/optimize/tests/test_zeros.py b/scipy/optimize/tests/test_zeros.py
index 15e054b9ca74..ddf77b2fa1ab 100644
--- a/scipy/optimize/tests/test_zeros.py
+++ b/scipy/optimize/tests/test_zeros.py
@@ -97,7 +97,7 @@ def test_newton_full_output(self):
                     x, r = zeros.newton(f1, x0, maxiter=iters, disp=True, **kwargs)
 
     def test_deriv_zero_warning(self):
-        func = lambda x: x**2
+        func = lambda x: x**2 - 2.0
         dfunc = lambda x: 2*x
         assert_warns(RuntimeWarning, cc.newton, func, 0.0, dfunc)
 
@@ -169,3 +169,36 @@ def f_2(x, *a):
     y = zeros.newton(f, z, args=coeffs, fprime=f_1, fprime2=f_2, tol=1e-6)
     # (-0.75000000000000078+1.1989578808281789j)
     assert_allclose(f(y, *coeffs), 0, atol=1e-6)
+
+
+@pytest.mark.filterwarnings("error")
+def test_gh8904_zeroder_at_root_fails():
+    """Test that Newton or Halley don't warn if zero derivative at root"""
+
+    # a function that has a zero derivative at it's root
+    def f_zeroder_root(x):
+        return x ** 3 - x ** 2
+
+    # should work with secant
+    r = zeros.newton(f_zeroder_root, x0=0)
+    assert_allclose(r, 0, atol=zeros._xtol, rtol=zeros._rtol)
+
+    # 1st derivative
+    def fder(x):
+        return 3 * x ** 2 - 2 * x
+
+    # should work with newton and halley
+    r = zeros.newton(f_zeroder_root, x0=0, fprime=fder)
+    assert_allclose(r, 0, atol=zeros._xtol, rtol=zeros._rtol)
+    r = zeros.newton(f_zeroder_root, x0=0, fprime=fder,
+                     fprime2=lambda x: 6 * x - 2)
+    assert_allclose(r, 0, atol=zeros._xtol, rtol=zeros._rtol)
+
+    # also test that if a root is found we do not raise RuntimeWarning even if
+    # the derivative is zero, EG: at x = 0.5, then fval = -0.125 and
+    # fder = -0.25 so the next guess is 0.5 - (-0.125/-0.5) = 0 which is the
+    # root, but if the solver continued with that guess, then it will calculate
+    # a zero derivative, so it should return the root w/o RuntimeWarning
+    r = zeros.newton(f_zeroder_root, x0=0.5, fprime=fder)
+    assert_allclose(r, 0, atol=zeros._xtol, rtol=zeros._rtol)
+    # doesn't apply to halley
diff --git a/scipy/optimize/zeros.py b/scipy/optimize/zeros.py
index 7272ff0fbd48..9c5bee352ac0 100644
--- a/scipy/optimize/zeros.py
+++ b/scipy/optimize/zeros.py
@@ -3,7 +3,7 @@
 import warnings
 
 from . import _zeros
-from numpy import finfo, sign, sqrt
+from numpy import finfo
 
 _iter = 100
 _xtol = 2e-12
@@ -191,14 +191,18 @@ def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50,
     if fprime is not None:
         # Newton-Raphson method
         for itr in range(maxiter):
+            # first evaluate fval
+            fval = func(p0, *args)
+            funcalls += 1
+            # If fval is 0, a root has been found, then terminate
+            if fval == 0:
+                return _results_select(full_output, (p0, funcalls, itr, _ECONVERGED))
             fder = fprime(p0, *args)
             funcalls += 1
             if fder == 0:
                 msg = "derivative was zero."
                 warnings.warn(msg, RuntimeWarning)
                 return _results_select(full_output, (p0, funcalls, itr + 1, _ECONVERR))
-            fval = func(p0, *args)
-            funcalls += 1
             newton_step = fval / fder
             if fprime2 is None:
                 # Newton step


Code review comments:
diff_hunk: @@ -169,3 +169,36 @@ def f_2(x, *a):
     y = zeros.newton(f, z, args=coeffs, fprime=f_1, fprime2=f_2, tol=1e-6)
     # (-0.75000000000000078+1.1989578808281789j)
     assert_allclose(f(y, *coeffs), 0, atol=1e-6)
+
+
+@pytest.mark.filterwarnings("error")
path: scipy/optimize/tests/test_zeros.py
text: @user1:
this line can be removed, this is the default behavior (see `pytest.ini` line 7 in the repo root).

@author:
Hi @user1, even though the `error` warning filter is set by default, it didn't seem to work when I use `python runtest.py -t scipy/optimize/tests/test_zeros.py` and more importantly I feel that it's more explicit to leave the decorator even if it is redundant.

@user1:
We want to control this globally. For one, we're switching this behavior between development versions and releases (we do want errors for the former, not for the latter).

@author:
Okay