Opened 12 years ago
Closed 7 years ago
#2607 closed defect (fixed)
find_minimum_on_interval() uses the wrong scipy function
Reported by: | AlexGhitza | Owned by: | jwmerrill |
---|---|---|---|
Priority: | major | Milestone: | sage-5.3 |
Component: | calculus | Keywords: | sd31, sd40.5 |
Cc: | robert.marik, kcrisman | Merged in: | sage-5.3.beta0 |
Authors: | Dan Drake, Andrey Novoseltsev, Andrzej Giniewicz, Volker Braun | Reviewers: | Karl-Dieter Crisman, Mike Hansen, Andrey Novoseltsev |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #13109 | Stopgaps: |
Description (last modified by )
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:
Attachments (5)
Change History (74)
comment:1 Changed 11 years ago by
- Component changed from calculus to numerical
comment:2 Changed 11 years ago by
- Component changed from numerical to calculus
comment:3 Changed 9 years ago by
- Report Upstream set to N/A
comment:4 Changed 9 years ago by
- Cc robert.marik added
comment:5 Changed 9 years ago by
Robert, that sounds like a great question for either sage-devel or the scipy list.
comment:6 Changed 9 years ago by
- Cc kcrisman added
comment:7 Changed 8 years ago by
- Keywords sd31 added
This Scipy tutorial page should be relevant. I will try to resolve this soon.
comment:8 Changed 8 years ago by
#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)
comment:9 Changed 8 years ago by
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.
comment:10 Changed 8 years ago by
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?
comment:11 Changed 8 years ago by
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/2*pi, rather than 3/2*pi, 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.
comment:12 Changed 8 years ago by
- Owner changed from was to jwmerrill
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.
comment:13 Changed 8 years ago by
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.
comment:14 Changed 8 years ago by
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...
comment:15 Changed 7 years ago by
I just ran into this bug with the following input:
find_maximum_on_interval(fast_float(8*e^(-x)*sin(x) - 1, x), 0, 8) (1.5791755355586754, 0.78539817769603915) find_maximum_on_interval(fast_float(8*e^(-x)*sin(x) - 1, x), 0, 9) (-0.9951835373923219, 7.0685835435476418)
...and was truly surprised that find_maximum_on_interval
is not monotonic (in the sense that a bigger interval should always give a (weakly) bigger maximum)!
At the VERY LEAST, we should fix the documentation to specify that this finds local extrema, and perhaps change the name of the function, too, since it does not always find the actual maximum value on the interval!
Note that one strange workaround is to simply plot the function; something like:
def find_maximum_on_interval(f, a, b): return plot(f, a, b).ymax()
seems like it would be pretty effective, despite being inelegant and crude!
comment:16 Changed 7 years ago by
- Keywords sd40.5 added
- Status changed from new to needs_review
Patch up, please review. I have just changed the documentation and added some suggestions for a workaround using plot(...).ymin()
. I think we should at least merge something like this right away and worry about fixing the algorithm later.
comment:17 Changed 7 years ago by
- Reviewers set to Karl-Dieter Crisman
- Status changed from needs_review to needs_work
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.
comment:18 Changed 7 years ago by
- Status changed from needs_work to needs_review
I fixed the versions of this function in symbolic.expression
to specify that they call the numerical.optimize
version. Please check the documentation; I looked at the bare html, but am working remotely and haven't viewed the result.
comment:19 Changed 7 years ago by
The whitespace patch removes trailing whitespace (and tabs! TABS!) from the relevant source files.
comment:20 Changed 7 years ago by
- Description modified (diff)
comment:21 Changed 7 years ago by
Needs work for sentence that doesn't end and :trac:
thing...
comment:22 Changed 7 years ago by
Add a "..." to fix doctest on OS X.
comment:23 follow-up: ↓ 31 Changed 7 years ago by
should trac2607.2.patch
be removed?
Paul
comment:24 Changed 7 years ago by
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
comment:25 follow-up: ↓ 26 Changed 7 years ago by
The patch doesn't change the behaviour, it only warns the user that it only returns some local minimum.
comment:26 in reply to: ↑ 25 Changed 7 years ago by
comment:27 Changed 7 years ago by
@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.
comment:28 Changed 7 years ago by
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.
comment:29 Changed 7 years ago by
Karl-Dieter, you mean the "stopgap" mechanism?
Paul
comment:30 Changed 7 years ago by
I think so. I've not yet used it, so I don't know if it would be appropriate here.
comment:31 in reply to: ↑ 23 Changed 7 years ago by
Replying to zimmerma:
should
trac2607.2.patch
be removed?
Yes. Or simply ignored.
Replying to zimmerma:
then I don't get any warning
The warning is in the documentation, which isn't the best...but it's better than right now, where it's impossible to figure out without reading a lot of code. There will, of course, be a followup ticket that actually fixes the functionality.
comment:32 follow-up: ↓ 33 Changed 7 years ago by
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
comment:33 in reply to: ↑ 32 Changed 7 years ago by
Replying to zimmerma:
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.
Fair enough. However, no one is offerring a fix to the code; it looks like we are going to change the documentation, and then later fix the actual function. This will require two tickets. We could use this ticket for the documentation and a new ticket for the code, or vice versa. We can change the description of this ticket if we want.
The exact way we do this doesn't seem very important to me, as long as we fix the documentation and later fix the code. If you would like to make a new ticket and move attachment:trac2607.patch and attachment:trac2607-whitespace.patch over, that's fine.
comment:34 Changed 7 years ago by
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.
comment:35 Changed 7 years ago by
Here's a relevant sage-support thread, showing some other problems with these functions: https://groups.google.com/forum/#!topic/sage-support/KCjW5QlB_sA
comment:36 Changed 7 years ago by
patchbot: apply trac2607.patch trac2607-whitespace.patch
comment:37 Changed 7 years ago by
- Status changed from needs_review to needs_work
- Work issues set to deprecation
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!
comment:38 Changed 7 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
- Work issues deprecation deleted
Apply trac2607.patch trac2607-whitespace.patch trac_2607_renaming.patch
comment:39 Changed 7 years ago by
According to the patchbot, these patches introduce some trailing whitespace :-) Otherwise, it looks good to me.
comment:40 Changed 7 years ago by
Given how many whitespaces were removed, I think these are fine to go...
comment:41 Changed 7 years ago by
- Reviewers changed from Karl-Dieter Crisman to Karl-Dieter Crisman, Mike Hansen
- Status changed from needs_review to positive_review
Okay, sounds good to me. Everything looks good.
comment:42 Changed 7 years ago by
- Status changed from positive_review to needs_work
Why 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.
comment:43 Changed 7 years ago by
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.)
comment:44 follow-up: ↓ 45 Changed 7 years ago by
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.
comment:45 in reply to: ↑ 44 Changed 7 years ago by
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?
comment:46 Changed 7 years ago by
- Dependencies set to #12950
In any case, this should be rebased to #12950.
comment:47 Changed 7 years ago by
- Status changed from needs_work to needs_review
The new whitespace patch now only affects the relevant functions and removes a couple tabs from expression.pyx
. It's far less invasive. All three relevant patches apply to 5.1.beta2.
comment:48 Changed 7 years ago by
- Dependencies changed from #12950 to sage-5.1.beta2
comment:49 Changed 7 years ago by
- Status changed from needs_review to positive_review
comment:50 Changed 7 years ago by
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.
comment:51 Changed 7 years ago by
- Description modified (diff)
Makes sense to me, new patch is attached. Please review!
comment:52 Changed 7 years ago by
- Status changed from positive_review to needs_work
comment:53 Changed 7 years ago by
- Status changed from needs_work to needs_review
comment:54 Changed 7 years ago by
- Status changed from needs_review to positive_review
I'd say it looks good now, applies cleanly, all tests passed, and new name is way easier to remember.
comment:55 Changed 7 years ago by
trac2607-whitespace.patch needs a proper commit message.
comment:56 Changed 7 years ago by
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 **********************************************************************
comment:57 Changed 7 years ago by
- Status changed from positive_review to needs_work
Changed 7 years ago by
patch trac_2607_renaming.2.patch rebased on top of trac2607-whitespace.2.patch
comment:58 Changed 7 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
I added the commit message to whitespace patch and rebased renaming patch based on it (keeping original authors/dates, I did not wanted to claim their work and it was minor rebase only to accommodate the commit message taken from comment by Dan).
I also uploaded trac2607-doctest-and-spacing.patch, which adjusts doctest for numerical noise, fixes spacing to be more consistent (3 places, spaces between arguments - i.e. "0,5" looks very close to "0.5" instead of "0, 5") and changed "f.find_local_minimum(" to "find_local_minimum(f" - because after all this test occurs in definition of find_local_minimum function, not method.
comment:59 Changed 7 years ago by
- Reviewers changed from Karl-Dieter Crisman, Mike Hansen to Karl-Dieter Crisman, Mike Hansen, Andrey Novoseltsev
- Status changed from needs_review to positive_review
Looks good and all tests pass.
comment:60 Changed 7 years ago by
- Milestone changed from sage-5.1 to sage-5.2
comment:61 Changed 7 years ago by
- Dependencies changed from sage-5.1.beta2 to #13109
Added patch to switch the deprecation to the new syntax.
comment:62 Changed 7 years ago by
- Description modified (diff)
- Status changed from positive_review to needs_work
comment:63 Changed 7 years ago by
- Status changed from needs_work to needs_review
comment:64 Changed 7 years ago by
- Status changed from needs_review to positive_review
Deprecation changes look good to me.
comment:65 Changed 7 years ago by
(and all tests still pass)
comment:66 Changed 7 years ago by
- Milestone changed from sage-5.2 to sage-pending
comment:67 Changed 7 years ago by
- Milestone changed from sage-pending to sage-5.2
comment:68 Changed 7 years ago by
- Milestone changed from sage-5.2 to sage-5.3
comment:69 Changed 7 years ago by
- Merged in set to sage-5.3.beta0
- Resolution set to fixed
- Status changed from positive_review to closed
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?