Repository: scipy/scipy
Pull request: BUG: optimize: fix max function call validation for minimize with Nelder-Mead
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Fix #14165
Related #13921

#### What does this implement/fix?
<!--Please explain your changes.-->
As reported in #14165, minimize with Nelder-Mead method has an issue of max function call validation.
The reason for this issue is that _minimize_neldermead check the max function call here once an optimization loop:
https://github.com/scipy/scipy/blob/fbcfba22f2e4f89510f0bc1e24591a6b91204378/scipy/optimize/optimize.py#L759

However, the objective function is evaluated in many places in the function:
https://github.com/scipy/scipy/blob/fbcfba22f2e4f89510f0bc1e24591a6b91204378/scipy/optimize/optimize.py#L749-L750
https://github.com/scipy/scipy/blob/fbcfba22f2e4f89510f0bc1e24591a6b91204378/scipy/optimize/optimize.py#L768
https://github.com/scipy/scipy/blob/fbcfba22f2e4f89510f0bc1e24591a6b91204378/scipy/optimize/optimize.py#L775
https://github.com/scipy/scipy/blob/fbcfba22f2e4f89510f0bc1e24591a6b91204378/scipy/optimize/optimize.py#L793
https://github.com/scipy/scipy/blob/fbcfba22f2e4f89510f0bc1e24591a6b91204378/scipy/optimize/optimize.py#L805
https://github.com/scipy/scipy/blob/fbcfba22f2e4f89510f0bc1e24591a6b91204378/scipy/optimize/optimize.py#L818

To fix the issue, I used the wrap function _wrap_function_with_validation which counts function evaluation and raises an exception when function evaluation reaches the maximum value, which is created in #13921

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
scipy.optimize.minimize method='Nelder-Mead': 'maxfev' is not checked before initial simplex is evaluated
**Is your feature request related to a problem? Please describe.**
In the instance that you have N variables, and hence an initial simplex of size N+1, but request maxfev < N +1, the initial simplex is assessed before checking maxfev, hence maxfev is violated.

I realize this example makes no actual sense if you wanted to minimize something, but in my case I was attempting to implement some new features in the optimize loop, and hence only wanted one function evaluation carried out. Since Nealder-Mead is often used for expensive objective functions, having to wait for N+1 itterations is a bit annoying, and also a bit confusing! 

**Describe the solution you'd like**.
I think the behavior is fine, I'd maybe just suggest adding a warning in so it's clear what's happening, because this just took me quite a while to figure out. Then again perhaps my use case is unique enough that no change is needed. Just a thought!

Hints from issue discussion:

Commit message: BUG: optimize: add a test to check iteration number of nelder-mead optimization

Patch to review:
diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py
index 3ee2381d6a60..5b3e30c3e2d2 100644
--- a/scipy/optimize/optimize.py
+++ b/scipy/optimize/optimize.py
@@ -479,7 +479,7 @@ class _MaxFuncCallError(RuntimeError):
     pass
 
 
-def _wrap_scalar_function_with_validation(function, args, maxfun):
+def _wrap_scalar_function_maxfun_validation(function, args, maxfun):
     # wraps a minimizer function to count number of evaluations
     # and to easily provide an args kwd.
     ncalls = [0]
