diff --git a/docs/source/explanation/aggregation_level.md b/docs/source/explanation/aggregation_level.md index f721a10c1..7cc781cd1 100644 --- a/docs/source/explanation/aggregation_level.md +++ b/docs/source/explanation/aggregation_level.md @@ -15,7 +15,7 @@ optionally `@om.mark.scalar`). That mark changes: 1. **Which specialized optimizers you can use** (for example pounders for least-squares, or BHHH for likelihood). 1. **How error penalties and derivatives are interpreted** when something goes wrong - (see {ref}`how-to-errors`). + (see {ref}`error-handling-explanation`). Any marked function can still be solved with a normal scalar optimizer; optimagic aggregates the vector output when needed (sum of squares for least-squares, sum of diff --git a/docs/source/explanation/error_handling.ipynb b/docs/source/explanation/error_handling.ipynb new file mode 100644 index 000000000..0d63f598b --- /dev/null +++ b/docs/source/explanation/error_handling.ipynb @@ -0,0 +1,279 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "(error-handling-explanation)=\n", + "\n", + "# Error handling during optimization\n", + "\n", + "This page explains what optimagic does when the objective function or its derivative\n", + "raises an error and `error_handling=\"continue\"` is set, why it works that way, and how it\n", + "relates to approaches in the literature. For practical advice, see {ref}`how-to-errors`.\n", + "\n", + "## The problem: hidden constraints\n", + "\n", + "Optimizers often try extreme parameter vectors at which the objective cannot be\n", + "evaluated: a simulation crashes, a matrix is not positive definite, or `np.exp`\n", + "overflows. In derivative-free and simulation-based optimization, regions where the\n", + "evaluation fails are called *hidden constraints* {cite}`eh-LeDigabel2024`. They are not\n", + "known in advance, and a failed evaluation only reveals that such a constraint is\n", + "violated, not by how much." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Design constraints\n", + "\n", + "Two requirements determine how optimagic handles errors:\n", + "\n", + "```{important}\n", + "1. **It has to work for solvers that cannot handle infinite values.** Many of the\n", + " optimizers that optimagic wraps, in particular gradient-based and model-based ones,\n", + " break or return nonsense when the objective returns `inf` or `nan`.\n", + "2. **It has to be implemented outside of the solver.** optimagic wraps optimizers from\n", + " many libraries and cannot change their internals. The only thing it controls is what\n", + " the objective function and its derivative return to the solver.\n", + "```\n", + "\n", + "These two requirements rule out the most common approaches from the literature (see\n", + "[below](error-handling-literature)) and lead to a replacement function with finite values and\n", + "a finite derivative." + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## The penalty function\n", + "\n", + "When an error occurs, optimagic replaces the objective by a penalty function and adjusts\n", + "the derivative accordingly. The penalty has two important properties:\n", + "\n", + "1. Its value is always worse than the objective value at the start parameters.\n", + "1. Its slope guides the optimizer back towards the start parameters. That is, if you are\n", + " minimizing, the direction of strongest decrease points towards the start parameters;\n", + " if you are maximizing, the direction of strongest increase does.\n", + "\n", + "Let $x_0$ be the internal start parameters and $f_0$ the scalar objective value at\n", + "$x_0$. For minimization, the penalty is\n", + "\n", + "$$\n", + "P(x) = c + s \\lVert x - x_0 \\rVert, \\qquad c = f_0 + |f_0| + C, \\qquad s = S,\n", + "$$\n", + "\n", + "with default values $C = 100$ and $S = 0.1$. For maximization, $c = f_0 - |f_0| - C$\n", + "and $s = -S$. You can set $c$ and $s$ via the `error_penalty` argument of `minimize` and\n", + "`maximize`.\n", + "\n", + "Least-squares and likelihood optimizers expect a vector of residuals or contributions\n", + "(see {ref}`aggregation_level`). For them, the penalty is split into $m$ identical\n", + "entries such that aggregating them gives back $P(x)$:\n", + "\n", + "- likelihood: $m$ contributions $P(x) / m$, which sum to $P(x)$;\n", + "- least squares: $m$ residuals $\\sqrt{P(x) / m}$, whose squares sum to $P(x)$.\n", + "\n", + "### Why this works\n", + "\n", + "Most local optimizers only accept steps that improve the objective. Since the penalty is\n", + "worse than the start value, and therefore worse than every iterate such an optimizer\n", + "accepts, a step into a region where errors occur is rejected: a line search backtracks,\n", + "and a trust-region method shrinks its radius. The optimizer does not simply return to a\n", + "previously evaluated point, which would lead to cyclical behavior. Instead, it takes a\n", + "shorter step in a similar direction. Since the penalty and its derivative are finite,\n", + "solvers that cannot handle infinite values keep working.\n", + "\n", + "### An example\n", + "\n", + "Consider the Jennrich–Sampson function from the Moré–Wild benchmark set, which is\n", + "numerically unstable for large parameter values. Started at $x_0 = (0.3, 0.4)$, the\n", + "penalty would look similar to the following:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import plotly.io as pio\n", + "from plotly import graph_objects as go\n", + "\n", + "pio.renderers.default = \"notebook_connected\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "def jennrich_sampson(x):\n", + " dim_out = 10\n", + " fvec = (\n", + " 2 * (1.0 + np.arange(1, dim_out + 1))\n", + " - np.exp(np.arange(1, dim_out + 1) * x[0])\n", + " - np.exp(np.arange(1, dim_out + 1) * x[1])\n", + " )\n", + " return fvec @ fvec\n", + "\n", + "\n", + "correct_params = np.array([0.2578252135686162, 0.2578252135686162])\n", + "correct_criterion = 124.3621823556148\n", + "\n", + "start_x = np.array([0.3, 0.4])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "def penalty(params):\n", + " start_params = np.array([0.3, 0.4])\n", + " # this is close to the actual value used by optimagic\n", + " constant = 8000\n", + " # the actual slope used by optimagic would be even smaller\n", + " slope = 10_000\n", + " diff = params - start_params\n", + " return constant + slope * np.linalg.norm(diff)" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "Now, let's plot the two functions. For better illustration, we assume that the jennrich_sampson function is only defined until it reaches a value of 100_000 and the penalty takes over from there. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "grid = np.linspace(0, 1)\n", + "params = [np.full(2, val) for val in grid]\n", + "values = np.array([jennrich_sampson(p) for p in params])\n", + "values = np.where(values <= 1e5, values, np.nan)\n", + "penalty_values = np.array([penalty(p) for p in params])\n", + "penalty_values = np.where(np.isfinite(values), np.nan, penalty_values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "fig = go.Figure()\n", + "fig.add_trace(go.Scatter(x=grid, y=values))\n", + "fig.add_trace(go.Scatter(x=grid, y=penalty_values))\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "We can see that the penalty is lower than the highest achieved value of `jennrich_sampson` but higher than the start values. It is also rather flat. Fortunately, that is all we need. " + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "(error-handling-literature)=\n", + "\n", + "## Related approaches in the literature\n", + "\n", + "Hidden constraints have been studied mainly in derivative-free, simulation-based and\n", + "Bayesian optimization. The main approaches are:\n", + "\n", + "1. **Extreme barrier.** Set the objective to $+\\infty$ when the evaluation fails. This\n", + " is the standard approach in direct-search methods such as mesh adaptive direct search\n", + " {cite}`eh-Audet2006`. It works because these methods only compare function values.\n", + " Nelder–Mead and most evolutionary methods tolerate it as well. It violates\n", + " requirement 1.\n", + "1. **Rejecting the step inside the solver.** Ipopt {cite}`eh-Waechter2005b` treats a\n", + " failed evaluation, or a non-finite function value, as a rejected trial step and cuts\n", + " back the step length. Trust-region methods can likewise shrink the radius. This is\n", + " arguably the cleanest solution, but it needs solver support and violates\n", + " requirement 2. optimagic's penalty can be seen as emulating this behavior from\n", + " outside the solver.\n", + "1. **Artificial function values.** SNOBFIT {cite}`eh-Huyer2008` assigns an artificial\n", + " value, based on the values at nearby points, to points where the evaluation failed.\n", + " This is the closest relative of optimagic's approach. SNOBFIT computes the\n", + " artificial values inside the solver; optimagic uses a fixed function of the distance\n", + " to the start parameters.\n", + "1. **Modeling the failure region.** Surrogate-based and Bayesian methods fit a classifier\n", + " or Gaussian process that predicts whether an evaluation succeeds and take it into\n", + " account when proposing new points {cite}`eh-Gramacy2011,eh-Gelbart2014,eh-Mueller2019`. This\n", + " pays off when evaluations are expensive, but requires a model inside the solver and\n", + " violates requirement 2.\n", + "\n", + "## Limitations\n", + "\n", + "- Model-based derivative-free optimizers such as pounders, DFO-LS or BOBYQA may include\n", + " penalized points in their interpolation set, which distorts their model of the\n", + " objective. Leaving failed points out of the model would be better, but this is not\n", + " possible from outside the solver.\n", + "- For global and population-based optimizers, the slope towards the start parameters\n", + " biases the search towards $x_0$. For solvers that support it, $+\\infty$ would be the\n", + " better choice.\n", + "- The boundary of the region where errors occur becomes an artificial discontinuity of\n", + " the objective. For optimizers that only accept improving steps this is harmless in\n", + " practice, but there is no convergence theory for it.\n", + "\n", + "Choosing the strategy per algorithm, for example returning $+\\infty$ to solvers that\n", + "handle it, is discussed in\n", + "[issue #717](https://github.com/optimagic-dev/optimagic/issues/717).\n", + "\n", + "## References\n", + "\n", + "```{eval-rst}\n", + ".. bibliography:: ../refs.bib\n", + " :filter: docname in docnames\n", + " :keyprefix: eh-\n", + "```" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/source/explanation/index.md b/docs/source/explanation/index.md index b06231439..3929a28d3 100644 --- a/docs/source/explanation/index.md +++ b/docs/source/explanation/index.md @@ -11,6 +11,7 @@ aggregation_level implementation_of_constraints internal_optimizers why_optimization_is_hard.ipynb +error_handling.ipynb explanation_of_numerical_optimizers tests_for_supported_optimizers numdiff_background diff --git a/docs/source/how_to/how_to_errors_during_optimization.ipynb b/docs/source/how_to/how_to_errors_during_optimization.ipynb index 407bbddb9..e3efdee80 100644 --- a/docs/source/how_to/how_to_errors_during_optimization.ipynb +++ b/docs/source/how_to/how_to_errors_during_optimization.ipynb @@ -154,81 +154,11 @@ "source": [ "### How does the magic work\n", "\n", - "When an error occurs and `error_handling` is set to `\"continue\"`, optimagic replaces your criterion with a dummy function (and adjusts the derivative accordingly). \n", - "\n", - "The dummy function has two important properties:\n", - "\n", - "1. Its value is always higher than criterion at start params. \n", - "2. Its slope guides the optimizer back towards the start parameters. I.e., if you are minimizing, the direction of strongest decrease is towards the start parameters; if you are maximizing, the direction of strongest increase is towards the start parameters. \n", - "\n", - "Therefore, when hitting an undefined area, an optimizer can take a few steps back until it is in better territory and then continue its work. \n", - "\n", - "Importantly, the optimizer will not simply go back to a previously evaluated point (which would just lead to cyclical behavior). It will just go back in the direction it originally came from.\n", - "\n", - "In the concrete example, the dummy function would look similar to the following:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "11", - "metadata": {}, - "outputs": [], - "source": [ - "def dummy(params):\n", - " start_params = np.array([0.3, 0.4])\n", - " # this is close to the actual value used by optimagic\n", - " constant = 8000\n", - " # the actual slope used by optimagic would be even smaller\n", - " slope = 10_000\n", - " diff = params - start_params\n", - " return constant + slope * np.linalg.norm(diff)" - ] - }, - { - "cell_type": "markdown", - "id": "12", - "metadata": {}, - "source": [ - "Now, let's plot the two functions. For better illustration, we assume that the jennrich_sampson function is only defined until it reaches a value of 100_000 and the dummy function takes over from there. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "13", - "metadata": {}, - "outputs": [], - "source": [ - "from plotly import graph_objects as go\n", - "\n", - "grid = np.linspace(0, 1)\n", - "params = [np.full(2, val) for val in grid]\n", - "values = np.array([jennrich_sampson(p) for p in params])\n", - "values = np.where(values <= 1e5, values, np.nan)\n", - "dummy_values = np.array([dummy(p) for p in params])\n", - "dummy_values = np.where(np.isfinite(values), np.nan, dummy_values)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "14", - "metadata": {}, - "outputs": [], - "source": [ - "fig = go.Figure()\n", - "fig.add_trace(go.Scatter(x=grid, y=values))\n", - "fig.add_trace(go.Scatter(x=grid, y=dummy_values))\n", - "fig.show()" - ] - }, - { - "cell_type": "markdown", - "id": "15", - "metadata": {}, - "source": [ - "We can see that the dummy function is lower than the highest achieved value of `jennrich_sampson` but higher than the start values. It is also rather flat. Fortunately, that is all we need. " + "When an error occurs and `error_handling` is set to `\"continue\"`, optimagic replaces your\n", + "criterion with a penalty function that is worse than the criterion at the start\n", + "parameters and whose slope guides the optimizer back towards them. For details, the\n", + "design constraints behind this approach, and related approaches in the literature, see\n", + "{ref}`error-handling-explanation`." ] } ], diff --git a/docs/source/refs.bib b/docs/source/refs.bib index 215f3471e..0dd26223c 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -1114,4 +1114,63 @@ @article{Ni2013 year = {2013} } +@article{LeDigabel2024, + author = {Le Digabel, S{\'e}bastien and Wild, Stefan M.}, + title = {A Taxonomy of Constraints in Black-Box Simulation-Based Optimization}, + journal = {Optimization and Engineering}, + year = {2024}, + volume = {25}, + pages = {1125--1143}, + url = {https://arxiv.org/abs/1505.07881} +} + +@article{Audet2006, + author = {Audet, Charles and Dennis, Jr., John E.}, + title = {Mesh Adaptive Direct Search Algorithms for Constrained Optimization}, + journal = {{SIAM} Journal on Optimization}, + year = {2006}, + volume = {17}, + number = {1}, + pages = {188--217}, + doi = {10.1137/040603371} +} + +@article{Huyer2008, + author = {Huyer, Waltraud and Neumaier, Arnold}, + title = {{SNOBFIT} -- Stable Noisy Optimization by Branch and Fit}, + journal = {{ACM} Transactions on Mathematical Software}, + year = {2008}, + volume = {35}, + number = {2}, + doi = {10.1145/1377612.1377613} +} + +@incollection{Gramacy2011, + author = {Gramacy, Robert B. and Lee, Herbert K. H.}, + title = {Optimization Under Unknown Constraints}, + booktitle = {Bayesian Statistics 9}, + publisher = {Oxford University Press}, + year = {2011}, + pages = {229--256} +} + +@inproceedings{Gelbart2014, + author = {Gelbart, Michael A. and Snoek, Jasper and Adams, Ryan P.}, + title = {Bayesian Optimization with Unknown Constraints}, + booktitle = {Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence}, + year = {2014}, + pages = {250--259} +} + +@article{Mueller2019, + author = {M{\"u}ller, Juliane and Day, Marcus}, + title = {Surrogate Optimization of Computationally Expensive Black-Box Problems with Hidden Constraints}, + journal = {{INFORMS} Journal on Computing}, + year = {2019}, + volume = {31}, + number = {4}, + pages = {689--702}, + doi = {10.1287/ijoc.2018.0864} +} + @Comment{jabref-meta: databaseType:bibtex;} diff --git a/src/optimagic/optimization/error_penalty.py b/src/optimagic/optimization/error_penalty.py index eab3463af..a652599fb 100644 --- a/src/optimagic/optimization/error_penalty.py +++ b/src/optimagic/optimization/error_penalty.py @@ -15,11 +15,27 @@ def _scalar_penalty( x: NDArray[np.float64], - constant: float | NDArray[np.float64], - slope: float | NDArray[np.float64], + constant: float, + slope: float, x0: NDArray[np.float64], dim_out: int | None = None, ) -> tuple[ScalarFunctionValue, NDArray[np.float64]]: # noqa: ARG001 + r"""Evaluate the scalar error penalty and its gradient. + + The penalty is :math:`P(x) = c + s \lVert x - x_0 \rVert` with gradient + :math:`\nabla P(x) = s (x - x_0) / \lVert x - x_0 \rVert`. + + Args: + x: Internal parameter vector of shape (n_params,). + constant: The constant :math:`c`. + slope: The slope :math:`s`. + x0: Internal start parameters of shape (n_params,). + dim_out: Unused; only present for a uniform signature across penalties. + + Returns: + The penalty value and its gradient of shape (n_params,). + + """ value = constant + slope * np.linalg.norm(x - x0) jac = slope * (x - x0) / np.linalg.norm(x - x0) return ScalarFunctionValue(value=value), jac @@ -27,11 +43,29 @@ def _scalar_penalty( def _likelihood_penalty( x: NDArray[np.float64], - constant: float | NDArray[np.float64], - slope: float | NDArray[np.float64], + constant: float, + slope: float, x0: NDArray[np.float64], dim_out: int, ) -> tuple[LikelihoodFunctionValue, NDArray[np.float64]]: # noqa: ARG001 + r"""Evaluate the likelihood error penalty and its Jacobian. + + The scalar penalty :math:`P(x)` (see :func:`_scalar_penalty`) is split into + :math:`m` identical contributions :math:`P(x) / m`, so they sum to + :math:`P(x)`. Each row of the Jacobian is :math:`\nabla P(x) / m`. + + Args: + x: Internal parameter vector of shape (n_params,). + constant: The constant :math:`c` of the scalar penalty. + slope: The slope :math:`s` of the scalar penalty. + x0: Internal start parameters of shape (n_params,). + dim_out: Number of contributions :math:`m`. + + Returns: + The contributions of shape (dim_out,) and the Jacobian of shape + (dim_out, n_params). + + """ factor = (constant + slope * np.linalg.norm(x - x0)) / dim_out contrib = np.ones(dim_out) * factor row = slope * (x - x0) / (dim_out * np.linalg.norm(x - x0)) @@ -41,11 +75,30 @@ def _likelihood_penalty( def _penalty_residuals( x: NDArray[np.float64], - constant: float | NDArray[np.float64], - slope: float | NDArray[np.float64], + constant: float, + slope: float, x0: NDArray[np.float64], dim_out: int, ) -> tuple[LeastSquaresFunctionValue, NDArray[np.float64]]: + r"""Evaluate the least-squares error penalty and its Jacobian. + + The scalar penalty :math:`P(x)` (see :func:`_scalar_penalty`) is split into + :math:`m` identical residuals :math:`\sqrt{P(x) / m}`, so their squares sum to + :math:`P(x)`. By the chain rule, each row of the Jacobian is + :math:`\nabla P(x) / (2 \sqrt{m P(x)})`. + + Args: + x: Internal parameter vector of shape (n_params,). + constant: The constant :math:`c` of the scalar penalty. + slope: The slope :math:`s` of the scalar penalty. + x0: Internal start parameters of shape (n_params,). + dim_out: Number of residuals :math:`m`. + + Returns: + The residuals of shape (dim_out,) and the Jacobian of shape + (dim_out, n_params). + + """ factor = np.sqrt((constant + slope * np.linalg.norm(x - x0)) / dim_out) contrib = np.ones(dim_out) * factor @@ -65,9 +118,29 @@ def get_error_penalty_function( solver_type: AggregationLevel, direction: Direction, ) -> Callable[[NDArray[np.float64]], tuple[SpecificFunctionValue, NDArray[np.float64]]]: + """Create the function that replaces the objective after an evaluation error. + + The penalty is worse than the start value and its slope points back towards + ``start_x``. Its aggregated scalar value is the same for every aggregation + level, so it does not depend on the solver type. + + Args: + start_x: Internal start parameters of shape (n_params,). + start_criterion: Objective function value at ``start_x``. + error_penalty: Optional dict with the keys "constant" and "slope". Missing + entries are filled with defaults derived from the scalar start value. + solver_type: Aggregation level the optimizer expects. + direction: Whether the problem is minimized or maximized. + + Returns: + Function that maps an internal parameter vector of shape (n_params,) to the + penalty value and its derivative, both at the aggregation level of + ``solver_type``. + + """ error_penalty = {} if error_penalty is None else error_penalty - first_value = start_criterion.internal_value(solver_type) + first_value = float(start_criterion.internal_value(AggregationLevel.SCALAR)) constant, slope = _process_error_penalty( error_penalty=error_penalty, @@ -84,8 +157,8 @@ def get_error_penalty_function( _penalty: Callable[ [ NDArray[np.float64], - float | NDArray[np.float64], - float | NDArray[np.float64], + float, + float, NDArray[np.float64], int, ], @@ -115,10 +188,31 @@ def penalty( def _process_error_penalty( error_penalty: dict[str, float] | None, - first_value: float | NDArray[np.float64], + first_value: float, direction: Direction, -) -> tuple[float | NDArray[np.float64], float | NDArray[np.float64]]: - """Add default options to error_penalty options.""" +) -> tuple[float, float]: + r"""Fill in the default constant and slope of the error penalty. + + With scalar start value :math:`f_0`, the defaults are + :math:`c = f_0 + |f_0| + C` and :math:`s = S` when minimizing, and + :math:`c = f_0 - |f_0| - C` and :math:`s = -S` when maximizing, where + :math:`C` and :math:`S` are ``CRITERION_PENALTY_CONSTANT`` and + ``CRITERION_PENALTY_SLOPE``. The constant is therefore always worse than + :math:`f_0`. + + Args: + error_penalty: Optional dict with the keys "constant" and "slope". Given + entries take precedence over the defaults. + first_value: Scalar objective value :math:`f_0` at the start parameters. + direction: Whether the problem is minimized or maximized. + + Returns: + The constant :math:`c` and the slope :math:`s`. + + Raises: + ValueError: If ``direction`` is not a valid Direction. + + """ if error_penalty is not None: error_penalty = error_penalty.copy() else: diff --git a/src/optimagic/optimization/optimize.py b/src/optimagic/optimization/optimize.py index b8a571126..ffb96c0b2 100644 --- a/src/optimagic/optimization/optimize.py +++ b/src/optimagic/optimization/optimize.py @@ -198,7 +198,7 @@ def maximize( optimization from there. For details see :ref:`how-to-errors`. error_penalty: A dictionary with the keys "slope" and "constant" that influences the magnitude of the penalty values. For maximization problems - both should be negative. For details see :ref:`how-to-errors`. + both should be negative. For details see :ref:`error-handling-explanation`. scaling: If None or False, the parameter space is not rescaled. If True, a heuristic is used to improve the conditioning of the optimization problem. To choose which heuristic is used and to customize the scaling, provide @@ -409,7 +409,7 @@ def minimize( optimization from there. For details see :ref:`how-to-errors`. error_penalty: A dictionary with the keys "slope" and "constant" that influences the magnitude of the penalty values. For minimization problems - both should be positive. For details see :ref:`how-to-errors`. + both should be positive. For details see :ref:`error-handling-explanation`. scaling: If None or False, the parameter space is not rescaled. If True, a heuristic is used to improve the conditioning of the optimization problem. To choose which heuristic is used and to customize the scaling, provide diff --git a/tests/optimagic/optimization/test_error_handling_end_to_end.py b/tests/optimagic/optimization/test_error_handling_end_to_end.py new file mode 100644 index 000000000..06a0ac52c --- /dev/null +++ b/tests/optimagic/optimization/test_error_handling_end_to_end.py @@ -0,0 +1,314 @@ +r"""End-to-end tests for ``error_handling="continue"``. + +There are two test problems, both maximum likelihood estimations with the linear +predictor :math:`\eta_i(x) = x_0 + x_1 t_i`: + +- Poisson regression for scalar and likelihood functions. The loglikelihood + contributions are :math:`y_i \eta_i - e^{\eta_i} - \log(y_i!)`. The data has a smaller + variance than the Poisson model implies, so the outer product of gradients + underestimates the Hessian and bhhh overshoots. +- Gaussian nonlinear regression for least-squares functions, with residuals + :math:`y_i - e^{\eta_i}`. + +Evaluating either objective fails if :math:`x_0 + x_1` exceeds its value at the optimum +by more than ``FAILURE_MARGIN``, similar to an overflow of the exponential. Starting at +``START_PARAMS``, all tested optimizers step into that region on their way to the +optimum. + +""" + +from dataclasses import dataclass, field + +import numpy as np +import pytest +import scipy.optimize +from numpy.testing import assert_array_almost_equal as aaae +from scipy.special import gammaln + +import optimagic as om +from optimagic.utilities import get_rng + +START_PARAMS = np.full(2, -1.0) +FAILURE_MARGIN = 0.01 + +BHHH_XFAIL = pytest.mark.xfail( + strict=True, + reason="bhhh steps uphill and accepts worse points; fixed in #719.", +) + + +def _generate_data(): + rng = get_rng(seed=0) + t = rng.uniform(size=200) + mean = np.exp(1 + t) + y_poisson = rng.binomial(10, mean / 10).astype(float) + y_gaussian = mean + rng.normal(scale=0.2, size=200) + return t, y_poisson, y_gaussian + + +T, Y_POISSON, Y_GAUSSIAN = _generate_data() + + +def _poisson_loglike(x): + eta = x[0] + x[1] * T + return Y_POISSON * eta - np.exp(eta) - gammaln(Y_POISSON + 1) + + +def _poisson_score(x): + """Return the derivative of `_poisson_loglike` of shape (n_obs, 2).""" + residuals = Y_POISSON - np.exp(x[0] + x[1] * T) + return residuals[:, None] * np.column_stack([np.ones_like(T), T]) + + +def _gaussian_residuals(x): + return Y_GAUSSIAN - np.exp(x[0] + x[1] * T) + + +def _gaussian_residuals_jac(x): + """Return the derivative of `_gaussian_residuals` of shape (n_obs, 2).""" + mean = np.exp(x[0] + x[1] * T) + return -mean[:, None] * np.column_stack([np.ones_like(T), T]) + + +POISSON_OPTIMUM = scipy.optimize.minimize( + lambda x: -_poisson_loglike(x).sum(), + x0=np.zeros(2), + jac=lambda x: -_poisson_score(x).sum(axis=0), + method="BFGS", + options={"gtol": 1e-10}, +).x + +GAUSSIAN_OPTIMUM = scipy.optimize.least_squares( + _gaussian_residuals, + x0=np.zeros(2), + jac=_gaussian_residuals_jac, + xtol=1e-15, + ftol=1e-15, + gtol=1e-15, +).x + + +@dataclass +class FailureCounter: + n_failures: int = field(default=0) + + +def _raise_in_failure_region(x, optimum, counter): + if x.sum() > optimum.sum() + FAILURE_MARGIN: + counter.n_failures += 1 + raise ValueError("Evaluation failed.") + + +def _poisson_functions(direction, counter): + """Return the loglikelihood contributions and their derivative of shape (n_obs, 2). + + For minimization, both are multiplied by -1. + + """ + sign = 1.0 if direction == "maximize" else -1.0 + + def contributions(x): + _raise_in_failure_region(x, POISSON_OPTIMUM, counter) + return sign * _poisson_loglike(x) + + def contributions_jac(x): + _raise_in_failure_region(x, POISSON_OPTIMUM, counter) + return sign * _poisson_score(x) + + return contributions, contributions_jac + + +def _scalar_problem(direction, counter): + contributions, contributions_jac = _poisson_functions(direction, counter) + fun = om.mark.scalar(lambda x: contributions(x).sum()) + jac = om.mark.scalar(lambda x: contributions_jac(x).sum(axis=0)) + return fun, jac + + +def _likelihood_problem(direction, counter): + contributions, contributions_jac = _poisson_functions(direction, counter) + return om.mark.likelihood(contributions), om.mark.likelihood(contributions_jac) + + +def _least_squares_problem(counter): + @om.mark.least_squares + def fun(x): + _raise_in_failure_region(x, GAUSSIAN_OPTIMUM, counter) + return _gaussian_residuals(x) + + @om.mark.least_squares + def jac(x): + _raise_in_failure_region(x, GAUSSIAN_OPTIMUM, counter) + return _gaussian_residuals_jac(x) + + return fun, jac + + +def _assert_poisson_recovered(res, counter): + sign = 1.0 if res.direction == "maximize" else -1.0 + optimal_value = sign * _poisson_loglike(POISSON_OPTIMUM).sum() + _assert_recovered(res, counter, POISSON_OPTIMUM, optimal_value) + + +def _assert_gaussian_recovered(res, counter): + optimal_value = (_gaussian_residuals(GAUSSIAN_OPTIMUM) ** 2).sum() + _assert_recovered(res, counter, GAUSSIAN_OPTIMUM, optimal_value) + + +def _assert_recovered(res, counter, optimum, optimal_value): + assert counter.n_failures > 0, "The failure region was never evaluated." + aaae(res.params, optimum, decimal=4) + # Penalty values are worse than the start value by at least 100, so an optimal + # value equal to the true one means that no penalty value leaked into the result. + assert np.isclose(res.fun, optimal_value, rtol=1e-8) + assert res.history is not None + history_values = [v for v in res.history.fun if v is not None] + best = max(history_values) if res.direction == "maximize" else min(history_values) + assert np.isclose(best, optimal_value, rtol=1e-8) + + +# ====================================================================================== +# scalar functions +# ====================================================================================== + + +def test_scalar_lbfgsb_minimize(): + counter = FailureCounter() + fun, _ = _scalar_problem("minimize", counter) + res = om.minimize( + fun, START_PARAMS, algorithm="scipy_lbfgsb", error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +def test_scalar_lbfgsb_maximize(): + counter = FailureCounter() + fun, _ = _scalar_problem("maximize", counter) + res = om.maximize( + fun, START_PARAMS, algorithm="scipy_lbfgsb", error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +def test_scalar_lbfgsb_minimize_with_jac(): + counter = FailureCounter() + fun, jac = _scalar_problem("minimize", counter) + res = om.minimize( + fun, START_PARAMS, algorithm="scipy_lbfgsb", jac=jac, error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +def test_scalar_lbfgsb_maximize_with_jac(): + counter = FailureCounter() + fun, jac = _scalar_problem("maximize", counter) + res = om.maximize( + fun, START_PARAMS, algorithm="scipy_lbfgsb", jac=jac, error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +def test_scalar_neldermead_minimize(): + counter = FailureCounter() + fun, _ = _scalar_problem("minimize", counter) + res = om.minimize( + fun, START_PARAMS, algorithm="scipy_neldermead", error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +def test_scalar_neldermead_maximize(): + counter = FailureCounter() + fun, _ = _scalar_problem("maximize", counter) + res = om.maximize( + fun, START_PARAMS, algorithm="scipy_neldermead", error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +# ====================================================================================== +# likelihood functions +# ====================================================================================== + + +def test_likelihood_lbfgsb_minimize(): + counter = FailureCounter() + fun, _ = _likelihood_problem("minimize", counter) + res = om.minimize( + fun, START_PARAMS, algorithm="scipy_lbfgsb", error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +def test_likelihood_lbfgsb_maximize(): + counter = FailureCounter() + fun, _ = _likelihood_problem("maximize", counter) + res = om.maximize( + fun, START_PARAMS, algorithm="scipy_lbfgsb", error_handling="continue" + ) + _assert_poisson_recovered(res, counter) + + +@BHHH_XFAIL +def test_likelihood_bhhh_minimize(): + counter = FailureCounter() + fun, _ = _likelihood_problem("minimize", counter) + res = om.minimize(fun, START_PARAMS, algorithm="bhhh", error_handling="continue") + _assert_poisson_recovered(res, counter) + + +@BHHH_XFAIL +def test_likelihood_bhhh_maximize(): + counter = FailureCounter() + fun, _ = _likelihood_problem("maximize", counter) + res = om.maximize(fun, START_PARAMS, algorithm="bhhh", error_handling="continue") + _assert_poisson_recovered(res, counter) + + +# ====================================================================================== +# least-squares functions +# ====================================================================================== + + +def test_least_squares_lbfgsb(): + counter = FailureCounter() + fun, _ = _least_squares_problem(counter) + res = om.minimize( + fun, START_PARAMS, algorithm="scipy_lbfgsb", error_handling="continue" + ) + _assert_gaussian_recovered(res, counter) + + +@BHHH_XFAIL +def test_least_squares_bhhh(): + counter = FailureCounter() + fun, _ = _least_squares_problem(counter) + res = om.minimize(fun, START_PARAMS, algorithm="bhhh", error_handling="continue") + _assert_gaussian_recovered(res, counter) + + +def test_least_squares_ls_trf(): + counter = FailureCounter() + fun, _ = _least_squares_problem(counter) + res = om.minimize( + fun, START_PARAMS, algorithm="scipy_ls_trf", error_handling="continue" + ) + _assert_gaussian_recovered(res, counter) + + +def test_least_squares_ls_trf_with_jac(): + counter = FailureCounter() + fun, jac = _least_squares_problem(counter) + res = om.minimize( + fun, START_PARAMS, algorithm="scipy_ls_trf", jac=jac, error_handling="continue" + ) + _assert_gaussian_recovered(res, counter) + + +def test_least_squares_pounders(): + counter = FailureCounter() + fun, _ = _least_squares_problem(counter) + res = om.minimize( + fun, START_PARAMS, algorithm="pounders", error_handling="continue" + ) + _assert_gaussian_recovered(res, counter) diff --git a/tests/optimagic/optimization/test_error_penalty.py b/tests/optimagic/optimization/test_error_penalty.py index 4fc4ac168..39d7e82e0 100644 --- a/tests/optimagic/optimization/test_error_penalty.py +++ b/tests/optimagic/optimization/test_error_penalty.py @@ -103,3 +103,101 @@ def test_penalty_aggregations_via_get_error_penalty(seed): assert np.isclose(scalar.value, contribs.value.sum()) # ty:ignore[unresolved-attribute] assert np.isclose(scalar.value, (root_contribs.value**2).sum()) # ty:ignore[unresolved-attribute] + + +SUPPORTED_COMBINATIONS = [ + (ScalarFunctionValue(-3.0), AggregationLevel.SCALAR, Direction.MINIMIZE), + (ScalarFunctionValue(-3.0), AggregationLevel.SCALAR, Direction.MAXIMIZE), + ( + LikelihoodFunctionValue(np.arange(-3.0, 2.0)), + AggregationLevel.SCALAR, + Direction.MINIMIZE, + ), + ( + LikelihoodFunctionValue(np.arange(-3.0, 2.0)), + AggregationLevel.SCALAR, + Direction.MAXIMIZE, + ), + ( + LikelihoodFunctionValue(np.arange(-3.0, 2.0)), + AggregationLevel.LIKELIHOOD, + Direction.MINIMIZE, + ), + ( + LikelihoodFunctionValue(np.arange(-3.0, 2.0)), + AggregationLevel.LIKELIHOOD, + Direction.MAXIMIZE, + ), + ( + LeastSquaresFunctionValue(np.arange(-2.0, 3.0)), + AggregationLevel.SCALAR, + Direction.MINIMIZE, + ), + ( + LeastSquaresFunctionValue(np.arange(-2.0, 3.0)), + AggregationLevel.LIKELIHOOD, + Direction.MINIMIZE, + ), + ( + LeastSquaresFunctionValue(np.arange(-2.0, 3.0)), + AggregationLevel.LEAST_SQUARES, + Direction.MINIMIZE, + ), +] + + +@pytest.mark.parametrize( + "start_criterion, solver_type, direction", SUPPORTED_COMBINATIONS +) +def test_default_penalty_matches_the_scalar_penalty_of_the_start_value( + start_criterion, solver_type, direction +): + """Without explicit options, every aggregation penalizes the same scalar value.""" + x = np.array([0.3, 0.7]) + x0 = np.array([0.1, 0.2]) + scalar_func = get_error_penalty_function( + start_x=x0, + start_criterion=ScalarFunctionValue( + start_criterion.internal_value(AggregationLevel.SCALAR) + ), + error_penalty=None, + solver_type=AggregationLevel.SCALAR, + direction=direction, + ) + func = get_error_penalty_function( + start_x=x0, + start_criterion=start_criterion, + error_penalty=None, + solver_type=solver_type, + direction=direction, + ) + expected, _ = scalar_func(x) + got, _ = func(x) + assert np.isclose( + got.internal_value(AggregationLevel.SCALAR), + expected.internal_value(AggregationLevel.SCALAR), + ) + + +@pytest.mark.parametrize( + "start_criterion, solver_type, direction", SUPPORTED_COMBINATIONS +) +def test_default_penalty_is_worse_than_the_start_value( + start_criterion, solver_type, direction +): + x0 = np.array([0.1, 0.2]) + func = get_error_penalty_function( + start_x=x0, + start_criterion=start_criterion, + error_penalty=None, + solver_type=solver_type, + direction=direction, + ) + start_value = start_criterion.internal_value(AggregationLevel.SCALAR) + for x in [x0 + 1e-8, np.array([0.3, 0.7]), np.array([-5.0, 10.0])]: + penalty, _ = func(x) + penalty_value = penalty.internal_value(AggregationLevel.SCALAR) + if direction == Direction.MINIMIZE: + assert penalty_value > start_value + else: + assert penalty_value < start_value