Repository: scipy/scipy
Pull request: ENH: optimize.root_scalar: allow newton without f', secant without x1
#### Reference issue
Closes gh-17570

#### What does this implement/fix?
`root_scalar` does not allow `method='newton'` to be used without an analytical derivative even though `scipy.optimize.root` does, and it does not allow `method='secant'` to be used without an explicit `x1` even though `scipy.optimize.newton` does. This PR removes `root_scalar`'s restrictions on use of these methods and updates the defaults method resolution accordingly.

#### Additional information
I don't think there are any BC concerns because it is [documented ](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root_scalar.html) that
![image](https://user-images.githubusercontent.com/6570539/210158382-4d6b25b3-b031-4609-9151-afd407069ad0.png)

Problem statement:
ENH: optimize.root_scalar: default to `newton` when only `x0` is provided
### Is your feature request related to a problem? Please describe.

When `method` is not specified, `optimize.root_scalar` attempts to select the best method given the information available. I was always surprised that the following:

```python3
from scipy import optimize

def f(x):
    return x**2 - 1

res = optimize.root_scalar(f, x0=2)
```

fails with the error:

```
Traceback (most recent call last):

  File "C:\Users\matth\AppData\Local\Temp\ipykernel_35340\1932188727.py", line 6, in <module>
    res = optimize.root_scalar(f, x0=2)

  File "C:\Users\matth\scipy\scipy\optimize\_root_scalar.py", line 280, in root_scalar
    raise ValueError('x1 must not be None for %s' % method)

ValueError: x1 must not be None for secant
```

I expected Newton's method to be used because `x0` is enough information to use Newton's method.

When I try to force Newton to be used,

```python3
res = optimize.root_scalar(f, x0=2, method='newton')
```

we see why `newton` wasn't selected automatically:

```
Traceback (most recent call last):

  File "C:\Users\matth\AppData\Local\Temp\ipykernel_35340\1942526262.py", line 1, in <module>
    res = optimize.root_scalar(f, x0=2, method='newton')

  File "C:\Users\matth\scipy\scipy\optimize\_root_scalar.py", line 289, in root_scalar
    raise ValueError('fprime must be specified for %s' % method)

ValueError: fprime must be specified for newton
```

But `scipy.optimize.root` does not require an analytical derivative, so I don't think `root_scalar` should. 



### Describe the solution you'd like.

Use numerical differentiation for gradient-based `root_scalar` methods when `fprime` is not provided, and default to a gradient-based method when only `x0` is provided.

### Describe alternatives you've considered.

_No response_

### Additional context (e.g. screenshots, GIFs)

I'll do it. I just thought I'd create an issue to see if there were others that supported the idea.

Hints from issue discussion:
+1 from me. I don't see why we would be required to provide fprime.

Commit message: ENH: optimize.root_scalar: allow newton without f', secant without x1

Patch to review:
diff --git a/scipy/optimize/_root_scalar.py b/scipy/optimize/_root_scalar.py
index 54ffd95c3dc2..6f75150a450c 100644
--- a/scipy/optimize/_root_scalar.py
+++ b/scipy/optimize/_root_scalar.py
@@ -9,6 +9,8 @@
 import numpy as np
 
 from . import _zeros_py as optzeros
+from ._differentiable_functions import ScalarFunction
+from ._hessian_update_strategy import HessianUpdateStrategy
 
 __all__ = ['root_scalar']
 
@@ -253,8 +255,10 @@ def root_scalar(f, args=(), method=None, bracket=None,
                     method = 'halley'
                 else:
                     method = 'newton'
-            else:
+            elif x1 is not None:
                 method = 'secant'
+            else:
+                method = 'newton'
     if not method:
         raise ValueError('Unable to select a solver as neither bracket '
                          'nor starting point provided.')
@@ -276,8 +280,6 @@ def root_scalar(f, args=(), method=None, bracket=None,
     elif meth in ['secant']:
         if x0 is None:
             raise ValueError('x0 must not be None for %s' % method)
-        if x1 is None:
-            raise ValueError('x1 must not be None for %s' % method)
         if 'xtol' in kwargs:
             kwargs['tol'] = kwargs.pop('xtol')
         r, sol = methodc(f, x0, args=args, fprime=None, fprime2=None,
@@ -286,7 +288,23 @@ def root_scalar(f, args=(), method=None, bracket=None,
         if x0 is None:
             raise ValueError('x0 must not be None for %s' % method)
         if not fprime:
-            raise ValueError('fprime must be specified for %s' % method)
+            # approximate fprime with finite differences
+
+            class NoUpdate(HessianUpdateStrategy):
+                # dummy implementation, since we don't need the Hessian
+                def initialize(self, *args, **kwargs):
+                    return None
+
+                def update(self, *args, **kwargsself):
+                    return None
+
+            func = ScalarFunction(f, x0, tuple(), grad='2-point',
+                                  hess=NoUpdate(), finite_diff_rel_step=None,
+                                  finite_diff_bounds=[-np.inf, np.inf])
+
+            def fprime(x):
+                return func.grad(x).reshape(np.shape(x))
+
         if 'xtol' in kwargs:
             kwargs['tol'] = kwargs.pop('xtol')
         r, sol = methodc(f, x0, args=args, fprime=fprime, fprime2=None,
diff --git a/scipy/optimize/tests/test_zeros.py b/scipy/optimize/tests/test_zeros.py
index 3fabf9723d6f..e97b6ee79cb3 100644
--- a/scipy/optimize/tests/test_zeros.py
+++ b/scipy/optimize/tests/test_zeros.py
@@ -276,6 +276,9 @@ def test_newton_by_name(self):
         for f, f_1, f_2 in [(f1, f1_1, f1_2), (f2, f2_1, f2_2)]:
             r = root_scalar(f, method='newton', x0=3, fprime=f_1, xtol=1e-6)
             assert_allclose(f(r.root), 0, atol=1e-6)
+        for f, f_1, f_2 in [(f1, f1_1, f1_2), (f2, f2_1, f2_2)]:
+            r = root_scalar(f, method='newton', x0=3, xtol=1e-6)  # without f'
+            assert_allclose(f(r.root), 0, atol=1e-6)
 
     def test_secant_by_name(self):
         r"""Invoke secant through root_scalar()"""
@@ -284,6 +287,9 @@ def test_secant_by_name(self):
             assert_allclose(f(r.root), 0, atol=1e-6)
             r = root_scalar(f, method='secant', x0=3, x1=5, xtol=1e-6)
             assert_allclose(f(r.root), 0, atol=1e-6)
+        for f, f_1, f_2 in [(f1, f1_1, f1_2), (f2, f2_1, f2_2)]:
+            r = root_scalar(f, method='secant', x0=3, xtol=1e-6)  # without x1
+            assert_allclose(f(r.root), 0, atol=1e-6)
 
     def test_halley_by_name(self):
         r"""Invoke halley through root_scalar()"""
@@ -293,13 +299,11 @@ def test_halley_by_name(self):
             assert_allclose(f(r.root), 0, atol=1e-6)
 
     def test_root_scalar_fail(self):
-        with pytest.raises(ValueError):
-            root_scalar(f1, method='secant', x0=3, xtol=1e-6)  # no x1
-        with pytest.raises(ValueError):
-            root_scalar(f1, method='newton', x0=3, xtol=1e-6)  # no fprime
-        with pytest.raises(ValueError):
+        message = 'fprime2 must be specified for halley'
+        with pytest.raises(ValueError, match=message):
             root_scalar(f1, method='halley', fprime=f1_1, x0=3, xtol=1e-6)  # no fprime2
-        with pytest.raises(ValueError):
+        message = 'fprime must be specified for halley'
+        with pytest.raises(ValueError, match=message):
             root_scalar(f1, method='halley', fprime2=f1_2, x0=3, xtol=1e-6)  # no fprime
 
     def test_array_newton(self):
@@ -463,6 +467,32 @@ def test_maxiter_int_check(self):
                     match="'float' object cannot be interpreted as an integer"):
                 method(f1, 0.0, 1.0, maxiter=72.45)
 
+    def test_gh17570_defaults(self):
+        # Previously, when fprime was not specified, root_scalar would default
+        # to secant. When x1 was not specified, secant failed.
+        # Check that without fprime, the default is secant if x1 is specified
+        # and newton otherwise.
+        res_newton_default = root_scalar(f1, method='newton', x0=3, xtol=1e-6)
+        res_secant_default = root_scalar(f1, method='secant', x0=3, x1=2,
+                                         xtol=1e-6)
+        # `newton` uses the secant method when `x1` and `x2` are specified
+        res_secant = newton(f1, x0=3, x1=2, tol=1e-6, full_output=True)[1]
+
+        # all three foun a root
+        assert_allclose(f1(res_newton_default.root), 0, atol=1e-6)
+        assert_allclose(f1(res_secant_default.root), 0, atol=1e-6)
+        assert_allclose(f1(res_secant.root), 0, atol=1e-6)
+
+        # Defaults are correct
+        assert (res_secant_default.root
+                == res_secant.root
+                != res_newton_default.iterations)
+        assert (res_secant_default.iterations
+                == res_secant_default.function_calls - 1  # true for secant
+                == res_secant.iterations
+                != res_newton_default.iterations
+                == res_newton_default.function_calls/2)  # newton 2-point diff
+
 
 def test_gh_5555():
     root = 0.1


Code review comments:
diff_hunk: @@ -286,7 +288,23 @@ def root_scalar(f, args=(), method=None, bracket=None,
         if x0 is None:
             raise ValueError('x0 must not be None for %s' % method)
         if not fprime:
-            raise ValueError('fprime must be specified for %s' % method)
+            # approximate fprime with finite differences
+
+            class NoUpdate(HessianUpdateStrategy):
+                # dummy implementation, since we don't need the Hessian
+                def initialize(self, *args, **kwargs):
+                    return None
+
+                def update(self, *args, **kwargsself):
+                    return None
+
+            func = ScalarFunction(f, x0, tuple(), grad='2-point',
path: scipy/optimize/_root_scalar.py
text: @author:
I can change this to use `approx_derivative`.