Base the default error penalty on the scalar start value

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

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
+}

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

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`."
    ]
   }
  ],