# Classification-Based Linear Surrogate Modeling of Constraints for AL-CMA-ES

Oskar Girardin, Nikolaus Hansen, Dimo Brockhoff, Anne Auger

► To cite this version:

Oskar Girardin, Nikolaus Hansen, Dimo Brockhoff, Anne Auger. Classification-Based Linear Surrogate Modeling of Constraints for AL-CMA-ES. GECCO 2025 - Genetic and Evolutionary Computation Conference, Jul 2025, Malaga, Spain. pp.728-736, 10.1145/3712256.3726435. hal-05406428

HAL Id: hal-05406428

https://hal.science/hal-05406428v1

Submitted on 19 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Classification-Based Linear Surrogate Modeling of Constraints for AL-CMA-ES

Oskar Girardin

Nikolaus Hansen

Inria, CMAP, CNRS, École Polytechnique

Institut Polytechnique de Paris

91120 Palaiseau, France

firstname.lastname@inria.fr

Dimo Brockhoff

Anne Auger

Inria, CMAP, CNRS, École Polytechnique

Institut Polytechnique de Paris

91120 Palaiseau, France

firstname.lastname@inria.fr

## Abstract

We introduce linear surrogate functions for modeling inequality constraints to solve constrained blackbox optimization problems with the Augmented Lagrangian CMA-ES. Each surrogate is constructed from a binary classifier that predicts the sign of the constraint value. The classifier, and consequently the resulting algorithm, is invariant under sign preserving transformations of the constraint values and can handle binary, flat, and deceptive constraints. Somewhat surprisingly, we find that adopting a sign-based classification model of the constraints allows to solve classes of constrained problems which can not be solved with the original Augmented Lagrangian method using the true constraint value.

## ACM Reference Format:

