diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 631cb3985..5366aa859 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -34,10 +34,18 @@ jobs: cache-environment: true create-args: | python=${{ matrix.python-version }} - - name: run pytest + - name: run pytest (for python 3.13) + shell: bash -l {0} + if: runner.os == 'Linux' && matrix.python-version == '3.13' + run: | + micromamba activate optimagic + pytest --cov-report=xml --cov=./ + - name: run pytest (for python < 3.13 with pip install pyensmallen) shell: bash -l {0} + if: runner.os == 'Linux' && matrix.python-version < '3.13' run: | micromamba activate optimagic + pip install pyensmallen pytest --cov-report=xml --cov=./ - name: Upload coverage report. if: runner.os == 'Linux' && matrix.python-version == '3.10' diff --git a/docs/source/algorithms.md b/docs/source/algorithms.md index bd8837b9a..01bbc91cb 100644 --- a/docs/source/algorithms.md +++ b/docs/source/algorithms.md @@ -4701,6 +4701,39 @@ package. To use it, you need to have - **seed**: Seed for the random number generator for reproducibility. ``` +## Optimizers from the Ensmallen C++ library + +optimagic supports some optimizers from the Ensmallen C++ library. Optimizers from this +library are made available in Python through the pyensmallen python wrapper. To use +optimizers from Ensmallen, you need to have +[pyensmallen](https://pypi.org/project/pyensmallen-experimental/) installed (pip install +pyensmallen_experimental). + +````{eval-rst} +.. dropdown:: ensmallen_lbfgs + + .. code-block:: + + "ensmallen_lbfgs" + + Minimize a scalar function using the “LBFGS” algorithm. + + L-BFGS is an optimization algorithm in the family of quasi-Newton methods that approximates the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm using a limited amount of computer memory. + + Detailed description of the algorithm is given in :cite:`Matthies1979`. + + - **limited_memory_storage_length** (int): Maximum number of saved gradients used to approximate the hessian matrix.. + - **stopping.maxiter** (int): Maximum number of iterations for the optimization (0 means no limit and may run indefinitely). + - **armijo_constant** (float): Controls the accuracy of the line search routine for determining the Armijo condition. default is 1e-4. + - **wolfe_condition** (float): Parameter for detecting the Wolfe condition. default is 0.9. + - **convergence.gtol_abs** (float): Stop when the absolute gradient norm is smaller than this. + - **convergence.ftol_rel** (float): Stop when the relative improvement between two iterations is below this. + - **max_line_search_trials** (int): The maximum number of trials for the line search (before giving up). default is 50. + - **min_step_for_line_search** (float): The minimum step of the line search. default is 1e-20. + - **max_step_for_line_search** (float): The maximum step of the line search. default is 1e20. + + + ## References ```{eval-rst} @@ -4709,3 +4742,4 @@ package. To use it, you need to have :filter: docname in docnames :style: unsrt ``` +```` diff --git a/docs/source/refs.bib b/docs/source/refs.bib index f8005d2e9..635291494 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -893,6 +893,19 @@ @book{Conn2009 URL = {https://epubs.siam.org/doi/abs/10.1137/1.9780898718768}, } + +@article{Matthies1979, + author = {H. Matthies and G. Strang}, + title = {The Solution of Nonlinear Finite Element Equations}, + journal = {International Journal for Numerical Methods in Engineering}, + volume = {14}, + number = {11}, + pages = {1613-1626}, + year = {1979}, + doi = {10.1002/nme.1620141104} +} + + @article{JAMES1975343, title = {Minuit - a system for function minimization and analysis of the parameter errors and correlations}, journal = {Computer Physics Communications}, diff --git a/pyproject.toml b/pyproject.toml index c74752252..dcad32bea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -334,6 +334,7 @@ ignore_errors = true [[tool.mypy.overrides]] module = [ + "pyensmallen_experimental", "pybaum", "scipy", "scipy.linalg", diff --git a/src/optimagic/algorithms.py b/src/optimagic/algorithms.py index f86792478..b681f59fe 100644 --- a/src/optimagic/algorithms.py +++ b/src/optimagic/algorithms.py @@ -53,6 +53,7 @@ NloptVAR, ) from optimagic.optimizers.pounders import Pounders +from optimagic.optimizers.pyensmallen_optimizers import EnsmallenLBFGS from optimagic.optimizers.pygmo_optimizers import ( PygmoBeeColony, PygmoCmaes, @@ -995,6 +996,7 @@ class GradientBasedLocalScalarAlgorithms(AlgoSelection): nlopt_slsqp: Type[NloptSLSQP] = NloptSLSQP nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS scipy_bfgs: Type[ScipyBFGS] = ScipyBFGS scipy_conjugate_gradient: Type[ScipyConjugateGradient] = ScipyConjugateGradient scipy_lbfgsb: Type[ScipyLBFGSB] = ScipyLBFGSB @@ -2186,6 +2188,7 @@ class GradientBasedLocalAlgorithms(AlgoSelection): nlopt_slsqp: Type[NloptSLSQP] = NloptSLSQP nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS scipy_bfgs: Type[ScipyBFGS] = ScipyBFGS scipy_conjugate_gradient: Type[ScipyConjugateGradient] = ScipyConjugateGradient scipy_lbfgsb: Type[ScipyLBFGSB] = ScipyLBFGSB @@ -2299,6 +2302,7 @@ class GradientBasedScalarAlgorithms(AlgoSelection): nlopt_slsqp: Type[NloptSLSQP] = NloptSLSQP nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS scipy_bfgs: Type[ScipyBFGS] = ScipyBFGS scipy_basinhopping: Type[ScipyBasinhopping] = ScipyBasinhopping scipy_conjugate_gradient: Type[ScipyConjugateGradient] = ScipyConjugateGradient @@ -3010,6 +3014,7 @@ class LocalScalarAlgorithms(AlgoSelection): nlopt_sbplx: Type[NloptSbplx] = NloptSbplx nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS scipy_bfgs: Type[ScipyBFGS] = ScipyBFGS scipy_cobyla: Type[ScipyCOBYLA] = ScipyCOBYLA scipy_conjugate_gradient: Type[ScipyConjugateGradient] = ScipyConjugateGradient @@ -3448,6 +3453,7 @@ class GradientBasedAlgorithms(AlgoSelection): nlopt_slsqp: Type[NloptSLSQP] = NloptSLSQP nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS scipy_bfgs: Type[ScipyBFGS] = ScipyBFGS scipy_basinhopping: Type[ScipyBasinhopping] = ScipyBasinhopping scipy_conjugate_gradient: Type[ScipyConjugateGradient] = ScipyConjugateGradient @@ -3670,6 +3676,7 @@ class LocalAlgorithms(AlgoSelection): nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR pounders: Type[Pounders] = Pounders + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS scipy_bfgs: Type[ScipyBFGS] = ScipyBFGS scipy_cobyla: Type[ScipyCOBYLA] = ScipyCOBYLA scipy_conjugate_gradient: Type[ScipyConjugateGradient] = ScipyConjugateGradient @@ -3907,6 +3914,7 @@ class ScalarAlgorithms(AlgoSelection): nlopt_sbplx: Type[NloptSbplx] = NloptSbplx nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS pygmo_bee_colony: Type[PygmoBeeColony] = PygmoBeeColony pygmo_cmaes: Type[PygmoCmaes] = PygmoCmaes pygmo_compass_search: Type[PygmoCompassSearch] = PygmoCompassSearch @@ -4111,6 +4119,7 @@ class Algorithms(AlgoSelection): nlopt_tnewton: Type[NloptTNewton] = NloptTNewton nlopt_var: Type[NloptVAR] = NloptVAR pounders: Type[Pounders] = Pounders + ensmallen_lbfgs: Type[EnsmallenLBFGS] = EnsmallenLBFGS pygmo_bee_colony: Type[PygmoBeeColony] = PygmoBeeColony pygmo_cmaes: Type[PygmoCmaes] = PygmoCmaes pygmo_compass_search: Type[PygmoCompassSearch] = PygmoCompassSearch diff --git a/src/optimagic/config.py b/src/optimagic/config.py index ce6cd4d60..9aaa16a88 100644 --- a/src/optimagic/config.py +++ b/src/optimagic/config.py @@ -39,7 +39,7 @@ def _is_installed(module_name: str) -> bool: IS_IMINUIT_INSTALLED = _is_installed("iminuit") IS_NEVERGRAD_INSTALLED = _is_installed("nevergrad") IS_BAYESOPT_INSTALLED = _is_installed("bayes_opt") - +IS_PYENSMALLEN_INSTALLED = _is_installed("pyensmallen") # ====================================================================================== # Check if pandas version is newer or equal to version 2.1.0 diff --git a/src/optimagic/optimizers/pyensmallen_optimizers.py b/src/optimagic/optimizers/pyensmallen_optimizers.py new file mode 100644 index 000000000..571d4dfc4 --- /dev/null +++ b/src/optimagic/optimizers/pyensmallen_optimizers.py @@ -0,0 +1,100 @@ +"""Implement ensmallen optimizers.""" + +from dataclasses import dataclass + +import numpy as np +from numpy.typing import NDArray + +from optimagic import mark +from optimagic.config import IS_PYENSMALLEN_INSTALLED +from optimagic.optimization.algo_options import ( + CONVERGENCE_FTOL_REL, + CONVERGENCE_GTOL_ABS, + LIMITED_MEMORY_STORAGE_LENGTH, + MAX_LINE_SEARCH_STEPS, + STOPPING_MAXITER, +) +from optimagic.optimization.algorithm import Algorithm, InternalOptimizeResult +from optimagic.optimization.internal_optimization_problem import ( + InternalOptimizationProblem, +) +from optimagic.typing import AggregationLevel, NonNegativeFloat, PositiveInt + +if IS_PYENSMALLEN_INSTALLED: + import pyensmallen as pye + +MIN_LINE_SEARCH_STEPS = 1e-20 +"""The minimum step of the line search.""" +MAX_LINE_SEARCH_TRIALS = 50 +"""The maximum number of trials for the line search (before giving up).""" +ARMIJO_CONSTANT = 1e-4 +"""Controls the accuracy of the line search routine for determining the Armijo +condition.""" +WOLFE_CONDITION = 0.9 +"""Parameter for detecting the Wolfe condition.""" + +STEP_SIZE = 0.001 +"""Step size for each iteration.""" +BATCH_SIZE = 32 +"""Step size for each iteration.""" +EXP_DECAY_RATE_FOR_FIRST_MOMENT = 0.9 +"""Exponential decay rate for the first moment estimates.""" +EXP_DECAY_RATE_FOR_WEIGHTED_INF_NORM = 0.999 +"""Exponential decay rate for the first moment estimates.""" + + +@mark.minimizer( + name="ensmallen_lbfgs", + solver_type=AggregationLevel.SCALAR, + is_available=IS_PYENSMALLEN_INSTALLED, + is_global=False, + needs_jac=True, + needs_hess=False, + supports_parallelism=False, + supports_bounds=False, + supports_linear_constraints=False, + supports_nonlinear_constraints=False, + disable_history=False, +) +@dataclass(frozen=True) +class EnsmallenLBFGS(Algorithm): + limited_memory_storage_length: PositiveInt = LIMITED_MEMORY_STORAGE_LENGTH + stopping_maxiter: PositiveInt = STOPPING_MAXITER + armijo_constant: NonNegativeFloat = ARMIJO_CONSTANT # needs review + wolfe_condition: NonNegativeFloat = WOLFE_CONDITION # needs review + convergence_gtol_abs: NonNegativeFloat = CONVERGENCE_GTOL_ABS + convergence_ftol_rel: NonNegativeFloat = CONVERGENCE_FTOL_REL + max_line_search_trials: PositiveInt = MAX_LINE_SEARCH_TRIALS + min_step_for_line_search: NonNegativeFloat = MIN_LINE_SEARCH_STEPS + max_step_for_line_search: NonNegativeFloat = MAX_LINE_SEARCH_STEPS + + def _solve_internal_problem( + self, problem: InternalOptimizationProblem, x0: NDArray[np.float64] + ) -> InternalOptimizeResult: + optimizer = pye.L_BFGS( + numBasis=self.limited_memory_storage_length, + maxIterations=self.stopping_maxiter, + armijoConstant=self.armijo_constant, + wolfe=self.wolfe_condition, + minGradientNorm=self.convergence_gtol_abs, + factr=self.convergence_ftol_rel, + maxLineSearchTrials=self.max_line_search_trials, + minStep=self.min_step_for_line_search, + maxStep=self.max_step_for_line_search, + ) + + def objective_function( + x: NDArray[np.float64], grad: NDArray[np.float64] + ) -> np.float64: + fun, jac = problem.fun_and_jac(x) + grad[:] = jac + return np.float64(fun) + + best_x = optimizer.optimize(objective_function, x0, report) + + res = InternalOptimizeResult( + x=best_x, + fun=problem.fun(best_x), + ) + + return res diff --git a/tests/optimagic/optimizers/test_pyensmallen_optimizers.py b/tests/optimagic/optimizers/test_pyensmallen_optimizers.py new file mode 100644 index 000000000..7a6bdfd17 --- /dev/null +++ b/tests/optimagic/optimizers/test_pyensmallen_optimizers.py @@ -0,0 +1,22 @@ +"""Tests for pyensmallen optimizers.""" + +import numpy as np +import pytest + +import optimagic as om +from optimagic.config import IS_PYENSMALLEN_INSTALLED +from optimagic.optimization.optimize import minimize + + +@pytest.mark.skipif(not IS_PYENSMALLEN_INSTALLED, reason="pyensmallen not installed.") +def test_stop_after_one_iteration(): + algo = om.algos.ensmallen_lbfgs(stopping_maxiter=1) + expected = np.array([0, 0.81742581, 1.63485163, 2.45227744, 3.26970326]) + res = minimize( + fun=lambda x: x @ x, + fun_and_jac=lambda x: (x @ x, 2 * x), + params=np.arange(5), + algorithm=algo, + ) + + assert np.allclose(res.x, expected)