@@ -709,8 +709,6 @@ def _minimize_neldermead(func, x0, args=(), callback=None,
     maxfun = maxfev
     retall = return_all
 
-    fcalls, func = _wrap_scalar_function(func, args)
-
     x0 = asfarray(x0).flatten()
 
     if adaptive:
@@ -786,8 +784,17 @@ def _minimize_neldermead(func, x0, args=(), callback=None,
     one2np1 = list(range(1, N + 1))
     fsim = np.empty((N + 1,), float)
 
-    for k in range(N + 1):
-        fsim[k] = func(sim[k])
+    fcalls, func = _wrap_scalar_function_maxfun_validation(func, args, maxfun)
+
+    try:
+        for k in range(N + 1):
+            fsim[k] = func(sim[k])
+    except _MaxFuncCallError:
+        pass
+    finally:
+        ind = np.argsort(fsim)
+        sim = np.take(sim, ind, 0)
+        fsim = np.take(fsim, ind, 0)
 
     ind = np.argsort(fsim)
     fsim = np.take(fsim, ind, 0)
@@ -797,74 +804,78 @@ def _minimize_neldermead(func, x0, args=(), callback=None,
     iterations = 1
 
     while (fcalls[0] < maxfun and iterations < maxiter):
-        if (np.max(np.ravel(np.abs(sim[1:] - sim[0]))) <= xatol and
-                np.max(np.abs(fsim[0] - fsim[1:])) <= fatol):
-            break
-
-        xbar = np.add.reduce(sim[:-1], 0) / N
-        xr = (1 + rho) * xbar - rho * sim[-1]
-        if bounds is not None:
-            xr = np.clip(xr, lower_bound, upper_bound)
-        fxr = func(xr)
-        doshrink = 0
+        try:
+            if (np.max(np.ravel(np.abs(sim[1:] - sim[0]))) <= xatol and
+                    np.max(np.abs(fsim[0] - fsim[1:])) <= fatol):
+                break
 
-        if fxr < fsim[0]:
-            xe = (1 + rho * chi) * xbar - rho * chi * sim[-1]
+            xbar = np.add.reduce(sim[:-1], 0) / N
+            xr = (1 + rho) * xbar - rho * sim[-1]
             if bounds is not None:
-                xe = np.clip(xe, lower_bound, upper_bound)
-            fxe = func(xe)
-
-            if fxe < fxr:
-                sim[-1] = xe
-                fsim[-1] = fxe
-            else:
-                sim[-1] = xr
-                fsim[-1] = fxr
-        else:  # fsim[0] <= fxr
-            if fxr < fsim[-2]:
-                sim[-1] = xr
-                fsim[-1] = fxr
-            else:  # fxr >= fsim[-2]
-                # Perform contraction
-                if fxr < fsim[-1]:
-                    xc = (1 + psi * rho) * xbar - psi * rho * sim[-1]
-                    if bounds is not None:
-                        xc = np.clip(xc, lower_bound, upper_bound)
-                    fxc = func(xc)
-
-                    if fxc <= fxr:
-                        sim[-1] = xc
-                        fsim[-1] = fxc
-                    else:
-                        doshrink = 1
+                xr = np.clip(xr, lower_bound, upper_bound)
+            fxr = func(xr)
+            doshrink = 0
+
+            if fxr < fsim[0]:
+                xe = (1 + rho * chi) * xbar - rho * chi * sim[-1]
+                if bounds is not None:
+                    xe = np.clip(xe, lower_bound, upper_bound)
+                fxe = func(xe)
+
+                if fxe < fxr:
+                    sim[-1] = xe
+                    fsim[-1] = fxe
                 else:
-                    # Perform an inside contraction
-                    xcc = (1 - psi) * xbar + psi * sim[-1]
-                    if bounds is not None:
-                        xcc = np.clip(xcc, lower_bound, upper_bound)
-                    fxcc = func(xcc)
-
-                    if fxcc < fsim[-1]:
-                        sim[-1] = xcc
-                        fsim[-1] = fxcc
-                    else:
-                        doshrink = 1
+                    sim[-1] = xr
+                    fsim[-1] = fxr
+            else:  # fsim[0] <= fxr
+                if fxr < fsim[-2]:
+                    sim[-1] = xr
+                    fsim[-1] = fxr
+                else:  # fxr >= fsim[-2]
+                    # Perform contraction
+                    if fxr < fsim[-1]:
+                        xc = (1 + psi * rho) * xbar - psi * rho * sim[-1]
+                        if bounds is not None:
+                            xc = np.clip(xc, lower_bound, upper_bound)
+                        fxc = func(xc)
 
-                if doshrink:
-                    for j in one2np1:
-                        sim[j] = sim[0] + sigma * (sim[j] - sim[0])
+                        if fxc <= fxr:
+                            sim[-1] = xc
+                            fsim[-1] = fxc
+                        else:
+                            doshrink = 1
+                    else:
+                        # Perform an inside contraction
+                        xcc = (1 - psi) * xbar + psi * sim[-1]
                         if bounds is not None:
-                            sim[j] = np.clip(sim[j], lower_bound, upper_bound)
-                        fsim[j] = func(sim[j])
+                            xcc = np.clip(xcc, lower_bound, upper_bound)
+                        fxcc = func(xcc)
 
-        ind = np.argsort(fsim)
-        sim = np.take(sim, ind, 0)
-        fsim = np.take(fsim, ind, 0)
-        if callback is not None:
-            callback(sim[0])
-        iterations += 1
-        if retall:
-            allvecs.append(sim[0])
+                        if fxcc < fsim[-1]:
+                            sim[-1] = xcc
+                            fsim[-1] = fxcc
+                        else:
+                            doshrink = 1
+
+                    if doshrink:
+                        for j in one2np1:
+                            sim[j] = sim[0] + sigma * (sim[j] - sim[0])
+                            if bounds is not None:
+                                sim[j] = np.clip(
+                                    sim[j], lower_bound, upper_bound)
+                            fsim[j] = func(sim[j])
+            iterations += 1
+        except _MaxFuncCallError:
+            pass
+        finally:
+            ind = np.argsort(fsim)
+            sim = np.take(sim, ind, 0)
+            fsim = np.take(fsim, ind, 0)
+            if callback is not None:
+                callback(sim[0])
+            if retall:
+                allvecs.append(sim[0])
 
     x = sim[0]
     fval = np.min(fsim)
@@ -3117,7 +3128,7 @@ def _minimize_powell(func, x0, args=(), callback=None, bounds=None,
 
     # we need to use a mutable object here that we can update in the
     # wrapper function
-    fcalls, func = _wrap_scalar_function_with_validation(func, args, maxfun)
+    fcalls, func = _wrap_scalar_function_maxfun_validation(func, args, maxfun)
 
     if direc is None:
         direc = eye(N, dtype=float)
diff --git a/scipy/optimize/tests/test_optimize.py b/scipy/optimize/tests/test_optimize.py
index 9f2f6154df60..d359adb047ee 100644
--- a/scipy/optimize/tests/test_optimize.py
+++ b/scipy/optimize/tests/test_optimize.py
@@ -533,7 +533,8 @@ def cost(x):
         return rng.random(1) * 1000  # never converged problem
 
     for imaxfev in [1, 10, 50]:
-        for method in ['Powell']:  # TODO: extend to more methods
+        # TODO: extend to more methods
+        for method in ['Powell', 'Nelder-Mead']:
             result = optimize.minimize(cost, rng.random(10),
                                        method=method,
                                        options={'maxfev': imaxfev})
@@ -546,7 +547,7 @@ def func_(x):
         return x
 
     fcalls, func = optimize.optimize.\
-        _wrap_scalar_function_with_validation(func_, np.asarray(1), 5)
+        _wrap_scalar_function_maxfun_validation(func_, np.asarray(1), 5)
 
     for i in range(5):
         func(np.asarray(i))
@@ -557,7 +558,7 @@ def func_(x):
         func(np.asarray(i))  # exceeded maximum function call
 
     fcalls, func = optimize.optimize.\
-        _wrap_scalar_function_with_validation(func_, np.asarray(1), 5)
+        _wrap_scalar_function_maxfun_validation(func_, np.asarray(1), 5)
 
     msg = "The user-provided objective function must return a scalar value."
     with assert_raises(ValueError, match=msg):
@@ -571,6 +572,14 @@ def test_obj_func_returns_scalar():
     with assert_raises(ValueError, match=match):
         optimize.minimize(lambda x: x, np.array([1, 1]), method='BFGS')
 
+
+def test_neldermead_iteration_num():
+    x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
+    res = optimize._minimize._minimize_neldermead(optimize.rosen, x0,
+                                                  xatol=1e-8)
+    assert res.nit == 339
+
+
 def test_neldermead_xatol_fatol():
     # gh4484
     # test we can call with fatol, xatol specified


Code review comments:
diff_hunk: @@ -786,8 +784,17 @@ def _minimize_neldermead(func, x0, args=(), callback=None,
     one2np1 = list(range(1, N + 1))
     fsim = np.empty((N + 1,), float)
path: scipy/optimize/optimize.py
text: @user1:
This should probably be:
`fsim = np.full((N + 1,), np.inf, dtype=float)`

@user1:
Otherwise we'll get undefined behaviour if the initial simplex evaluation doesn't run to completion.

diff_hunk: @@ -571,6 +572,14 @@ def test_obj_func_returns_scalar():
     with assert_raises(ValueError, match=match):
         optimize.minimize(lambda x: x, np.array([1, 1]), method='BFGS')
 
+
+def test_neldermead_iteration_num():
+    x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
+    res = optimize._minimize._minimize_neldermead(optimize.rosen, x0,
+                                                  xatol=1e-8)
+    assert res.nit == 339
path: scipy/optimize/tests/test_optimize.py
text: @user1:
Hmm, I don't see why we have to hard code the number of iterations in a test.

@author:
I'm on neutral (ref: https://github.com/scipy/scipy/pull/14805#discussion_r722177925), what do you think? @user2 

@user1:
@author, if you change this to `assert res.nit <= 339`, I'd be happy for it to be merged.

@user2:
Indeed, I agree with Andrew here. Would make the test a bit more robust.

@author:
Done.

@user2:
Thanks!

@user1, I will let you merge then as you wanted to do it 😃 (or squash with putting yourself as co-author IMHO)