Oskar Girardin, Nikolaus Hansen, Dimo Brockhoff, and Anne Auger. 2025. Classification-Based Linear Surrogate Modeling of Constraints for AL-CMA-ES. In *Genetic and Evolutionary Computation Conference (GECCO '25), July 14–18, 2025, Malaga, Spain*. ACM, New York, NY, USA, 9 pages. https://doi.org/10.1145/3712256.3726435

## 1 Introduction

Many real-world numerical optimization problems come with constraints. We consider in this paper the design of constraint handling mechanisms. The CMA-ES algorithm [10] is an adaptive stochastic search algorithm often recognized as state-of-the-art to address difficult numerical optimization problems [13, 16]. In the unconstrained setting, CMA-ES is affine-invariant, i.e., invariant to affine transformations of the search space, and invariant to strictly increasing transformations of the objective function [12]. Such invariance properties ensure generalization of performance and are certainly one reason for the success of CMA-ES.

The default handling of constraints within CMA-ES is an augmented Lagrangian approach which was introduced in the context of Evolution Strategies for plus selection [1] and later on introduced for comma selection and then incorporated within CMA-ES [2, 8]. In the augmented Lagrangian CMA-ES (AL-CMA-ES), candidate solutions are evaluated on an augmented Lagrangian function that balances the objective function and the constraints using Lagrangian multipliers and penalty parameters. These parameters are adapted and become new state variables of the algorithm. The

invariance to strictly increasing transformations of the objective function from CMA-ES is lost in AL-CMA-ES while the algorithm is still invariant to the affine rewriting of $f$ [4].

In this context, we build upon AL-CMA-ES and propose a new method that views $m$ constraints as a set of classification problems with two classes per constraint, feasible and infeasible. In the proposed ALLS-CMA-ES, we learn a linear classifier for each constraint. We use these classifiers to replace the original constraints with their linear surrogates in the AL-CMA-ES. The resulting algorithm, ALLS-CMA-ES, is invariant to sign-preserving transformations of the constraints. We evaluate and compare ALLS-CMA-ES on a set of test problems and compare it with AL-CMA-ES, COBYLA, SLSQP and another classical solver. We illustrate in particular the impact of applying sign-preserving transformation on the constraints. While our ultimate goal is to handle nonlinear constraints, we evaluate our approach on constraints with linear decision boundaries. Our tests still provide a relatively general assessment of the local convergence behavior, because the decision boundary of a nonlinear constraint is usually locally linear in the limit (at the optimum).

This paper is organized as follows. In Section 2 we present the problem statement and AL-CMA-ES. In Section 3 we show limitations of current approaches when applying sign-preserving transformations to the constraints. In Section 4 we present the ALLS-CMA-ES whose performance is evaluated in Section 5. Section 6 gives a summary and discusses some limitations. Additional material can be found at https://oskargirardin.github.io/Classification-Based-Linear-Surrogate-for-AL-CMA-ES.

## Symbols and Notations

$n \in \mathbb{N}^*$ - search space dimension

$g_k : \mathbb{R}^n \rightarrow \mathbb{R}^k$ - $k$-th constraint with linear boundary given by $\beta_k^\top x + \alpha_k = 0$

$g = (g_1, \dots, g_m)$ - vector valued constraints function

$s_k : \mathbb{R}^n \rightarrow \mathbb{R}$, $x \mapsto \operatorname{sgn}(g_k(x))$ - see Equations (10) and (11)

$\operatorname{sgn} : \mathbb{R} \rightarrow \mathbb{R}$, $y \mapsto 1$ if $y > 0$ else $-1$ - sign function

$\| \cdot \|$ - Euclidean norm

$\| \cdot \|_A$ - norm such that $\|x\|_A = x^\top A x$ and $A$ is a symmetric positive definite matrix

$\operatorname{iqr}$ - interquartile range of a set of values

$\operatorname{ker}$ - kernel (or null space) of a matrix

### CMA-ES.

$m \in \mathbb{R}^n$ - mean of the sample distribution

$P \in \mathbb{N}^*$ - population size

$x \in \mathbb{R}^n$ - single individual

$m \in \mathbb{N}$ - number of constraints

<figure><img src="image_2.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. GECCO '25, July 14–18, 2025, Malaga, Spain

© 2025 Copyright held by the owner/author(s).

ACM ISBN 979-8-4007-1465-8/2025/07

https://doi.org/10.1145/3712256.3726435
$C \in S_{++}^d$ - covariance matrix of CMA-ES

$\sigma \in \mathbb{R}_+$ - step-size

$C_\sigma = \sigma^2 C \in S_{++}^d$ - sample covariance matrix

Augmented Lagrangian.

$\lambda_k \in \mathbb{R}_+$ - Lagrangian coefficient for constraint $k$

$\lambda \in \mathbb{R}_+$ - vector of Lagrangian coefficients

$\mu_k \in \mathbb{R}_+$ - penalty factor for constraint $k$

$\mu \in \mathbb{R}_+$ - vector of penalty factors

$\chi \in \mathbb{R}_{>1}$ - controls the magnitude of changes in $\mu$

$d_\lambda \in \mathbb{R}_+$ - damping factor for changes in $\lambda$

$F_{n+5} \in \mathbb{R}^{m \times (n+5)}$ - feasibility sign matrix

$P_k^g \in [0, 1]$ - feasibility ratio of $k$-th constraint

Surrogate.

$\hat{g}_k : \mathbb{R}^n \to \mathbb{R}$ - linear surrogate for constraint $k$

$\hat{\beta}_k \in \mathbb{R}^n$ - coefficients of linear surrogate for constraint $k$

where $\|\hat{\beta}_k\| = 1$

$\hat{\alpha}_k \in \mathbb{R}$ - intercept of linear surrogate for constraint $k$

$\mathbf{B}_k \in \mathbb{R}^{n \times n}$ - points that characterize the surrogate constraint boundary with index $k$

$\bar{\mathbf{b}}_k \in \mathbb{R}^n$ - mean of points in $\mathbf{B}_k$

$\mathcal{P}_k = \{\mathbf{x} \in \mathbb{R}^n : \hat{\beta}_k^\top \mathbf{x} + \hat{\alpha}_k = 0\}$ - set of points that are included in the surrogate constraint boundary with index $k$

$\delta \in \mathbb{R}^+$ - norm of interval during interval halving

## 2 Background

We desire to minimize an objective function $f: \mathbb{R}^n \to \mathbb{R}$ such that all inequality constraints $g_k: \mathbb{R}^n \to \mathbb{R}, k = 1, \dots, m$ are satisfied. We call a point $\mathbf{x} \in \mathbb{R}^n$ feasible if all its $g_k$ values are nonpositive. To define the optimization goal, we consider the infimum $f$-value that has a feasible sublevel set with nonzero volume (i.e., with positive spatial suboptimality [9]), called essential infimum, formalized as

$$f^* = \underset{\{\mathbf{x} \in \mathbb{R}^n \mid g_k(\mathbf{x}) \le 0, k = 1, \dots, m\}}{\operatorname{essinf}} f(\mathbf{x}). \quad (1)$$

For continuous functions, the essential infimum coincides with their infimum and, if it exists, with their minimum $f$-value. The optimization goal is to find a feasible solution with an $f$-value close to $f^*$. We consider a black-box context and do not use gradients of $f$ or $g_k$ or assume that they exist. In particular, we investigate problems where the constraint is discontinuous (and hence nondifferentiable) at the optimum.

### 2.1 Related Literature

Augmented Lagrangian constraint handling has been investigated before in the context of Evolution Strategies. The approach has been used for the (1 + 1)-ES [1] and tested on linearly constrained ellipsoid functions. The method was generalized to an Evolution Strategy with covariance matrix adaptation [2] and multiple (linear) constraints [3]. The authors analyze the stability of the underlying Markov Chain and provide results that suggest linear convergence.

More recently, linear surrogate models for the constraints were studied within AL-CMA-ES [8] using linear regression. This allowed to handle constraints with a flat feasible part (as in ReLu).

This paper constructs a linear model of the constraints by using only the feasibility classification rather than the actual constraint values. Hence, in contrast to previous works, the algorithm is invariant to sign-preserving transformations of the constraints.

### 2.2 AL-CMA-ES

The CMA-ES [10] is an adaptive randomized optimization algorithm introduced originally for unconstrained black-box optimization. The algorithm iteratively adapts a multivariate Gaussian distribution. At iteration $t$, this distribution is parameterized by a mean vector $\mathbf{m}^{(t)} \in \mathbb{R}^n$ representing the current favorite solution (or incumbent) and a sample covariance matrix $C_\sigma^{(t)}$ that encodes a metric based on which candidate solutions are sampled. The sample covariance matrix is decomposed as $C_\sigma^{(t)} = [\sigma^{(t)}]^2 C^{(t)}$ into a step-size and a symmetric positive definite covariance matrix. The state variables of CMA-ES naturally include the parameters of the multivariate Gaussian distribution, namely the mean vector $\mathbf{m}^{(t)}$, step-size $\sigma^{(t)}$ and covariance matrix $C^{(t)}$ but also two additional path vectors representing directions of momentum. We denote all state-variables as $\theta^{(t)}$. The update of the state variable at each iteration of CMA-ES reads in an abstract way as

$$\theta^{(t+1)} = CM\mathcal{A}_f(\theta^{(t)}, U^{(t+1)})$$

where $U^{(t+1)}$ contains all random input, sampled at iteration $t$. The update can be decomposed into three main steps. First, given the current state $\theta^{(t)}$ and some random input $U^{(t+1)}$, a so-called ask function [7] samples $P$ candidate solutions according to the multivariate Gaussian distribution with mean $\mathbf{m}^{(t)}$ and covariance matrix $C_\sigma^{(t)}$,

$$(\mathbf{x}_1^{(t+1)}, \dots, \mathbf{x}_P^{(t+1)}) = \text{Sample}(\theta^{(t)}, U^{(t+1)}). \quad (2)$$

These candidate solutions are then evaluated on $f$ and ranked according to their $f$-values

$$f(\mathbf{x}_{\pi(1)}^{(t+1)}) \le \dots \le f(\mathbf{x}_{\pi(P)}^{(t+1)}) \quad (3)$$

where $\pi$ is a permutation such that $\pi^{-1}(k)$ is the $f$-rank of the $k$-th candidate solution. From the current state $\theta^{(t)}$ and the ranked candidate solutions, the new state is computed as

$$\theta^{(t+1)} = CM\mathcal{A}_{\text{update}}\left(\theta^{(t)}, \mathbf{x}_{\pi(1)}^{(t+1)}, \dots, \mathbf{x}_{\pi(P)}^{(t+1)}\right). \quad (4)$$

We refer to [12] for details on how the update of $\theta^{(t)}$ is performed.

The default constraint handling technique for CMA-ES [11] is an augmented Lagrangian approach. The augmented Lagrangian function aggregates the objective function and the constraints,

$$\mathcal{L}_g(\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}) := f(\mathbf{x}) + \sum_{k=1}^{m} \begin{cases} \lambda_k g_k(\mathbf{x}) + \frac{\mu_k}{2} g_k^2(\mathbf{x}) & \text{if } \mu_k g_k(\mathbf{x}) \ge -\lambda_k \\ -\frac{\lambda_k^2}{2\mu_k} & \text{otherwise} \end{cases} \quad (5)$$

where $\lambda_k \in \mathbb{R}^m$ are the Lagrange multipliers and $\mu_k \in \mathbb{R}^m$ are penalty parameters. The $f$-landscape remains unchanged when $g_k(\mathbf{x}) \le -\lambda_k/\mu_k$ for all constraints $k$. This approach was proposed for the (1+1)-ES [1] and adopted for nonelitist Evolution Strategies and CMA-ES [2, 3].
**Algorithm 1 AL(LS)-CMA-ES outer optimization loop**

**Require:** $\mathcal{L}_g$ and $\hat{g}$ (if model)

1 initialize all CMA-ES parameters $\theta^{(0)}$

2 $\lambda^{(0)}, \mu^{(0)} = \text{init\_al()}$ ▶ initialize state variables for Algorithm 2 (AL update) and $\mathcal{L}_g$

3 **if model then**

4 set $\hat{\beta}_k = 0, \hat{\alpha}_k = 0 \forall k = 1, \dots, m$

5 set $g$ as $\hat{g}$

6 **while** not happy **do**

7 sample population $\mathbf{x}_1^{(t)}, \dots, \mathbf{x}_P^{(t)} \sim N(\mathbf{m}^{(t)}, \mathbf{C}_\sigma^{(t)})$

8 set $\pi$ such that it satisfies the $\mathcal{L}_g(\cdot, \lambda^{(t)}, \mu^{(t)})$ ranking of (6)

9 $\theta^{(t+1)} = CM\mathcal{A}_{\text{update}}(\theta^{(t)}, \mathbf{x}_{\pi(1)}^{(t+1)}, \dots, \mathbf{x}_{\pi(P)}^{(t+1)})$

10 **if** model **then**

11 **for** $k = 1, \dots, m$ **do**

12 $\hat{\beta}_k, \hat{\alpha}_k \leftarrow \text{model.update}(\mathbf{m}^{(t+1)}, \mathbf{C}_\sigma^{(t+1)}, s_k, \hat{\beta}_k, \hat{\alpha}_k)$

13 $\lambda^{(t+1)}, \mu^{(t+1)} \leftarrow \text{update\_al}(g(\mathbf{m}^{(t+1)}), \lambda^{(t)}, \mu^{(t)}) \triangleright \text{Algorithm 2}$

14 **return** $\mathbf{m}$

The AL-CMA-ES algorithm builds on CMA-ES with all its state variables $\theta^{(t)}$ and adds the state variables $\lambda^{(t)}$ and $\mu^{(t)}$ for the adaptation of the augmented Lagrangian function. Given the state $(\theta^{(t)}, \lambda^{(t)}, \mu^{(t)})$ at iteration $t$, the first step of AL-CMA-ES is the same as in CMA-ES and consists in sampling candidate solutions according to (2). Then, the augmented Lagrangian function is used to determine the permutation $\pi$ by

$$ \mathcal{L}_g \left( \mathbf{x}_{\pi(1)}^{(t+1)}, \lambda^{(t)}, \mu^{(t)} \right) \le \dots \le \mathcal{L}_g \left( \mathbf{x}_{\pi(P)}^{(t+1)}, \lambda^{(t)}, \mu^{(t)} \right). \quad (6) $$

With this ranking permutation, the state variables of CMA-ES are updated according to (4). In the last step, AL-CMA-ES updates the parameters $\lambda^{(t)}$ and $\mu^{(t)}$ according to Algorithm 2. The update depends on the constraint values of the new mean vector, $g(\mathbf{m}^{(t+1)})$, on the previous values of $(\lambda^{(t)}, \mu^{(t)})$, and on the feasibility ratios of the last $n+5$ mean values for each constraint, later denoted as $p_k^g = \frac{1}{t+2-\max(0,t-n-3)} \sum_{i=\max(0,t-n-3)}^{t+1} 1[g_k(\mathbf{m}^{(i)}) \le 0]$. Let $F_{n+5}^{(t+1)}$ denote the $m \times (n+5)$ binary matrix $[s_k(\mathbf{m}^{(j)})]_{k,j}$, $k = 1, \dots, m, j = t-n-3, \dots, t+1$, indicating the feasibility sign values of the $n+5$ previous sample means. The update of $(\lambda^{(t)}, \mu^{(t)})$ can be written with a transition function $\mathcal{AL}$ which is independent of any other time-changing variables as

$$ (\lambda^{(t+1)}, \mu^{(t+1)}, F_{n+5}^{(t+1)}) = \mathcal{AL}(g(\mathbf{m}^{(t+1)}), F_{n+5}^{(t)}, \lambda^{(t)}, \mu^{(t)}). \quad (7) $$

This will be particularly useful to formalize invariance properties.

The pseudocode for the AL-CMA-ES is given in Algorithm 1. The AL update in Algorithm 2 describes the default setting (algorithm == 3) from GitHub/CMA-ES/pycma [11]. We used the code with commit tag v4.1.0 for this paper. In contrast to previous updates [4, 8], after initialization, our update\_al procedure in Algorithm 2 does not depend on the objective function $f$ and uses predominantly feasibility ratios. We use the default termination criteria but disable tolstagnation.

**Algorithm 2 AL update**

1 **procedure** setup\_al():

2 $\chi = 2^{\frac{1}{5n}}, d_\lambda = 5$

3 **for all** $k = 1, \dots, m$ **do**

4 $\lambda_k^{(0)} = \mu_k^{(0)} \leftarrow 0$

5 **procedure** init\_al():

6 **for all** $k = 1, \dots, m$ **do**

7 **if** $\mu_k^{(t)} = 0$ and more than 10% of the sampled solutions are infeasible **then**

8 $\lambda_k^{(t)} = \frac{\text{iqr}(f)}{(n \times \text{iqr}(g_k))}$

9 $\mu_k^{(t)} = \frac{2 \text{iqr}(f)}{5n(\text{iqr}(g_k) + \text{iqr}(g_k^2))}$

