Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add AGS global solver #194

Merged
merged 23 commits into from
Jul 26, 2018
Merged

Conversation

sovrasov
Copy link
Contributor

@sovrasov sovrasov commented Jul 15, 2018

Related to #192
This PR adds AGS global solver. Unlike most other global methods presented in NLOpt it can handle nonlinear inequality constraints.

@sovrasov
Copy link
Contributor Author

PR is done, but may be I missed some tricky integration moments.

@@ -31,7 +31,7 @@ Better yet, run some algorithm for a really long time until the minimum *f*<sub>
Global optimization
-------------------

All of the global-optimization algorithms currently require you to specify bound constraints on all the optimization parameters. Of these algorithms, only ISRES and ORIG_DIRECT support nonlinear inequality constraints, and only ISRES supports nonlinear equality constraints. (However, any of them can be applied to nonlinearly constrained problems by combining them with the [augmented Lagrangian method](#Augmented_Lagrangian_algorithm.md) below.)
All of the global-optimization algorithms currently require you to specify bound constraints on all the optimization parameters. Of these algorithms, only ISRES, AGS and ORIG_DIRECT support nonlinear inequality constraints, and only ISRES supports nonlinear equality constraints. (However, any of them can be applied to nonlinearly constrained problems by combining them with the [augmented Lagrangian method](#Augmented_Lagrangian_algorithm.md) below.)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add serial comma after AGS

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Didn't know that, used rules from Russian :)

AGS can handle arbitrary objectives and nonlinear inequality constraints. Also bound constraints are required for this method. To guarantee convergence, objectives and constraints should satisfy the Lipschitz condition on the specified hyperrectangle.
AGS is derivative-free and employs the Hilbert curve to reduce the source problem to the univariate one. The algorithm divides the univariate space into intervals, generating new points by using posterior probabilities. On each trial AGS tries to evaluate the constraints consequently one by one. If some constraint is violated at this point, the next ones won't be evaluated. If all constraints are preserved, i.e. the trial point is feasible, AGS will evaluate the objective. Thus, some of constraints (except the first one) and objective can be partially undefined inside the search hyperrectangle. Current implementation of AGS doesn't support vector constraints.

Limitations of the machine arithmetic don't allow to build a tight approximation for Hilbert when the space dimension is greater than 5, so this implementation of AGS is restricted in that sense.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If it is just slightly less reliable for dimensions > 5, it is probably worth allowing it to be used there so that people can experiment… or does the code not work at all in higher dimensions?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, increased the limit to 10, added a warning that sometimes the method can stop early for high-dimensional problems.

@@ -437,7 +437,8 @@ static int inequality_ok(nlopt_algorithm algorithm) {
|| AUGLAG_ALG(algorithm)
|| algorithm == NLOPT_GN_ISRES
|| algorithm == NLOPT_GN_ORIG_DIRECT
|| algorithm == NLOPT_GN_ORIG_DIRECT_L);
|| algorithm == NLOPT_GN_ORIG_DIRECT_L
|| algorithm == NLOPT_GN_AGS);
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

indenting

@stevengj
Copy link
Owner

Overall, this looks great, modulo minor comments, though I haven't tried out your code myself yet.

@jschueller
Copy link
Collaborator

I got these warnings compiling with -Wall, you might want to fix them:

[ 82%] Building CXX object CMakeFiles/nlopt.dir/src/algs/ags/ags.cc.o
/tmp/nlopt/src/algs/ags/evolvent.cc: In function ‘void mapd(double, int, double*, int, int)’:
/tmp/nlopt/src/algs/ags/evolvent.cc:142:5: warning: ‘is’ may be used uninitialized in this function [-Wmaybe-uninitialized]                                                                                                                 
     if (is == (nexp - 1)) i = -1;                                                                                                                                                                                                          
     ^~                                                                                                                                                                                                                                     
In file included from /tmp/nlopt/src/algs/ags/solver.hpp:10,                                                                                                                                                                                
                 from /tmp/nlopt/src/algs/ags/solver.cc:8:                                                                                                                                                                                  
/tmp/nlopt/src/algs/ags/solver.cc: In member function ‘void ags::NLPSolver::SetProblem(std::shared_ptr<ags::IGOProblem<double> >)’:                                                                                                         
/tmp/nlopt/src/algs/ags/solver.cc:76:54: warning: comparison of integer expressions of different signedness: ‘int’ and ‘unsigned int’ [-Wsign-compare]                                                                                      
   NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= solverMaxConstraints,                                                                                                                                                              
                     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~                                                                                                                                                               
/tmp/nlopt/src/algs/ags/data_types.hpp:14:43: note: in definition of macro ‘NLP_SOLVER_ASSERT’                                                                                                                                              
 #define NLP_SOLVER_ASSERT(expr, msg) if(!(expr)) NLP_SOLVER_ERROR(msg)                                                                                                                                                                     
                                           ^~~~                                                                                                                                                                                             
/tmp/nlopt/src/algs/ags/solver.cc: In member function ‘void ags::NLPSolver::SetProblem(const std::vector<std::function<double(const double*)> >&, const std::vector<double>&, const std::vector<double>&)’:                                 
/tmp/nlopt/src/algs/ags/solver.cc:88:54: warning: comparison of integer expressions of different signedness: ‘int’ and ‘unsigned int’ [-Wsign-compare]                                                                                      
   NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= solverMaxConstraints,                                                                                                                                                              
                     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~                                                                                                                                                               
/tmp/nlopt/src/algs/ags/data_types.hpp:14:43: note: in definition of macro ‘NLP_SOLVER_ASSERT’                                                                                                                                              
 #define NLP_SOLVER_ASSERT(expr, msg) if(!(expr)) NLP_SOLVER_ERROR(msg)                                                                                                                                                                     
                                           ^~~~                                                                                                                                                                                             
/tmp/nlopt/src/algs/ags/solver.cc: In member function ‘void ags::NLPSolver::EstimateOptimum()’:                                                                                                                                             
/tmp/nlopt/src/algs/ags/solver.cc:324:54: warning: suggest parentheses around ‘&&’ within ‘||’ [-Wparentheses]                                                                                                                              
         mOptimumEstimation.idx == mNextPoints[i].idx &&                                                                                                                                                                                    
         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~                                                                                                                                                                                    
         mOptimumEstimation.g[mOptimumEstimation.idx] > mNextPoints[i].g[mNextPoints[i].idx])                                                                                                                                               
         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                                                                                                                                                
/tmp/nlopt/src/algs/ags/solver.cc: In member function ‘void ags::NLPSolver::UpdateH(double, int)’:                                                                                                                                          
/tmp/nlopt/src/algs/ags/solver.cc:338:70: warning: suggest parentheses around ‘&&’ within ‘||’ [-Wparentheses]                                                                                                                              
   if (newValue > mHEstimations[index] || mHEstimations[index] == 1.0 && newValue > zeroHLevel)                                                                                                                                             
/tmp/nlopt/src/algs/ags/ags.cc: In function ‘int ags_minimize(unsigned int, nlopt_func, void*, unsigned int, nlopt_constraint*, double*, double*, const double*, const double*, nlopt_stopping*)’:                                          
/tmp/nlopt/src/algs/ags/ags.cc:89:22: warning: comparison of integer expressions of different signedness: ‘int’ and ‘unsigned int’ [-Wsign-compare]                                                                                         
     if (optPoint.idx != m)                                                                                                                                                                                                                 
         ~~~~~~~~~~~~~^~~~                                                                                                                                                                                                                  
/tmp/nlopt/src/algs/ags/ags.cc:94:9: warning: comparison of integer expressions of different signedness: ‘unsigned int’ and ‘int’ [-Wsign-compare]                                                                                          
   if (m == optPoint.idx)                                                                                                                                                                                                                   
       ~~^~~~~~~~~~~~~~~ 

@jschueller
Copy link
Collaborator

jschueller commented Jul 18, 2018

If I modify test/t_python.py to use AGS:

-opt = nlopt.opt(nlopt.LD_MMA, 2)
-opt.set_lower_bounds([-float('inf'), 0])
+opt = nlopt.opt(nlopt.GN_AGS, 2)
+opt.set_lower_bounds([-800, 0])
+opt.set_upper_bounds([800, 800])

I've got:

$ PYTHONPATH=$PWD/install/lib/python3.6/site-packages/ python3 ../test/t_python.py 
optimum at  0.1953125 400.09765625
minimum value =  20.002441257256574
result code =  5
nevals =  0

and same bounds with LD_MMA:

$ PYTHONPATH=$PWD/install/lib/python3.6/site-packages/ python3 ../test/t_python.py 
optimum at  0.333333374789 0.296296204114
minimum value =  0.5443309692765502
result code =  4
nevals =  11

Does that look ok?

@sovrasov
Copy link
Contributor Author

I've fixed the warnings.
About the test problem:
For global methods that assume Lipschitz continuity (like DIRECT and AGS) search domain should be chosen carefully. AGS acts on the unite hypercube and uses a runtime evaluation of Lipschitz constants. Constraints have approximately the same constant as x^3 on [-800, 800], and this is 3*800^2, really huge value. The method stucks because it uses one global estimation of L-constant for the whole search domain and expects very high variation speed for constraints.

If I reduce the bounds to

opt.set_lower_bounds([-8, 0])
opt.set_upper_bounds([8, 8])

output will be OK

optimum at  0.333984375 0.3134765625
minimum value =  0.559889777099
result code =  5
nevals =  46

The result is slightly imprecise because of the Hilbert curve, but usually the global solution is coarse and can be used as a starting point for a fast local optimizer.
Also I've fixed numevals for the objective, now output is nonzero.

Result of GN_ORIG_DIRECT with the restricted domain:

optimum at  0.33165675964 0.300157496317
minimum value =  0.547866312449
result code =  4
nevals =  0

Evals counter is seems to be broken.

@jschueller
Copy link
Collaborator

jschueller commented Jul 18, 2018

Yes, the eval counter is always 0 with AGS maybe your ++*stop->nevals_p is misplaced ?
Is it initialized to 0 in the common code at each optimization, then each solver must increment it each time the function is called.

@sovrasov
Copy link
Contributor Author

sovrasov commented Jul 18, 2018

I fixed the counter in the last commit by adding ++*stop->nevals_p but it's broken for DIRECT.
Although, from C counter always works.

@jschueller
Copy link
Collaborator

Yes, it's ok now.
What do you mean it's broken for DIRECT, using GN_DIRECT on the same test I've got nevals>0.

@jschueller
Copy link
Collaborator

tst.cc should be moved to /tests and added to the testsuite

@sovrasov
Copy link
Contributor Author

I've meant GN_ORIG_DIRECT.
tst.cc is rather sample than a test. May be rename it to example.cc? but I can add a standalone test executable.
In CMake some tests are disabled now. Should I rather add one more test executable than enable the existing common test?

@jschueller
Copy link
Collaborator

jschueller commented Jul 19, 2018

It's fine to add it to the common test as you did.

@jschueller
Copy link
Collaborator

I dont see AGS in inequality_ok function (src/api/options.c:~432), did you test inequality constraints ?

@sovrasov
Copy link
Contributor Author

It's on the line 441, after NLOPT_GN_ORIG_DIRECT_L:

|| algorithm == NLOPT_GN_ORIG_DIRECT_L
|| algorithm == NLOPT_GN_AGS);

@jschueller
Copy link
Collaborator

I looked at the wrong branch, sorry

@jschueller
Copy link
Collaborator

Why do you need c++11 again ? You're using c++, but you don't even seem to instanciate objects: all your functions could be marked as static and attributes are used as class attributes. I may be wrong but plain C might as well be used.

@sovrasov
Copy link
Contributor Author

Yes, one can use plain pointers in public interface and write sufficient wrappers for nlopt functions.
But solver itself uses a number of c++11 features: range-based loops, new methods, new constants and functions from STL, autos, and in-class initialization.
I can spend some time to downgrade my code, but in this case the code from my repo and from nlopt will diverge. I have plans to upgrade this implementation of AGS in future to make it converge faster on simple problems (like DIRECT-L faster than DIRECT). Support of 2 different versions of the same code is much more costly.
Also now on modern Linux distros it's hard to find a C++ compiler without C++11 support. On Windows one can have VS10 on Win10 and it's a problem.

@jschueller
Copy link
Collaborator

Ok, dont change your code, cxx11 is fine by me. It's just that it did not look very c++-ish at first glance, so I had to ask.

@sovrasov
Copy link
Contributor Author

What am I supposed to do next to get this merged?

@stevengj stevengj merged commit eab47e7 into stevengj:master Jul 26, 2018
@stevengj
Copy link
Owner

Whoops, I just noticed you are missing a license and readme file. I assume it is the same as in https://github.com/sovrasov/glob_search_nlp_solver and will copy the license file from there.

@sovrasov
Copy link
Contributor Author

Thanks for merge! I"ll add them soon in the next PR.

@sovrasov sovrasov deleted the vs/ags_global_nlp_solver branch July 26, 2018 14:39
stevengj added a commit that referenced this pull request Jul 26, 2018
@stevengj
Copy link
Owner

I just added them in 174ed18 … submit a PR to edit them as needed.

@sovrasov
Copy link
Contributor Author

Your version is ok for me.

@stevengj
Copy link
Owner

I just pushed a fix for the numevals counter in ORIG_DIRECT


ags::SolverParameters params;
params.r = ags_r;
params.itersLimit = stop->maxeval != 0 ? stop->maxeval : 5000;
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is wrong — maxeval == 0 denotes no maximum, not 5000 iterations. If you must set it to a specific value, use INT_MAX.

@stevengj
Copy link
Owner

@sovrasov, it doesn't look like you support the forced_stop or maxtime stopping criteria. Do you think you can submit a PR to use them?

For example, see how CRS does it:

nlopt/src/algs/crs/crs.c

Lines 135 to 138 in 24fc75f

if (nlopt_stop_forced(d->stop)) return NLOPT_FORCED_STOP;
if (d->p[0] < worst->k[0]) break;
if (nlopt_stop_evals(d->stop)) return NLOPT_MAXEVAL_REACHED;
if (nlopt_stop_time(d->stop)) return NLOPT_MAXTIME_REACHED;

forced_stop is especially important for languages like Python, because it is how we handle exceptions thrown by the objective function.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants