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

[34f90a52-c114-47db-a93b-3f83978622c0]

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.

[e0d4dd3a-58d6-47e6-b509-f9cf71bc0f0a]

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.

[e0d4dd3a-58d6-47e6-b509-f9cf71bc0f0a]

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

[56048686-9665-4c5e-973e-6c3add3aa805]

comment:25

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

[zimmermann6]

comment:26

Replying to @sagetrac-dsm:

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

then I don't get any warning (see comment [comment:24]).

Paul

[kcrisman]

comment:27

@sagetrac-dsm: Correct, and we HAVE to open a new ticket to get something better eventually. But that turns out to be hard with the current tools.

Hmm, maybe Paul is suggesting we should use the new "system warning that you won't get what you think/mathematically correct result" in Sage whose name I forget? Paul, the warning is just in the documentation, not the function itself.

I also have a tiny change to this so that the documentation looks better coming up.

[kcrisman]

comment:28

Actually, if we need to do more along the lines of warnings, I'll just put it here. The point is the tilde, as opposed to having the big long thing show up in the doc (while still telling people that it's in another place, which I imagine was the point of not having the tilde before).

-Uses :func:`sage.numerical.optimize.find_minimum_on_interval`
+Uses the :func:`~sage.numerical.optimize.find_minimum_on_interval`
+function in the numerical optimization module of Sage.

[zimmermann6]

comment:29

Karl-Dieter, you mean the "stopgap" mechanism?

Paul

[kcrisman]

comment:30

I think so. I've not yet used it, so I don't know if it would be appropriate here.

[zimmermann6]

comment:32

I'm not in favour of giving a positive review, since the proposed patch does not solve the problem described in the description of that ticket.

Paul

[benjaminfjones]

comment:34

It seems reasonable to me to fix the docs now, add stopgap, and open a new ticket to address the global / local optimization issue; especially given how old this ticket is. In this case the added documentation that refers to "See :trac:2607" should be updated to point to the new ticket.

[kini]

comment:36

patchbot: apply trac2607.patch trac2607-whitespace.patch

[novoselt]

comment:37

I talked to Dan and plan to write the following patch:

• rename this functions to clearly mark that they compute local extrema;
• keep old names with deprecation warning and explanation of their behaviour.

Will do it on top of the current patches, so no changes please!

[novoselt]

Work Issues: deprecation

[novoselt]

Description changed:

--- 
+++ 
@@ -10,4 +10,7 @@

 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 trac2607.patch and trac2607-whitespace.patch
