adding jitted scalar maximization routine, first build

[jstac]

Hi @mmcky, in this PR I'm adding an optimize subpackage and a first function within that subpackage called maximize_scalar.

There's an example of usage within the documentation for that function.

The code and algorithm come from SciPy's fminbound. I've stripped out some tests and jitted the function. We require numba version 0.38 and above because the routine takes a jitted function as its first argument (the scalar function to be maximized).

The intention is to add further scalar maximization routines and at least one multivariate routine, as well as perhaps root finders. Essentially, the structure will mimic scipy.optimize and the feature set will no doubt be smaller, but every function should be jitted and accept jitted functions to act upon.

The purpose of functions in this subpackage is not necessarily to maximize speed within the optimization routine itself, but rather to ensure that we can jit compile outer functions like Bellman operators, which include a maximization routine as one step of the algorithm.

Thanks to @natashawatkins for helping me pull this together.

@sglyon @albop @cc7768 @chrishyland FYI

[coveralls]

Coverage decreased (-0.07%) to 95.107% when pulling 7f24600 on add_optim_subpackage into 649cc55 on master.

[mmcky]

this will be set in setup.py
I will update this tomorrow morning.

[mmcky]

thanks @jstac this look great - I will review more thoroughly tomorrow morning. We can open issues to track which other optimize routines we wish to add in and tag them for contributions. This will be a source of good projects.

[oyamad]

Very nice!

[oyamad]

It will be helpful to add *args (to pass to func).

[natashawatkins]

@oyamad thanks I have added this, although unfortunately you need to be quite careful to pass in a tuple and not a scalar, as I can't figure out how to check the type inside a jitted function

[oyamad]

Sorry I don't see what you mean. Can you elaborate?

And any difference between adding args=() and *args?

[natashawatkins]

In scipy, if something in args is passed that is not a tuple (ie. a scalar), the function will convert it to a tuple. I don't seem to be able to get isinstance to work inside a jitted function. If you only want to set one fixed argument, you need to pass args=(y,) which is somewhat annoying.

I guess we could use *args - I was just following scipy's style.

[oyamad]

I see thanks.

Following exactly scipy's style makes sense, while *args looks more Pythonic, allowing passing e.g. y=5. I can't say which is better...

[oyamad]

I think a status flag (ierr in scipy.fminbound) should also be returned, which in scipy.fminbound is 0 if converged and 1 if maxiter is reached.

[oyamad]

The minimum version number in

QuantEcon.py/setup.py

Line 116 in 649cc55

'numba>=0.36.2',

should be increased to 0.38.

[QBatista]

Did you mean "maximizer"?

[natashawatkins]

Thanks @QBatista !

[mmcky]

It looks like numba=0.38 is the default version in the latest Anaconda 5.2 distribution.
https://docs.anaconda.com/anaconda/packages/py3.6_linux-64. Continuum have changed the way packages are reported from what I can tell 0.38 is the default (rather than just an available package in the anaconda repository).

[albop]

This is very interesting!
If you give me a bit of time (end of this week), I'd like to look into it a bit.
In particular, I'm curious about how this function behaves in a loop, and whether the way function arguments are passed makes any difference.

[jstac]

@mmcky Great news. (BTW, lots of people are having problems with the library because they're not at the latest version of Numba. I guess there's only so much we can do on this front...)

@oyamad Thanks for your feedback.

@albop It would be great to have your input. @chrishyland and his friend will be working on these and subsequent additions close to full time in July, so they can build on your results.

@chrishyland Please loop your friend Paul into this discussion.

[jstac]

@albop Did you still want to look at this before we merge?

[albop]

yes @jstac, 2h more please ;-)

…

On Sun, Jun 24, 2018 at 10:50 PM John Stachurski ***@***.***> wrote: @albop <https://github.com/albop> Did you still want to look at this before we merge? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#416 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAQ5KYU4aZ6iaNfKCuhFTqmBRe_nVEJqks5t__uYgaJpZM4UciQO> .

[thomassargent30]

Dear John, maybe we should put a news note on the quantecon page noting the need to have the latest version of numba

…

