OneDimensionOptimizationTools::bracketMinimum - clarification on line 93

[Keren Halabi]

Hi dear team,

I am encountering an issue when using BrentOneDimension for optimization of a function.

Say I have a parameter with the constraint [0.05,0.95].

When I call BrentOneDimension::doInit --> which in turn calls OneDimensionOptimizationTools::bracketMinimum with a=0.05 and b=0.95, the guess of the 3rd point in line 93 of OneDimensionOptimizationTools.cpp ends up being 2.406. which exceeds the bounds of my parameter.

I looked up golden ratio bracketing and the formula in line 93 didn't fully make sense to me:

line 93 is:

bracket.c.x = bracket.b.x + NumConstants::GOLDEN_RATIO_PHI() * (bracket.b.x - bracket.a.x);

Where NumConstants::GOLDEN_RATIO_PHI() = (1+sqrt(5))/2

This computation yields in my case:

bracket.c.x = 0.95 + (1+sqrt(5))/2*(0.95-0.05)=2.40623059 >> 0.95

However, when I follow equation (3) described here I get:

b-c=r(b-a) -> c=b-r(b-c)

where r= (-1+sqrt(5))/2

And then accordingly change line 93 to:

bracket.c.x = bracket.b.x - NumConstants::GOLDEN_RATIO_R() * (bracket.b.x - bracket.a.x);

I get:

bracket.c.x = 0.95 - (-1+sqrt(5))/2*(0.95-0.05)=0.39376941

The latter value of bracket.c.x is within the interval, thus not causing constraintException upon setting the value of the parameter.

I therefore don't fully understand something in the function bracketMinimum - maybe you could clarify?

Many thanks!

Keren

  

[Julien Y. Dutheil]

Dear Keren,

The algorithm used in this function is the one described in "Numerical Recipes in C". Maybe there are different variants. Anyway, whatever algorithm you use for optimizing, none of them are meant to be used with constraints (apart from the BFGS). They all assume a function with parameters defined on ]-inf, +inf[. There are therefore two possibilities:

- either you reparameterize your function. This is what you should do in theory, but in practice, it might slow down your optimization. Note that there is a ReparametrizationWrapper in bpp-core that proposes an automatic reparametrization that should work in the general case (but better reparametrization might be found depending on the function).

- or you use "AutoParameters" (this can be done automatically by using the setConstraintPolicy function of the optimizer), a dirty trick that consists in pushing the parameter values back to the closest valid bound. This is quite efficient and works in most cases, but can have you stuck in some local minima for more complex functions.

Hope this helps,

Julien.

[Keren Halabi]

Thank you for the quick reply, Julien,

Your answer is indeed helpful. I see that your bracketMinimum function follows the algorithm described in numerical recipes, page 352, for "look outward" bracketing. That is, it finds a third minimal point while expanding the initial interval provided to it, [a,b].

Just to make sure I follow your second suggestion - in the case of frequency parameter, which is constrained to be within the bounds [0,1],  whenever the assigned value is greater than 1, this solution would result in repeated computation of the function at value 1.

I saw another bracketing option in the book, which conducts "look inward" bracketing, that is, finds the third point within the interval [a,b] without modifying it. Could a bracketing achieved in this manner yield a more efficient local search with Brent, when using a constrained parameter? I have written an implementation of this bracketing and I could add it to OneDimensionOptimizationTools class, if you wish to add it to Bio++. Please let me know and I will pull request if needed.

If you think there is a better approach that "inward bracketing", I would love to hear your input.

Many thanks again!

Keren

[Julien Y. Dutheil]

Dear Keren,

Just to make sure I follow your second suggestion - in the case of frequency parameter, which is constrained to be within the bounds [0,1],  whenever the assigned value is greater than 1, this solution would result in repeated computation of the function at value 1.

Yes, which is not optimal for bracketing. Parameter transformation is a more elegant solution. 

I saw another bracketing option in the book, which conducts "look inward" bracketing, that is, finds the third point within the interval [a,b] without modifying it. Could a bracketing achieved in this manner yield a more efficient local search with Brent, when using a constrained parameter? I have written an implementation of this bracketing and I could add it to OneDimensionOptimizationTools class, if you wish to add it to Bio++. Please let me know and I will pull request if needed.

That could be an interesting option to have indeed, in particular for bounded parameters like frequencies. If you make a pull request I will look into it.

Cheers,

Julien.

[Keren Halabi]

Thank you Julien,

Before I pull request, I would like to consult on the best incorporation of this feature into bpp-core. 

Because I didn't want to modify any existing code written by the team, I have implemented an alternative class of BreneOneDimension that uses inward bracketing, but aside that is completely similar to the original BreneOneDimension class. A more sensible solution might be to add a non-mandatory argument to the function BrentOneDimension::doInit that dictates whether inward or outward bracketing should be used. Naturally, the default would be outward, as currently done. However, this approach would require that I modify and existing function written by you.

Which approach would you prefer? if there is another approach that didn't occur to me and that you would prefer, please let me know.

Many thanks!

Keren

[Julien Y. Dutheil]

Dear Keren,

Sorry for the late reply! I think it is fine to modify code that we did, git will sort that out :) But that's why it is also important to make specific branches for each purpose, to make things easier to check.

As for your suggestion: I think it would be better to let the doInit method  as is, and then have an inner option like

BrentOneDimention::BRACKET_OUTWARD and BrentOneDimention::BRACKET_INWARD (can be an arbitrary 'short' number, or an explicit 'string' if you prefer) with a corresponding bracketing_ field set to be equal to BRACKET_OUTWARD, for backward compatibility. Then you create two new methods setBracketing() and getBracketing() which permit modifying/assessing this tag. The doInit() function then just read the bracketing option. I believe this is the most backward compatible way in terms of library interface.

Cheers,

Julien.

[Keren Halabi]

Hi Julien,

It's not a problem. Thank you for the elaborated answer. I have pulled request (#19).

Thanks!

Keren