+**Apply:**
+* [attachment: trac2607.patch](https://github.com/sagemath/sage-prod/files/10639961/trac2607.patch.gz)
+* [attachment: trac2607-whitespace.patch](https://github.com/sagemath/sage/files/ticket2607/trac2607-whitespace.patch.gz)
+* [attachment: trac_2607_renaming.patch](https://github.com/sagemath/sage/files/ticket2607/trac_2607_renaming.patch.gz)

[novoselt]

Changed author from Dan Drake to Dan Drake, Andrey Novoseltsev

[novoselt]

Changed work issues from deprecation to none

[novoselt]

comment:38

Apply trac2607.patch trac2607-whitespace.patch trac_2607_renaming.patch

[mwhansen]

comment:39

According to the patchbot, these patches introduce some trailing whitespace :-) Otherwise, it looks good to me.

[novoselt]

comment:40

Given how many whitespaces were removed, I think these are fine to go...

[mwhansen]

comment:41

Okay, sounds good to me. Everything looks good.

[mwhansen]

Changed reviewer from Karl-Dieter Crisman to Karl-Dieter Crisman, Mike Hansen

[jdemeyer]

comment:42

Why attachment: trac2607-whitespace.patch? These kind of patches are very annoying for rebasing. Removing whitespace from code which you change is fine (and even encouraged!), simply removing whitespace all over the place is bad.

[novoselt]

comment:43

Hi Jeroen, this patch removes some surprising TABs and I already rebased another patch on top of this whitespace removal, so maybe we can leave it... (Although in general I agree that meddling with whitespaces only complicates life and perhaps patchbot plugin was not such a great idea.)

[kini]

comment:44

What's wrong with the patchbot plugin? It only checks whether you added new trailing whitespace. I'm pretty sure everyone agrees you should at least try to avoid doing that.

[novoselt]

comment:45

Replying to @kini:

What's wrong with the patchbot plugin? It only checks whether you added new trailing whitespace. I'm pretty sure everyone agrees you should at least try to avoid doing that.

The problem is that some people get too eager to remove whitespace or insist on patchbot not showing any whitespace mistakes while others don't want to work on new patches if whitespaces are the only issue, and I cannot blame them. I am pretty sure that nobody adds whitespaces on purpose, but having some is not such a problem. Perhaps limiting the overall line length would be better. By the way - will the switch to git help somehow with whitespaces and patches that have conflicts because of them only?

[jdemeyer]

Dependencies: #12950

[jdemeyer]

comment:46

In any case, this should be rebased to #12950.

[novoselt]

Changed dependencies from #12950 to sage-5.1.beta2

[f8667504-6eb9-4678-970d-c85e19c6bba2]

comment:50

If in this ticket find_maximum_on_interval is renamed to find_local_maximum_on_interval, why not use the change to name it just find_local_maximum and keep it consistent with find_root? After all they both work same:

find_root(f, a, b)
find_local_maximum(f, a, b)

find_local_maximum_on_interval is 31 characters long, I think that such long names become hard to type. I know, it works on interval, but find_root also looks for root on [a,b] and isn't called find_root_on_interval.

[novoselt]

comment:51

Makes sense to me, new patch is attached. Please review!

[novoselt]

Description changed:

--- 
+++ 
@@ -13,4 +13,4 @@
 **Apply:**
 * [attachment: trac2607.patch](https://github.com/sagemath/sage-prod/files/10639961/trac2607.patch.gz)
 * [attachment: trac2607-whitespace.patch](https://github.com/sagemath/sage/files/ticket2607/trac2607-whitespace.patch.gz)
-* [attachment: trac_2607_renaming.patch](https://github.com/sagemath/sage/files/ticket2607/trac_2607_renaming.patch.gz)
+* [attachment: trac_2607_renaming.2.patch](https://github.com/sagemath/sage/files/ticket2607/trac_2607_renaming.2.patch.gz)

[f8667504-6eb9-4678-970d-c85e19c6bba2]

comment:54

I'd say it looks good now, applies cleanly, all tests passed, and new name is way easier to remember.

[jdemeyer]

comment:55

attachment: trac2607-whitespace.patch needs a proper commit message.

[jdemeyer]

comment:56

More importantly, this fails on OS X 10.4 ppc with numerical noise:

sage -t  --long -force_lib devel/sage/sage/numerical/optimize.py
**********************************************************************
File "/Users/buildbot/build/sage/moufang-1/moufang_full/build/sage-5.1.beta5/devel/sage-main/sage/numerical/optimize.py", line 135:
    sage: find_local_maximum(fast_float(8*e^(-x)*sin(x) - 1, x), 0, 8)
Expected:
    (1.5791755355586754, 0.78539817769603...)
Got:
    (1.5791755355586754, 0.7853981777050254)
GLPK Simplex Optimizer, v4.44
6 rows, 3 columns, 8 non-zeros
Preprocessing...
2 rows, 2 columns, 4 non-zeros
Scaling...
 A: min|aij| =  2.400e+01  max|aij| =  5.000e+01  ratio =  2.083e+00
GM: min|aij| =  8.128e-01  max|aij| =  1.230e+00  ratio =  1.514e+00
EQ: min|aij| =  6.606e-01  max|aij| =  1.000e+00  ratio =  1.514e+00
Constructing initial basis...
Size of triangular part = 2
*     0: obj =  -5.100000000e+01  infeas =  0.000e+00 (0)
*     1: obj =  -5.225000000e+01  infeas =  0.000e+00 (0)
OPTIMAL SOLUTION FOUND
**********************************************************************

[f8667504-6eb9-4678-970d-c85e19c6bba2]

patch trac2607-whitespace.patch with commit message

[f8667504-6eb9-4678-970d-c85e19c6bba2]

Attachment: trac2607-whitespace.2.patch.gz

patch trac_2607_renaming.2.patch rebased on top of trac2607-whitespace.2.patch