Nelder-Mead stops prematurely if function returns Inf

[MaximilianWollner]

Issue
Optimizing a function which returns Inf at one of the vertices causes the Nelder-Mead solver to stop prematurely without giving a proper reason for its termination. I did not encounter this issue in previous versions of Optim, where Inf was a valid function output.

Example

using Optim, LinearAlgebra

function test(x)
    if norm(x) > 1
        return Inf
    else
        return sum(x)
    end
end

res = optimize(x -> begin; y = test(x); println(y); y; end, [0.7, 0.3], NelderMead(), Optim.Options(show_trace=true))

foreach((key, value) -> println("$key: $value"), keys(res.stopped_by), values(res.stopped_by))
println("Termination code: $(res.termination_code)")

Versions
Julia v1.12.4
Optim v2.0.0

[pkofod]

Interesting, thanks. Will try to look into it.

[pkofod]

SInce you were looking into it already. What is printed in v1.12.4?

[MaximilianWollner]

Here is the output, slightly adjusted for better readability:

Function output at [0.7, 0.3]: 1.0
Function output at [1.0749999999999997, 0.3]: Inf
Function output at [0.7, 0.475]: 1.1749999999999998
Function output at [0.7, 0.38749999999999996]: 1.0875

Iter     Function value    √(Σ(yᵢ-ȳ)²)/n 
------   --------------    --------------
     0     1.000000e+00              NaN
 * time: 1.5974044799804688e-5
 
x_converged: false
f_converged: false
g_converged: false
f_limit_reached: false
g_limit_reached: false
h_limit_reached: false
time_limit: false
callback: false
f_increased: false
ls_failed: false
iterations: false

Termination code: NotImplemented

The optimization returns:

Status: failure

 * Candidate solution
    Final objective value:     1.000000e+00

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    0
    f(x) calls:    4

I am fairly confident I have encountered Infs in the first iteration in previous versions of Optim without the solver stopping prematurely . Setting the penalty to a large number, e.g., 1e10 instead of Inf, leads to the expected behavior.

[pkofod]

Hm, it seems to me that (as you say I think) the failure indeed happened on the last version of v1.x.y prior to v2. In any event, we can try to do better here just like for line searches where we backtrack until finite values are observed instead of panicking at alpha=1 => phi(alpha)=Inf

[MaximilianWollner]

As far as I understand, the Nelder-Mead algorithm never makes use of the values of the objective function in any algebraic operation, but only compares said values, see Nelder & Mead (1965, Fig. 1):

This might be a bit naive, but would it be enough to change the strict ineqalities to greater or equal? The issue seems to be that any newly constructed point could potentially also return Inf, which is actually okay as long as we end up with a contraction assuming the starting point returns a finite value.

[pkofod]

how will you guarantee that the contracted value is finite? I suppose you are envisioning a case where going "too far" away from x0 puts you into a region where the returned value is inf, but as long as you're close enough it's OK? In practice, these Inf regions can surely be isolated regions of "unfortunate" parameter combinations. So you can, in some cases, even contract into such a region. However, my reference above to backtracking indeed referred to the possibility that we could add some kind of "get out of bad regions" heuristic that contracts the point if this happens. Alternatively, the user can input a full starting simplex already although the interface is a bit funny.

[MaximilianWollner]

As far as I can tell, we do not need to guarantee that the contracted values are finite. I was thinking about the following scenarios assuming that the starting point returns a finite value:

1. If all $y_i$ are finite, we are good.
2. If $y_l$ is finite and any other $y_i$ is infinite, then $y_h$ is infinite:
   1. In case $y^\ast$ is infinite, then $y^\ast \geq y_i \forall i$ and we contract:
      1. If $y^{\ast\ast}$ is finite, we are good.
      2. If $y^{\ast\ast}$ is infinite, then $y^{\ast\ast} \geq y_h = \infty$ and we move all $P_i$ toward $P_l$ which returned a finite value.
   2. In case $y^\ast$ is finite, then $y^\ast < y_h = \infty$ and we are good.
3. In case all $y_i$ return infinite values, one could shrink the whole simplex towards the starting point until a vertex returns a finite number. Then, we are in case 1 or 2.

Maybe I am overlooking something, but case 3 should eventually lead to case 1 or 2. And iterating case 2 should get rid of all remaining vertices with infinite values.

As long as we are returning Infs, the algorithm is neither reflecting nor expanding, only contracting towards a ball around the starting point, in which all vertices return finite values.

To implement this approach, one would only need to replace the three rightmost strict inequalities with $\geq$ since Inf $\geq$ Inf returns 'true'.

[pkofod]

Are you thinking about this as working for the starting points only? What if this happens during the optimization?

[MaximilianWollner]

As far as I can tell, this should work throughout the optimization, not only for the initial simplex.