10 **procedure** update\_al((g₁, ..., gₘ), $\lambda^{(t)}, \mu^{(t)}$):

11 $\text{init\_al()}$

12 **for** $k = 1, \dots, m$ **do**

13 **if** $\mu_k^{(t)} = 0$ **then**

14 $\lambda_k^{(t+1)} = \mu_k^{(t+1)} \leftarrow 0$

15 **continue**

16 $\lambda_k^{(t+1)} = \max[0, \lambda_k^{(t)} + \frac{\mu_k^{(t)}}{d_\lambda} g_k]$

17 $p_k^g \leftarrow$ feasibility ratio of last $n+5$ passed $g_k$ values

18 **if** $g_k > 0$ **or** $[g_k > \frac{-\lambda_k^{(t+1)}}{\mu_k^{(t)}} \text{ and } p_k^g > 0.95]$ **then**

19 $\mu_k^{(t+1)} \leftarrow \mu_k^{(t)} \chi^{1/4}$

▸ increase

20 **if** $g_k \le 0$ **and** $p_k^g < 0.9$ **then**

21 $\mu_k^{(t+1)} \leftarrow \mu_k^{(t)} \chi^{-1/4}$

▸ decrease

22 **return** $\lambda^{(t+1)}, \mu^{(t+1)}$

## 3 Limitation Cases of Constraints Handling

In this paper, we show results on variations of the ellipsoid problem,

$$ \begin{align} \text{f}_\text{elli} : \mathbb{R}^n &\rightarrow \mathbb{R}, \quad \mathbf{x} \mapsto (\mathbf{x} - \mathbf{e})^\top \mathbf{A}_a (\mathbf{x} - \mathbf{e}) \tag{8} \\ \mathbf{g}_k : \mathbb{R}^n &\rightarrow \mathbb{R}, \quad \mathbf{x} \mapsto x_k \quad k = 1, \dots, m \quad (m \le n) \nonumber \end{align} $$

where $\mathbf{e} \in \mathbb{R}^n$, $e_k = 1$ if $k \le m$ and 0 otherwise, and $A_a = \text{diag}[a_{n-i}, i = 1, \dots, n]$ where $a \in \{1, 10^6\}$ is the condition number. The unconstrained minimizer of $f_\text{elli}$ is $\mathbf{e}$, the constrained minimizer is $\mathbf{0}$. The $k$-th constraint only affects the $k$-th component of $\mathbf{x}$. We also use the rotated ellipsoid $\mathbf{x} \mapsto f_\text{elli}(\mathbf{R}\mathbf{x})$ subject to $g_1(\mathbf{R}\mathbf{x}), \dots, g_k(\mathbf{R}\mathbf{x}) \le 0$ where $\mathbf{R}$ is an $n \times n$ orthogonal matrix with its rows uniformly distributed on the $n-1$-dimensional hypersphere, fixed for each $(n, m)$ in all experiments.

Constraint handling algorithms often utilize the information encoded in the constraint value as a guide towards feasible solutions. This makes them sensitive to constraints that have non-informative, noisy, or misleading values. To assess such limitations of existing algorithms, we defined the $g$-value transformations shown in Figure 1 (their expressions are given in the caption). While the transformations change the constraint value, they leave the sign unaffected and thus do not change the feasible set.

In this section, we show results of four existing algorithms (in comparison to ALLS-CMA-ES as outlined later in Section 4) on the sphere function with one constraint, that is, problem (8) with

$$ n = 4, m = 1, a = 1 \tag{9} $$
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Sign-preserving transformations of a constraint function. The Step function equals sgn as given in (11), *LinearStep*: $y \mapsto y$ if $y \le 0$ else $y + 1$, *ReLu*: $y \mapsto y$ if $y > 0$ else 0, *Noisy*: $y \mapsto y \exp(0.1N(0, 1))$, *Deceiving*: $y \mapsto \operatorname{sgn}(y) \exp(-y \operatorname{sgn}(y))$.</figcaption></figure>

where we replace the single constraint $g_1: x \mapsto x_1$ with $x \mapsto T(x_1)$ where $T$ is one of the transformations presented in Figure 1.

*Experimental setup and algorithms.* The tested algorithms are 1) COBYLA [15], 2) SLSQP [14], 3) A trust-region interior point method described in [5] and 4) AL-CMA-ES described in Algorithm 2. For the first three algorithms we used their Python implementation in the library `scipy` [19] (version 1.11.1). For any given run, we sample a feasible initial point for all algorithms and use the default stopping conditions. The initial distribution is $x_0 = (U(-1, 0), U(0, 1), U(0, 1), U(0, 1))^\top$ where $U(a, b)$ is the uniform distribution on the interval $[a, b]$.

Figure 2 shows 11 runs of the algorithms in 4-D on the constrained problem described in (9) with the different constraint transformations depicted in Figure 1. Shown is the distance to the optimum of each evaluated point versus the number of constraint evaluations. The ALLS-CMA-ES is invariant under the presented transformations, hence, it behaves the same in all transformed scenarios up to stochastic variations. We discuss our observations for each transformation in turn.

*No Transformation.* All five algorithms lower the distance to the optimum with increasing iterations reaching a distance of $10^{-6}$,

except for COBYLA where runs stop at around $10^{-4}$. COBYLA, trust-region and SLSQP need typically far fewer evaluations compared to AL-CMA-ES. For example, to reach a distance of $10^{-4}$, they need less than 100 evaluations, while AL-CMA-ES needs between about 600 and 1000 evaluations and ALLS-CMA-ES needs typically around 1000 evaluations.

*LinearStep.* SLSQP and trust-region are not visibly affected by the transformation whereas COBYLA fails to get values below $10^{-1}$. The AL-CMA-ES needs about twice as long as without transformation to get close to the optimum.

*ReLu.* COBYLA and SLSQP are not visibly affected by the transformation whereas the trust-region algorithm shows a highly variable behaviour after around 50 iterations. Most of the runs stagnate at a distance of $10^{-2}$ or $10^{-3}$. For AL-CMA-ES, some runs take about 50% more evaluations than without transformation to reach a distance to the optimum of below $10^{-6}$.

*Step.* COBYLA, SLSQP and trust-region are strongly affected by the *Step* transformation and show only a small initial improvement after which they stall. The AL-CMA-ES requires almost about twice the number of evaluations compared to no transformation.

