Repository: scipy/scipy
Pull request: MAINT: optimize.minimize_scalar: enforce output shape consistency
#### Reference issue
Closes gh-16196

#### What does this implement/fix?
gh-16196 reported that the output shape of `minimize_scalar`'s result object `x` attribute was consistent when an objective function returned an array. This PR ensures that `x` attribute is always a NumPy scalar.

#### Additional information
The user is only supposed to pass in a 1D array for `bounds`/`bracket`, so ensuring 0D output would not be wrong. But we should vectorize this function in the future, and while we're fixing this bug, it would better for the policy be that the shape of `res.x` (the values of `x` that zero all outputs of `f(x)`) should match that of `res.fun == f(res.x)`, which is controlled by the user. 

This will also be less disruptive to users. Those who use `brent` and `bounded` have already accommodated for `res.x` being an array if `f` returns an array, e.g. `res.x[0]`. If we were to enforce that `res.x` should be a scalar, they will begin to get an error. Users of `golden` will begin getting the array that they should have been getting all along, but size 1 NumPy arrays tend to behave just like scalars in computations in most cases. And this way, we don't change things _again_ when we want to vectorize things.


Problem statement:
BUG: OptimizeResult from optimize.minimize_scalar changes 'x' type if convergence is found in first step.
### Describe your issue.

Executing an optimization of the same function through `optimize.minimize_scalar`, but with different bounds or tolerance values, is causing the type of the `x` attribute of the `OptimizeResult` object to change between executions.

In this example, the function returns after the first iteration. The `x` value is a scalar.

 ```
 Func-count     x          f(x)          Procedure
    1     3.2482e-09  -3.8548e+08        initial

Optimization terminated successfully;
The returned value satisfies the termination criteria
(using xtol =  1e-05 )
     fun: array([-3.8547964e+08])
 message: 'Solution found.'
    nfev: 1
     nit: 1
  status: 0
 success: True
       x: 3.2481988817183283e-09
```

In this example, the function is called with an explicit tolerance value that forces it to iterate further, and returns after the 8th iteration. The `x` value is now a numpy array (as it should be, according to the documentation).

``` 
 Func-count     x          f(x)          Procedure
    1     3.2482e-09  -3.8548e+08        initial
    2    4.34094e-09 -1.14571e+07        golden
    3    2.57285e-09      -543583        golden
    4    3.46288e-09 -9.51544e+08        parabolic
    5    3.74413e-09 -1.55055e+09        parabolic
    6    3.97209e-09 -9.13931e+08        golden
    7    3.71079e-09 -1.53441e+09        parabolic
    8    3.77746e-09 -1.53008e+09        parabolic

Optimization terminated successfully;
The returned value satisfies the termination criteria
(using xtol =  1e-10 )
     fun: array([-1.55054734e+09])
 message: 'Solution found.'
    nfev: 8
     nit: 8
  status: 0
 success: True
       x: array([3.74412799e-09])
```
I understand that multiple optimization approaches in SciPy have different output formats, but the same optimizer should not change output types depending on iterations.

Looking at the code, I could not pinpoint the location where this `x` value changes to an nd array, but it consistently returns as a scalar if only one iteration is performed before termination. 



### Reproducing Code Example

```python
import numpy as np
from scipy import stats, optimize

testArray = np.array([3.98870996e-09, 3.98826001e-09, 3.67428003e-09, 3.76826001e-09,
       3.61183008e-09, 3.77692993e-09, 3.95031006e-09, 2.96020996e-09,
       3.35727002e-09, 3.52889990e-09, 3.38438989e-09, 3.95433008e-09,
       3.37065991e-09, 3.19607007e-09, 3.86244995e-09, 3.35210010e-09,
       3.60525000e-09, 3.87412988e-09, 3.79013989e-09, 3.75905005e-09,
       3.89471997e-09, 3.40269995e-09, 3.72817993e-09, 3.72587988e-09,
       3.89586011e-09, 3.75660010e-09, 3.77033008e-09, 3.41007007e-09,
       3.39072998e-09, 4.01754004e-09, 3.57187012e-09, 3.65805005e-09,
       3.54221997e-09, 3.86550000e-09, 3.43247998e-09, 2.99565991e-09,
       3.63393994e-09, 3.74633008e-09, 3.83310010e-09, 3.70584009e-09,
       3.64812988e-09, 4.07268994e-09, 3.72121997e-09, 3.58945996e-09,
       3.82662988e-09, 3.98991992e-09, 3.66976001e-09, 3.42454004e-09,
       3.54378003e-09, 3.76230005e-09])

kdeModelForces = stats.gaussian_kde(testArray, bw_method="silverman")

result = optimize.minimize_scalar(lambda x: -kdeModelForces(x), 
                                  bounds=[testArray.min()*0.5,
                                          testArray.max()*1.5], 
                                  method="Bounded", 
                                  options={'disp':3}, 
                                  tol=None)

print(result)

print(" ------- ")

result = optimize.minimize_scalar(lambda x: -kdeModelForces(x), 
                                  bounds=[testArray.min()*0.5,
                                          testArray.max()*1.5], 
                                  method="Bounded",
                                  options={'disp':3},
                                  tol=1.0e-10)

print(result)
```


### Error message

```shell
---
```


### SciPy/NumPy/Python version information

`1.8.0 1.22.3 sys.version_info(major=3, minor=9, micro=12, releaselevel='final', serial=0)`

