find_minimum_on_interval() uses the wrong scipy function

[aghitza]

This was reported by Dean Moore on sage-support. Consider:

sage: f(x) = -x*sin(x^2)
sage: f.find_minimum_on_interval(-2.5, 2)
(-1.3076194129914434, 1.35521114057)
sage: f.find_minimum_on_interval(-2.5, -1)
(-2.1827697846777219, -2.19450274985)

So find_minimum_on_interval() returns a local minimum as opposed to the global one. (The same issue applies to find_maximum_on_interval.) This is due to the fact that the function wraps scipy.optimize.fminbound, which is only a local optimizer (Carl Witty pointed this out). We should instead use one of the global optimizers, i.e. scipy.optimize.anneal or scipy.optimize.brute.

Apply:

• attachment: trac2607.patch
• attachment: trac2607-whitespace.2.patch
• attachment: trac_2607_renaming.3.patch
• attachment: trac2607-doctest-and-spacing.patch
• attachment: trac_2607_update_deprecation.patch

Depends on #13109

CC: @robert-marik @kcrisman

Component: calculus

Keywords: sd31, sd40.5

Author: Dan Drake, Andrey Novoseltsev, Andrzej Giniewicz, Volker Braun

Reviewer: Karl-Dieter Crisman, Mike Hansen, Andrey Novoseltsev

Merged: sage-5.3.beta0

Issue created by migration from https://trac.sagemath.org/ticket/2607

[robert-marik]

comment:3

It seems that there are two global optimizers in scipy: scipy.optimize.anneal and scipy.optimize.brute. Does anybody more experienced in numerics know, which one is better for including into Sage?

[jasongrout]

comment:5

Robert, that sounds like a great question for either sage-devel or the scipy list.

[kcrisman]

comment:7

This Scipy tutorial page should be relevant. I will try to resolve this soon.

[kcrisman]

Changed keywords from none to sd31

[kcrisman]

comment:8

#5960 was a dup. Here are the examples from there.

sage: h(x) =  -sin(x) - 2*sin(2*x)
sage: h.find_minimum_on_interval(0, 2*pi)
(-1.3271810224585345, 3.8298351449342838)
But there is another local minimum at h(0.8666760871050464) = -2.73581510406


sage: find_minimum_on_interval(x*(x-1)*(x+1), -2, 2)
(-0.38490017945975047, 0.57735026913115706)
The minimum on this interval is the endpoint h(-2) = 6.


sage: find_minimum_on_interval((x-2)*(x-1)*x*(x+1) - x, -2, 2)
(-0.43749999999999994, -0.49999999973911674)

but
sage: find_minimum_on_interval((x-2)*(x-1)*x*(x+1) - x, 0, 2)
(-2.6642135623730949, 1.7071067879138031)

[kcrisman]

comment:9

The brent algorithm will also not work.

Triple (a,b,c) where (a<b<c) and func(b) < func(a),func(c). If bracket consists of two numbers (a,c) then they are assumed to be a starting interval for a downhill bracket search (see bracket); it doesn’t always mean that the obtained solution will satisfy a<=x<=c.

Which is not the same as constrained minimization, for us.

[jicama]

comment:10

It seems like finding a local minimum might be the best you can hope for with a general function. Wouldn't finding an absolute minimum (on an interval) be intractable unless you can exploit some special structure of the function?

[kcrisman]

comment:11

Answering the question first...

Sure, but it would be nice if we at least got the 'right' answer for 'easy' functions. That's all I'm looking for here, not things like finding a minimum on a set of measure zero...

---

I think we can fix Brent to use this. Compare:

sage: optimize.fminbound(h._fast_float_(x),0,6,full_output=True)
(3.8298366870225147, -1.327181022449951, 0, 10)
sage: optimize.fminbound(h._fast_float_(x),0,3,full_output=True)
(0.86667541098916612, -2.7358151040622416, 0, 9)
sage: optimize.brent(h._fast_float_(x),brack=(0,6),full_output=True)
(0.86667608708813437, -2.73581510406422, 11, 12)

This shows that brent does give the 'right' answer in this case. So when does it give a 'wrong' answer?

sage: j(x) = sin(x)
sage: optimize.brent(j._fast_float_(x),brack=(0,6),full_output=True)
(10.995574367047061, -0.99999999999999689, 10, 11)
sage: 3.5*pi.n()
10.9955742875643

Well, of course - there IS no calculus-style minimum of sin between 0 and pi! Only a minimum relative to the interval itself. Interesting that it goes all the way to 7/2pi, rather than 3/2pi, but oh well!

So the fix is to switch to Brent, and then if it gives an answer outside the interval, pick the 'lower' endpoint. This would need lots of testing with well-behaved functions to make sure they actually work correctly.

[jicama]

comment:12

But how will you explain to users which functions are 'easy', and when they should expect to get the 'right' answer? I think it's better design to just change the contract of this function to admit that it is only looking for local minima.

[kcrisman]

comment:13

It doesn't matter. Or, at worst, we add some documentation to clarify that.

The reason it doesn't matter is that this is still better than the other. Unless you can produce some (natural) examples where optimize.brent does the same.

The Scipy documentation is not 100% clear on what is done, and it's conceivable they are the same. It's certainly possible that in fact using optimize.brent in the way I'm suggesting would be just as 'bad' as the previous one, or even essentially equivalent. But it would be nice to have an explicit example of this before we resort to that.

[kcrisman]

comment:14

Here's another example where brent does better.

sage: j(x) = sin(x^8)-.1*x
sage: optimize.brent(j._fast_float_(x),brack=(0,2),full_output=True)
(2.0000389609484905, -1.2000038913452364, 22, 23)
sage: optimize.fminbound(j._fast_float_(x),0,2,full_output=True)
(1.5288339777087034, -1.152883200877608, 0, 16)

Of course, this does cause problems for my supposed algorithm to then go back to an endpoint, since it's just outside of it...

[kcrisman]

comment:17

Suggestions in person about how to further enhance the messages in documentation. Thank you so much for doing this - don't forget to open a followup ticket.

[kcrisman]

Reviewer: Karl-Dieter Crisman

[kcrisman]

comment:21

Needs work for sentence that doesn't end and :trac: thing...

[zimmermann6]

comment:23

should trac2607.2.patch be removed?

Paul

[zimmermann6]

comment:24

after applying both patches on Sage 5.0, I still get:

sage: f(x) = -x*sin(x^2)                 
sage: f.find_minimum_on_interval(-2.5, 2)
(-1.3076194129914434, 1.3552111405712108)

Did I something wrong?

Paul

[sagetrac-dsm]

comment:25

The patch doesn't change the behaviour, it only warns the user that it only returns some local minimum.

[vbraun]

Changed dependencies from sage-5.1.beta2 to #13109

[vbraun]

comment:61

Added patch to switch the deprecation to the new syntax.

[vbraun]

Updated patch

[vbraun]

comment:62

Attachment: trac_2607_update_deprecation.patch.gz

[vbraun]

Changed author from Dan Drake, Andrey Novoseltsev, Andrzej Giniewicz to Dan Drake, Andrey Novoseltsev, Andrzej Giniewicz, Volker Braun

[jhpalmieri]

comment:64

Deprecation changes look good to me.

[sagetrac-aginiewicz]

comment:65

(and all tests still pass)

[sagetrac-aginiewicz]

comment:67

As I understand that ticket moved to sage-pending together with #13109 because it depends on it. But #13109 is back in 5.2 with positive review, so let's move this one back too.

[jdemeyer]

Merged: sage-5.3.beta0