Error handling during optimization

This page explains what optimagic does when the objective function or its derivative raises an error and error_handling="continue" is set, why it works that way, and how it relates to approaches in the literature. For practical advice, see How to handle errors during optimization.

The problem: hidden constraints

Optimizers often try extreme parameter vectors at which the objective cannot be evaluated: a simulation crashes, a matrix is not positive definite, or np.exp overflows. In derivative-free and simulation-based optimization, regions where the evaluation fails are called hidden constraints [LDW24]. They are not known in advance, and a failed evaluation only reveals that such a constraint is violated, not by how much.

Design constraints

Two requirements determine how optimagic handles errors:

Important

1. It has to work for solvers that cannot handle infinite values. Many of the optimizers that optimagic wraps, in particular gradient-based and model-based ones, break or return nonsense when the objective returns inf or nan.

2. It has to be implemented outside of the solver. optimagic wraps optimizers from many libraries and cannot change their internals. The only thing it controls is what the objective function and its derivative return to the solver.

These two requirements rule out the most common approaches from the literature (see below) and lead to a replacement function with finite values and a finite derivative.

The penalty function

When an error occurs, optimagic replaces the objective by a penalty function and adjusts the derivative accordingly. The penalty has two important properties:

1. Its value is always worse than the objective value at the start parameters.

2. Its slope guides the optimizer back towards the start parameters. That is, if you are minimizing, the direction of strongest decrease points towards the start parameters; if you are maximizing, the direction of strongest increase does.

Let \(x_0\) be the internal start parameters and \(f_0\) the scalar objective value at \(x_0\). For minimization, the penalty is

\[ P(x) = c + s \lVert x - x_0 \rVert, \qquad c = f_0 + |f_0| + C, \qquad s = S, \]

with default values \(C = 100\) and \(S = 0.1\). For maximization, \(c = f_0 - |f_0| - C\) and \(s = -S\). You can set \(c\) and \(s\) via the error_penalty argument of minimize and maximize.

Least-squares and likelihood optimizers expect a vector of residuals or contributions (see Problem types (AggregationLevel)). For them, the penalty is split into \(m\) identical entries such that aggregating them gives back \(P(x)\):

• likelihood: \(m\) contributions \(P(x) / m\), which sum to \(P(x)\);

• least squares: \(m\) residuals \(\sqrt{P(x) / m}\), whose squares sum to \(P(x)\).

Why this works

Most local optimizers only accept steps that improve the objective. Since the penalty is worse than the start value, and therefore worse than every iterate such an optimizer accepts, a step into a region where errors occur is rejected: a line search backtracks, and a trust-region method shrinks its radius. The optimizer does not simply return to a previously evaluated point, which would lead to cyclical behavior. Instead, it takes a shorter step in a similar direction. Since the penalty and its derivative are finite, solvers that cannot handle infinite values keep working.

An example

Consider the Jennrich–Sampson function from the Moré–Wild benchmark set, which is numerically unstable for large parameter values. Started at \(x_0 = (0.3, 0.4)\), the penalty would look similar to the following:

import numpy as np
import plotly.io as pio
from plotly import graph_objects as go

pio.renderers.default = "notebook_connected"

def jennrich_sampson(x):
    dim_out = 10
    fvec = (
        2 * (1.0 + np.arange(1, dim_out + 1))
        - np.exp(np.arange(1, dim_out + 1) * x[0])
        - np.exp(np.arange(1, dim_out + 1) * x[1])
    )
    return fvec @ fvec


correct_params = np.array([0.2578252135686162, 0.2578252135686162])
correct_criterion = 124.3621823556148

start_x = np.array([0.3, 0.4])

def penalty(params):
    start_params = np.array([0.3, 0.4])
    # this is close to the actual value used by optimagic
    constant = 8000
    # the actual slope used by optimagic would be even smaller
    slope = 10_000
    diff = params - start_params
    return constant + slope * np.linalg.norm(diff)

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.

