Function value convergence test for `brent_max`

[QBatista]

Without a function value convergence test, brent_max can fail on simple functions, for instance:

>>> import quantecon as qe
>>> import numpy as np
>>> from numba import njit
>>> @njit
... def f(x): return 1.
... 
>>> qe.optimize.brent_max(f, -np.inf, np.inf)
(nan, 1.0, (0, 1))

Adding a function value convergence test should solve this issue.

[chrishyland]

Hi @QBatista
Sorry for the slow response but I just finished most of my exams for the semester this week. Will have more time to help out.

Thank you for pointing that out. Just a few quick (naive) questions which would be helpful to enlighten me. What did you mean by function value convergence test? It would be helpful for you to direct me to some resources on this. Currently we checked whether if the iterations have gone on too long. In the cases of the above example, did you want the function to raise an exception instead?

Thank you.

[QBatista]

@chrishyland By function value convergence test, I mean a termination condition that depends on the sequence of function evaluations. An example of an optimization method that uses such a test is Nelder-Mead, where during each iteration, the algorithm is terminated if after evaluating the function at each of the current simplex points, the highest and lowest function values are sufficiently close. In the example I gave above, all the function evaluations would be the same so checking convergence of the sequence would allow termination at a local maximum.

After checking SciPy's source code for brent, it seems that they are not using this type of tests, so it might be that it is actually not a good idea or is not feasible for this particular method. Nevertheless, their implementation works on the example above:

>>> from scipy.optimize import minimize_scalar
>>> from numba import njit
>>> @njit
... def f(x): return 1
... 
>>> minimize_scalar(f)
     fun: 1
    nfev: 40
     nit: 36
 success: True
       x: 2.6180339603380443

Do you have an idea of what is causing this difference? If not, I can try to look into it.

[chrishyland]

Hi @QBatista
Thanks for the response and explanation. Yep, the brent_max was just copied over from Scipy and then jitted with some things removed. First guess for me would be that the get_bracket_info function wasn't ported over and could be a potential source of issue? Or it could just be the case there was a typo/logical mistake when porting the code over. I'll be able to work on it in my downtime throughout the week if you're busy.

On another note, I am somewhat confused on how for a function that returns any constant like in the example you gave, the solution is always 2.618... Why is this even the solution of the function?

[chrishyland]

Hi @QBatista ,

Just an update after playing around with the code this morning. An issue that was occurring was that the while loop to actually execute the parabolic interpolation/golden section search is not actually going through. In particular, on line while (np.abs(xf - xm) > (tol2 - 0.5 * (b - a))): is evaluating as False right away and hence not executing unlike the Scipy function which runs. On further inspection of this, it is because np.abs(xf - xm) and (tol2 - 0.5 * (b - a)) are both nan's unlike the Scipy function which gives back numbers like 0.309 and -1.309 respectively. I will have to look more into why are we getting nans tonight after work. Will investigate this further.

[chrishyland]

Hi @QBatista

Another update. The issue is the line
fulc = a + golden_mean * (b-a) near the top of the brent_max function. It's interesting because when you set the lower bound and upper bound a,b to infinity, the aforementioned line has issues. If you set the bounds to be finite numbers, then you start getting stuff that's kind of strange. When you change the bounds of brent_max, the maximizer returned increases as the bounds are relaxed (bounds are bigger). An example would be that

qe.optimize.brent_max(f, -100,100)
qe.optimize.brent_max(f, -200,200)

would return different results as the maximizer.

However for scipy's function, it consistently returns 2.618 as the optimal solution regardless of the bounds.

Any suggestions on how do you want the behavior of the function to look like? Should we be getting the same answer as the Scipy function?

[QBatista]

@chrishyland Thanks for looking into this! I think you're correct -- the difference is that SciPy's brent function uses automatic bracketing internally. When no bracket is passed, a bracket is automatically created using (0, 1) as a starting point. It's possible to replicate SciPy's brent function as follows:

qe.optimize.brent_max(f, optimize.bracket(f)[0], optimize.bracket(f)[2])

Also, since the function is constant, any x is a maximizer and minimizer so it's fine if the solution differs from SciPy as long as it's not nan.
The following actually fails and raises a runtime warning pointing to the line corresponding to the one you mentioned:

from scipy import optimize
optimize.fminbound(f, -np.inf, np.inf)

It makes sense that you're getting nan because when b=np.inf and a=-np.inf, a + golden_mean * (b-a) is indeterminate. However, the runtime error is not particularly helpful -- I think raising an error when a and b are not finite would be a better option. I guess we could also look into adding automatic bracketing if it's useful.

@jstac @natashawatkins Do you think a version of brent_max where the lower bound a and upper bound b are endogenously determined so that you are not required to pass them as arguments be useful?

[chrishyland]

Thank you for the explanation @QBatista . I was wondering about why those particular values were the maximizer for a function that is just a constant, your explanation cleared up alot.

I'm more than happy to keep working on resolving this issue, just let me know or if you have any other things you want me to take a look at, just let me know.