IntervalOptimisation surprising result

Specific Domains Optimization (Mathematical)
1

I need to find the maximiser and the maximum of functions that have a narrow bump. I had some issues with numerics, so I decided to try my luck with intervals.

Here is a minimal working (?) example:

using IntervalArithmetic, IntervalOptimisation
f(x) = exp(1/x)/(x^2*(1 + exp(1/x))^2)
maximise(f, 0..∞)

This results to

([0.439228, ∞], Interval{Float64}[[0, 0.000655297], [0.000655296, 0.00132092], [0.00132091, 0.00198653], [0.415145, 0.415866], [0.414436, 0.415146], [0.413728, 0.414437], [0.412322, 0.413031], [0.417947, 0.418656], [0.41303, 0.413729], [0.416551, 0.41725]  …  [0.389124, 0.389822], [0.442803, 0.443524], [0.388426, 0.389125], [0.443523, 0.444255], [0.387739, 0.388427], [0.444254, 0.444964], [0.387031, 0.38774], [0.444963, 0.445683], [0.445682, 0.446403], [0.386333, 0.387032]])

The interval maximum that it gives has the actual maximum as lower bound, but the upper bound is infinity. However, the main problem is that it gives me 88 maximasers.

For reference the function looks like that
Capture

I also tried to invert the argument, hoping that this would make things better

using IntervalArithmetic, IntervalOptimisation
f(x) = 2*exp(x)*x^2/(1 + exp(x))^2
maximise(f, 0..∞)

but I waited 10 minutes and killed the process.

Does anyone have any insight on what’s happening here and how I can get the maximiser?

EDIT: I realised that but restricting the upper bound of the starting interval it gives me the maximum. So I have calculated a crude upper bound. The problem is that I get a huge list of maximisers (it grows longer with lower tolerance). Here’s a minimal working example:

using IntervalArithmetic, IntervalOptimisation
f(x) = exp(x)*x^2/(1 + exp(x))^2
maximise(f, 0..4)

This gives me

([0.439228, 0.440029], Interval{Float64}[[2.39909, 2.4001], [2.3981, 2.3991], [2.39614, 2.39714], [2.40401, 2.40501], [2.39223, 2.39323], [2.40597, 2.40698], [2.40697, 2.40797], [2.39713, 2.39811], [2.40205, 2.40304], [2.39516, 2.39615]  …  [2.31803, 2.31903], [2.45531, 2.45582], [2.31705, 2.31804], [2.48043, 2.48136], [2.45581, 2.45633], [2.31607, 2.31706], [2.48135, 2.48229], [2.45928, 2.45979], [2.45978, 2.4603], [2.48228, 2.48322]])

EDIT: Calling wrong package fixed.

2

It probably stalls because of dependency. I can’t test your MWE (edit: needed to load IntervalOptimisation.jl), but if the cosh function is defined in the interval arithmetic library, you can try:

g(x) = x^2/(1 + cosh(x))

edit: we get

julia> maximise(g, 0..4, tol=1e-6)
([0.878457, 0.878459], Interval{Float64}[[2.39934, 2.39935], [2.39934, 2.39935], [2.39933, 2.39934], [2.39933, 2.39934], [2.39932, 2.39933], [2.39931, 2.39932], [2.3994, 2.39941], [2.39928, 2.39929], [2.39943, 2.39944], [2.39945, 2.39946]  …  [2.39755, 2.39756], [2.39755, 2.39756], [2.39756, 2.39757], [2.39754, 2.39756], [2.39762, 2.39763], [2.39757, 2.39758], [2.39761, 2.39762], [2.39761, 2.39762], [2.39756, 2.39757], [2.40065, 2.40066]])

We do get a lot of intervals, but if you take the union, you get a pretty good enclosure of the optimal solution.

2 Likes
3

Thanks for the link and the suggestion about using cosh. I will try it out.

I was surprised by the number of intervals. In some cases, they don’t contain the maximum!

I got a working solution by taking the union of the intervals and passing it back to maximise with lower tolerance. Far from optimal, but thankfully this code needs to run only once.

1 Like
4

Great.
Can you post the example where the union of intervals doesn’t contain the minimizer? That would be a severe bug (and it defeats the purpose of using interval arithmetic).

5

Take a look at my first minimal working example. The maximiser of the function is approximately 0.41678, but the first interval is [0, 0.000655297].

EDIT: I guess I should have led with that :slight_smile:

6

Another approach is to use ApproxFun:

julia> using ApproxFun

julia> f(x) = x<0.008 ? zero(x) : 2*exp(1/x)/(x^2*(1 + exp(1/x))^2)
f (generic function with 1 method)

julia> fun = Fun(f, 0..100)
Fun(Chebyshev(0..100),[0.03434119405627297, -0.06697215513709168, 0.06259849738569075, -0.05623338233565475, 0.04846593045117235, -0.03980574884125622, 0.03068712490678202, -0.021473860801793494, 0.012464556120606804, -0.003898171386952791  …  -6.0202358506331044e-15, 5.316927453868914e-15, -4.644722106927901e-15, 3.999892864842147e-15, -3.381572365873664e-15, 2.7849901204635152e-15, -2.2079777242667298e-15, 1.6418073497947994e-15, -1.0888642418271921e-15, 5.423179266772493e-16])

julia> findmax(fun)
(0.8784576797813086, 0.41677827980041116)

Note: In your original post, you left out the factor of 2 in defining f, but it was included in the plotted function.

7

Ohhh I see what you mean. That’s not how it works :wink:
See the documentation:

They return an Interval that is guaranteed to contain the global
minimum (maximum), and a Vector of Intervals or IntervalBoxes
whose union contains all the minimisers.
8

Ah, I see. Thanks for clarifying.

By the way I figured out that it is faster to use the derivative and find its root.

1 Like
9

Thanks for the suggestion. I haven’t used ApproxFun before. I will take a look, but for now I solved my problem by finding the root of the derivative. :slight_smile:

10

Definitely! If I remember correctly, IntervalOptimisation.jl is a 0th-order method, that is it evaluates interval ranges but doesn’t use derivatives (yet). So it can be quite slow and overestimation due to dependency decreases linearly with the size of the interval.
In some cases, setting the gradient to 0 is the way to go :slight_smile:

2 Likes