*Noisy.* SLSQP and trust-region are strongly affected by the *Noisy* transformation and do not notably improve the initial solution. In contrast, COBYLA performs virtually indistinguishable compared to no transformation. The AL-CMA-ES converges much slower than without transformation and does not reach a final distance to the optimum below $10^{-3}$.

*Deceiving.* All algorithms but ALLS-CMA-ES are strongly affected by this transformation. Consistent with the deceiving nature of the constraint, COBYLA, SLSQP and trust-region *increase* the distance to the optimum. For AL-CMA-ES, some runs converge within 2000 evaluations to the optimum whereas in some runs, the distance to the optimum increases progressively.

*Interpretation.* Most tested algorithms are sensitive to at least some of the sign-preserving transformations of the constraint *value* of a single constraint with a linear feasibility decision boundary even on the 4-dimensional sphere function. Remarkably, different algorithms are sensitive to different transformations. COBYLA performs well on the *Noisy* transformation which is challenging for the other algorithms. The AL-CMA-ES algorithm fails only on the *Deceiving* and *Noisy* transformations and stands out by solving the *Step* transformation. However, on a 20-dimensional rotated ellipsoid with three or more constraints (see below, Figure 4), AL-CMA-ES is also affected by the other transformations. As illustrated before and as we will prove below, ALLS-CMA-ES is by design unaffected by any sign-preserving transformation of the constraint values.

# 4 ALLS-CMA-ES

To address the limitations observed in the previous section, we interpret the inequality constraint, $g_k$, as a classifier with two classes, "feasible" and "infeasible". To create a surrogate model for the constraint $k$, we build and iteratively update a binary classifier with a linear decision boundary, which in turn represents a linear constraint function, $\hat{g}_k$. We consider a linear decision boundary to be
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2: Distance to the optimum versus g-evaluations of COBYLA, SLSQP, trust-region, AL-CMA-ES, and ALLS-CMA-ES on the 4-dimensional sphere function defined in (9) with different transformations of a single linear constraint as shown in Figure 1. The run with the median final f-value is annotated and plotted thick. Without transformation, all algorithms drop below $10^{-3}$, however of the 5 transformations, they typically drop below $10^{-3}$ (only) on 2, 2, 1, 3, and 5 of them, respectively.</figcaption></figure>

reasonably general, because when approaching a local (or global) optimum, the decision boundary of any smooth nonlinear constraint is locally linear in the limit. Even a value jump at the boundary, like with the *Step* or *LinearStep* transitions from Figure 1, does not change the contour set of the decision boundary in x-space.

The classifier is based on sign information only, namely a sequence of pairs $(x, s_k(x)) \in \mathbb{R}^n \times \{-1, 1\}$, where

$$s_k : x \in \mathbb{R}^n \mapsto \operatorname{sgn}(g_k(x)) \qquad (10)$$

$$\operatorname{sgn} : y \in \mathbb{R} \mapsto \begin{cases} 1 & \text{if } y > 0 \text{ (i.e., infeasible)} \\ -1 & \text{otherwise} \end{cases} . \qquad (11)$$

Thereby, we disregard the information given by the exact constraint value $g_k(x)$ for constructing the surrogate model $\hat{g}_k$.

The surrogate values are computed as

$$\hat{g}_k : \mathbb{R}^n \to \mathbb{R}, \quad x \mapsto \hat{\beta}_k^\top x + \hat{\alpha}_k \qquad (12)$$

with $\hat{\beta}_k \in \mathbb{R}^n$, $\hat{\alpha}_k \in \mathbb{R}$ and $\|\hat{\beta}_k\| = 1$. Internally, the linear model is represented by $n$ points on the $n-1$-dimensional decision boundary, denoted as $B_k \in \mathbb{R}^{n \times n}$. The updates of $B_k$, $\hat{\beta}_k$, $\hat{\alpha}_k$ are given in Algorithms 3 and 4. The initial $B_k$ is created at the cost of a single line search when the mean passes the feasibility boundary for the first time: Algorithm 4 is applied with $a = m^{(t)}$, $b = m^{(t+1)}$, $N_{\text{expanded}} = 0$ to find the first element of $B_k$ which is assigned to be “the newest”. The mean shift direction is considered as the gradient of $\hat{g}_k$ thereby defining the decision boundary. The $n-1$ remaining elements of $B_k$ are sampled from the current sample distribution (2) and projected onto the decision boundary. Further updates replace the oldest point in $B_k$ with a new point determined by a line search starting from two points close to a projection of $m^{(t+1)}$ onto the decision boundary. The mapping $B_k \mapsto (\hat{\beta}_k, \hat{\alpha}_k)$ for (12) is given in lines 15 and 16 of Algorithm 3. This approach resembles active learning [6], as the line search “actively” determines the samples
**Algorithm 3 Update of k-th surrogate**

**Require:** $m, C_σ, s_k$ and previous model parameters $\hat{\beta}_k, \hat{\alpha}_k, B_k$ and largest final precision $\delta_{max}$ among points in $B_k$

1. Compute the projection of the mean onto the hyperplane $P_k = \{x \in \mathbb{R}^n : \hat{\beta}_k^\top x + \hat{\alpha}_k = 0\}$:

$$m_p = \arg \min_{x \in P_k} \|x - m\|_{C_\sigma^{-1}} = m - \frac{\hat{\beta}_k^\top m + \hat{\alpha}_k}{\hat{\beta}_k^\top C_\sigma \hat{\beta}_k} C_\sigma \hat{\beta}_k$$

2. $v \leftarrow (m_p - m)/(\|m_p - m\|)$

3. **if** # updates $\ge n$ **then**

4. $x_0, x_1 \leftarrow m_p \pm 2 \times \delta_{max} \times v$

5. **else**

6. $x_0, x_1 \leftarrow m, m_p$

► $B_k$ still contains "low-quality" points

7. $i, j \leftarrow 0, x_2 \leftarrow x_{1-i}$

8. **while** $x_0$ and $x_1$ have the same feasibility do

9. $j \leftarrow j + 1$

► exploring iterations

10. $i \leftarrow 1 - i$

► i remains in $\{0, 1\}$

11. $x_2 \leftarrow x_i$

► save previous $x_i$

12. $x_i \leftarrow x_i + (x_i - x_{1-i})$

► expand further

13. $x_0, x_1 \leftarrow \text{IntervalHalving}(x_i, x_2, s_k, j)$

