ENH: Only use Halley's adjustment to Newton if close enough.

scipy/optimize/tests/test_zeros.py

@@ -7,9 +7,10 @@
                            assert_allclose,
                            assert_equal)
 import numpy as np
+from numpy import finfo, power
 
 from scipy.optimize import zeros as cc
-from scipy.optimize import zeros
+from scipy.optimize import zeros, newton
 
 # Import testing parameters
 from scipy.optimize._tstutils import functions, fstrings
@@ -357,3 +358,29 @@ def fder2(x):
     r = zeros.newton(f_zeroder_root, x0=[0.5]*10, fprime=fder)
     assert_allclose(r, 0, atol=zeros._xtol, rtol=zeros._rtol)
     # doesn't apply to halley
+
+
+def test_gh_8881():
+    r"""Test that Halley's method realizes that the 2nd order adjustment
+    is too big and drops off to the 1st order adjustment."""
+    n = 9
+
+    def f(x):
+        return power(x, 1.0/n) - power(n, 1.0/n)
+
+    def fp(x):
+        return power(x, (1.0-n)/n)/n
+
+    def fpp(x):
+        return power(x, (1.0-2*n)/n) * (1.0/n) * (1.0-n)/n
+
+    x0 = 0.1
+    # The root is at x=9.
+    # The function has positive slope, x0 < root.
+    # Newton succeeds in 8 iterations
+    rt, r = newton(f, x0, fprime=fp, full_output=True)
+    assert(r.converged)
+    # Before the Issue 8881/PR 8882, halley would send x in the wrong direction.
+    # Check that it now succeeds.
+    rt, r = newton(f, x0, fprime=fp, fprime2=fpp, full_output=True)
+    assert(r.converged)

scipy/optimize/zeros.py

@@ -264,14 +264,19 @@ class of similar problems can be solved together.
                 warnings.warn(msg, RuntimeWarning)
                 return _results_select(full_output, (p0, funcalls, itr + 1, _ECONVERR))
             newton_step = fval / fder
-            if fprime2 is None:
-                # Newton step
-                p = p0 - newton_step
-            else:
+            if fprime2:
                 fder2 = fprime2(p0, *args)
                 funcalls += 1
-                # Halley's method
-                p = p0 - newton_step / (1.0 - 0.5 * newton_step * fder2 / fder)
+                # Halley's method:
+                #   newton_step /= (1.0 - 0.5 * newton_step * fder2 / fder)
+                # Only do it if denominator stays close enough to 1
+                # Rationale:  If 1-adj < 0, then Halley sends x in the
+                # opposite direction to Newton.  Doesn't happen if x is close
+                # enough to root.
+                adj = newton_step * fder2 / fder / 2
+                if abs(adj) < 1:
+                    newton_step /= 1.0 - adj
+            p = p0 - newton_step
             if abs(p - p0) < tol:
                 return _results_select(full_output, (p, funcalls, itr + 1, _ECONVERGED))
             p0 = p