icecube / pisa

Monte Carlo-based data analysis
http://icecube.github.io/pisa/
Apache License 2.0
19 stars 47 forks source link

Support all gradient-free algorithms from NLOPT #683

Closed atrettin closed 2 years ago

atrettin commented 2 years ago

Support for all gradient-free global and local optimization routines from the NLOPT package

With this PR I add a new option to the box of algorithms available in PISA that allows the user to use any of the global (https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#global-optimization) and local (https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#local-derivative-free-optimization) algorithms. The reason why I only chose to support the gradient free methods is that NLOPT doesn't apply a finite-difference method automatically and expects every function to supply its gradient itself. Most gradient based algorithms are also included in scipy anyway and for our purposes the gradient-free methods and global searches are more interesting.

NLOPT can be dropped in place of scipy and iminuit by writing a dictionary with "method": "nlopt" and choosing the algorithm by its name of the form NLOPT_{G,L}{N}_XXXX. For example, to use the Nelder-Mead algorithm, you would write

nlopt_settings = {
    "method": "nlopt",
    "method_kwargs": {
        "algorithm": "NLOPT_LN_COBYLA",
        "ftol_abs": 1e-5,
        "ftol_rel": 1e-5,
        # other options that can be set here: 
        # xtol_abs, xtol_rel, stopval, maxeval, maxtime
        # after maxtime seconds, stop and return best result so far
        "maxtime": 60
    },
    "local_fit_kwargs": None  # no further nesting available
}

and then run the fit with

best_fit_info = ana.fit_recursively(
    data_dist,
    dm,
    "chi2",
    None,
    **nlopt_settings
)

. Of course, you can also nest the nlopt_settings dictionary in any of the fit_octants, fit_ranges and so on by passing it as local_fit_kwargs.

Adding constraints

Adding inequality constraints to algorithms that support it is possible by writing a lambda function in a string that expects to get the current parameters as a ParamSet and returns a float. The result will satisfy that the passed function stays negative (to be consistent with scipy). The string will be passed to eval to build the callable function. For example, a silly way to bound delta_index > 0.1 would be:

    "method_kwargs": {
        "algorithm": "NLOPT_LN_COBYLA",
        "ftol_abs": 1e-5,
        "ftol_rel": 1e-5,
        "maxtime": 30,
        "ineq_constraints": [
            # be sure to convert parameters to their magnitude
            "lambda params: params.delta_index.m - 0.1"
        ]
    }

Adding inequality constraints to algorithms that don't support it can be done by either nesting the local fit in the constrained_fit method or to use NLOPT's AUGLAG method that adds a penalty for constraint violations internally. For example, we could do this to fulfill the same constraint with the PRAXIS algorithm:

"method_kwargs": {
    "algorithm": "NLOPT_AUGLAG",
    "ineq_constraints":[
        "lambda params: params.delta_index.m - 0.1"
    ],
    "local_optimizer": {
        # supports all the same options as above
        "algorithm": "NLOPT_LN_PRAXIS",
        "ftol_abs": 1e-5,
        "ftol_rel": 1e-5,
    }
}

Using global searches with local subsidiary minimizers

Some global searches, like evolutionary strategies, use local subsidiary minimizers. These can be defined just as above by passing a dictionary with the settings to the local_optimizer keyword. Note that, again, only gradient-free methods are supported. Here is an example for the "Multi-Level single linkage" (MLSL) algorithm, using PRAXIS as the local optimizer:

"method_kwargs": {
    "algorithm": "NLOPT_G_MLSL_LDS",
    "local_optimizer": {
        "algorithm": "NLOPT_LN_PRAXIS",
        "ftol_abs": 1e-5,
        "ftol_rel": 1e-5,
    }
}

For some evolutionary strategies such as ISRES, the population option can also be set.

"method_kwargs": {
    "algorithm": "NLOPT_GN_ISRES",
    "population": 100,
}
atrettin commented 2 years ago

I forgot to mention: Algorithm-specific parameters (https://nlopt.readthedocs.io/en/latest/NLopt_Reference/#algorithm-specific-parameters) are passed as a dictionary to the keyword algorithm_params.

philippeller commented 2 years ago

Could the above useful description be added to the code as docstring? Then this will be handy and also render nicely in the online pisa docs. Other than that this PR looks good to me 👍

atrettin commented 2 years ago

I added a lot of documentation and formatted the docstrings in such a way that they will actually be rendered nicely when building the docs as HTML.

philippeller commented 2 years ago

I am happy to merge this, it adds nlopt support and a lot of documentation, while not changing any of the existing code