minimise error because not minimum found

[Suavesito-Olimpiada]

While working on a project I found that minimise throwed an error because I didn't find a minimum. The MWE is the next

julia> using IntervalArithmetic, IntervalOptimisation
julia> cost(s, a) = 1 + (s^2 + a^2) * abs(s-a) + abs(a)
cost (generic function with 2 methods)
julia> x = -1..1
[-1, 1]
julia> minimise(X -> cost(X[1], X[2]), x×x)
ArgumentError: reducing over an empty collection is not allowed

IntervalOptimisation.jl version is v0.4.3, versioninfo() is

Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 4

[Suavesito-Olimpiada]

Same happens in master.

[Suavesito-Olimpiada]

I think is something numerical happening here, minimise suceeds if I move one axis a little bit, so mid(x×x) is not exactly [0,0] (the minimum). The code is here

julia> using IntervalArithmetic, IntervalOptimisation
julia> cost(s, a) = 1 + (s^2 + a^2) * abs(s-a) + abs(a)
cost (generic function with 2 methods)
julia> x = -1..1
[-1, 1]
julia> minimise(X -> cost(X[1]-1e-10, X[2]), x×x)
([1, 1.00001], IntervalBox{2, Float64}[[-0.00131321, -0.000359587] × [-0.000359588, 0.000564693], [0.000564692, 0.00150353] × [-0.000359588, 0.000564693], [-0.000359588, 0.000564693]²])

[dpsanders]

Ah interesting. It should be easy to work out what's going on by debugging or printing out exactly what's being calculated, then.

[Suavesito-Olimpiada]

Yeah, I been doing that. If I come up with a solution I'll post a PR. :)

[Suavesito-Olimpiada]

The easiest way is to change the src/optimise.jl:40 from

m = sup(f(Interval.(mid.(X))))   # evaluate at midpoint of current interval

to

m = sup(f(Interval.(mid.(X, 0.45))))   # evaluate at midpoint of current interval

but there is not needed it to be 0.45, it can be changed to any other value.

The most robust solution is to change the filter_elements! so it uses <= intead of < for HeapedVector and cutoff instead of cutoff-1 for SortedArray.

Which one you think is better? I have PR's for both.

cc @dpsanders

[Suavesito-Olimpiada]

The problem surged originally because the partition that was computed had the same Interval's as result, and the correct one was thrown away because both seemed to have the correct global_minimum (1.0). As I understand this cannot be solved as such because we have repeated variables.

An even smaller example that fails is the next

julia> using IntervalArithmetic, IntervalOptimisation
julia> cost(s, a) = abs(s-a) + abs(a)
cost (generic function with 1 methods)
julia> x = -1..1
[-1, 1]
julia> minimise(X -> cost(X[1], X[2]), x×x)
ArgumentError: reducing over an empty collection is not allowed

[dpsanders]

Sorry I still don't see what exactly is going on.

[Suavesito-Olimpiada]

The problem arise because of several things. First the (unique) real global_min is hit in the first run of the while loop in minimise.

In the second pass of the loop, it happens that the bisections X1 and X2 return the same Interval and the true minimum as X_min. When they are push!ed into the HeapedVector (or SortedVector) they are ordered in a way that makes that the interval not containing the true minimum is popfirst! in the nex run of the loop.

In the third run of the loop the interval not containing the true global_min is popfirst! (so that working.data just contains one element, the interval with the true min). When filter_elemts! happens, because the algorithm has the real global_min AND filter_elements! use < instead of <=, the interval with the true global minimum is thrown away (similar for SortedVector).

Finally the algorithm continues but never finds an interval containing the true global minimum with the desired acuracy. When the loop finished minimizers is empty and minimum fails.

The robust solution is to use <= instead of < so if the algorithms hits the real global_minimum (as floating-point precision permits) it will not throw the interval containing the true minimum.

[Suavesito-Olimpiada]

The problem is that this solution can make worst the cases where you have repeated variables in the expression. But I think is better to have this instead of crashing.