From 9b11ed16896bd48b8f148a93b569b03bb79a7e3d Mon Sep 17 00:00:00 2001 From: Hans-Martin von Gaudecker Date: Fri, 16 Oct 2026 09:26:36 +0200 Subject: [PATCH 1/6] Base the default error penalty on the scalar start value For least-squares and likelihood solvers the default penalty constant was built from the start criterion's residual or contribution vector, so it became an array: a least-squares penalty raised InvalidFunctionError and a likelihood penalty landed at the wrong level. The default now derives from the scalar start value for every aggregation level. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01Hs8mfXPYyEa24KaAp6BwCW --- src/optimagic/optimization/error_penalty.py | 2 +- .../optimization/test_error_penalty.py | 40 +++++++++++++++++++ 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/src/optimagic/optimization/error_penalty.py b/src/optimagic/optimization/error_penalty.py index eab3463af..0b9ad3a35 100644 --- a/src/optimagic/optimization/error_penalty.py +++ b/src/optimagic/optimization/error_penalty.py @@ -67,7 +67,7 @@ def get_error_penalty_function( ) -> Callable[[NDArray[np.float64]], tuple[SpecificFunctionValue, NDArray[np.float64]]]: error_penalty = {} if error_penalty is None else error_penalty - first_value = start_criterion.internal_value(solver_type) + first_value = start_criterion.internal_value(AggregationLevel.SCALAR) constant, slope = _process_error_penalty( error_penalty=error_penalty, diff --git a/tests/optimagic/optimization/test_error_penalty.py b/tests/optimagic/optimization/test_error_penalty.py index 4fc4ac168..cdf9005b2 100644 --- a/tests/optimagic/optimization/test_error_penalty.py +++ b/tests/optimagic/optimization/test_error_penalty.py @@ -103,3 +103,43 @@ 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] + + +@pytest.mark.parametrize( + "start_criterion, solver_type", + [ + (LikelihoodFunctionValue(np.arange(1.0, 6.0)), AggregationLevel.LIKELIHOOD), + ( + LeastSquaresFunctionValue(np.arange(1.0, 6.0)), + AggregationLevel.LEAST_SQUARES, + ), + ], +) +def test_default_penalty_matches_the_scalar_penalty_of_the_start_value( + start_criterion, solver_type +): + """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.MINIMIZE, + ) + func = get_error_penalty_function( + start_x=x0, + start_criterion=start_criterion, + error_penalty=None, + solver_type=solver_type, + direction=Direction.MINIMIZE, + ) + expected, _ = scalar_func(x) + got, _ = func(x) + assert np.isclose( + got.internal_value(AggregationLevel.SCALAR), + expected.internal_value(AggregationLevel.SCALAR), + ) From da9a9a1513153e8c7abcfa32c9ecf6b14639e4c3 Mon Sep 17 00:00:00 2001 From: Janos Gabler Date: Mon, 19 Oct 2026 15:49:39 +0200 Subject: [PATCH 2/6] Narrow error penalty types and document the penalty functions The default penalty constant is now always computed from the scalar start value, so constant and slope are floats for every aggregation level. Co-Authored-By: Claude Opus 5.5 --- src/optimagic/optimization/error_penalty.py | 118 ++++++++++++++++++-- 1 file changed, 106 insertions(+), 12 deletions(-) diff --git a/src/optimagic/optimization/error_penalty.py b/src/optimagic/optimization/error_penalty.py index 0b9ad3a35..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(AggregationLevel.SCALAR) + 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: From 2d704d32b46d865df69b47e8dcac64a6aee11fbf Mon Sep 17 00:00:00 2001 From: Janos Gabler Date: Mon, 19 Oct 2026 15:49:40 +0200 Subject: [PATCH 3/6] Move error handling background into an explanation page Split the explanation of the error penalty out of the how-to guide, state the two design constraints (finite values, implemented outside the solver) and add background on hidden constraints from the literature. Co-Authored-By: Claude Opus 5.5 --- docs/source/explanation/aggregation_level.md | 2 +- docs/source/explanation/error_handling.ipynb | 279 ++++++++++++++++++ docs/source/explanation/index.md | 1 + .../how_to_errors_during_optimization.ipynb | 80 +---- docs/source/refs.bib | 59 ++++ src/optimagic/optimization/optimize.py | 4 +- 6 files changed, 347 insertions(+), 78 deletions(-) create mode 100644 docs/source/explanation/error_handling.ipynb 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/optimize.py b/src/optimagic/optimization/optimize.py index cfefb8e6f..7cdbb5bb5 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 From 4a299b53bb537f2d697352598d4eacc9ed4d9834 Mon Sep 17 00:00:00 2001 From: Janos Gabler Date: Mon, 19 Oct 2026 16:01:24 +0200 Subject: [PATCH 4/6] Add end-to-end tests for error_handling="continue" Cover all supported combinations of function type, solver type and direction, with numerical and closed-form derivatives. Extend the penalty unit tests to maximization and check that the default penalty is worse than the start value. bhhh cases are marked xfail because bhhh steps uphill when the derivative is consistent with the objective. Co-Authored-By: Claude Opus 5.5 --- .../test_error_handling_end_to_end.py | 166 ++++++++++++++++++ .../optimization/test_error_penalty.py | 80 +++++++-- 2 files changed, 235 insertions(+), 11 deletions(-) create mode 100644 tests/optimagic/optimization/test_error_handling_end_to_end.py 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..706774941 --- /dev/null +++ b/tests/optimagic/optimization/test_error_handling_end_to_end.py @@ -0,0 +1,166 @@ +"""End-to-end tests for ``error_handling="continue"``. + +The test problem has the residuals :math:`r(x) = e^x - e` with optimum +:math:`x^* = (1, 1)` and optimal value 0. Starting at :math:`x_0 = (-2, -2)`, the first +steps of all tested optimizers overshoot into the region :math:`x_i > 1.5`, where the +objective raises an error. The tests cover all supported combinations of function +type, solver type and direction. + +""" + +from dataclasses import dataclass, field + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal as aaae + +import optimagic as om +from optimagic.exceptions import InvalidFunctionError + +START_PARAMS = np.full(2, -2.0) +OPTIMUM = np.ones(2) +FAILURE_THRESHOLD = 1.5 + +BHHH_XFAIL = pytest.mark.xfail( + strict=True, + reason="bhhh steps uphill if the derivative is consistent with the objective.", +) + + +@dataclass +class FailureCounter: + n_failures: int = field(default=0) + + +def _raise_in_failure_region(x, counter): + if np.any(x > FAILURE_THRESHOLD): + counter.n_failures += 1 + raise ValueError("Evaluation failed.") + + +def _get_fun_and_jac(fun_type, direction, counter): + """Return the objective and its derivative for the test problem. + + Args: + fun_type: One of "scalar", "likelihood" or "least_squares". + direction: One of "minimize" or "maximize". + counter: Counts the evaluations in the failure region. + + Returns: + The marked objective and its derivative. The derivative has shape (2,) for + scalar and shape (2, 2) for likelihood and least-squares functions. + + """ + sign = -1.0 if direction == "maximize" else 1.0 + + def residuals(x): + _raise_in_failure_region(x, counter) + return np.exp(x) - np.e + + if fun_type == "scalar": + + @om.mark.scalar + def fun(x): + r = residuals(x) + return sign * 0.5 * r @ r + + @om.mark.scalar + def jac(x): + return sign * residuals(x) * np.exp(x) + + elif fun_type == "likelihood": + + @om.mark.likelihood + def fun(x): + return sign * 0.5 * residuals(x) ** 2 + + @om.mark.likelihood + def jac(x): + return sign * np.diag(residuals(x) * np.exp(x)) + + elif fun_type == "least_squares": + fun = om.mark.least_squares(residuals) + + @om.mark.least_squares + def jac(x): + _raise_in_failure_region(x, counter) + return np.diag(np.exp(x)) + + return fun, jac + + +CASES = [ + # scalar functions + ("scalar", "scipy_lbfgsb", "minimize", False), + ("scalar", "scipy_lbfgsb", "maximize", False), + ("scalar", "scipy_lbfgsb", "minimize", True), + ("scalar", "scipy_lbfgsb", "maximize", True), + ("scalar", "scipy_neldermead", "minimize", False), + ("scalar", "scipy_neldermead", "maximize", False), + # likelihood functions + ("likelihood", "scipy_lbfgsb", "minimize", False), + ("likelihood", "scipy_lbfgsb", "maximize", False), + pytest.param("likelihood", "bhhh", "minimize", False, marks=BHHH_XFAIL), + pytest.param("likelihood", "bhhh", "maximize", False, marks=BHHH_XFAIL), + # least-squares functions + ("least_squares", "scipy_lbfgsb", "minimize", False), + pytest.param("least_squares", "bhhh", "minimize", False, marks=BHHH_XFAIL), + ("least_squares", "scipy_ls_trf", "minimize", False), + ("least_squares", "scipy_ls_trf", "minimize", True), + ("least_squares", "pounders", "minimize", False), +] + + +@pytest.mark.parametrize("fun_type, algorithm, direction, with_jac", CASES) +def test_optimizer_recovers_from_errors(fun_type, algorithm, direction, with_jac): + counter = FailureCounter() + fun, jac = _get_fun_and_jac(fun_type, direction, counter) + optimize = om.maximize if direction == "maximize" else om.minimize + + res = optimize( + fun=fun, + jac=jac if with_jac else None, + params=START_PARAMS, + algorithm=algorithm, + error_handling="continue", + ) + + assert counter.n_failures > 0, "The failure region was never evaluated." + aaae(res.params, OPTIMUM, decimal=4) + # Penalty values are at least 100 away from the start value, so an optimal value + # close to 0 means that no penalty value leaked into the result. + assert abs(res.fun) < 1e-6 + 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 direction == "maximize" else min(history_values) + assert abs(best) < 1e-6 + + +def test_least_squares_function_cannot_be_maximized(): + fun, _ = _get_fun_and_jac("least_squares", "minimize", FailureCounter()) + with pytest.raises(InvalidFunctionError, match="cannot be maximized"): + om.maximize( + fun=fun, + params=START_PARAMS, + algorithm="scipy_lbfgsb", + error_handling="continue", + ) + + +@pytest.mark.parametrize( + "fun_type, algorithm", + [ + ("scalar", "bhhh"), + ("scalar", "scipy_ls_trf"), + ("likelihood", "scipy_ls_trf"), + ], +) +def test_unsupported_solver_type_raises(fun_type, algorithm): + fun, _ = _get_fun_and_jac(fun_type, "minimize", FailureCounter()) + with pytest.raises(InvalidFunctionError): + om.minimize( + fun=fun, + params=START_PARAMS, + algorithm=algorithm, + error_handling="continue", + ) diff --git a/tests/optimagic/optimization/test_error_penalty.py b/tests/optimagic/optimization/test_error_penalty.py index cdf9005b2..39d7e82e0 100644 --- a/tests/optimagic/optimization/test_error_penalty.py +++ b/tests/optimagic/optimization/test_error_penalty.py @@ -105,18 +105,52 @@ def test_penalty_aggregations_via_get_error_penalty(seed): 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", - [ - (LikelihoodFunctionValue(np.arange(1.0, 6.0)), AggregationLevel.LIKELIHOOD), - ( - LeastSquaresFunctionValue(np.arange(1.0, 6.0)), - AggregationLevel.LEAST_SQUARES, - ), - ], + "start_criterion, solver_type, direction", SUPPORTED_COMBINATIONS ) def test_default_penalty_matches_the_scalar_penalty_of_the_start_value( - start_criterion, solver_type + start_criterion, solver_type, direction ): """Without explicit options, every aggregation penalizes the same scalar value.""" x = np.array([0.3, 0.7]) @@ -128,14 +162,14 @@ def test_default_penalty_matches_the_scalar_penalty_of_the_start_value( ), error_penalty=None, solver_type=AggregationLevel.SCALAR, - direction=Direction.MINIMIZE, + direction=direction, ) func = get_error_penalty_function( start_x=x0, start_criterion=start_criterion, error_penalty=None, solver_type=solver_type, - direction=Direction.MINIMIZE, + direction=direction, ) expected, _ = scalar_func(x) got, _ = func(x) @@ -143,3 +177,27 @@ def test_default_penalty_matches_the_scalar_penalty_of_the_start_value( 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 From 6632eea121430ce01e9d4abdd32aac30e0079b23 Mon Sep 17 00:00:00 2001 From: Janos Gabler Date: Mon, 19 Oct 2026 19:23:39 +0200 Subject: [PATCH 5/6] Use actual likelihood problems in the error handling end-to-end tests Replace the toy problem by a Poisson regression for scalar and likelihood functions and a Gaussian nonlinear regression for least-squares functions, so that bhhh is tested on problems it is designed for. Co-Authored-By: Claude Opus 5.5 --- .../test_error_handling_end_to_end.py | 165 ++++++++++++------ 1 file changed, 116 insertions(+), 49 deletions(-) diff --git a/tests/optimagic/optimization/test_error_handling_end_to_end.py b/tests/optimagic/optimization/test_error_handling_end_to_end.py index 706774941..c0e595907 100644 --- a/tests/optimagic/optimization/test_error_handling_end_to_end.py +++ b/tests/optimagic/optimization/test_error_handling_end_to_end.py @@ -1,10 +1,20 @@ -"""End-to-end tests for ``error_handling="continue"``. +r"""End-to-end tests for ``error_handling="continue"``. -The test problem has the residuals :math:`r(x) = e^x - e` with optimum -:math:`x^* = (1, 1)` and optimal value 0. Starting at :math:`x_0 = (-2, -2)`, the first -steps of all tested optimizers overshoot into the region :math:`x_i > 1.5`, where the -objective raises an error. The tests cover all supported combinations of function -type, solver type and direction. +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. The tests cover all supported combinations of function type, solver type and +direction. """ @@ -12,34 +22,87 @@ 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.exceptions import InvalidFunctionError +from optimagic.utilities import get_rng -START_PARAMS = np.full(2, -2.0) -OPTIMUM = np.ones(2) -FAILURE_THRESHOLD = 1.5 +START_PARAMS = np.full(2, -1.0) +FAILURE_MARGIN = 0.01 BHHH_XFAIL = pytest.mark.xfail( strict=True, - reason="bhhh steps uphill if the derivative is consistent with the objective.", + 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, counter): - if np.any(x > FAILURE_THRESHOLD): +def _raise_in_failure_region(x, optimum, counter): + if x.sum() > optimum.sum() + FAILURE_MARGIN: counter.n_failures += 1 raise ValueError("Evaluation failed.") -def _get_fun_and_jac(fun_type, direction, counter): - """Return the objective and its derivative for the test problem. +def _get_problem(fun_type, direction, counter): + """Return the objective, its derivative, the optimum and the optimal value. Args: fun_type: One of "scalar", "likelihood" or "least_squares". @@ -47,46 +110,50 @@ def _get_fun_and_jac(fun_type, direction, counter): counter: Counts the evaluations in the failure region. Returns: - The marked objective and its derivative. The derivative has shape (2,) for - scalar and shape (2, 2) for likelihood and least-squares functions. + The marked objective, its marked derivative, the optimal parameters of shape + (2,) and the aggregated objective value at the optimum. The derivative has + shape (2,) for scalar and shape (n_obs, 2) for likelihood and least-squares + functions. """ - sign = -1.0 if direction == "maximize" else 1.0 + # The loglikelihood is maximized; minimization uses its negative. + sign = 1.0 if direction == "maximize" else -1.0 - def residuals(x): - _raise_in_failure_region(x, counter) - return np.exp(x) - np.e + if fun_type in ("scalar", "likelihood"): + optimum = POISSON_OPTIMUM - if fun_type == "scalar": + def contributions(x): + _raise_in_failure_region(x, optimum, counter) + return sign * _poisson_loglike(x) - @om.mark.scalar - def fun(x): - r = residuals(x) - return sign * 0.5 * r @ r + def contributions_jac(x): + _raise_in_failure_region(x, optimum, counter) + return sign * _poisson_score(x) - @om.mark.scalar - def jac(x): - return sign * residuals(x) * np.exp(x) + if fun_type == "scalar": + fun = om.mark.scalar(lambda x: contributions(x).sum()) + jac = om.mark.scalar(lambda x: contributions_jac(x).sum(axis=0)) + else: + fun = om.mark.likelihood(contributions) + jac = om.mark.likelihood(contributions_jac) + optimal_value = sign * _poisson_loglike(optimum).sum() - elif fun_type == "likelihood": + elif fun_type == "least_squares": + optimum = GAUSSIAN_OPTIMUM - @om.mark.likelihood + @om.mark.least_squares def fun(x): - return sign * 0.5 * residuals(x) ** 2 - - @om.mark.likelihood - def jac(x): - return sign * np.diag(residuals(x) * np.exp(x)) - - elif fun_type == "least_squares": - fun = om.mark.least_squares(residuals) + _raise_in_failure_region(x, optimum, counter) + return _gaussian_residuals(x) @om.mark.least_squares def jac(x): - _raise_in_failure_region(x, counter) - return np.diag(np.exp(x)) + _raise_in_failure_region(x, optimum, counter) + return _gaussian_residuals_jac(x) + + optimal_value = (_gaussian_residuals(optimum) ** 2).sum() - return fun, jac + return fun, jac, optimum, optimal_value CASES = [ @@ -114,7 +181,7 @@ def jac(x): @pytest.mark.parametrize("fun_type, algorithm, direction, with_jac", CASES) def test_optimizer_recovers_from_errors(fun_type, algorithm, direction, with_jac): counter = FailureCounter() - fun, jac = _get_fun_and_jac(fun_type, direction, counter) + fun, jac, optimum, optimal_value = _get_problem(fun_type, direction, counter) optimize = om.maximize if direction == "maximize" else om.minimize res = optimize( @@ -126,18 +193,18 @@ def test_optimizer_recovers_from_errors(fun_type, algorithm, direction, with_jac ) assert counter.n_failures > 0, "The failure region was never evaluated." - aaae(res.params, OPTIMUM, decimal=4) - # Penalty values are at least 100 away from the start value, so an optimal value - # close to 0 means that no penalty value leaked into the result. - assert abs(res.fun) < 1e-6 + 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 direction == "maximize" else min(history_values) - assert abs(best) < 1e-6 + assert np.isclose(best, optimal_value, rtol=1e-8) def test_least_squares_function_cannot_be_maximized(): - fun, _ = _get_fun_and_jac("least_squares", "minimize", FailureCounter()) + fun, *_ = _get_problem("least_squares", "minimize", FailureCounter()) with pytest.raises(InvalidFunctionError, match="cannot be maximized"): om.maximize( fun=fun, @@ -156,7 +223,7 @@ def test_least_squares_function_cannot_be_maximized(): ], ) def test_unsupported_solver_type_raises(fun_type, algorithm): - fun, _ = _get_fun_and_jac(fun_type, "minimize", FailureCounter()) + fun, *_ = _get_problem(fun_type, "minimize", FailureCounter()) with pytest.raises(InvalidFunctionError): om.minimize( fun=fun, From 3dbe012b9f09159c4611e6757ab7f15d566d75b6 Mon Sep 17 00:00:00 2001 From: Janos Gabler Date: Mon, 19 Oct 2026 19:52:43 +0200 Subject: [PATCH 6/6] Replace parametrized error handling tests with individual test cases Drop the tests for maximizing least-squares problems (already covered in test_function_formats_ls.py) and for incompatible solver types (moved to #720). Co-Authored-By: Claude Opus 5.5 --- .../test_error_handling_end_to_end.py | 315 +++++++++++------- 1 file changed, 198 insertions(+), 117 deletions(-) diff --git a/tests/optimagic/optimization/test_error_handling_end_to_end.py b/tests/optimagic/optimization/test_error_handling_end_to_end.py index c0e595907..06a0ac52c 100644 --- a/tests/optimagic/optimization/test_error_handling_end_to_end.py +++ b/tests/optimagic/optimization/test_error_handling_end_to_end.py @@ -13,8 +13,7 @@ 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. The tests cover all supported combinations of function type, solver type and -direction. +optimum. """ @@ -27,7 +26,6 @@ from scipy.special import gammaln import optimagic as om -from optimagic.exceptions import InvalidFunctionError from optimagic.utilities import get_rng START_PARAMS = np.full(2, -1.0) @@ -101,97 +99,63 @@ def _raise_in_failure_region(x, optimum, counter): raise ValueError("Evaluation failed.") -def _get_problem(fun_type, direction, counter): - """Return the objective, its derivative, the optimum and the optimal value. +def _poisson_functions(direction, counter): + """Return the loglikelihood contributions and their derivative of shape (n_obs, 2). - Args: - fun_type: One of "scalar", "likelihood" or "least_squares". - direction: One of "minimize" or "maximize". - counter: Counts the evaluations in the failure region. - - Returns: - The marked objective, its marked derivative, the optimal parameters of shape - (2,) and the aggregated objective value at the optimum. The derivative has - shape (2,) for scalar and shape (n_obs, 2) for likelihood and least-squares - functions. + For minimization, both are multiplied by -1. """ - # The loglikelihood is maximized; minimization uses its negative. sign = 1.0 if direction == "maximize" else -1.0 - if fun_type in ("scalar", "likelihood"): - optimum = POISSON_OPTIMUM - - def contributions(x): - _raise_in_failure_region(x, optimum, counter) - return sign * _poisson_loglike(x) - - def contributions_jac(x): - _raise_in_failure_region(x, optimum, counter) - return sign * _poisson_score(x) - - if fun_type == "scalar": - fun = om.mark.scalar(lambda x: contributions(x).sum()) - jac = om.mark.scalar(lambda x: contributions_jac(x).sum(axis=0)) - else: - fun = om.mark.likelihood(contributions) - jac = om.mark.likelihood(contributions_jac) - optimal_value = sign * _poisson_loglike(optimum).sum() - - elif fun_type == "least_squares": - optimum = GAUSSIAN_OPTIMUM - - @om.mark.least_squares - def fun(x): - _raise_in_failure_region(x, optimum, counter) - return _gaussian_residuals(x) - - @om.mark.least_squares - def jac(x): - _raise_in_failure_region(x, optimum, counter) - return _gaussian_residuals_jac(x) - - optimal_value = (_gaussian_residuals(optimum) ** 2).sum() - - return fun, jac, optimum, optimal_value - - -CASES = [ - # scalar functions - ("scalar", "scipy_lbfgsb", "minimize", False), - ("scalar", "scipy_lbfgsb", "maximize", False), - ("scalar", "scipy_lbfgsb", "minimize", True), - ("scalar", "scipy_lbfgsb", "maximize", True), - ("scalar", "scipy_neldermead", "minimize", False), - ("scalar", "scipy_neldermead", "maximize", False), - # likelihood functions - ("likelihood", "scipy_lbfgsb", "minimize", False), - ("likelihood", "scipy_lbfgsb", "maximize", False), - pytest.param("likelihood", "bhhh", "minimize", False, marks=BHHH_XFAIL), - pytest.param("likelihood", "bhhh", "maximize", False, marks=BHHH_XFAIL), - # least-squares functions - ("least_squares", "scipy_lbfgsb", "minimize", False), - pytest.param("least_squares", "bhhh", "minimize", False, marks=BHHH_XFAIL), - ("least_squares", "scipy_ls_trf", "minimize", False), - ("least_squares", "scipy_ls_trf", "minimize", True), - ("least_squares", "pounders", "minimize", False), -] - - -@pytest.mark.parametrize("fun_type, algorithm, direction, with_jac", CASES) -def test_optimizer_recovers_from_errors(fun_type, algorithm, direction, with_jac): - counter = FailureCounter() - fun, jac, optimum, optimal_value = _get_problem(fun_type, direction, counter) - optimize = om.maximize if direction == "maximize" else om.minimize - - res = optimize( - fun=fun, - jac=jac if with_jac else None, - params=START_PARAMS, - algorithm=algorithm, - error_handling="continue", - ) + 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 @@ -199,35 +163,152 @@ def test_optimizer_recovers_from_errors(fun_type, algorithm, direction, with_jac 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 direction == "maximize" else min(history_values) + best = max(history_values) if res.direction == "maximize" else min(history_values) assert np.isclose(best, optimal_value, rtol=1e-8) -def test_least_squares_function_cannot_be_maximized(): - fun, *_ = _get_problem("least_squares", "minimize", FailureCounter()) - with pytest.raises(InvalidFunctionError, match="cannot be maximized"): - om.maximize( - fun=fun, - params=START_PARAMS, - algorithm="scipy_lbfgsb", - error_handling="continue", - ) - - -@pytest.mark.parametrize( - "fun_type, algorithm", - [ - ("scalar", "bhhh"), - ("scalar", "scipy_ls_trf"), - ("likelihood", "scipy_ls_trf"), - ], -) -def test_unsupported_solver_type_raises(fun_type, algorithm): - fun, *_ = _get_problem(fun_type, "minimize", FailureCounter()) - with pytest.raises(InvalidFunctionError): - om.minimize( - fun=fun, - params=START_PARAMS, - algorithm=algorithm, - error_handling="continue", - ) +# ====================================================================================== +# 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)