14. Replace oldest point and its associated precision in $B_k$ with $(x_0 + x_1)/2$ and $\|x_1 - x_0\|$

15. $\hat{\beta} \leftarrow \ker(B_k - \bar{b}_k 1_n)$

► computed with SVD

16. $\hat{\alpha}_k \leftarrow -\hat{\beta}_k^\top \bar{b}_k$

17. **return** $\hat{\beta}_k, \hat{\alpha}_k$

**Algorithm 4 IntervalHalving**

**Require:** Two points $a, b \in \mathbb{R}^n$ of opposite feasibility, feasibility function $s$, additional iterations $N_{expanded}$

1. **for** $N_{expanded} + 6$ **do**

2. **if** $\|a - b\| / (\|a\| + \|b\|) \le 10^{-15}$ **then**

3. **break**

► numerical safeguard

4. $c = (a + b)/2$

5. **if** $s(c) = s(a)$ **then**

6. $a \leftarrow c$

7. **else**

8. $b \leftarrow c$

9. **return** $a, b$

► $a$ and $b$ still have opposite feasibility

used to build the model. Our first attempts to use logistic regression [18] or support vector machines [17] to build the linear model were discouraging. We mainly attributed this to model fluctuations from the lack of data close enough to the decision boundary.

The linear surrogates are coupled with AL-CMA-ES as shown in Algorithm 1. We set $g$ as $\hat{g}$ in Equations (5)–(7) and hence use $\mathcal{L}_{\hat{g}(t)}$ instead of $\mathcal{L}_{g(t)}$ to determine the ranking $\pi$. After updating $\theta$ in (4), we update the model parameters,

$$\hat{g}_k^{(t+1)} = \mathcal{SL}_{s_k} \left( m^{(t+1)}, C_\sigma^{(t+1)}, \hat{g}_k^{(t)} \right) \quad (13)$$

where $\mathcal{SL}$ denotes the transition function for the model according to Algorithm 3. Finally, we update $\mu_k$ and $\lambda_k$ with (5) using $\hat{g}_k$ in place of $g_k$.

## 4.1 Invariance properties

We show that ALLS-CMA-ES exhibits two important invariance properties: invariance to sign-preserving transformations of the constraints $g$ and invariance to affine transformations of the objective function $f$. Informally, an algorithm is invariant to a class of transformations of the problem if it performs the same on any two problems from the invariance class. Formally, given a stochastic algorithm optimizing a constrained problem $x \in \mathbb{R}^n \mapsto (f(x), g(x)) \in \mathbb{R} \times \mathbb{R}^m$, we say that an algorithm with state $\psi^{(t)}$ and update function

$$\psi^{(t+1)} = \mathcal{A}_{(f,g)}(\psi^{(t)}, U^{(t+1)})$$

(where $\{U^{(t)}, t \in \mathbb{N}\}$ represents the independent random input at each iteration) is invariant to a problem transformation $T(f, g)$ if we can find a bijective state-space transformation $\Phi$ such that

$$\Phi^{-1} \left( \mathcal{A}_{T(f,g)}(\Phi(\psi^{(t)}), U^{(t+1)}) \right) = \mathcal{A}_{(f,g)}(\psi^{(t)}, U^{(t+1)}) . \quad (14)$$

The previous equation shows that up to changing the states via the bijective transformation $\Phi$, running the algorithm on $(f, g)$ is the same as running it on the transformed problem $T(f, g)$.

When the state-space transformation $\Phi$ is the identity (i.e., we do not need to apply any change of state) and hence does not depend on $T$, we say that the algorithm is unconditionally invariant. The CMA-ES algorithm is unconditionally invariant to any strictly increasing transformation of the objective function, i.e. to any $\gamma \circ f$ where $\gamma : \text{Im}(f) \to \mathbb{R}$ is strictly increasing. Yet, this invariance is lost in the augmented Lagrangian CMA-ES [4] since its ranking (6) depends on the weighted sum of $f$ and the constraints $g$ in (5). While in ALLS-CMA-ES we do not recover this invariance, we however obtain invariance to sign-preserving transformations of the constraints as by design, ALLS-CMA-ES uses only the sign of the constraints. Hence, a sign-preserving transformation of the constraints (which does not change the feasible set) has no effect on the update of the state variables. More formally, we call a function $\varphi : \mathbb{R} \to \mathbb{R}$ sign-preserving if:

$$\varphi(x) > 0 \iff x > 0 . \quad (15)$$

**PROPOSITION 4.1.** *The ALLS-CMA-ES algorithm is unconditionally invariant to sign-preserving transformations of the constraints. In other words, assuming the same initialization and random input, the algorithm will produce the same sequence of states when optimizing a constrained problem $(f, g_1, \dots, g_m)$ or any $(f, \varphi_1(g_1), \dots, \varphi_m(g_m))$ where $\varphi_i$ are sign-preserving transformations.*

**PROOF.** We observe that in the equations defining the update of the state of ALLS-CMA-ES, namely (2), (4), (5)–(7) with $\hat{g}$ instead of $g$, and (13), the constraints information enters in the algorithm via the functions $s_k(x) = \text{sgn}(g_k(x))$ or via $\hat{g}_k$ which only depend on $s_k$ (see Algorithm 3). Yet since for any sign preserving transformation $\varphi_k$, $\text{sgn}(\varphi_k(g_k(x))) = \text{sgn}(g_k(x))$, the different updates are invariant to sign-preserving transformations. $\square$

While the AL-CMA-ES algorithm is not invariant to strictly increasing transformations of the objective function (in contrast to CMA-ES), it is invariant to affine rewriting of $f$, i.e. there exists a state-space transformation such that, up to this transformation, the algorithm sequence is the same when optimizing $f$ and $\tau_1 f + \tau_2$ for
<figure><img src="image_6.png" /><figcaption>Figure 3: Typical single run of ALLS-CMA-ES on the 20-dimensional ellipsoid function from (8) with $m = 1$ constraint and $a = 10^6$. Left: shown are the step size $\sigma$, the axis ratio (ratio of the axis lengths) of the covariance matrix $C$, the smallest and largest standard deviations of the sampling distribution and the Euclidean distance of the mean to the optimum. Right: The top plot shows the Euclidean distance of the surrogate model parameters $\hat{\beta}_k$ and $\hat{\alpha}_k$ to their respective true values, the estimation error of the Lagrange coefficient $|\lambda_k - \lambda_k^*|$, the ratio between the second smallest and largest singular value of the empirical covariance matrix of the points in $B_1$ (the smallest singular value is always zero), denoted as diversity, and the penalty factor $\mu_1$. The bottom plot shows the number of line search iterations, the number of exploring iterations thereof, and the classification accuracy of the surrogate model for the current population.</figcaption></figure>

