Repository: scipy/scipy
Pull request: MAINT: optimize.bracket: don't fail silently
#### Reference issue
Closes gh-14858

#### What does this implement/fix?
gh-14858 reported that `scipy.optimize.bracket` may silently return an invalid bracket, causing scalar minimization functions to fail. This PR resolves the issue by checking for failure before returning. It also begins to improve the function documentation and tests.

_Update: `bracket` is used in `brent`, `golden`, and `minimize_scalar` with methods `'brent'` and `'golden'`. It's probably appropriate for `brent` and `golden` to raise errors when `bracket` fails, but `minimize_scalar` should terminate gracefully with `success=False`. This requires plumbing the status information from `bracket` through the optimizers, which could this PR quite a bit more complex than intended - see a6aa40db0ca4e2a76d4aa369346700a48d49f851. I ended up refactoring it that makes it much simpler, albeit a bit unorthodox._

#### Additional information
In main, `bracket` is not tested outside of `minimize_scalar`, and the documentation is not explicit about the function's behavior. This is intended to be an _improvement_, but more work is certainly needed in followup PRs. For instance, this function still needs:

- A reference
- A description of what the function does and how to select the initial points. Besides expanding the bracket by the golden ratio, I don't really know the other things it tries.
- Unit tests (separate from `minimize_scalar`)
- (Possibly) Improvements to mitigate difficulties caused be erroneous documentation (https://github.com/scipy/scipy/issues/5899#issuecomment-1368781013)
 
For this PR, let's stick within the existing scope of fixing gh-14858 and these doc improvements so as not to delay progress. 

Problem statement:
BUG: scipy.optimize.bracket sometimes fails silently
### Describe your issue.

Sometimes scipy.optimize.bracket fails to find a bracket (e.g. see example below).  Thats probably ok.  There are some weird functions out there.  But at least one of the following should probably be true

- scipy.optimize.bracket should throw an exception if it fails
- functions that use scipy.optimize.bracket (e.g. brent) should be aware that it might return an invalid bracket

Otherwise you can have people using functions like brent and getting answers that don't make sense.

### Reproducing Code Example

```python
import scipy.optimize
import numpy

def f(x):
    return x**2 if ((-1<x)&(x<1)) else 100.0
brack=scipy.optimize.bracket(f,-1,1)[:3]

# this isn't a valid bracket
print(brack)  # <-- (-1, 1, 4.2360679999999995) 
# f(-1)=f(1)=f(4.23)=100, so thats not a valid bracket.
# valid bracket should have f(a) > f(b) < f(c)

# now if you ask brent to find a bracket for you,
# it calls scipy.optimize.bracket, gets an invalid bracket (namely -1,1,4.23)
# tries to use this invalid bracket, and yields a nonsense answer 
print(scipy.optimize.brent(f,brack=(-1,1)))  # <--4.23, not the argmin

# though if you explicitly tell brent to use this brack
# it will notice that it is a  Bad Bracket
print(scipy.optimize.brent(f,brack=brack)) # <-- raises ValueError (as it should)
```

BTW I'm not entirely sure what the right definition of "failure" is.   If f(xa)=f(xb)<f(xc) then it seems to me that this is actually still a useable bracket.  But my 1d optimization intuitions have failed me before.

### Error message

```shell
No problematic error messages.  The problem is the lack of error messages :).
```


### SciPy/NumPy/Python version information

1.7.1 1.19.5 sys.version_info(major=3, minor=9, micro=6, releaselevel='final', serial=0)

Hints from issue discussion:
I think the problem is that `optimize.bracket` does not check the results whether the results satisfies " f(xa) > f(xb) < f(xc)" or not here.
https://github.com/scipy/scipy/blob/47bb6febaa10658c72962b9615d5d5aa2513fa3a/scipy/optimize/optimize.py#L2661-L2667


I came up with two options to fix the issue:
[Option 1]
Throwing exception here if the results do not satisfy " f(xa) > f(xb) < f(xc)" like "optimize.brent" and "optimize.golden".
But this might break someone's code (backward compatibility).

[Option 2]
Update the doc that the objective function criteria from " f(xa) > f(xb) < f(xc)" to " f(xa) >= f(xb) <= f(xc)".
Both sound good to me!  If we go with option 2 we probably ought to also update brent and golden to be consistent (i.e. with option 2, pathological cases like one I noticed won't throw an error in bracket but brent and golden would throw an error after noticing that bracket failed).

Commit message: MAINT: optimize.bracket: correct error condition, add Raises section to docs

Patch to review:
diff --git a/scipy/optimize/_minimize.py b/scipy/optimize/_minimize.py
index 04bb46e4f2dc..b2909a93b79a 100644
--- a/scipy/optimize/_minimize.py
+++ b/scipy/optimize/_minimize.py
@@ -19,7 +19,8 @@
                         _minimize_bfgs, _minimize_newtoncg,
                         _minimize_scalar_brent, _minimize_scalar_bounded,
                         _minimize_scalar_golden, MemoizeJac, OptimizeResult,
-                        _wrap_callback)
+                        _wrap_callback, BracketError,
+                        _recover_from_bracket_error)
 from ._trustregion_dogleg import _minimize_dogleg
 from ._trustregion_ncg import _minimize_trust_ncg
 from ._trustregion_krylov import _minimize_trust_krylov