grid = np.linspace(0, 1)
params = [np.full(2, val) for val in grid]
values = np.array([jennrich_sampson(p) for p in params])
values = np.where(values <= 1e5, values, np.nan)
penalty_values = np.array([penalty(p) for p in params])
penalty_values = np.where(np.isfinite(values), np.nan, penalty_values)

fig = go.Figure()
fig.add_trace(go.Scatter(x=grid, y=values))
fig.add_trace(go.Scatter(x=grid, y=penalty_values))
fig.show()

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.

Related approaches in the literature

Hidden constraints have been studied mainly in derivative-free, simulation-based and Bayesian optimization. The main approaches are:

1. Extreme barrier. Set the objective to \(+\infty\) when the evaluation fails. This is the standard approach in direct-search methods such as mesh adaptive direct search [AD06]. It works because these methods only compare function values. Nelder–Mead and most evolutionary methods tolerate it as well. It violates requirement 1.

2. Rejecting the step inside the solver. Ipopt [WB05] treats a failed evaluation, or a non-finite function value, as a rejected trial step and cuts back the step length. Trust-region methods can likewise shrink the radius. This is arguably the cleanest solution, but it needs solver support and violates requirement 2. optimagic’s penalty can be seen as emulating this behavior from outside the solver.

3. Artificial function values. SNOBFIT [HN08] assigns an artificial value, based on the values at nearby points, to points where the evaluation failed. This is the closest relative of optimagic’s approach. SNOBFIT computes the artificial values inside the solver; optimagic uses a fixed function of the distance to the start parameters.

4. Modeling the failure region. Surrogate-based and Bayesian methods fit a classifier or Gaussian process that predicts whether an evaluation succeeds and take it into account when proposing new points [GSA14, GL11, MullerD19]. This pays off when evaluations are expensive, but requires a model inside the solver and violates requirement 2.

Limitations

• Model-based derivative-free optimizers such as pounders, DFO-LS or BOBYQA may include penalized points in their interpolation set, which distorts their model of the objective. Leaving failed points out of the model would be better, but this is not possible from outside the solver.

• For global and population-based optimizers, the slope towards the start parameters biases the search towards \(x_0\). For solvers that support it, \(+\infty\) would be the better choice.

• The boundary of the region where errors occur becomes an artificial discontinuity of the objective. For optimizers that only accept improving steps this is harmless in practice, but there is no convergence theory for it.

Choosing the strategy per algorithm, for example returning \(+\infty\) to solvers that handle it, is discussed in issue #717.

References

[AD06]

Charles Audet and John E. Dennis, Jr. Mesh adaptive direct search algorithms for constrained optimization. SIAM Journal on Optimization, 17(1):188–217, 2006. doi:10.1137/040603371.

[GSA14]

Michael A. Gelbart, Jasper Snoek, and Ryan P. Adams. Bayesian optimization with unknown constraints. In Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, 250–259. 2014.

[GL11]

Robert B. Gramacy and Herbert K. H. Lee. Optimization under unknown constraints. In Bayesian Statistics 9, pages 229–256. Oxford University Press, 2011.

[HN08]

Waltraud Huyer and Arnold Neumaier. SNOBFIT – stable noisy optimization by branch and fit. ACM Transactions on Mathematical Software, 2008. doi:10.1145/1377612.1377613.

[LDW24]

Sébastien Le Digabel and Stefan M. Wild. A taxonomy of constraints in black-box simulation-based optimization. Optimization and Engineering, 25:1125–1143, 2024. URL: https://arxiv.org/abs/1505.07881.

[MullerD19]

Juliane Müller and Marcus Day. Surrogate optimization of computationally expensive black-box problems with hidden constraints. INFORMS Journal on Computing, 31(4):689–702, 2019. doi:10.1287/ijoc.2018.0864.

[WB05]

Andreas Wächter and Lorenz T. Biegler. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1):25–57, apr 2005. doi:10.1007/s10107-004-0559-y.