On Mon, Jun 11, 2018, 4:01 PM John Stachurski ***@***.***> wrote: @mmcky <https://github.com/mmcky> Great news. (BTW, lots of people are having problems with the library because they're not at the latest version of Numba. I guess there's only so much we can do on this front...) @oyamad <https://github.com/oyamad> Thanks for your feedback. @albop <https://github.com/albop> It would be great to have your input. @chrishyland <https://github.com/chrishyland> and his friend will be working on these and subsequent additions close to full time in July, so they can build on your results. @chrishyland <https://github.com/chrishyland> Please loop your friend Paul into this discussion. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#416 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AD22TrqOBDyDq5O_LvvWUFMKUDnORJcyks5t7ngZgaJpZM4UciQO> .

[albop]

@jstac: I've read the code carefully and did some timings. It all sounds good to me. In particular, the overhead of using this function vs the manual inclusion of it's logic, seems small or in-existent in the examples I have tried.
I have only two remarks related to the API:

• The method name bears no reference to the algorithm that is used (brent I guess). But then if you want to include another algorithm, how would it be called?
• Even in a simplified optimization function, it seems risky to me to not return the convergence information (number of iterations, tolerances, etc.). There is often some "gridpoint that doesn't converge" and it's useful to know which and why. There would be at least two options:
  • return a jitted class with val, x0, iter etc... as fields
  • return a tuple val, x0, infos where infos would be a structure (or jitted class) with the informations. If the compiler optimizes across function, one could even hope the last one is not computed if you do val, x0, = maximize_scalar( )

[albop]

optimize subpackage is a very good idea.

[jstac]

Thanks all for your input.
@oyamad @albop I've added status_flag as a return value. In particular, the function now returns a tuple with three elements, the last of which is status_flag. It takes the value 0 if the function executes normally and 1 if the maximum number of iterations is attained.

@albop It's a good point about the name. Perhaps we can leave it as is and later add a method keyword argument that defaults to brent, if necessary.

[jstac]

@thomassargent30 That's a good point. A few people had problems with this at the recent workshops.

@natashawatkins Would you mind to add a news item when this is merged that emphasizes the need to update to the latest version of Numba, preferably by obtaining the latest version of Anaconda..

[albop]

@jstac : status_flag is really useful. I would still like to have some more information about convergence or even a trace. I would be perfectly happy if there was a way to return it some way, optionally, without too much performance cost (and I can do a PR).
method keyword is a good solution user-wise. I wonder whether the time of the string-check would be significant in an intensive loop, but it's not clear to me.
In a nutshell, I'm happy with the PR as it is.

[jstac]

Thanks @albop .

If we return the number of function calls that brings us up to parity with fminbound. I've implemented this via

fval, fx, info = maximize_scalar(f, a, b)

where info is a tuple containing status_flag and num_iter. I've modified the docstring accordingly.

[jstac]

I'm open to changing the name to, say, brent_max if you think it's better.

[albop]

Don't want to nitpick, but don't you want to return fval, fx, status, info like fminbound then ? The advantage is that for all the possible minimization routines, status, always has the same simple type, while info can hold more information and be method-specific ?
As for the method name, I like brent-max, but don't feel too strongly about it. It's really up to you. If there is to be more than one method, then it would be natural to make the general method an alias to the specific one.

[oyamad]

• Better to order the first two returns as xf, fval as in fminbound (instead of fval, xf)?
• Regarding the Numba version issue, see add minimum version number for numba support #409.

[jstac]

I switched the order of xf and fval as requested by @oyamad and changed the function name to brent_max. The return value is unchanged otherwise, but I think this simplifies the function call and gives us enough flexibility going forward.

Unless there are other comments I think this is ready to merge @mmcky .

[oyamad]

The optimize subpackage should be registered in docs/qe_apidoc.py; refer to c7c9fdf for example.

[mmcky]

@jstac I will add this to the docs.

[mmcky]

@jstac documentation is ready to go: https://quanteconpy.readthedocs.io/en/add_optim_subpackage/

[mmcky]

Thanks everyone for this PR. I will merge tonight if there are no other comments.

[jstac]

Hi @mmcky, thanks for updating the docs. Could you please merge this now?

[jstac]

Hi @mmcky, I suspect that you're otherwise occupied with more important things, so just this once I'll go ahead and merge. I hope that's OK.