@@ -924,14 +925,16 @@ def minimize_scalar(fun, bracket=None, bounds=None, args=(),
     if meth == '_custom':
         res = method(fun, args=args, bracket=bracket, bounds=bounds, **options)
     elif meth == 'brent':
-        res = _minimize_scalar_brent(fun, bracket, args, **options)
+        res = _recover_from_bracket_error(_minimize_scalar_brent,
+                                          fun, bracket, args, **options)
     elif meth == 'bounded':
         if bounds is None:
             raise ValueError('The `bounds` parameter is mandatory for '
                              'method `bounded`.')
         res = _minimize_scalar_bounded(fun, bounds, args, **options)
     elif meth == 'golden':
-        res = _minimize_scalar_golden(fun, bracket, args, **options)
+        res = _recover_from_bracket_error(_minimize_scalar_golden,
+                                          fun, bracket, args, **options)
     else:
         raise ValueError('Unknown solver %s' % method)
 
diff --git a/scipy/optimize/_optimize.py b/scipy/optimize/_optimize.py
index 21b2ace1af58..3a813d6cb9b1 100644
--- a/scipy/optimize/_optimize.py
+++ b/scipy/optimize/_optimize.py
@@ -2877,20 +2877,19 @@ def _minimize_scalar_golden(func, brack=None, args=(),
 
 def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
     """
-    Bracket the minimum of the function.
+    Bracket the minimum of a function.
 
     Given a function and distinct initial points, search in the
     downhill direction (as defined by the initial points) and return
-    new points xa, xb, xc that bracket the minimum of the function
-    f(xa) > f(xb) < f(xc). It doesn't always mean that obtained
-    solution will satisfy xa<=x<=xb.
+    three points that bracket the minimum of the function.
 
     Parameters
     ----------
     func : callable f(x,*args)
         Objective function to minimize.
     xa, xb : float, optional
-        Bracketing interval. Defaults `xa` to 0.0, and `xb` to 1.0.
+        Initial points. Defaults `xa` to 0.0, and `xb` to 1.0.
+        A local minimum need not be contained within this interval.
     args : tuple, optional
         Additional arguments (if present), passed to `func`.
     grow_limit : float, optional
@@ -2898,15 +2897,29 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
     maxiter : int, optional
         Maximum number of iterations to perform. Defaults to 1000.
 
+    Raises
+    ------
+    BracketError
+        If no valid bracket is found before the algorithm terminates.
+        See notes for conditions of a valid bracket.
+
     Returns
     -------
     xa, xb, xc : float
-        Bracket.
+        Final points of the bracket.
     fa, fb, fc : float
-        Objective function values in bracket.
+        Objective function values at the bracket points.
     funcalls : int
         Number of function evaluations made.
 
+    Notes
+    -----
+    The algorithm attempts to find three strictly ordered points (i.e.
+    :math:`x_a < x_b < x_c` or :math:`x_c < x_b < x_a`) satisfying
+    :math:`f(x_b) ≤ f(x_a)` and :math:`f(x_b) ≤ f(x_c)`, where one of the
+    inequalities must be satistfied strictly and all :math:`x_i` must be
+    finite.
+
     Examples
     --------
     This function can find a downward convex region of a function:
@@ -2918,7 +2931,7 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
     ...     return 10*x**2 + 3*x + 5
     >>> x = np.linspace(-2, 2)
     >>> y = f(x)
-    >>> init_xa, init_xb = 0, 1
+    >>> init_xa, init_xb = 0.1, 1
     >>> xa, xb, xc, fa, fb, fc, funcalls = bracket(f, xa=init_xa, xb=init_xb)
     >>> plt.axvline(x=init_xa, color="k", linestyle="--")
     >>> plt.axvline(x=init_xb, color="k", linestyle="--")
@@ -2928,9 +2941,17 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
     >>> plt.plot(xc, fc, "bx")
     >>> plt.show()
 
+    Note that both initial points were to the right of the minimum, and the
+    third point was found in the "downhill" direction: the direction
+    in which the function appeared to be decreasing (to the left).
+    The final points are strictly ordered, and the function value
+    at the middle point is less than the function values at the endpoints;
+    it follows that a minimum must lie within the bracket.
+
     """
     _gold = 1.618034  # golden ratio: (1.0+sqrt(5.0))/2.0
     _verysmall_num = 1e-21
+    xa, xb = np.asarray([xa, xb])
     fa = func(*(xa,) + args)
     fb = func(*(xb,) + args)
     if (fa < fb):                      # Switch so fa > fb
@@ -2950,8 +2971,11 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
             denom = 2.0 * val
         w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom
         wlim = xb + grow_limit * (xc - xb)
+        msg = ("No valid bracket was found before the iteration limit was "
+               "reached. Consider trying different initial points or "
+               "increasing `maxiter`.")
         if iter > maxiter:
-            raise RuntimeError("Too many iterations.")
+            raise RuntimeError(msg)
         iter += 1
         if (w - xc) * (xb - w) > 0.0:
             fw = func(*((w,) + args))
@@ -2961,11 +2985,11 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
                 xb = w
                 fa = fb
                 fb = fw
-                return xa, xb, xc, fa, fb, fc, funcalls
+                break
             elif (fw > fb):
                 xc = w
                 fc = fw
-                return xa, xb, xc, fa, fb, fc, funcalls
+                break
             w = xc + _gold * (xc - xb)
             fw = func(*((w,) + args))
             funcalls += 1
@@ -2994,9 +3018,59 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
         fa = fb
         fb = fc
         fc = fw
+
+    # three conditions for a valid bracket
+    cond1 = (fb < fc and fb <= fa) or (fb < fa and fb <= fc)
+    cond2 = (xa < xb < xc or xc < xb < xa)
+    cond3 = np.isfinite(xa) and np.isfinite(xb) and np.isfinite(xc)
+    msg = ("The algorithm terminated without finding a valid bracket. "
+           "Consider trying different initial points.")
+    if not (cond1 and cond2 and cond3):
+        e = BracketError(msg)
+        e.data = (xa, xb, xc, fa, fb, fc, funcalls)
+        raise e
+
     return xa, xb, xc, fa, fb, fc, funcalls
 
 
+class BracketError(RuntimeError):
+    pass
+
+
+def _recover_from_bracket_error(solver, fun, bracket, args, **options):
+    # `bracket` was originally written without checking whether the resulting
+    # bracket is valid. `brent` and `golden` built on top of it without
+    # checking the returned bracket for validity, and their output can be
+    # incorrect without warning/error if the original bracket is invalid.
+    # gh-14858 noticed the problem, and the following is the desired
+    # behavior:
+    # - `scipy.optimize.bracket`, `scipy.optimize.brent`, and
+    #   `scipy.optimize.golden` should raise an error if the bracket is
+    #   invalid, as opposed to silently returning garbage
+    # - `scipy.optimize.minimize_scalar` should return with `success=False`
+    #   and other information
+    # The changes that would be required to achieve this the traditional
+    # way (`return`ing all the required information from bracket all the way
+    # up to `minimizer_scalar`) are extensive and invasive. (See a6aa40d.)
+    # We can achieve the same thing by raising the error in `bracket`, but
+    # storing the information needed by `minimize_scalar` in the error object,
+    # and intercepting it here.
+    try:
+        res = solver(fun, bracket, args, **options)
+    except BracketError as e:
+        msg = str(e)
+        xa, xb, xc, fa, fb, fc, funcalls = e.data
+        xs, fs = [xa, xb, xc], [fa, fb, fc]
+        if np.any(np.isnan([xs, fs])):
+            x, fun = np.nan, np.nan
+        else:
+            imin = np.argmin(xs)
+            x, fun = xs[imin], fs[imin]
+        return OptimizeResult(fun=fun, nfev=funcalls, x=x,
+                              nit=0, success=False, message=msg)
+    return res
+
+
 def _line_for_search(x0, alpha, lower_bound, upper_bound):
     """
     Given a parameter vector ``x0`` with length ``n`` and a direction
@@ -3089,7 +3163,9 @@ def myfunc(alpha):
         return ((fval, p, xi) if fval is not None else (func(p), p, xi))
     elif lower_bound is None and upper_bound is None:
         # non-bounded minimization
-        alpha_min, fret, _, _ = brent(myfunc, full_output=1, tol=tol)
+        res = _recover_from_bracket_error(_minimize_scalar_brent,
+                                          myfunc, None, tuple(), xtol=tol)
+        alpha_min, fret = res.x, res.fun
         xi = alpha_min * xi
         return squeeze(fret), p + xi, xi
     else:
diff --git a/scipy/optimize/tests/test_optimize.py b/scipy/optimize/tests/test_optimize.py
index 6505423b7324..559a99d4588e 100644
--- a/scipy/optimize/tests/test_optimize.py
+++ b/scipy/optimize/tests/test_optimize.py
@@ -1971,6 +1971,47 @@ def f(x):
         options['disp'] = False
         optimize.minimize_scalar(f, method=method, options=options, **kwargs)
 
+class TestBracket:
+
+    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
+    def test_errors_and_status_false(self):
+        # Check that `bracket` raises the errors it is supposed to
+        def f(x):  # gh-14858
+            return x**2 if ((-1 < x) & (x < 1)) else 100.0
+
+        message = "The algorithm terminated without finding a valid bracket."
+        with pytest.raises(RuntimeError, match=message):
+            optimize.bracket(f, -1, 1)
+        with pytest.raises(RuntimeError, match=message):
+            optimize.bracket(f, -1, np.inf)
+        with pytest.raises(RuntimeError, match=message):
+            optimize.brent(f, brack=(-1, 1))
+        with pytest.raises(RuntimeError, match=message):
+            optimize.golden(f, brack=(-1, 1))
+
+        def f(x):  # gh-5899
+            return -5 * x**5 + 4 * x**4 - 12 * x**3 + 11 * x**2 - 2 * x + 1
+
+        message = "No valid bracket was found before the iteration limit..."
+        with pytest.raises(RuntimeError, match=message):
+            optimize.bracket(f, -0.5, 0.5, maxiter=10)
+
+    @pytest.mark.parametrize('method', ('brent', 'golden'))
+    def test_minimize_scalar_success_false(self, method):
+        # Check that status information from `bracket` gets to minimize_scalar
+        def f(x):  # gh-14858
+            return x**2 if ((-1 < x) & (x < 1)) else 100.0
+
+        message = "The algorithm terminated without finding a valid bracket."
+
+        res = optimize.minimize_scalar(f, bracket=(-1, 1), method=method)
+        assert not res.success
+        assert message in res.message
+        assert res.nfev == 3
+        assert res.nit == 0
+        assert res.fun == 100
+
+
 def test_brent_negative_tolerance():
     assert_raises(ValueError, optimize.brent, np.cos, tol=-.01)
 
diff --git a/tools/unicode-check.py b/tools/unicode-check.py
index 2a1e5b09f3fa..4a223163126d 100755
--- a/tools/unicode-check.py
+++ b/tools/unicode-check.py
@@ -11,7 +11,7 @@
 # allow in the source code.
 latin1_letters = set(chr(cp) for cp in range(192, 256))
 box_drawing_chars = set(chr(cp) for cp in range(0x2500, 0x2580))
-extra_symbols = set(['®', 'ő', 'λ', 'π', 'ω', '∫', '≠', '≥', 'μ'])
+extra_symbols = set(['®', 'ő', 'λ', 'π', 'ω', '∫', '≠', '≥', '≤', 'μ'])
 allowed = latin1_letters | box_drawing_chars | extra_symbols
 
 


Code review comments:
diff_hunk: @@ -2928,9 +2941,17 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
     >>> plt.plot(xc, fc, "bx")
     >>> plt.show()
 
+    Note that both initial points were to the right of the minimum, and the
+    third point was found in the "downhill" direction: the direction
+    in which the function appeared to be decreasing (to the left).
+    The final points are strictly ordered, and the function value
+    at the middle point is less than the function values at the endpoints;
+    it follows that a minimum must lie within the bracket.
+
     """
     _gold = 1.618034  # golden ratio: (1.0+sqrt(5.0))/2.0
     _verysmall_num = 1e-21
+    xa, xb = np.asarray([xa, xb])
path: scipy/optimize/_optimize.py
text: @user1:
Do we need this line as `xa` and `xb` are advertised as being floats?

@author:
~~I can guess why I did it, but if I don't remember and I didn't test it, it must not have been important. I'll remove it.~~ My guess was right. This makes `xa` and `xb` NumPy floats (without explicitly specifying their type in case they are already NumPy floats). The example in gh-5899 used to result in an `OverflowError`; now it finds that the function is unbounded below. I think I didn't test it because this was intended (originally) to be a quick bug fix, not a complete overhaul.

@user1:
Ah ok, neat trick

diff_hunk: @@ -19,7 +19,8 @@
                         _minimize_bfgs, _minimize_newtoncg,
                         _minimize_scalar_brent, _minimize_scalar_bounded,
                         _minimize_scalar_golden, MemoizeJac, OptimizeResult,
-                        _wrap_callback)
+                        _wrap_callback, BracketError,
+                        _recover_from_bracket_error)
path: scipy/optimize/_minimize.py
text: @author:
```suggestion
                        _wrap_callback, _recover_from_bracket_error)
```
We can commit this with `[skip ci]` when this is ready to be merged.

