Numerical optimization

Using simple examples, this tutorial shows how to do an optimization with optimagic. More details on the topics covered here can be found in the how to guides.

import numpy as np
import pandas as pd

import optimagic as om

Basic usage of minimize

The basic usage of optimagic.minimize is very similar to scipy.optimize.minimize

def sphere(params):
    return params @ params

lbfgsb_res = om.minimize(
    fun=sphere,
    params=np.arange(5),
    algorithm="scipy_lbfgsb",
)

lbfgsb_res.params.round(5)

array([ 0., -0., -0., -0., -0.])

params do not have to be vectors

In optimagic, params can by arbitrary pytrees. Examples are (nested) dictionaries of numbers, arrays and pandas objects. This is very useful if you have many parameters!

def dict_sphere(params):
    return params["a"] ** 2 + params["b"] ** 2 + (params["c"] ** 2).sum()

nm_res = om.minimize(
    fun=dict_sphere,
    params={"a": 0, "b": 1, "c": pd.Series([2, 3, 4])},
    algorithm="scipy_neldermead",
)

nm_res.params

{'a': np.float64(3.885293415171424e-09),
 'b': np.float64(-9.144916295625622e-09),
 'c': 0   -6.348281e-09
 1   -4.599364e-09
 2   -2.724007e-09
 dtype: float64}

You can compare optimizers

In practice, it is super hard to pick the right optimizer for your problem. With optimagic, you can simply try a few and compare their results!

results = {"lbfgsb": lbfgsb_res, "nelder_mead": nm_res}
fig = om.criterion_plot(results, max_evaluations=300)
fig.show(renderer="png")

You can also zoom in on the history of specific parameters. This can be super helpful to diagnose problems in the optimization.

fig = om.params_plot(
    nm_res,
    max_evaluations=300,
    # optionally select a subset of parameters to plot
    selector=lambda params: params["c"],
)
fig.show(renderer="png")

There are many optimizers

By default, optimagic comes with optimizers from scipy, including global optimizers and least-squares optimizers. But we also have wrappers for algorithms from NlOpt, Pygmo, as well as several optimizers from individual packages like fides, ipopt, pybobyqa and dfols.

To use optimizers that are not from scipy, follow our installation guide for optional dependencies. To see which optimizers we have, check out the full list.

If you are missing your favorite optimizer in the list, let us know with an issue

Amazing autocomplete

Assume you need a gradient-free optimizer that supports bounds on the parameters. Moreover, you have a fixed computational budget, so you want to set stopping options.

In most optimizer libraries, you would have to spend a few minutes with the docs to find an optimizer that fits your needs and the stopping options it supports. In optimagic, all of this is discoverable in your editor!

If you type om.algos., your editor will show you all available optimizers and a list of categories you can use to filter the results. In our case, we select GradientFree and Bounded, and we could do that in any order we want.

After selecting one of the displayed algorithms, in our case scipy_neldermead, the editor shows all tuning parameters of that optimizer. If you start to type stopping, you will see all stopping criteria that are available.

Adding bounds

As any optimizer library, optimagic lets you specify bounds for the parameters.

bounds = om.Bounds(lower=np.arange(5) - 2, upper=np.array([10, 10, 10, np.inf, np.inf]))

res = om.minimize(
    fun=sphere,
    params=np.arange(5),
    algorithm="scipy_lbfgsb",
    bounds=bounds,
)

res.params.round(5)

array([0., 0., 0., 1., 2.])

Fixing parameters

On top of bounds, you can also fix one or more parameters during the optimization.

res = om.minimize(
    fun=sphere,
    params=np.arange(5),
    algorithm="scipy_lbfgsb",
    constraints=om.FixedConstraint(selector=lambda params: params[[1, 3]]),
)

res.params.round(5)

array([0., 1., 0., 3., 0.])

Other constraints

As an example, let’s impose the constraint that the first three parameters are valid probabilities, i.e. they are between zero and one and sum to one:

res = om.minimize(
    fun=sphere,
    params=np.array([0.1, 0.5, 0.4, 4, 5]),
    algorithm="scipy_lbfgsb",
    constraints=om.ProbabilityConstraint(selector=lambda params: params[:3]),
)

res.params.round(5)

array([ 0.33334,  0.33333,  0.33333, -0.     ,  0.     ])

For a full overview of the constraints we support and the corresponding syntaxes, check out the documentation.

Note that "scipy_lbfgsb" is not a constrained optimizer. If you want to know how we achieve this, check out the explanations.

There is also maximize