Hints from issue discussion:
@seberg I came up with many ways to convert a python scalar, NumPy scalar, or size 1 NumPy array to a NumPy scalar of the appropriate `dtype`. The simplest is probably:
```python3
x = np.atleast_1d(x)[0]
```
or
```python3
x = np.squeeze(x)[()]
```
but these take 10x as long on my machine as
```python3
x = np.asarray(x).squeeze()[()]
```
What is the preferred way?

---
Related: I was surprised by this:
```python3
np.float32(np.array([1])).shape  # (1,)
np.float64(np.array([1])).shape  # ()
```

Commit message: MAINT: optimize.minimize_scalar: always return `x` as NumPy scalar

Patch to review:
diff --git a/scipy/optimize/_minimize.py b/scipy/optimize/_minimize.py
index 563fc2cfbea8..5b62b25769d9 100644
--- a/scipy/optimize/_minimize.py
+++ b/scipy/optimize/_minimize.py
@@ -921,19 +921,22 @@ def minimize_scalar(fun, bracket=None, bounds=None, args=(),
         options['disp'] = 2 * int(disp)
 
     if meth == '_custom':
-        return method(fun, args=args, bracket=bracket, bounds=bounds, **options)
+        res = method(fun, args=args, bracket=bracket, bounds=bounds, **options)
     elif meth == 'brent':
-        return _minimize_scalar_brent(fun, bracket, args, **options)
+        res = _minimize_scalar_brent(fun, bracket, args, **options)
     elif meth == 'bounded':
         if bounds is None:
             raise ValueError('The `bounds` parameter is mandatory for '
                              'method `bounded`.')
-        return _minimize_scalar_bounded(fun, bounds, args, **options)
+        res = _minimize_scalar_bounded(fun, bounds, args, **options)
     elif meth == 'golden':
-        return _minimize_scalar_golden(fun, bracket, args, **options)
+        res = _minimize_scalar_golden(fun, bracket, args, **options)
     else:
         raise ValueError('Unknown solver %s' % method)
 
+    res.x = np.asarray(res.x).squeeze()[()]
+    return res
+
 
 def _remove_from_bounds(bounds, i_fixed):
     """Removes fixed variables from a `Bounds` instance"""
diff --git a/scipy/optimize/tests/test_optimize.py b/scipy/optimize/tests/test_optimize.py
index 4422a636f422..e5963037f039 100644
--- a/scipy/optimize/tests/test_optimize.py
+++ b/scipy/optimize/tests/test_optimize.py
@@ -1878,6 +1878,25 @@ def test_minimize_unbounded_method_with_bounds_gh10911(self, method):
         with pytest.raises(ValueError, match=msg):
             optimize.minimize_scalar(np.sin, method=method, bounds=(1, 2))
 
+    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
+    @pytest.mark.parametrize("method", MINIMIZE_SCALAR_METHODS)
+    @pytest.mark.parametrize("tol", [1, 1e-6])
+    def test_minimize_scalar_output_dimensionality_gh16196(self, method, tol):
+        # gh-16196 reported that the output shape of `minimize_scalar` was not
+        # consistent when an objective function returned an array.
+        # Check that the output is always a NumPy scalar.
+        def f(x):
+            return np.array(x**4).reshape(1)
+
+        a, b = -0.1, 0.2
+        kwargs = (dict(bracket=(a, b)) if method != "bounded"
+                  else dict(bounds=(a, b)))
+        kwargs.update(dict(method=method, tol=tol))
+
+        res = optimize.minimize_scalar(f, **kwargs)
+        assert res.x.shape == tuple()
+        assert res.fun.shape == f(res.x).shape
+
 
 def test_brent_negative_tolerance():
     assert_raises(ValueError, optimize.brent, np.cos, tol=-.01)


Code review comments:
diff_hunk: @@ -921,19 +921,22 @@ def minimize_scalar(fun, bracket=None, bounds=None, args=(),
         options['disp'] = 2 * int(disp)
 
     if meth == '_custom':
-        return method(fun, args=args, bracket=bracket, bounds=bounds, **options)
+        res = method(fun, args=args, bracket=bracket, bounds=bounds, **options)
     elif meth == 'brent':
-        return _minimize_scalar_brent(fun, bracket, args, **options)
+        res = _minimize_scalar_brent(fun, bracket, args, **options)
     elif meth == 'bounded':
         if bounds is None:
             raise ValueError('The `bounds` parameter is mandatory for '
                              'method `bounded`.')
-        return _minimize_scalar_bounded(fun, bounds, args, **options)
+        res = _minimize_scalar_bounded(fun, bounds, args, **options)
     elif meth == 'golden':
-        return _minimize_scalar_golden(fun, bracket, args, **options)
+        res = _minimize_scalar_golden(fun, bracket, args, **options)
     else:
         raise ValueError('Unknown solver %s' % method)
 
+    res.x = np.asarray(res.x).squeeze()[()]
path: scipy/optimize/_minimize.py
text: @author:
The user is only supposed to pass in a 1D array for `bounds`/`bracket`, so ensuring 0D output would not be wrong. But we should vectorize this function in the future, and while we're fixing this bug, it would better for the policy be that the shape of `res.x` (the values that zero all outputs of `f`) should match that of `res.fun`, which is controlled by the user. 

This will also be less disruptive to users. Those who use `brent` and `bounded` have already accommodated for `res.x` being an array if `f` returns an array, e.g. `res.x[0]`. If we were to enforce that `res.x` should be a scalar, they will begin to get an error. Users of `golden` will begin getting the array that they should have been getting all along, but size 1 NumPy arrays tend to behave just like scalars in computations in most cases. And this way, we don't change things _again_ when we want to vectorize things.