any $\tau_1 > 0$ and $\tau_2 \in \mathbb{R}$ [4]. We prove in the next proposition that this invariance is preserved in ALLS-CMA-ES.

**PROPOSITION 4.2.** The ALLS-CMA-ES algorithm is invariant to affine rewriting of the objective function. Consider $f$ and $\tau_1 f + \tau_2$ for any $\tau_1 > 0$ and $\tau_2 \in \mathbb{R}$, the equivariance relation (14) is satisfied with the state-space transformation $\lambda' = \tau_1 \lambda$ and $\mu' = \tau_1 \mu$ and all other state variables unchanged.

**PROOF.** We consider the optimization of $f$ and $\tau_1 f + \tau_2$ and denote the state variable of the algorithm optimizing $\tau_1 f + \tau_2$ with a prime. We first prove that the selection of the candidate solutions is the same provided $\lambda^{(t)} = \tau_1 \lambda^{(t)}$ and $\mu^{(t)} = \tau_1 \mu^{(t)}$. We observe that

$$ \tau_1 \mathcal{L}_{(f,g)}(x, \lambda, \mu) + \tau_2 = \mathcal{L}_{(\tau_1 f + \tau_2, \hat{g})}(x, \tau_1 \lambda, \tau_1 \mu) $$

which implies that the selection of candidate solutions will be the same on $(f, \hat{g}^{(t)})$ using $(\lambda^{(t)}, \mu^{(t)})$ and $(\tau_1 f + \tau_2, \hat{g}^{(t)})$ using $(\lambda^{(t)}, \mu^{(t)})$. Hence the update of the CMA-ES state in (4) will be the same on $f$ or $\tau_1 f + \tau_2$ provided $\theta^{(t)} = \theta^{(t)}'$. Equation (13) tells us that the surrogate building is invariant (since we input the same state variables). Last, we need to verify the equivariance relation for the update (7) with $\hat{g}$ instead of $g$. Remark that the matrix $F_{n+5}^{(t)}$ is invariant since it contains the information whether previous means are feasible or not and the means are unchanged and the surrogate as well. Now if $\lambda^{(t)}' = \tau_1 \lambda^{(t)}$ and $\mu^{(t)}' = \tau_1 \mu^{(t)}$ with all other state variables identical, we remark that the AL-CMA-ES update satisfies $\lambda^{(t+1)}' = \tau_1 \lambda^{(t+1)}$ and $\mu^{(t+1)}' = \tau_1 \mu^{(t+1)}$ (this can be easily verified looking at each update). □

# 5 Empirical Results

To give insights into the functioning of ALLS-CMA-ES, Figure 3 shows its internal state variables when run on the 20-dimensional ellipsoid function $f_{elli}$ (8) with condition number $a = 10^6$ and a single constraint. All shown states, except for the axis ratio, tend to decrease with increasing iteration. The distance to the optimum $||m - x^*||$ reaches $10^{-5}$ after around 1200 iterations and stagnates shortly after (which we interpret as due to “numerical noise”). During the adaptation phase of the covariance matrix over the first 900 iterations, the distance to the optimum changes only by a factor of roughly $\exp(-0.05/n)$ per iteration, whereas after iteration 1000 the final convergence factor is roughly $\exp(-1/n)$.

The ratio between largest and smallest standard deviation of the sampling distribution steadily increases and reaches around $10^3$ reflecting the square root of the condition number of the ellipsoid function. After $C$ is adapted, both standard deviations decrease quickly and with the same rate as the distance to the optimum.

The Lagrange coefficient $\lambda_1$ (right subfigure) takes about 500 iterations to get moderately close to its correct value $2 \times 10^6$. After that, the distribution mean (left subfigure) assumes also feasible values (blue color). The final relative precision of $\lambda_1$ is below $10^{-8}$. On the 20-dimensional ellipsoid, the comparatively slow adaptation of $\lambda_1$ is not a performance limiting factor as the adaptation of the covariance matrix takes about twice as many iterations.

The surrogate model parameters $\hat{\alpha}$ and $\hat{\beta}$ linearly converge towards their true values reaching a distance of $10^{-12}$ after about 250 iterations (right subfigure). From iteration 500, the evolution of $\hat{\alpha}$ is similar to the evolution of the standard deviations (left subfigure) with another convergence phase after iteration 800, reaching a final error of about $10^{-22}$. The typical number of total line search
<figure><img src="image_7.png" /><figcaption>Figure 4: Distance to the optimum for ALLS-CMA-ES and AL-CMA-ES versus the number of $f$ evaluations plus the average number of constraints evaluations per iteration with 3 and 8 constraints from 7 runs (4-D) and 5 runs (20-D) on the rotated ellipsoid function (8). From left to right: *No Transformation, LinearStep, ReLu*, and *Step* transformations from Figure 1. Top: $n = 4, m = 3$, middle: $n = 20, m = 3$, bottom: $n = 20, m = 8$. Crosses depict the final distance of a run, for long runs they are drawn at the right axis limit.</figcaption></figure>

iterations drops significantly after the first convergence phase at around iteration 250 and assumes a median value of (only) three. The number of exploring iterations (line 9 in Algorithm 3) averages to the low value of 0.8 throughout the run.

Figure 4 shows the evolution of the distance to the optimum of ALLS-CMA-ES in comparison to AL-CMA-ES on the rotated ill-conditioned ellipsoid (8) with $a = 10^6$ and varying $n, m$, and transformation—from left to right without transformation and with *LinearStep*, *ReLu* and *Step* transformation for all constraints (see Figure 1). The initial point is drawn from the multivariate normal distribution $N(0, 2I)$ until we find a feasible point. Runs of ALLS-CMA-ES are stopped beyond the displayed range after 1000, 4000 and 5000 iterations, respectively for each row. Both algorithms are invariant under rotation of the search space. As above, we interpret the observed stagnation of the distance to the optimum below $10^{-5}$ as “numerical noise”.