If you ever forgot to switch back the sign of your criterion function after doing a maximization with scipy.optimize.minimize, there is good news:

def upside_down_sphere(params):
    return -params @ params

res = om.maximize(
    fun=upside_down_sphere,
    params=np.arange(5),
    algorithm="scipy_bfgs",
)

res.params.round(5)

array([ 0., -0., -0.,  0., -0.])

optimagic got your back.

Speeding up your optimization with derivatives

You can speed up your optimization by providing closed form derivatives. Those derivatives can be hand-coded or calculated with JAX!

def sphere_gradient(params):
    return 2 * params


res = om.minimize(
    fun=sphere,
    params=np.arange(5),
    algorithm="scipy_lbfgsb",
    jac=sphere_gradient,
)
res.params.round(5)

array([ 0., -0., -0., -0., -0.])

Alternatively, you can let optimagic calculate numerical derivatives with parallelized finite differences. This is very handy if you do not want to invest the time to derive the derivatives of your criterion function.

res = om.minimize(
    fun=sphere,
    params=np.arange(5),
    algorithm="scipy_lbfgsb",
    numdiff_options=om.NumdiffOptions(n_cores=6),
)

res.params.round(5)

array([ 0.,  0.,  0., -0., -0.])

For more details and examples check out how-to speed up your optimization with derivatives

Turn local optimizers global with multistart

Multistart optimization requires finite soft bounds on all parameters. Those bounds will be used for sampling but not enforced during optimization.

bounds = om.Bounds(soft_lower=np.full(10, -5), soft_upper=np.full(10, 15))

res = om.minimize(
    fun=sphere,
    params=np.arange(10),
    algorithm="scipy_neldermead",
    bounds=bounds,
    multistart=om.MultistartOptions(convergence_max_discoveries=5),
)
res.params.round(5)

array([-0., -0., -0., -0., -0., -0., -0.,  0., -0., -0.])

And plot the history of all local optimizations

fig = om.criterion_plot(res)
fig.show(renderer="png")

Exploit the structure of your optimization problem

Many estimation problems have a least-squares structure. If so, specialized optimizers that exploit this structure can be much faster than standard optimizers. The sphere function from above is the simplest possible least-squarse problem you could imagine: the least-squares residuals are just the params.

To use least-squares optimizers in optimagic, you need to declare mark your function with a decorator and return the least-squares residuals instead of the aggregated function value.

More details can be found here

@om.mark.least_squares
def ls_sphere(params):
    return params

res = om.minimize(
    fun=ls_sphere,
    params=np.arange(5),
    algorithm="pounders",
)
res.params.round(5)

array([ 0.,  0., -0.,  0., -0.])

Of course, any least-squares problem can also be solved with a standard optimizer.

There are also specialized optimizers for likelihood functions.

Using and reading persistent logging

For long-running and difficult optimizations, it can be worthwhile to store the progress in a persistent log file. You can do this by providing a path to the logging argument:

res = om.minimize(
    fun=sphere,
    params=np.arange(5),
    algorithm="scipy_lbfgsb",
    logging="my_log.db",
    log_options={"if_database_exists": "replace"},
)

/home/docs/checkouts/readthedocs.org/user_builds/optimagic/conda/stable/lib/python3.10/site-packages/optimagic/deprecations.py:504: FutureWarning:

Usage of the parameter log_options is deprecated and will be removed in a future version. Provide a LogOptions instance for the parameter `logging`, if you need to configure the logging.

/home/docs/checkouts/readthedocs.org/user_builds/optimagic/conda/stable/lib/python3.10/site-packages/optimagic/deprecations.py:512: FutureWarning:


Using log_options={'if_database_exists': 'replace'} to create an instance of SQLiteLogOptions. This mechanism will be removed in the future.

You can read the entries in the log file (while the optimization is still running or after it has finished) as follows:

reader = om.OptimizeLogReader("my_log.db")
reader.read_history().keys()

/home/docs/checkouts/readthedocs.org/user_builds/optimagic/conda/stable/lib/python3.10/site-packages/optimagic/logging/read_log.py:25: FutureWarning:

OptimizeLogReader is deprecated and will be removed in a future version. Please use optimagic.logging.SQLiteLogReader instead.

dict_keys(['params', 'fun', 'time'])

For more information on what you can do with the log file and LogReader object, check out the logging tutorial

The persistent log file is always instantly synchronized when the optimizer tries a new parameter vector. This is very handy if an optimization has to be aborted and you want to extract the current status. It can be displayed in criterion_plot and params_plot, even while the optimization is running.