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 optimagic as om
import pandas as pd
Basic usage of minimize#
def sphere(params):
return params @ params
res = om.minimize(
fun=sphere,
params=np.arange(5),
algorithm="scipy_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.
def dict_sphere(params):
return params["a"] ** 2 + params["b"] ** 2 + (params["c"] ** 2).sum()
res = om.minimize(
fun=dict_sphere,
params={"a": 0, "b": 1, "c": pd.Series([2, 3, 4])},
algorithm="scipy_powell",
)
res.params
{'a': -6.821706323446569e-25,
'b': 2.220446049250313e-16,
'c': 0 8.881784e-16
1 8.881784e-16
2 1.776357e-15
dtype: float64}
The result contains all you need to know#
res = om.minimize(
fun=dict_sphere,
params={"a": 0, "b": 1, "c": pd.Series([2, 3, 4])},
algorithm="scipy_neldermead",
)
res
Minimize with 5 free parameters terminated successfully after 805 criterion evaluations and 507 iterations.
The value of criterion improved from 30.0 to 1.6760003634613059e-16.
The scipy_neldermead algorithm reported: Optimization terminated successfully.
Independent of the convergence criteria used by scipy_neldermead, the strength of convergence can be assessed by the following criteria:
one_step five_steps
relative_criterion_change 1.968e-15*** 2.746e-15***
relative_params_change 9.834e-08* 1.525e-07*
absolute_criterion_change 1.968e-16*** 2.746e-16***
absolute_params_change 9.834e-09** 1.525e-08*
(***: change <= 1e-10, **: change <= 1e-8, *: change <= 1e-5. Change refers to a change between accepted steps. The first column only considers the last step. The second column considers the last five steps.)
You can visualize the convergence#
fig = om.criterion_plot(res, max_evaluations=300)
fig.show(renderer="png")
fig = om.params_plot(
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#
If you install some optional dependencies, you can choose from a large (and growing) set of optimization algorithms – all with the same interface!
For example, we wrap optimizers from scipy.optimize, nlopt, cyipopt, pygmo, fides, tao and others.
We also have some optimizers that are not part of other packages. Examples are a parallel Nelder-Mead algorithm, The BHHH algorithm and a parallel Pounders algorithm.
See the full list [here](../how_to_guides/optimization/how_to_specify_algorithm_and_algo_options
You can add bounds#
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.])
You can fix parameters#
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.])
Or impose 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.33333, 0.33333, 0.33334, -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.
You can provide closed form derivatives#
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.])
Or use parallelized numerical derivatives#
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.])
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 criterion 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.
Customize your optimizer#
Most algorithms have a few optional arguments. Examples are convergence criteria or tuning parameters. You can find an overview of supported arguments here.
algo_options = {
"convergence.ftol_rel": 1e-9,
"stopping.maxiter": 100_000,
}
res = om.minimize(
fun=sphere,
params=np.arange(5),
algorithm="scipy_lbfgsb",
algo_options=algo_options,
)
res.params.round(5)
array([ 0., -0., -0., -0., -0.])