Figure 4 shows results for $n = 4$ with $m = 3$ in the top row and for $n = 20$ with $m = 3$ (middle) and $m = 8$ (bottom row). Without transformation of the constraints (left column), ALLS-CMA-ES performs similar to AL-CMA-ES—it is slightly slower in dimension 4 and slightly faster in dimension 20. The ALLS-CMA-ES is, by design, unaffected by any transformation. Its dependency on the number of constraints is here negligible too. In contrast to the setting in Figure 2, AL-CMA-ES performs poorly on all transformations.

# 6 Summary and Limitations

We presented ALLS-CMA-ES which builds on AL-CMA-ES and considers the feasibility decision boundary of a constraint for defining a classification problem. A linear classifier is built based on repeated line searches and a linear surrogate model of each constraint is then used in AL-CMA-ES. We prove that the resulting algorithm gains invariance to sign-preserving transformations of the constraints while retaining the invariance to affine rewriting of the objective function from AL-CMA-ES. We also introduce a new update of the Lagrange parameters of AL-CMA-ES which is mainly based on observed feasibility ratios and which is independent of the objective function $f$. We illustrate the performance of ALLS-CMA-ES on different test problems with constraints with linear decision boundaries and compare it to AL-CMA-ES and three other numerical solvers from classical optimization. In the studied scenarios, ALLS-CMA-ES is by far the least susceptible to failure.

Our systematic experimental assessment was limited to orthogonal constraints with linear decision boundaries. First experiments with nonlinear decision boundaries were encouraging (not shown), however they also indicated that building a linear model for nonlinear decision boundaries reliably and efficiently requires some algorithmic modifications.

The presented algorithm generates the candidate solutions for the evaluation of the constraints *independently for each constraint*. This is a disadvantage when the cost to evaluate all constraints at the same time is considerably smaller than the combined cost of evaluating each constraint individually.
## References

[1] Dirk V Arnold and Jeremy Porter. 2015. Towards an augmented Lagrangian constraint handling approach for the (1+1)-ES. In *Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation*. MIT Press, Cambridge, MA, USA, 249–256.

[2] Asma Atamna, Anne Auger, and Nikolaus Hansen. 2016. Augmented Lagrangian constraint handling for CMA-ES—Case of a single linear constraint. In *International Conference on Parallel Problem Solving from Nature*. Springer, Berlin, Heidelberg, Germany, 181–191.

[3] Asma Atamna, Anne Auger, and Nikolaus Hansen. 2017. Linearly convergent evolution strategies via augmented Lagrangian constraint handling. In *Foundations of Genetic Algorithms (FOGA 2017)*. ACM, New York, NY, USA, 149–161.

[4] Asma Atamna, Anne Auger, and Nikolaus Hansen. 2020. On Invariance and Linear Convergence of Evolution Strategies with Augmented Lagrangian Constraint Handling. *Theoretical Computer Science* 832 (Sept. 2020), 68–97. doi:10.1016/j.tcs.2018.10.006

[5] Richard H Byrd, Mary E Hribar, and Jorge Nocedal. 1999. An interior point algorithm for large-scale nonlinear programming. *SIAM Journal on Optimization* 9, 4 (1999), 877–900.

[6] David A Cohn, Zoubin Ghahramani, and Michael I Jordan. 1996. Active learning with statistical models. *Journal of artificial intelligence research* 4 (1996), 129–145.

[7] Yann Collette, Nikolaus Hansen, Gilles Pujol, Daniel Salazar Aponte, and Rodolphe Le Riche. 2013. *Object-Oriented Programming of Optimizers – Examples in Scilab*. John Wiley & Sons, Ltd, Hoboken, NJ, USA, Chapter 14, 499–538. doi:10.1002/9781118600153.ch14 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118600153.ch14

[8] Paul Dufossé and Nikolaus Hansen. 2021. Augmented Lagrangian, penalty techniques and surrogate modeling for constrained optimization with CMA-ES. In *Genetic and Evolutionary Computation Conference (GECCO 2021)*. ACM, New York, NY, USA, 519–527.

[9] Tobias Glasmachers. 2020. Global convergence of the (1+ 1) evolution strategy to a critical point. *Evolutionary computation* 28, 1 (2020), 27–53.

[10] Nikolaus Hansen. 2016. The CMA Evolution Strategy: A Tutorial. CoRR abs/1604.00772 (2016), 39 pages. arXiv:1604.00772 http://arxiv.org/abs/1604.00772

[11] Nikolaus Hansen, Youhei Akimoto, and Petr Baudis. 2019. CMA-ES/pycma on Github. Zenodo. DOI:10.5281/zenodo.2559634. doi:10.5281/zenodo.2559634

[12] Nikolaus Hansen and Anne Auger. 2013. Principled design of continuous stochastic search: From theory to practice. In *Theory and principled methods for the design of metaheuristics*. Springer, Berlin, Heidelberg, Germany, 145–180.

[13] Nikolaus Hansen, Anne Auger, Raymond Ros, Steffen Finck, and Petr Pošík. 2010. Comparing Results of 31 Algorithms from the Black-Box Optimization Benchmarking BBOB-2009. In *Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation* (Portland Oregon USA, 2010-07-07). ACM, New York, NY, USA, 1689–1696. doi:10.1145/1830761.1830790

[14] Dieter Kraft. 1994. Algorithm 733: TOMP-Fortran modules for optimal control calculations. *ACM Transactions on Mathematical Software (TOMS)* 20, 3 (1994), 262–281.

[15] Michael JD Powell. 1994. *A direct search optimization method that models the objective and constraint functions by linear interpolation*. Springer, Berlin, Heidelberg, Germany.

[16] Luis Miguel Rios and Nikolaos V Sahinidis. 2013. Derivative-free optimization: a review of algorithms and comparison of software implementations. *Journal of Global Optimization* 56, 3 (2013), 1247–1293.

[17] Bernhard Scholkopf and Alexander J Smola. 2018. *Learning with kernels: support vector machines, regularization, optimization, and beyond*. MIT Press, Cambridge, MA, USA.

[18] Shai Shalev-Shwartz and Shai Ben-David. 2014. *Understanding Machine Learning: From Theory to Algorithms* (1 ed.). Cambridge University Press, Cambridge, MA, USA. doi:10.1017/CBO9781107298019

[19] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. 2020. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods* 17 (2020), 261–272. doi:10.1038/s41592-019-0686-2