diff_hunk: @@ -2994,9 +3018,59 @@ def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
         fa = fb
         fb = fc
         fc = fw
+
+    # three conditions for a valid bracket
+    cond1 = (fb < fc and fb <= fa) or (fb < fa and fb <= fc)
+    cond2 = (xa < xb < xc or xc < xb < xa)
+    cond3 = np.isfinite(xa) and np.isfinite(xb) and np.isfinite(xc)
+    msg = ("The algorithm terminated without finding a valid bracket. "
+           "Consider trying different initial points.")
+    if not (cond1 and cond2 and cond3):
+        e = BracketError(msg)
+        e.data = (xa, xb, xc, fa, fb, fc, funcalls)
+        raise e
+
     return xa, xb, xc, fa, fb, fc, funcalls
 
 
+class BracketError(RuntimeError):
+    pass
+
+
+def _recover_from_bracket_error(solver, fun, bracket, args, **options):
+    # `bracket` was originally written without checking whether the resulting
+    # bracket is valid. `brent` and `golden` built on top of it without
+    # checking the returned bracket for validity, and their output can be
+    # incorrect without warning/error if the original bracket is invalid.
+    # gh-14858 noticed the problem, and the following is the desired
+    # behavior:
+    # - `scipy.optimize.bracket`, `scipy.optimize.brent`, and
+    #   `scipy.optimize.golden` should raise an error if the bracket is
+    #   invalid, as opposed to silently returning garbage
+    # - `scipy.optimize.minimize_scalar` should return with `success=False`
+    #   and other information
+    # The changes that would be required to achieve this the traditional
+    # way (`return`ing all the required information from bracket all the way
+    # up to `minimizer_scalar`) are extensive and invasive. (See a6aa40d.)
+    # We can achieve the same thing by raising the error in `bracket`, but
+    # storing the information needed by `minimize_scalar` in the error object,
+    # and intercepting it here.
+    try:
+        res = solver(fun, bracket, args, **options)
+    except BracketError as e:
+        msg = str(e)
+        xa, xb, xc, fa, fb, fc, funcalls = e.data
+        xs, fs = [xa, xb, xc], [fa, fb, fc]
+        if np.any(np.isnan([xs, fs])):
+            x, fun = np.nan, np.nan
+        else:
+            imin = np.argmin(xs)
path: scipy/optimize/_optimize.py
text: @user1:
Possible mistake?
```suggestion
            imin = np.argmin(fs)
```

@author:
Unfortunately, I'm having trouble coming up with a realistic test case that exercises this with `minimize_scalar`. Sure, we could simulate calling `_recover_from_backet_error` with an artificial `solve` argument that raises a `BracketError` with out-of-order `xs` and `fs` and show that this picks the right one, but I don't think it's worth it.