# Direct-search methods in the year 2025: Theoretical guarantees and algorithmic paradigms

K.J. Dzahini, F. Rinaldi, C.W. Royer, D. Zeffiro

► To cite this version:

K.J. Dzahini, F. Rinaldi, C.W. Royer, D. Zeffiro. Direct-search methods in the year 2025: Theoretical guarantees and algorithmic paradigms. EURO Journal on Computational Optimization, 2025, 13, pp.100110. 10.1016/j.ejco.2025.100110. hal-05235664

HAL Id: hal-05235664

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

Submitted on 12 Feb 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
# Direct-search methods in the year 2025: Theoretical guarantees and algorithmic paradigms

K.J. Dzahini$^{a, b, 1}$, F. Rinaldi$^{b, c, *}$, C.W. Royer$^{c, d, 2}$, D. Zeffiro$^{b}$

## ABSTRACT

Optimizing a function without using derivatives is a challenging paradigm, that precludes from using classical algorithms from nonlinear optimization, and may thus seem intractable other than by using heuristics. Nevertheless, the field of derivative-free optimization has succeeded in producing algorithms that do not rely on derivatives and yet are endowed with convergence guarantees. One class of such methods, called direct-search methods, is particularly popular thanks to its simplicity of implementation, even though its theoretical underpinnings are not always easy to grasp.

In this work, we survey contemporary direct-search algorithms from a theoretical viewpoint, with the aim of highlighting the key theoretical features of these methods. We provide a basic introduction to the main classes of direct-search methods, including line-search techniques that have received little attention in earlier surveys. We also put a particular emphasis on probabilistic direct-search techniques and their application to noisy problems, a topic that has undergone significant algorithmic development in recent years. Finally, we complement existing surveys by reviewing the main theoretical advances for solving constrained and multiobjective optimization using direct-search algorithms.

## 1. Introduction

Direct-search methods, whose name was coined by Hooke and Jeeves [61], form one of the main classes of derivative-free optimization (DFO) algorithms. Those methods work in a very easy and intuitive way: they sample the objective function at a finite number of points at each iteration and decide which actions to take next solely based on those function values, either by function value comparison through ranking or based on numerical values, and without fundamentally relying upon any explicit or implicit derivative approximation or model building. This simple principle partly explains why direct-search approaches have been highly popular since the early days of numerical optimization. Despite the existence of tools that allow to easily use first- or second-order derivative information in an algorithm (i.e., tools that compute derivatives automatically), there is still a need for methods that rely solely on function values in simulation-based problems where derivative information can hardly be obtained [9]. As a result, direct-search methods remain a popular choice of practitioners. Although classical techniques such as Nelder-Mead's simplex method [81] or evolutionary algorithms [16] are not built on strong theoretical foundations, modern direct-search algorithms are typically equipped with convergence guarantees, that certify their ability to converge to a (possibly approximate) solution of the problem at hand asymp-

https://doi.org/10.1016/j.ejco.2025.100110

<footnote>$^{a}$ Argonne National Laboratory, 9700 S. Cass Avenue, Lemont, IL 60439, USA</footnote>

<footnote>$^{b}$ Dipartimento di Matematica “Tullio Levi-Civita”, Università di Padova, Padua, Italy</footnote>

<footnote>$^{c}$ LAMSADE, CNRS, Université Paris Dauphine-PSL, Place du Maréchal de Lattre de Tassigny, 75016 Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: kdzahini@anl.gov (K.J. Dzahini), rinaldi@math.unipd.it (F. Rinaldi), clement.royer@lamsade.dauphine.fr (C.W. Royer), zeffiro@math.unipd.it (D. Zeffiro).</footnote>

<footnote n="1">¹ Funding of this author was supported by the CAMPA and ComPASS-4, projects of the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research and Office of High Energy Physics, Scientific Discovery through Advanced Computing (SciDAC) program under Contract No. DE-AC02-06CH11357.</footnote>

<footnote n="2">² Funding for this author's research was partially provided by Agence Nationale de la Recherche through program ANR-19-P3IA-0001 (PRAIRIE 3IA Institute) and by CNRS under the grant IAE BONUS.</footnote>
totically (global convergence) or in finite time (complexity bounds). Still, those guarantees are somewhat scattered in the literature, in research papers as well as surveys and books on derivative-free optimization.

Indeed, direct-search methods have been the subject of several dedicated surveys [4,65,73,90]. The 1996 survey of Wright [90] was arguably the first on direct-search schemes, with an emphasis on simplicial methods that may lack convergence guarantees. The later survey of Lewis et al. [73] discussed in particular pattern search methods along with the convergence results that relied on simple decrease. The landmark survey of Kolda et al. [65] followed and provided a thorough investigation of convergence results on direct-search schemes based on simple and sufficient decrease, highlighting developments in the (linearly) constrained setting but focusing mainly on smooth objectives. To the best of our knowledge, the most recent survey article is due to Audet [4], who described the developments of direct-search methods with an emphasis on mesh adaptive direct-search (MADS) variants and their use in applications. More broadly, direct-search methods have been featured in textbooks and reviews on derivative-free optimization. Conn et al. [32] covered the basics of direct-search methods with convergence results. Audet and Hare [9] focused on direct-search schemes based on simple decrease, and provided a thorough introduction to pattern search and MADS. Finally, in the recent overview of derivative-free algorithms, Larson et al. [70] provided complexity results for direct-search methods in comparison with zeroth-order algorithms, especially in the convex setting.

This survey aims at complementing existing work by reviewing the main direct-search techniques proposed in the literature and their theoretical guarantees. First, more emphasis is put on line-search techniques, that have been somewhat overlooked in previous surveys. Secondly, connections are drawn between the directional and the mesh-based direct-search techniques, two classes of methods that tend to be discussed separately. We believe that readers should be aware of the two classes of algorithms, along with their specifics, that are highlighted throughout our work. Third, this survey begins with a thorough treatment of direct-search methods in the unconstrained setting, in both a smooth and a nonsmooth setting (note that we do not leverage convexity in this survey's results). The underlying goal is to ease the transition from this setting to more complicated setups, such as noisy or constrained optimization. We believe that a reader interested in using direct-search methods would likely face such challenges, and that a good understanding of an algorithm's behavior in the simplest possible setting provides valuable insight regarding its extension to other settings.

The rest of the survey is structured as follows. Section 2 provides a generic direct-search template for unconstrained problems, along with the key mathematical concepts used to establish convergence results. Section 3 reviews the main classes of algorithms in the unconstrained setting and represents a starting point for any of the subsequent sections. Section 4 describes algorithms that tolerate noise in the function values. Section 5 is concerned with handling constrained formulations. Section 6 presents the main results in a multiobjective setting. The final comments are given in Section 7.

## 2. Main algorithm and definitions

**Notations** Throughout, the vectors are written in lowercase boldface (e.g., $v \in \mathbb{R}^n, n \ge 2$) while matrices are written in uppercase boldface (e.g., $S \in \mathbb{R}^{n \times p}$). The set of column vectors of a matrix $D$ is denoted by $\mathcal{D}$. Other sets are denoted by calligraphic or upper Greek letters (e.g., $\mathcal{K}$, $\Omega$), except for the sets $\mathbb{N}$ and $\mathbb{R}$ of natural and real numbers. The set of nonnegative real numbers are denoted by $\mathbb{R}_{\ge 0}$. Given two integers $a, b$, we write $[a : b]$ for the set of all integer values $a \le z \le b$. Sequences indexed by $\mathbb{N}$ are denoted by $\{a_k\}_{k \in \mathbb{N}}$ or $\{a_k\}$ in absence of ambiguity. Random objects are denoted using underlined letters (e.g., $z, v, S$). The expected value operator is denoted by $\mathbb{E}[\cdot]$. When necessary, we specify the variables on which the expected value is computed (e.g. $\mathbb{E}_{\zeta}[\cdot]$ for the expected value over $\zeta$). Given a point $x \in \mathbb{R}^n$ and a set $\Omega \subset \mathbb{R}^n$, we use $\pi_{\Omega}(\cdot)$ for the projection operator on $\Omega$. The notation $B_{\delta}(x)$ refers to the ball of radius $\delta$ centered at $x$. Finally, given a quantity $A$, the notation $\mathcal{O}(A)$ means a constant times $A$, where the constant does not depend on $A$.

### 2.1. The direct-search framework

In this section, we introduce the main ideas behind direct-search schemes by considering an unconstrained minimization problem

$$ \min_{x \in \mathbb{R}^n} f(x), \qquad (1) $$

where $f: \mathbb{R}^n \to \mathbb{R}$ is the objective function of our problem. In this survey, we describe theoretical guarantees that require some regularity on $f$ (continuity, smoothness). However, we emphasize that a direct-search method only requires access to the values of $f$ to proceed.

Algorithm 1 describes a generic direct-search framework, that is used as a baseline for developing more involved variants in the subsequent sections. The method initially selects a point in the variable space (also called the *incumbent solution* in the direct-search literature), as well as an initial value for the so-called *stepsize*. At every iteration, the function is evaluated at trial points, obtained from the current point by a displacement along certain directions called *poll directions*. If one of those points provides a decrease with respect to the current function value as measured by a *forcing function*, this point can be accepted as the new point, and the stepsize can be increased (or kept constant). Otherwise, the next point is set to be the current point, and the stepsize is decreased.

Each step of Algorithm 1 highlights the key ingredients that define a given direct-search method. First, the choice of poll directions $D_k$ is instrumental to deriving convergence results. These directions can be chosen in a deterministic or randomized fashion, fixed throughout the iterations or varying at every iteration. In the rest of Section 2, we review the most classical polling strategies.

Secondly, the condition for success

$$ f(x_k + s) < f(x_k) - \rho(\alpha_k), \qquad (2) $$
**Algorithm 1:** Basic direct-search method.

1: **Inputs:** Starting point/incumbent solution $x_0 \in \mathbb{R}^n$, initial stepsize $\alpha_0 > 0$, stepsize update parameters $0 < \theta < 1 \le \gamma$, forcing function $\rho: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$.

2: **for** $k = 0, 1, 2, \dots$ **do**

3: Select a set $\mathcal{D}_k$ of poll directions.

4: if $f(x_k + s) < f(x_k) - \rho(\alpha_k)$ holds for some $s \in \{\alpha_k d : d \in \mathcal{D}_k\}$ **then**

5: Declare the iteration as successful, set $x_{k+1} = x_k + s$ and $\alpha_{k+1} = \gamma\alpha_k$.

6: **else**

7: Declare the iteration as unsuccessful, set $x_{k+1} = x_k$ and $\alpha_{k+1} = \theta\alpha_k$.

8: **end if**

9: **end for**

characterizes a major divide between direct-search methods, depending on whether or not any possible improvement in the objective function is accepted. Algorithm 1 (or an instance thereof) is said to be based on a sufficient decrease condition when $\rho(\alpha) > 0$ if $\alpha > 0$. Otherwise, the algorithm is said to be based on simple decrease.

Thirdly, in order to check the decrease condition (2), Algorithm 1 requires access to the function value at both $x_k + s$ and $x_k$. When exact values of $f$ are available, one only needs to evaluate $f$ once at a given point, hence values can be saved between iterations to avoid extra calculations. On the contrary, repeated evaluations (also termed replications in the literature) are particularly useful in the context of stochastic optimization, where only noisy evaluations of the objective are available. Such algorithms are an important part of this survey and are reviewed in Section 4.

Fourthly, the set of points that are considered at every iteration may be broader than those obtained through the polling directions and the current stepsize. In particular, direct-search methods of line-search type, which are discussed in Section 3.3, evaluate additional points beyond those considered in polling in order to perform extrapolation.

Finally, the decision made at every iteration of Algorithm 1 may not require all points in $\{\alpha_k d : d \in \mathcal{D}_k\}$ to be checked in case of success (although this remains necessary to declare an iteration as unsuccessful). Algorithm 1 (or an instance thereof) is said to use opportunistic polling if Step 2 stops as soon as it finds a vector $s$ satisfying (2). Otherwise, the algorithm is said to use complete polling.

**Remark 2.1.** Practical direct-search methods typically include an optional search step, that allows for using problem-specific heuristics as well as additional strategies for optimization [9, Section 7.2]. Nevertheless, theoretical guarantees of direct-search techniques do not rely on the properties of this search step. Since our survey focuses on these guarantees, we present direct-search methods without a search step, but emphasize that search steps are often crucial to the practical performance of these algorithms.

## 2.2. Main mathematical concepts

In this section, we gather key algebraic definitions associated with the choice of directions in a direct-search method, for future reference in the rest of the survey. We also provide standard analytical definitions associated with the objective function (and the constraint functions when applicable).

### 2.2.1. Geometry of polling sets

One of the main features of a direct-search method is its strategy to choose the poll directions. Those directions play a key role in guaranteeing convergence, provided they possess favorable geometric properties. The first property below is an intrinsic property of a direction set.

**Definition 2.1.** ([57, Section 2.1]) For a vector $v \in \mathbb{R}^n \setminus \{0\}$ and a set of directions $\mathcal{D}$, the cosine measure of $\mathcal{D}$ with respect to $v$ is defined as

$$ \mathrm{cm}(\mathcal{D}, v) := \max_{d \in \mathcal{D}} \frac{\langle d, v \rangle}{\|d\| \|v\|} $$

The cosine measure of $\mathcal{D}$ is then

$$ \mathrm{cm}(\mathcal{D}) := \min_{v \in \mathbb{R}^n \setminus \{0\}} \mathrm{cm}(\mathcal{D}, v). $$

The cosine measure of a set $\mathcal{D}$ captures how much of the space is covered by the cone spanned using the vectors in $\mathcal{D}$. The best scenario corresponds to spanning the entire space, and corresponds to the following notion.

**Definition 2.2.** A set $\mathcal{D}$ of $n_D$ vectors in $\mathbb{R}^n$ is called a positive spanning set if its cosine measure is positive, or equivalently, if every vector in $\mathbb{R}^n$ can be written as a nonnegative linear combination of vectors in $\mathcal{D}$.

Positive spanning sets are the typical choice of poll directions for direct-search methods. In addition to relying on positive spanning sets, certain direct-search methods ensure that the points considered for evaluation lie on carefully designed grids. We adopt here the terminology of meshes and frames [9], and provide the definitions independently of any algorithm.
**Definition 2.3.** Let $\mathbf{D} := \mathbf{GY} \in \mathbb{R}^{n \times n_D}$ be the matrix representation of a set $\mathbf{D}$, where $\mathbf{G} \in \mathbb{R}^{n \times n}$ is invertible and the columns of $\mathbf{Y} \in \mathbb{Z}^{n \times n_D}$ form a positive spanning set for $\mathbb{R}^n$. For any vector $\mathbf{x} \in \mathbb{R}^n$ and any $\alpha^m > 0$, the mesh of coarseness $\alpha^m$ corresponding to $\mathbf{D}$ and centered at $\mathbf{x}$, is defined by

$$ \mathcal{M} = \{ \mathbf{x} + \alpha^m \mathbf{D} \mathbf{z} : \mathbf{z} \in \mathbb{N}^{n_D} \} \subset \mathbb{R}^n. \quad (3) $$

The frame $\mathcal{F} \subset \mathcal{M}$ of extent $\alpha^p \ge \alpha^m$ generated by $\mathbf{D}$, centered at $\mathbf{x}$, is defined by

$$ \mathcal{F} = \{ \mathbf{y} \in \mathcal{M} : \|\mathbf{y} - \mathbf{x}\|_{\infty} \le \alpha^p d_{\max} \}, \quad (4) $$

where $d_{\max} := \max \{ \|\mathbf{d}'\|_{\infty} : \mathbf{d}' \in \mathbf{D} \}$.

One interesting feature of meshes and frames is that they lead to asymptotic density, a concept that has proven fundamental in analyzing direct-search methods, especially in the nonsmooth setting [7,9].

**Definition 2.4.** An infinite sequence $\{\mathbb{D}_k\}_{k \in \mathcal{K}}$ of vector sets defines an asymptotically dense set of directions if the set

$$ \mathbb{D}_{\infty} := \{ \mathbf{d} / \|\mathbf{d}\|_{\infty} : \mathbf{d} \in \mathbb{D}_k \}_{k \in \mathcal{K}} $$

is dense in the unit sphere $\mathbb{S}^{n-1}$ of $\mathbb{R}^n$, i.e. $\forall \mathbf{u} \in \mathbb{S}^{n-1}, \forall \epsilon > 0, \exists \mathbf{d}_{\infty} \in \mathbb{D}_{\infty}$ such that $\|\mathbf{u} - \mathbf{d}_{\infty}\|_2 < \epsilon$.

### 2.2.2. Function classes and stationary points

By design, a direct-search method can be applied to any function for which the objective value (or approximations thereof, see Section 4) can be queried. Still, in order to establish convergence and convergence rates, one typically assumes a certain level of regularity in the objective function. For this reason, we assume continuity of the objective function, though we note that direct-search methods have also been applied to extended-value and discontinuous problems [14,24,35]. In this survey, we consider both Lipschitz continuous functions or functions with Lipschitz continuous derivative, in the sense of the two definitions below.

**Definition 2.5.** A function $F : \mathbb{R}^n \to \mathbb{R}^{n'}$ with $n' \in \{1, n\}$ is said to be *L-Lipschitz continuous* for $L > 0$ if

$$ \forall (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^n \times \mathbb{R}^n, \quad \|F(\mathbf{x}) - F(\mathbf{y})\| \le L \|\mathbf{x} - \mathbf{y}\|. $$

**Definition 2.6.** A function $f : \mathbb{R}^n \to \mathbb{R}$ is said to be *L-smooth* for $L > 0$ if $f$ is continuously differentiable and $\nabla f : \mathbb{R}^n \to \mathbb{R}^n$ is *L-Lipschitz continuous*.

When optimizing a smooth function without constraints, direct-search methods generally aim at converging towards a point with zero gradient, called a *stationary point*. We present results of that form in the upcoming sections. Existing analyzes of direct-search methods also provide guarantees on nonsmooth functions using generalized derivatives. The next two definitions from Clarke differential calculus [31] correspond to the most classical notions of generalized derivatives used in the direct-search literature.

**Definition 2.7.** The Clarke directional derivative of $f$ at $\mathbf{x} \in \mathbb{R}^n$ in the direction $\mathbf{v} \in \mathbb{R}^n$ is

$$ f^{\circ}(\mathbf{x}; \mathbf{v}) = \limsup_{\substack{\mathbf{y} \to \mathbf{x} \\ t \to 0}} \frac{f(\mathbf{y} + t\mathbf{v}) - f(\mathbf{y})}{t}. \quad (5) $$

The Clarke subdifferential of $f$ at $\mathbf{x}$ is then defined as

$$ \partial f(\mathbf{x}) = \{ \mathbf{g} \in \mathbb{R}^n \mid \mathbf{g}^{\top} \mathbf{v} \le f^{\circ}(\mathbf{x}; \mathbf{v}) \quad \forall \mathbf{v} \in \mathbb{R}^n \}. \quad (6) $$

**Definition 2.8.** The point $\mathbf{x}^*$ is Clarke stationary for $f$ if $f^{\circ}(\mathbf{x}^*; \mathbf{v}) \ge 0$ for every $\mathbf{v} \in \mathbb{R}^n$ or, equivalently, if $\mathbf{0} \in \partial f(\mathbf{x}^*)$.

The notion of Clarke stationarity is weaker than many others used in variational analysis. As a result, it has been used extensively to derive guarantees for direct-search schemes [5]. More recently, an alternate notion of stationarity due to Goldstein [55] has been employed to obtain complexity results for direct-search methods.

**Definition 2.9.** Given $\delta > 0$, the point $\mathbf{x}^*$ is called $\delta$-Goldstein stationary for $f$ if

$$ \inf \left\{ \|\mathbf{g}\| \mid \mathbf{g} \in \operatorname{conv} \left( \cup_{y \in B_{\delta}(\mathbf{x}^*)} \partial f(y) \right) \right\} = 0 $$

where $\operatorname{conv}(S)$ denotes the set of convex combinations of vectors in $S$.
Note that this survey also considers approximate definitions of stationary points in order to derive complexity results. For instance, in the smooth setting, rather than a point with zero gradient, we might be interested in converging towards a point with small enough gradient norm. In the nonsmooth setting, the notion of an approximate stationary point is more intricate, and we provide the definition below.

**Definition 2.10.** Given $\delta > 0$ and $\epsilon > 0$, a vector $x \in \mathbb{R}^n$ is called a $(\delta, \epsilon)$-Goldstein stationary point for a function $f$ if

$$ \inf \left\{ \|g\| \mid g \in \operatorname{conv} \cup_{y \in B_{\delta}(x)} \partial f(y) \} \le \epsilon. \quad (7) $$

The concept of approximate stationary points are also discussed in Section 3.1.

## 2.3. Basic convergence ingredients

In this section, we provide the main assumptions and technical lemmas that are used while analyzing a direct-search method. As in the rest of Section 2, we focus on the unconstrained, deterministic setting, and provide results for both smooth and nonsmooth objectives.

### 2.3.1. Main assumptions

We begin by listing the key assumptions that are made on the objective function.

**Assumption 2.1.** The objective $f$ is lower bounded, i.e., there exists $f^* \in \mathbb{R}$ such that $f(x) \ge f^* \forall x \in \mathbb{R}^n$.

Assumption 2.1 is typically used together with a guarantee of decrease at every iteration of the algorithm to establish convergence.

**Assumption 2.2.** The sublevel set $L_0 := \{x \in \mathbb{R}^n \mid f(x) \le f(x_0)\}$ is bounded, where $x_0 \in \mathbb{R}^n$ is the initial point of Algorithm 1.

Assumption 2.2 certifies that the iterates of Algorithm 1 remain in a bounded set, which is particularly helpful to guarantee convergence of a subsequence of the iterates.

In addition to assumptions on the objective function, an analysis of Algorithm 1 typically requires conditions on the polling directions, that are stated below.

**Assumption 2.3.** For $d \in \mathbb{D}_k$, $k \in \mathbb{N}_0$, we have $\|d\| \in [d_{\min}, d_{\max}]$ where $d_{\min}$ and $d_{\max}$ are positive constants.

The upper bound ensures that the information obtained in unsuccessful steps is relevant for the directional derivative, while the lower bound ensures that the sufficient decrease condition can indeed be satisfied for good descent directions.

The last classical assumption concerns the forcing function.

**Assumption 2.4.** The forcing function $\rho: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$, is continuous, nondecreasing and $\rho(\alpha) = o(\alpha)$ as $\alpha \to 0$.

Note that using simple decrease (i.e., $\rho \equiv 0$) leads to immediate satisfaction of Assumption 2.4. Popular choices for $\rho$ in direct search based on sufficient decrease include $\rho(\alpha) = c \alpha^q$ with $c > 0$ and $q > 1$.

### 2.3.2. Key lemmas

A remarkable feature of direct-search schemes is that their analysis relies upon unsuccessful iterations, at which the iterate does not change. On those iterations, one can relate a stationarity criterion of interest to the stepsize parameter. In a smooth setting, the following lemma is the first step towards establishing theoretical guarantees of a direct-search technique.

**Lemma 2.1.** Let $f$ be an $L$-smooth function, and consider that Algorithm 1 is applied to the minimization of $f$. Suppose that the $k$th iteration of Algorithm 1 is unsuccessful. Then, for any $d \in \mathbb{D}_k$,

$$ -\langle \nabla f(\mathbf{x}_k), \mathbf{d} \rangle \le \alpha_k \frac{L \|\mathbf{d}\|^2}{2} + \frac{\rho(\alpha_k)}{\alpha_k}. \quad (8) $$

**Proof.** Since iteration $k$ is unsuccessful, for any $d \in \mathbb{D}_k$, the decrease condition is not satisfied, and thus

$$ f(\mathbf{x}_k + \alpha_k \mathbf{d}) - f(\mathbf{x}_k) \ge -\rho(\alpha_k). \quad (9) $$

Since $f$ is $L$-smooth by assumption, using for example [9, Lemma 9.4] yields

$$ f(\mathbf{x}_k + \alpha_k \mathbf{d}) \le f(\mathbf{x}_k) + \alpha_k \langle \nabla f(\mathbf{x}_k), \mathbf{d} \rangle + \frac{L \alpha_k^2}{2} \|\mathbf{d}\|^2. $$
Combining this inequality with (9) and dividing by $\alpha_k > 0$ yields the desired result. ☐

**Lemma 2.1** highlights the importance of choosing a good set of polling directions, so as to relate the left-hand side of (8) with a convergence criterion. The use of positive spanning sets is particularly suited in that regard, as it leads to the following result.

**Proposition 2.1.** Let the assumptions of Lemma 2.1 hold, in particular that the kth iteration of Algorithm 1 is unsuccessful. Suppose further that Assumption 2.3 holds, and that $\mathbb{D}_k$ is a positive spanning set with cosine measure $\tau > 0$. Then,

$$ \|\nabla f(\mathbf{x}_k)\| \le \frac{1}{\tau} \left[ \alpha_k \frac{L d_{\max}}{2} + \frac{\rho(\alpha_k)}{\alpha_k d_{\min}} \right]. \qquad (10) $$

The result of Proposition 2.1 can be used to guarantee convergence of the gradient norm to zero for a subsequence of unsuccessful iterates. We illustrate such a result in Section 3.1.

In absence of smoothness, the result of Lemma 2.1 no longer holds. However, one can still guarantee that

$$ \frac{f(\mathbf{x}_k + \alpha_k \mathbf{d}) - f(\mathbf{x}_k)}{\alpha_k} \ge - \frac{\rho(\alpha_k)}{\alpha_k}, \qquad (11) $$

for any $\mathbf{d} \in \mathbb{D}_k$. The left-hand side of (11) resembles the formula of the Clarke derivative, and is used together with the notion of refining subsequence to yield asymptotic guarantees on the Clarke directional derivatives. We define the former concept below.

**Definition 2.11.** Consider a sequence of iterates $\{\mathbf{x}_k\}_{k \in \mathbb{N}}$, a sequence of stepsize parameters $\{\alpha_k\}_{k \in \mathbb{N}}$ and a sequence of polling sets $\{\mathbb{D}_k\}_{k \in \mathbb{N}}$ associated with a direct-search algorithm. Let $\mathcal{K} \subset \mathbb{N}$ be an infinite index subsequence such that the kth iteration is unsuccessful for every $k \in \mathcal{K}$ and the subsequence $\{\mathbf{x}_k\}_{k \in \mathcal{K}}$ converges towards $\hat{\mathbf{x}} \in \mathbb{R}^n$.

* The subsequence $\{\mathbf{x}_k\}_{k \in \mathcal{K}}$ is called a refining subsequence if $\lim_{k \in \mathcal{K}} \alpha_k = 0$, in which case $\hat{\mathbf{x}}$ is called a refined point.

* If $\{\mathbf{x}_k\}_{k \in \mathcal{K}}$ is a refining subsequence, a direction $\mathbf{d} \in \mathbb{R}^n$ is called a refining direction if there exists an infinite subset $\mathcal{L} \subset \mathcal{K}$ and a sequence $\{\mathbf{d}_\ell\}_{\ell \in \mathcal{L}}$ such that $\mathbf{d}_\ell \in \mathbb{D}_\ell$ for any $\ell \in \mathcal{L}$ and $\lim_{\ell \in \mathcal{L}} \frac{\mathbf{d}_\ell}{\|\mathbf{d}_\ell\|} = \frac{\mathbf{d}}{\|\mathbf{d}\|}$.

By reasoning on refining subsequences, one can directly establish asymptotic results on (Clarke) directional derivatives. This property is crucial in analyzing algorithms based on simple decrease such as those described in Section 3.2.

**Lemma 2.2.** Suppose that $f$ is Lipschitz continuous, and consider that Algorithm 1 is applied to the minimization of $f$ with a forcing function that satisfies Assumption 2.4. Let $\mathcal{K} \subset \mathbb{N}$ be an infinite index subsequence such that $\{\mathbf{x}_k\}_{k \in \mathcal{K}}$ is a refining subsequence with refined point $\hat{\mathbf{x}}$. Then, for any refining direction $\mathbf{d}$, it holds that

$$ f^{\circ}(\hat{\mathbf{x}}; \mathbf{d}) \ge 0. \qquad (12) $$

**Proof.** Let $\{\mathbf{d}_\ell\}_{\ell \in \mathcal{L}}$ be a sequence of directions such that $\mathbf{d}_\ell \in \mathbb{D}_\ell$ for every $\ell \in \mathcal{L}$ and $\lim_{\ell \in \mathcal{L}} \frac{\mathbf{d}_\ell}{\|\mathbf{d}_\ell\|} = \frac{\mathbf{d}}{\|\mathbf{d}\|}$. Without loss of generality, suppose that all iterations in $\mathcal{L}$ denote unsuccessful iterations. Applying (11) at iteration $\ell \in \mathcal{L}$ to the direction $\mathbf{d}_\ell$ gives

$$ \frac{f(\mathbf{x}_\ell + \alpha_\ell \mathbf{d}_\ell) - f(\mathbf{x}_\ell)}{\alpha_\ell} \ge - \frac{\rho(\alpha_\ell)}{\alpha_\ell}. $$

Taking the limit on both sides as $\ell$ goes to infinity, we have $\mathbf{d}_\ell \to \mathbf{d}^*$, $\alpha_\ell \to 0$ as well as $\mathbf{x}_\ell \to \hat{\mathbf{x}}$. It follows that

$$ f^{\circ}(\hat{\mathbf{x}}; \mathbf{d}) \ge \lim_{\ell \to \infty} \frac{f(\mathbf{x}_\ell + \alpha_\ell \mathbf{d}_\ell) - f(\mathbf{x}_\ell)}{\alpha_\ell} \ge 0. \qquad \square $$

This lemma is used in the subsequent analyzes of direct-search algorithms in the nonsmooth setting.

## 3. Algorithms with guarantees for unconstrained optimization

In this section, we review the main instances of Algorithm 1 that possess theoretical guarantees. We begin by reviewing directional algorithms, arguably the simplest class of such instances. We then present a simple mesh-based algorithm, and present its convergence guarantees in the unconstrained setting, that follow from standard analyzes of MADS-type algorithms. Finally, we review the theory behind line-search algorithms.

### 3.1. Directional algorithms

Directional direct-search methods directly follow the template of Algorithm 1. Those schemes differ in their choice of polling directions, polling strategy, and forcing function. Under standard assumptions, a directional method generally comes with global
convergence guarantees [65]. Variants with sufficient decrease are also amenable to a complexity analysis [89], and we stress out that sufficient decrease is currently the only paradigm in which complexity guarantees have been obtained.

We first review the main results for deterministic directional direct-search methods on smooth problems in Section 3.1.1. We then describe the corresponding results for probabilistic directional direct-search schemes in Section 3.1.2. Finally, we cover existing guarantees for deterministic methods on nonsmooth problems in Section 3.1.3.

### 3.1.1. Deterministic convergence results in the smooth case

We begin the analysis of Algorithm 1 by focusing on the smooth setting. In that context, the polling sets are typically required to be positive spanning sets with positive cosine measure. We formalize this requirement through the following assumption.

**Assumption 3.1.** The polling sets {$D_k$} used in Algorithm 1 satisfy $cm(D_k) \ge \tau$ for every $k \in \mathbb{N}$, where $\tau > 0$.

A direct application of Proposition 2.1 under this assumption yields the following convergence result.

**Theorem 3.1.** Suppose that Algorithm 1 is applied to an L-smooth function $f$ under Assumptions 2.1 and 2.4. Suppose also that the sequence {$D_k$} satisfies Assumptions 2.3 and 3.1. Let $\mathcal{K}$ be an index set of unsuccessful iterations such that $\lim_{k \in \mathcal{K}} \alpha_k \to 0$. Then,

$$ \lim_{k \in \mathcal{K}} \| \nabla f(\mathbf{x}_k) \| = 0. \qquad (13) $$

**Proof.** All the assumptions made in Proposition 2.1 hold whence taking the limit when $k \to \infty$, $k \in \mathcal{K}$ in (10) yields (13) since $\rho(\alpha_k)/\alpha_k \to 0$ per Assumption 2.4. $\square$

Theorem 3.1 is a weak convergence result, in that it only guarantees that a subsequence of gradient norms is going to zero provided there exists a subsequence of stepsizes that converges to 0. Note that neither this result nor that of Lemma 2.1 do leverage the use of sufficient decrease. A first step towards improving the results of Theorem 3.1 consists of exploiting sufficient decrease to guarantee that the entire sequence of stepsizes converges to 0.

**Lemma 3.1.** Suppose that Algorithm 1 is applied to an L-smooth function $f$ under Assumptions 2.1, 2.3 and 2.4. Suppose further that the algorithm uses sufficient decrease. Then, $\lim_{k \to \infty} \alpha_k = 0$.

Combining Lemma 3.1 with Theorem 3.1 does not guarantee that $\lim_{k \to \infty} \| \nabla f(\mathbf{x}_k) \| = 0$, since the result of Theorem 3.1 only concerns subsequences of unsuccessful iterations. A popular approach to obtain these results requires complete polling to be performed [32,65], which corresponds to our assumptions in Theorem 3.2 below. Note that similar results can also be obtained in the context of locally strongly convex objectives [43,65], or in the context of line-search methods [79], to be discussed in Section 3.3.

**Theorem 3.2.** Let the assumptions of Lemma 3.1 hold. In addition, suppose that Algorithm 1 is applied using complete polling and a sequence of polling sets satisfying Assumption 3.1. Then,

$$ \lim_{k \to \infty} \| \nabla f(\mathbf{x}_k) \| = 0. $$

*Complexity guarantees* In addition to global convergence results, the sufficient decrease condition also endows direct-search schemes with complexity bounds, or non-asymptotic results, that quantify how many iterations or function evaluations are needed in the worst case to reach an approximate optimality criterion, such as $\| \nabla f(\mathbf{x}_k) \| \le \epsilon$ (for nonconvex objectives) or $f(\mathbf{x}_k) - \min_{x \in \mathbb{R}^n} f(\mathbf{x}) \le \epsilon$ (for convex and strongly convex objectives). The first bounds of this nature for direct-search schemes were established by Vicente [89] in the nonconvex setting, and later extended to the convex and strongly convex setting by Dodangeh and Vicente [41]. A unified and simplified argument was presented in [68], assuming however no expansion in successful steps and a suitable initialization of the algorithm.

One thing these results have in common is that they exhibit additional dependencies on the problem dimension compared to those for derivative-based algorithms, a common trait among derivative-free algorithms [28]. In this survey, we focus on analyzing the nonconvex setting, and give the corresponding result below.

**Theorem 3.3.** Let the assumptions of Theorem 3.2 hold. Suppose further that $\rho(\alpha) = c \alpha^2$ for some $c > 0$, and that $|D_k| \le m$ for every $k \in \mathbb{N}$. Then, for any $\epsilon > 0$, Algorithm 1 reaches an iterate $\mathbf{x}_k$ such that $\| \nabla f(\mathbf{x}_k) \| \le \epsilon$ in at most

$$ \mathcal{O}(\tau^{-2}\epsilon^{-2}) \text{ iterations and } \mathcal{O}(m\tau^{-2}\epsilon^{-2}) \text{ function evaluations,} \qquad (14) $$

where $\tau$ is the bound on the cosine measures given by Assumption 3.1.

As stated, the result of Theorem 3.3 does not exhibit explicit dependencies on the problem dimension. However, if one aims at minimizing the complexity bounds as functions of $m$ and $\tau$, it has been shown that the best choice consists of selecting $D_k$ as
an orthonormal basis and its negative, in which case the results hold with $m = 2n$ and $\tau = \frac{1}{\sqrt{n}}$ (see, e.g., [42,89]). As a result, the iteration complexity bounds exhibit a dependency in $n$, while the evaluation complexity bounds have a dependency in $n^2$. In the next section, we describe variants of directional direct-search algorithms that attempt to reduce that dependency through randomization.

**Remark 3.1.** Under convexity assumptions on $f$, it is possible to improve the dependency on $\epsilon$ from $\epsilon^{-2}$ to $\epsilon^{-1}$ in the convex case and $\log(\epsilon^{-1})$ in the strongly convex case, thereby following similar improvements in the complexity of gradient descent techniques [41]. Note however that the dependencies on $m$ and $\tau$, which are intrinsic to the use of a direct-search scheme, are identical to that in the nonconvex setting.

### 3.1.2. Probabilistic convergence results in the smooth case

An important subclass of directional direct-search algorithms, that has grown in importance over the past decade, relies on assuming probabilistic properties of the polling sets rather than relying on deterministic properties such as Assumption 3.1. The associated schemes employ a small number of randomly chosen directions to effectively reduce the number of points considered at every iteration compared to a deterministic approach from $\mathcal{O}(1)$ to $\mathcal{O}(n)$. We note that such randomized strategies have been applied to zeroth-order algorithms that attempt to build a gradient approximation, yielding improved theoretical guarantees compared to deterministic variants [19,82].

In the context of directional direct-search methods, Gratton et al. [57] established theoretical guarantees for randomized approaches thanks to the notion of probabilistic descent, that considers the polling sets as random objects.

**Assumption 3.2.** The random polling set sequence $\{\mathbb{D}_k\}$ satisfies

$$ \mathbb{P}(\text{cm}(\mathbb{D}_k, -\nabla f(\mathbf{x}_k)) \ge \tau \mid \mathcal{F}_{k-1}) \ge \beta \quad \text{for all } k, \qquad (15) $$

where $\tau > 0$, $\beta \in (0, 1)$, $\mathcal{F}_{k-1}$ is the $\sigma$-algebra generated by $\mathbb{D}_0, \dots, \mathbb{D}_{k-1}$ and $\mathcal{F}_{-1} = \emptyset$.

Note that Assumption 3.2 is clearly satisfied when Assumption 3.1 holds, i.e., when the polling sets are positive spanning sets. We note that since the objective function is deterministic and $\mathbb{D}_0, \dots, \mathbb{D}_{k-1}$ are $\mathcal{F}_{k-1}$-measurable by construction, then $\nabla f(\mathbf{x}_k)$ is also $\mathcal{F}_{k-1}$-measurable and hence deterministic with respect to $\mathcal{F}_{k-1}$ since its distribution before the construction of $\mathbb{D}_k$ is only determined by $\mathbb{D}_0, \dots, \mathbb{D}_{k-1}$. When the $\mathbb{D}_i, i = 0, 1, \dots, k$, are mutually independent, then in particular $\mathbb{D}_k$ is independent of $\mathcal{F}_{k-1}$ so that for all $k$, requiring $\mathbb{P}(\text{cm}(\mathbb{D}_k, \nu) \ge \tau) \ge \beta$ to hold for any arbitrarily fixed $\nu \in \mathbb{R}^n$ is enough to satisfy (15). When the elements in these independent sets $\mathbb{D}_i$ (which do not need to be positive spanning sets) are independent uniformly distributed random vectors on the unit sphere of $\mathbb{R}^n$, the latter inequality is shown in [57, Lemma B.2] to be satisfied, thereby providing a demonstration of the satisfiability of (15) and its realisticness. We finally emphasize that a more formal and rigorous demonstration of how the above inequality leads to (15) is provided by [57, Proposition B.3].

Under Assumption 3.2, one can derive a high-probability function evaluation complexity for Algorithm 1, that implies almost-sure global convergence of the algorithmic process. Note that expected value results can also be obtained by that process [57].

**Theorem 3.4.** Suppose that Algorithm 1 is applied to an L-smooth function $f$ with $\rho(\alpha) = c \alpha^2$ for some $c > 0$ and $\gamma = \theta^{-1}$. Suppose that Assumption 2.3 holds for all realizations of $\{\mathbb{D}_k\}$, and that Assumption 3.2 holds with $\beta > 1/2$ and $|\mathbb{D}_k| \le m$ for every $k \in \mathbb{N}$. Then, for any $\epsilon > 0$, Algorithm 1 reaches an iterate $\mathbf{x}_k$ such that $\|\nabla f(\mathbf{x}_k)\| \le \epsilon$ in at most

$$ \mathcal{O}(\tau^{-2}\epsilon^{-2}) \text{ iterations and } \mathcal{O}(m\tau^{-2}\epsilon^{-2}) \text{ function evaluations,} \qquad (16) $$

with probability at least $1 - e^{\mathcal{O}(\tau^{-2}\epsilon^{-2})}$.

When every $\mathbb{D}_k$ consists of a direction uniformly selected in the unit sphere and its negative, probabilistic descent (Assumption 3.2) is satisfied with $m = 2$ and $\tau = \mathcal{O}(1/\sqrt{n})$. It follows that the algorithm reaches an approximate stationary point in at most $\mathcal{O}(n\epsilon^{-2})$ iterations and function evaluations with probability $1 - e^{\mathcal{O}(n\epsilon^{-2})}$. The complexity bounds, that are probabilistic, improve over the deterministic result of Theorem 3.3 by a factor of $n$.

Similar results were also established for random search schemes in Bergou et al. [19], in which a stepsize dependent on the Lipschitz constant of the gradient was used. Direct-search methods employing a random direction and its negative at every iteration can be viewed as adaptive versions of those schemes, that maintain an estimate of the Lipschitz constant of the gradient through the stepsize parameter. The analysis of direct-search methods using probabilistic descent directions was recently extended by Roberts and Royer [87] to allow for unbounded sets of descent directions, highlighting the connection with the theory of random matrices. In that work, each descent direction is defined through combinations of columns of a random matrix $\mathbf{P}_k \in \mathbb{R}^{n \times r}$, that represents a subspace of possible directions at iteration $k$. The complexity result from Gratton et al. [57] is recovered, with an additional dependence on the probabilistic bounds for the extreme singular values of $\mathbf{P}_k$. Such probabilistic bounds are non trivial to compute and have been the subject of recent studies [48].
### 3.1.3. Deterministic convergence results in the nonsmooth case

We now turn to the nonsmooth setting, where we merely assume that $f$ is Lipschitz continuous. Convergence results for deterministic methods based on sufficient decrease require asymptotic density of the set of directions.

**Theorem 3.5.** Suppose that Algorithm 1 is applied to a Lipschitz continuous function $f$. Let $\mathcal{K}$ be a sequence of unsuccessful steps such that $\{x_k\}_{k \in \mathcal{K}}$ is a refining subsequence converging to $\hat{x}$. If the corresponding set of refining directions is dense in the unit sphere, then $\hat{x}$ is a Clarke stationary point.

**Proof.** By Lemma 2.2, for any refining direction $d$ associated with $\{x_k\}_{k \in \mathcal{K}}$, we have $f^\circ(\hat{x}; d) \ge 0$. Since the set of refining directions is dense in the unit sphere, for any nonzero vector $u \in \mathbb{R}^n$, there exists a sequence of refining directions $\{d_i\}_{i \in \mathbb{N}}$ such that $d_i \to \frac{u}{\|u\|}$ as $i \to \infty$. As a result, we obtain in the limit that $f^\circ(\hat{x}; \frac{u}{\|u\|}) \ge 0$. Since the Clarke directional derivative is positively homogeneous [31, Proposition 2.1.1], we also have $f^\circ(\hat{x}; u) \ge 0$, from which the desired conclusion follows. $\square$

Under additional assumptions, one can obtain a stronger result akin to Theorem 3.2 for nonsmooth objectives. While this result follows from well established techniques, we are not aware of any existing reference in the literature.

**Theorem 3.6.** Let the assumptions of Theorem 3.5 hold. Suppose further that Assumptions 2.1, 2.4 and 2.3 hold. For any $k \in \mathbb{N}$, define $\tilde{D}_k = \{d/\|d\| \mid d \in \tilde{D}_k\}$. If $\cup_{k \in \mathcal{K}} \tilde{D}_k$ is dense in the unit sphere for every infinite subsequence $\mathcal{K}$, then every limit point of $\{x_k\}_{k \in \mathbb{N}}$ is Clarke stationary.

Unlike in the smooth setting, complexity bounds for direct-search algorithms remain rather unexplored in the nonsmooth setting. Garmanjani and Vicente [53] provided complexity results for a special class of nonsmooth objectives that can be approximated by a family of smooth functions $\{\tilde{f}(\cdot, \mu)\}_{\mu>0}$ such that $\tilde{f}(\cdot, \mu) \to f$ as $\mu \to 0$. Examples of such functions include Gaussian or uniform smoothing, for which one can show that $\tilde{f}(\cdot, \mu)$ is $O(\frac{\sqrt{n}}{\mu})$-smooth [82]. We present below a simplified version of the analysis in Garmanjani and Vicente [53] based on a fixed smoothing parameter.

**Theorem 3.7.** Let $f$ be a Lipschitz continuous function satisfying Assumption 2.1. Given $\epsilon > 0$, suppose that we apply Algorithm 1 to the smoothed function $\tilde{f}(\cdot, \epsilon)$ under Assumption 2.2 with $\rho(\alpha) = \epsilon \alpha^{3/2}$ and polling sets satisfying Assumptions 2.3, 3.1 with $\tau = \frac{1}{\sqrt{n}}$.

Then, the algorithm takes at most $O(n^{3/2}\epsilon^{-3})$ iterations and $O(n^{5/2}\epsilon^{-3})$ function evaluations of the smoothed objective to obtain a point such that $\|\nabla_x \tilde{f}(x, \epsilon)\| < \epsilon$.

Very recently, Lin et al. [74] obtained a function evaluation complexity of $O(n^{3/2}\epsilon^{-5})$ for zeroth-order, smoothing-based methods to achieve a similar outcome to that of Theorem 3.7. The former result has a lower dependence on $n$ but a higher dependence from $\epsilon$ compared to that of the latter. Note however that the algorithm of Lin et al. [74] only uses true function values rather than values for smoothed approximations.

### 3.2. Mesh-based algorithms

In contrast with several directional direct-search techniques using sufficient decrease, mesh-based algorithms rely on simple decrease, and enforce convergence thanks to a careful generation of search directions. Pioneering methods of that form include generalized pattern search (GPS) [23,88] and frame-based methods [33]. We focus our presentation on the mesh adaptive direct-search (MADS) paradigm [7], that is nowadays dominant in the landscape of mesh-based algorithms, thanks in part to the practical success of the NOMAD software [15].

Algorithm 2, inspired by [9, Algorithm 8.1], describes a basic mesh-based approach for unconstrained optimization. A key feature of these methods is the use of two stepsize parameters called mesh size parameter denoted by $\alpha_k^m$ and poll size parameter denoted by $\alpha_k^p$, in order to produce a rich set of trial points to be considered by the algorithm. The algorithm defines directions based on a single set $\mathbb{D}$ formed by the columns of $\mathbf{D} := \mathbf{G}\mathbf{Y} \in \mathbb{R}^{n \times n_D}$, where $\mathbf{G} \in \mathbb{R}^{n \times n}$ is invertible and the columns of $\mathbf{Y} \in \mathbb{Z}^{n \times n_D}$ form a positive spanning set. At iteration $k$, the mesh $\mathcal{M}_k$, and the frame $\mathcal{F}_k$ centered at $x_k \in \mathbb{R}^n$ are defined respectively as

$$ \mathcal{M}_k := \bigcup_{x \in \mathcal{V}_k} \{x + \alpha_k^m \mathbf{D} z : z \in \mathbb{N}^{n_D}\} \subset \mathbb{R}^n \quad \text{and} \quad \mathcal{F}_k := \{x \in \mathcal{M}_k : \|x - x_k\|_{\infty} \le \alpha_k^p d_{\max}\}, \quad (17) $$

where $\mathcal{V}_k$ is the set containing all the points where the objective function has been evaluated by the start of iteration $k$, $\alpha_k^m \le \alpha_k^p$ for all $k$, and $d_{\max} = \max\{\|d'\|_{\infty} : d' \in \mathbb{D}\}$, so that $\mathcal{F}_k \subseteq \mathcal{M}_k$. We note that the relation (17) justifies that the poll size parameter is sometimes referred to as frame size parameter in the literature. Note that the mesh is conceptual in the sense that it is actually never explicitly constructed nor stored. Trial points are generated so as to belong to the mesh, by defining $\mathcal{M}_k$ according to $\mathcal{V}_k$ and using the directions in $\mathbb{D}$.

The goal of a mesh-based iteration is to locate an improved mesh point by evaluating the objective at a finite number of points lying on the frame $\mathcal{F}_k$. Algorithm 3 describes how such directions can be chosen in order to lie on the frame. In practice, we point
**Algorithm 2: Mesh-based direct-search algorithm.**

1: **Inputs:** $\mathbf{x}_0 \in \mathbb{R}^n$, $\alpha_k^p > 0$ and $\theta \in (0, 1) \cap \mathbb{Q}$.

2: **for** $k = 0, 1, 2, \dots, \text{do}$

3: Select a positive spanning set $\mathbb{D}_k^p$ (e.g., by calling Algorithm 3) so as to satisfy $\mathbf{x}_k + \alpha_k^m \mathbf{d} \in \mathcal{F}_k \subseteq \mathcal{M}_k$ for all $\mathbf{d} \in \mathbb{D}_k^p$, where $\alpha_k^m = \min\{\alpha_k^p, (\alpha_k^p)^2\}$.

4: if $f(\mathbf{x}_k + \mathbf{s}) < f(\mathbf{x}_k)$ holds for some $\mathbf{s} \in \{\alpha_k^m \mathbf{d} : \mathbf{d} \in \mathbb{D}_k^p\}$ **then**

5: Declare the iteration as successful, set $\mathbf{x}_{k+1} = \mathbf{x}_k + \mathbf{s}$ and $\alpha_{k+1}^p = \theta^{-1}\alpha_k^p$.

6: **else**

7: Declare the iteration as unsuccessful, set $\mathbf{x}_{k+1} = \mathbf{x}_k$ and $\alpha_{k+1}^p = \theta\alpha_k^p$.

8: **end if**

9: **end for**

out that additional points on $\mathcal{M}_k \setminus \mathcal{F}_k$ are often considered during a search step. When iteration $k$ of Algorithm 2 is unsuccessful, the frame $\mathcal{F}_k$ is called a *minimal frame* and $\mathbf{x}_k$ is called a *minimal frame center* [7, Section 2]. In that case, we have $\alpha_{k+1}^p = \theta\alpha_k^p < \alpha_k^p$ and $\alpha_{k+1}^m < \alpha_k^m$, but allowing the poll size parameter $\alpha_k^m$ to decrease more rapidly than $\alpha_k^p$ leads to strong guarantees on the trial points that are considered by the method [9].

**Algorithm 3: Creating the set $\mathbb{D}_k^p$ of poll directions.**

1: **Inputs:** $\mathbf{v}_k \in \mathbb{R}^n$ with $\|\mathbf{v}_k\| = 1$ and $\alpha_k^p \ge \alpha_k^m > 0$.

2: Use $\mathbf{v}_k$ to create the Householder matrix $\mathbf{H}_k = \mathbf{I}_n - 2\mathbf{v}_k \mathbf{v}_k^\top = [\mathbf{h}_{1,k} \ \mathbf{h}_{2,k} \ \dots \ \mathbf{h}_{n,k}] \in \mathbb{R}^{n \times n}$.

3: Define $\mathbb{B}_k = \{\mathbf{h}_{1,k} \ \mathbf{h}_{2,k} \ \dots \ \mathbf{h}_{n,k}\}$ with $\mathbf{h}_{j,k} = \text{round}\left(\frac{\alpha_k^p}{\alpha_k^p + \|\mathbf{h}_{j,k}\|_{\infty}}\right) \in \mathbb{Z}^n, j = 1, 2, \dots, n$.

4: Set the poll set as $\mathbb{D}_k^p = \mathbb{B}_k \cup (-\mathbb{B}_k)$.

A first step towards establishing convergence of Algorithm 2 consists of observing that the coarseness level of the mesh vanishes asymptotically.

**Theorem 3.8.** ([9, Theorem 8.1]) Suppose that Algorithm 2 is applied to a function $f$ satisfying Assumption 2.2. Then, the stepsize parameter sequences satisfy

$$ \liminf_{k \to \infty} \alpha_k^m = \liminf_{k \to \infty} \alpha_k^p = 0. \qquad (18) $$

The proof of Theorem 3.8 relies on the properties of the meshes rather than sufficient decrease, as well as the use of a rational parameter to update both $\alpha_k^m$ and $\alpha_k^p$ [9, Lemmas 7.1-7.3]. In particular, the structure of the matrix $\mathbf{D}$ used in Algorithm 2 guarantees that the distance between distinct mesh points is lower bounded by a function of the mesh size parameter $\alpha_k^m$.

The second argument in a mesh-based convergence theory relates to the Clarke generalized directional derivative of the objective function $f$ at refined points in refining directions, in the sense of Definition 2.11. Crucially, the notion of the refining sequence uses $\{\alpha_k^m\}$ as the stepsize parameters but $\{\mathbb{D}_k^p\}$ as polling sets. By combining Theorem 3.8 with Assumption 2.2, one can then establish that there exists at least one refining subsequence with a refined point [9, Theorem 7.6]. Meanwhile, the existence of a refining direction follows from the compactness of the closed unit ball along with Bolzano-Weierstrass's theorem [9, Section 8.3]. Overall, the following convergence result can be derived.

**Theorem 3.9.** ([9, Theorem 8.3]) Suppose that Algorithm 2 is applied to a Lipschitz continuous function $f$ satisfying Assumption 2.2. Let $\{\mathbf{x}_k\}_{k \in \mathcal{K}}$ be a refining subsequence for the algorithm with refined point $\hat{\mathbf{x}} \in \mathbb{R}^n$, and let $\mathbf{d} \in \mathbb{R}^n$ be a refining direction for $\hat{\mathbf{x}}$. Then $f(\hat{\mathbf{x}}; \mathbf{d}) \ge 0$.

Theorem 3.9 is a direct consequence of Lemma 2.2, but only appeared in the more general, constrained setting throughout the mesh-based direct-search literature [7, Theorem 3.12], [9, Theorem 8.3]). As such, the result of Theorem 3.9 is relatively weak, in the sense that we have no information regarding the refining directions. A key to strengthening Theorem 3.9 lies in the fact that the mesh size parameter $\alpha_k^m$ goes to zero faster than the frame size parameter $\alpha_k^p$. By carefully designing the poll directions, it then becomes possible to guarantee asymptotic density of those directions.

**Theorem 3.10.** ([9, Theorem 8.6]) The sequence of sets $\mathbb{D}_k^p$ produced by Algorithm 3 forms an asymptotically dense set of directions if the directions $\mathbf{v}_k$ used to construct the Householder matrices $\mathbf{H}_k$ form a dense sequence $\{\mathbf{v}_k\}_{k \in \mathbb{N}}$ in the unit sphere $\mathbb{S}^{n-1}$.

**Corollary 3.1.** Let the assumptions of Theorem 3.9 hold, and suppose that the set of refining directions for $\hat{\mathbf{x}}$ is dense in the unit sphere. Then $\hat{\mathbf{x}}$ is a Clarke stationary point.

Remarkably, a simpler version of the result leads to convergence to a stationary point in the presence of smoothness, which we again emphasize is not the usual setup of mesh-based analyzes.
**Corollary 3.2.** Under the assumptions of Theorem 3.9, suppose further that $f$ is smooth near a refined point $\hat{x} \in \mathbb{R}^n$ and that the set of refining directions contains a positive spanning set. Then, $\hat{x}$ is a stationary point for $f$.

By design, a mesh-based algorithm aims at generating an asymptotically dense set of polling directions, and thus fulfills convergence conditions in both a nonsmooth and a smooth setting. In that sense, mesh-based algorithms enjoy asymptotic convergence guarantees for a broad class of functions. Non-asymptotic results, such as complexity results or convergence rates, have not been derived for mesh-based methods that rely on simple decrease. Note that combining a mesh-based approach with sufficient decrease does lead to complexity guarantees, however this combination has only been adopted in the context of stochastic optimization (see Section 4.2), and remains nonstandard in the deterministic case.

### 3.3. Line-search algorithms

Line-search-based derivative-free methods combine a sufficient decrease condition with a thorough search-direction exploration procedure. Those line-search approaches date back 40 years, with a first example of derivative-free line-search method proposed by De Leone, Gaudioso and Grippo [60], and later analyzed by Grippo, Lampariello and Lucidi [60]. Those strategies were then suitably adapted by Lucidi and Sciandrone [79] to define new globally convergent direct-search schemes that make use of an extrapolation step.

Algorithm 5 provides an example of exploration procedure that performs an extrapolation step once that a given trial point satisfies the acceptance condition has been found. Such a step ensures improved practical performances by accepting larger step sizes than those obtained using the basic strategy of Algorithm 1. As it can be easily seen by observing the scheme, the extrapolation indeed enables us to further explore the direction $d$ by suitably increasing the stepsize (see Steps 5-7) as long as sufficient decrease is guaranteed.

<caption>**Algorithm 4:** Line-search-based DS scheme.</caption>

1: **Inputs:** $x_0 \in \mathbb{R}^n$, $m \in \mathbb{N}$, $\{\alpha'_k\}_{j \in [1:m]}$, $\theta \in (0, 1), \gamma > 1$, a forcing function $\rho: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$.
2: **for** $k=0, 1, \dots, \text{do}$
3: &nbsp;&nbsp;&nbsp;&nbsp;Select a set of poll directions $D_k = \{d'_k\}_{j \in [1:m]}$.
4: &nbsp;&nbsp;&nbsp;&nbsp;Set $j(k) = \text{mod}(k, m)$, $\alpha'_{k+1} = \alpha'_k$ for $i \in [1:m] \setminus \{j(k)\}$.
5: &nbsp;&nbsp;&nbsp;&nbsp;Compute $x_{k+1}, \alpha_{k+1}^{j(k)}$ with Linesearch_extrapolation($x_k, d_k^{(k)}, \alpha_k^{(k)}, \theta, \gamma, \rho$).
6: **end for**

<caption>**Algorithm 5:** Linesearch_extrapolation($x, d, \alpha, \theta, \gamma, \rho$).</caption>

1: if $f(x + \alpha d) \ge f(x) - \rho(\beta)$ then
2: &nbsp;&nbsp;&nbsp;&nbsp;Return $x, \theta, \alpha$.
3: **end if**
4: **while** $f(x + \alpha d) < f(x) - \rho(\alpha)$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;$$\alpha \leftarrow \gamma\alpha$$.
6: **end while**
7: **Return** $(x + \frac{\alpha}{\gamma}d, \alpha/\gamma)$.

A full algorithm that adapts Algorithm 1 so as to allow for line searches is given in Algorithm 4. When using line-search techniques, it becomes necessary to maintain a distinct stepsize for each direction in the poll set $D_k$, and to ensure some consistency in the choice of $D_k$. More precisely, each sequence of directions $\{d'_k\}_{k \in \mathbb{N}}$ must satisfy some specific condition like in [79] (e.g., $\{d'_k\}_{k \in \mathbb{N}}$ bounded and such that every limit point $\{\tilde{d}\}_{j \in [1:m]}$ of $\{d'_k\}_{j \in [1:m]}$ is a positive spanning set), and distinct sequences of step sizes $\{\alpha'_k\}_{k \in \mathbb{N}}$ for $j \in [1:m]$ must be defined. For the sequence of directions $\{d'_k\}_{k \in \mathbb{N}}$, a natural approach is to use the same direction throughout all iterations. Note there is no distinction between successful and unsuccessful steps, since every iteration is guaranteed to produce some point that violates the sufficient decrease test.

Provided the function $f$ is smooth in the sense of Definition 2.6 and all polling sets correspond to the same positive spanning set, one can obtain similar results to that of deterministic directional direct-search methods [79]. Thanks to the use of the sufficient decrease, any limit point of the sequence generated by Algorithm 4 is a stationary point.

**Theorem 3.11.** Suppose that we apply Algorithm 5 to an L-smooth function $f$ under Assumptions 2.2, 2.3, 2.4 and 3.1. Suppose further that sufficient decrease is used. Then,

$$ \lim_{k \to \infty} \| \nabla f(x_k) \| = 0. \qquad (19) $$

In addition to global convergence results, worst-case complexity bounds have recently been derived for line-search techniques in the smooth setting. Brilli et al. [26] report iteration complexity bounds and evaluation complexity bounds of $O(n\epsilon^{-2})$ and $O(n^2\epsilon^{-2})$,
respectively, for a particular version of a line-search method. Interestingly, the most efficient version in practice is endowed with complexity bounds that are worse by a factor of $n$.

To establish convergence of Algorithm 4 on nonsmooth problems, the positive spanning set assumption is replaced by an asymptotic density assumption on the polling directions, akin to that made for mesh-based techniques. In particular, the choice $\mathbb{D}_k = \{g_k\}$, where $g_k$ is a random direction in the unit sphere, satisfies this assumption, and leads to the following result that corresponds to a special case of the analysis in Diouane et al. [40].

**Theorem 3.12.** Suppose that we apply Algorithm 5 to a Lipschitz continuous function $f$, under Assumptions 2.1 and 2.4, with sufficient decrease being used. Let $\{x_k\}_{k \in \mathcal{K}}$ be a refining subsequence for the algorithm with refined point $\hat{x}$ such that $\mathbb{D}_k = \{g_k\}$ for every $k \in \mathcal{K}$, where $\mathcal{K} \subset \mathbb{N}$ denotes a subsequence of unsuccessful iterations. If any subsequence of $\{\mathbb{D}_k\}_{k \in \mathcal{K}}$ is dense in the unit sphere, then $\hat{x}$ is a Clarke stationary point.

Finally, we comment on the use of non-monotone acceptance rules, an important feature of line-search techniques that can significantly enhance both robustness and efficiency of those methods, while also applying to other directional direct-search techniques. A non-monotone counterpart to the sufficient decrease condition can be defined as follows [59].

**Definition 3.1.** Let $M$ be a positive integer. A non-monotone sufficient decrease condition for Algorithm 1 or Algorithm 5 consists of replacing (2) by

$$f(\mathbf{x}_k + \alpha_k \mathbf{d}_k) < f_k - \rho(\alpha_k), \quad f(\mathbf{x}_k) \le f_k \le \max_{j \in \{0, \min(k, M)\}} \{f(\mathbf{x}_{k-j})\}.$$

Grippo and Rinaldi [59] incorporate those non-monotone strategies into a wide range of algorithms that make use of simplex gradients, thus promoting an extensive integration of direct-search and model-based approaches. Standard direct-search techniques have also been combined with the non-monotone principle [54]. Other non-monotone line-search derivative-free algorithms have been proposed by Diniz-Ehrhardt et al. [39] and García-Palomares and Rodríguez [52]. In all cases, a global convergence result similar to that of Theorems 3.1 or 3.12 is established.

# 4. Handling noisy and stochastic values

This section focuses on direct-search methods for unconstrained optimization problems where the objective function cannot be evaluated exactly due to noise in the calculation. Standard direct-search methods such as those described in Section 4 have long been observed to be robust to noise up to a certain magnitude associated with the stepsize [65]. However, provably handling noise in function evaluations requires modifying the basic direct-search framework to establish (probabilistic) convergence guarantees. Our goal in this section is to provide the reader with key approaches in this very active area of research.

## 4.1. Direct-search methods with controlled noise

In this section, we consider Problem (1) under a controlled noise model. At any point $\mathbf{x} \in \mathbb{R}^d$, any query for the true function value $f(\mathbf{x})$ returns a noisy evaluation $f(\mathbf{x}, \zeta)^3$ where $\zeta$ captures the problem noise, and $f(\mathbf{x}, \zeta)$ represents the noisy function values that can be queried by the algorithm. A controlled noise model assumes that such that $|f(\mathbf{x}, \zeta) - f(\mathbf{x})|$ is uniformly bounded by a fixed value called noise level. This noise level typically prevents the method from converging below a certain accuracy, yet it is often possible to guarantee convergence up to that threshold.

The controlled noise model has been analyzed in various ways in the direct-search literature. Kelley's implicit filtering method is a directional direct-search algorithm equipped with global convergence guarantees that depend on the noise level [64]. The thesis of Pauwels follows a similar approach to derive complexity guarantees [84, Chapter 3]. In line-search techniques, Lucidi and Scian-drone [78] consider a line-search algorithm similar to those of Section 3.3 under the assumption that the noise level at every point considered at iteration $k$ is bounded by a deterministic value $\Delta_k > 0$. Provided the sequence $\left\{ \frac{\Delta_k}{\min_{j \in \{1, \dots, \lfloor \Delta_k \rfloor\}} \alpha_{k,j}} \right\}$ (where $\alpha_{k,j}$ denotes the stepsize at iteration $k$ for direction $j$) converges to 0 as $k \to \infty$, convergence to stationarity points is established on smooth problems.

In the bilevel setting, Diouane et al. [40] describe directional and mesh-based approaches motivated by bilevel optimization to handle controlled noise of the form $|f(\mathbf{x}) - f(\mathbf{x})| \le \epsilon_f$ for some noise level $\epsilon_f > 0$. By using a lower bound $\alpha_{\min}$ on the stepsize, and assuming that the method stops at the first iteration $k$ satisfying $\alpha_k = \alpha_{\min}$, it is possible to provide complexity guarantees for the method. In the smooth case, denoting by $\mathbf{x}_k$ the current iterate, one has

$$\|\nabla f(\mathbf{x}_k)\| = \mathcal{O}\left(\frac{\epsilon_f}{\alpha_{\min}} + \alpha_{\min}\right)$$

<footnote n="3">³ We recall from Section 2 that we underline quantities or objects that depend on some form of noise, although that noise may be deterministic in nature.</footnote>
after at most $\mathcal{O}(\alpha_{\min}^{-2})$ iterations. By setting $\alpha_{\min} = \mathcal{O}(\epsilon_f^{1/2})$, one guarantees a gradient accuracy of order $\mathcal{O}(\epsilon_f^{1/2})$ in $\mathcal{O}(\epsilon_f^{-1})$ iterations. In the nonsmooth case, the method converges asymptotically to a $(\alpha_{\min}, \mathcal{O}(\epsilon_f/\alpha_{\min} + \alpha_{\min}))$-Goldstein stationary point in the sense of Definition 2.10. Setting again $\alpha_{\min} = \mathcal{O}(\epsilon_f^{1/2})$ guarantees convergence to $(\mathcal{O}(\sqrt{\epsilon_f}), \mathcal{O}(\sqrt{\epsilon_f}))$-Goldstein stationary points. Note that the noise is stochastic in nature, but the uniform bound is deterministic. As a result, the aforementioned complexity results hold for every realization of the method.

Finally, we mention another approach that has been used to deal with controlled noise in stochastic optimization. Berahas et al. [17,18] proposed line-search methods based on a relaxed Armijo condition that explicitly involves the noise level value. In a smooth setting, these methods can converge to an approximate stationary point with gradient norm of order $\mathcal{O}(\sqrt{\epsilon_f})$ where $\epsilon_f > 0$ is the noise level. Although these techniques typically employ a gradient approximation, a similar approach could be applied in direct search by using a forcing function of the form $\rho(\underline{\alpha} - \mathcal{O}(\epsilon_f))$.

## 4.2. Direct-search methods for general random noise

In this section, we consider more general forms of noise on function values that are not necessarily bounded. A typical formulation of the problem in that setting is

$$ \min_{\mathbf{x} \in \mathbb{R}^n} f(\mathbf{x}) = \mathbb{E}_{\zeta} [f(\mathbf{x}, \zeta)], \qquad (20) $$

and we again assume that the algorithm can only evaluate noisy function values of the form $f(\mathbf{x}, \zeta)$. Problem (20) has attracted significant attention in the derivative-free optimization community beyond the direct-search paradigm. Under smoothness assumptions on $f$ and $f(\cdot, \zeta)$, line-search techniques that employ a stochastic estimation of the gradient $\nabla f(\mathbf{x})$ have been investigated [18,21,62,83]. Interestingly, recent proposals of that form, termed “step search”, modify the search direction along with the stepsize, thus bringing them closer to direct-search methods in spirit [21]. Extensions of the direct-search schemes reviewed in Section 4 have also been proposed to tackle Problem (20) [3,13,29,46,47,49,86]. Audet et al. [13] introduced StoMADS, which we view as a stochastic variant of Algorithm 2 for the purpose of this survey. At every iteration, this method uses estimates $\underline{f}_k^0 \approx f(\mathbf{x}_k)$ and $\underline{f}_k^s \approx f(\mathbf{x}_k + \mathbf{s})$ that are computed using the noisy function $f(\cdot, \zeta)$ or function values computed for different values of the noise parameter $\zeta$. One drawback of using such estimates is that the simple decrease condition $\underline{f}_k^s < \underline{f}_k^0$ does not necessarily imply $f(\mathbf{x}_k + \mathbf{s}) < f(\mathbf{x}_k)$. To circumvent this issue, Audet et al. [13] used so-called $\epsilon_f$-accurate estimates, inspired by derivative-free trust-region algorithms [30] and proved that a sufficient decrease in such estimates guarantees an improvement in $f$. This notion is stated below for a realization of the algorithm at hand (hence all quantities or objects considered are deterministic).

**Definition 4.1.** ([13, Definition 1] and [46, Definition 1].) Consider a realization of Algorithm 2 applied to Problem (20), and let $\mathbf{x}_k$ be the incumbent at iteration $k$. Then, for any $\mathbf{s} \in \mathbb{R}^n$ and any $\epsilon_f > 0$, a value $\underline{f}_k^s$ is called an $\epsilon_f$-accurate estimate of $f(\mathbf{x} + \mathbf{s})$ for $\alpha_k^p$ if

$$ |\underline{f}_k^s - f(\mathbf{x}_k + \mathbf{s})| \le \epsilon_f (\alpha_k^p)^2. $$

Similar to the analysis in Section 3.1.2, probabilistic satisfaction of the accuracy property can be considered.

**Definition 4.2.** A sequence $\{\underline{f}_k^0, \underline{f}_k^s\}$ of random estimates produced by the stochastic version of Algorithm 2 is called $\beta$-probabilistically $\epsilon_f$-accurate with respect to the corresponding sequence $\{\underline{\mathbf{x}}_k, \underline{\mathbf{s}}_k, \underline{\alpha}_k^p\}$, if

$$ \mathbb{P} \left( \left\{ \left| \underline{f}_k^0 - f(\underline{\mathbf{x}}_k) \right| \le \epsilon_f (\underline{\alpha}_k^p)^2 \right\} \cap \left\{ \left| \underline{f}_k^s - f(\underline{\mathbf{x}}_k + \underline{\mathbf{s}}_k) \right| \le \epsilon_f (\underline{\alpha}_k^p)^2 \right\} \mid \mathcal{F}_{k-1} \right) \ge \beta. \qquad (21) $$

The convergence analysis of StoMADS uses the following key assumptions on the nature of the stochastic information.

**Assumption 4.1.** The sequence $\{\underline{f}_k^0, \underline{f}_k^s\}$ of random estimates is $\beta$-probabilistically $\epsilon_f$-accurate for some $\beta \in (1/2, 1)$. Moreover, for any $k \ge 0$, one has

$$ \mathbb{E} \left[ \left| \underline{f}_k^s - f(\underline{\mathbf{x}}_k + \underline{\mathbf{s}}_k) \right|^2 \mid \mathcal{F}_{k-1} \right] \le \kappa_f^2 (\underline{\alpha}_k^p)^4 \quad \text{and} \quad \mathbb{E} \left[ \left| \underline{f}_k^0 - f(\underline{\mathbf{x}}_k) \right|^2 \mid \mathcal{F}_{k-1} \right] \le \kappa_f^2 (\underline{\alpha}_k^p)^4 \qquad (22) $$

almost surely for some $\kappa_f \ge 0$.

Random estimates satisfying Assumption 4.1 can easily be obtained when the function estimates have bounded variance, i.e. there exists $V > 0$ such that $\mathbb{V}_{\zeta} [f(\mathbf{x}, \zeta)] \le V$ for all $\mathbf{x} \in \mathbb{R}^n$. Under this assumption, one constructs the estimates by averaging independent samples as follows:

$$ \underline{f}_k^0 := \frac{1}{p_k} \sum_{i=1}^{p_k} f(\underline{\mathbf{x}}_k, \underline{\zeta}_i^0) \quad \text{and} \quad \underline{f}_k^s := \frac{1}{p_k} \sum_{i=1}^{p_k} f(\underline{\mathbf{x}}_k + \underline{\mathbf{s}}_k, \underline{\zeta}_i^s). \qquad (23) $$
Provided $p_k \ge \frac{V}{\varepsilon_f^2(1-\sqrt{\beta})(\alpha_k^p)^4}$, the estimates (23) satisfy Assumption 4.1 with $\kappa_f^2 \ge \varepsilon_f^2(1-\sqrt{\beta})$.

Drawing inspiration again on derivative-free trust-region methods [20,30,83], the analysis of StoMADS provides a condition on the probability $\beta$ for guaranteeing convergence. This condition depends on an auxiliary parameter $\nu$ corresponding to the Lyapunov function $\nu(f(x_k) - \min_{x \in \mathbb{R}^n} f(x)) + (1 - \nu)\{\alpha_k^p\}^2$. The requirements on $\nu$ (and thus $\beta$) arise naturally from the analysis, and as such are difficult to draw intuition from. On the other hand, we note that the bounds (24) can be estimated for specific parameter values [30, Remark 4.13].

**Theorem 4.1.** [13, Theorem 1 and Corollary 1] Under Assumptions 2.1 and 4.1, if $\nu \in (0, 1)$ and $\beta \in (1/2, 1)$ satisfy

$$ \frac{\nu}{1-\nu} \ge \frac{2(\theta^{-4}-1)}{\varepsilon_f(\gamma-2)} \quad \text{and} \quad \frac{\beta}{\sqrt{1-\beta}} \ge \frac{4\nu\kappa_f}{(1-\nu)(1-\theta^2)}, \qquad (24) $$

the sequence $\{\alpha_k^m\}_{k \in \mathbb{N}}$ of random mesh size parameters satisfies $\sum_{k=0}^{\infty} \alpha_k^m < \infty$ almost surely. Consequently, $\lim_{k \to +\infty} \alpha_k^m = 0$ almost surely.

The last result of Theorem 4.1 shows that on an event of probability one, the meshes used by StoMADS become infinitely fine. Such a zeroth-order convergence result is stronger than that of Algorithm 2 (recall Theorem 3.8), and more generally of mesh-based techniques that do not rely upon sufficient decrease. In fact, sufficient decrease appears necessary to obtain convergence and complexity results not only for direct-search algorithms, but also for other derivative-free techniques [30].

The existence of a convergent refining subsequence was established for StoMADS under the assumption that all the iterates generated by the method belong to a compact set [13, Theorem 2]. Due to possibly inaccurate estimates, the algorithm can accept points that increase the actual function value, and may thus venture outside the initial level set. As a result, Assumption 2.2 cannot be used to establish convergence. Still, a compactness assumption yields the desired result.

**Theorem 4.2.** Under the assumptions of Theorem 4.1, suppose that all the iterates generated by StoMADS lie in a compact set for all realization. Then, there exists a refining subsequence $\{x_k\}_{k \in \mathcal{K}}$ defined by a random index set $\mathcal{K}$ that converges almost surely to a refined point $\hat{x}$.

Finally, the analysis of StoMADS provides a stationarity result for the nonsmooth setting, akin to Theorem 3.9.

**Theorem 4.3.** Let the assumptions of Theorem 4.2 hold. Then, there exists an almost-sure event $\mathcal{E}$ such that, conditioned on $\mathcal{E}$, for any refined point $\hat{x} \in \mathbb{R}^n$ and any refining direction $\hat{d} \in \mathbb{R}^n$ for $\hat{x}$, one has

$$ f^{\circ}(\hat{x}; \hat{d}) \ge 0. $$

The argument and analysis of StoMADS was adapted to directional direct-search methods by Dzahini [46], giving rise to the Stochastic Directional Direct-Search (SDDS) algorithm, which can be viewed as the extension of Algorithm 1 to handle noisy function values. The SDDS framework which, unlike StoMADS, does not require that trial points belong to a mesh, relies on a careful update parameter for the stepsize $\alpha_k$, as well as sufficient decrease [46, Section 2.1] based on a forcing function of the form $\rho(\alpha) = c \alpha^q$, $\alpha > 0$. Using martingale-based analysis [20] and assuming $\beta$-probabilistically $\varepsilon_f$-accurate estimates, one can show that the expected number of iterations required by SDDS algorithms to drive $\|\nabla f(x_k)\|_2$ below a given threshold $\varepsilon \in (0, 1)$ is of order $\mathcal{O}\left(\frac{\varepsilon^{-\frac{q}{\min(q-1, 1)}}}{(2\beta - 1)}\right)$ [46, Theorem 4]. More recently, Rinaldi et al. [86] proved convergence of a stochastic directional direct-search scheme for nonsmooth objectives. The convergence analysis uses a weak tail bound condition on the noise, a weaker assumption than prior results, that still leads to almost-sure convergence.

**Sample complexity** Sample complexity refers to bounding the number of evaluations of the noisy function estimates (also called samples) required to achieve a desired convergence level. To the best of our knowledge, there are no bounds available on the number of samples required to achieve a given accuracy for direct-search methods with adaptive sampling, despite recent developments in the general stochastic optimization setting [63]. On the other hand, bounds on the number of samples per iteration have been obtained in the direct-search literature. A first bound of that nature was given by Audet et al. [13], where estimates (23) required $\mathcal{O}(\alpha_k^{-4})$ samples in order to obtain convergence guarantees. A similar bound in $\mathcal{O}(\alpha_k^{-4})$ can be showed for SDDS [46]. An improvement to $\mathcal{O}(\alpha_k^{-2q})$ samples per iteration was proved by Rinaldi et al. [86] based on using a forcing function of the form $\rho(\alpha) = c \alpha^q$ with $c > 0$. Assuming a common random number generator can be used to sample the function estimates, and that $f(\cdot, \cdot)$ is Lipschitz continuous with a Lipschitz constant that is uniform across all $\zeta$, the sampling cost improves by a factor of $\alpha_k^2$ at every iteration [86]. Another notable result is the $\mathcal{O}(\log T)^{\frac{2}{3}} T^{\frac{2}{3}})$ regret bound with respect to a sample budget $T$ for stochastic direct-search applied to smooth strongly convex objectives proved by Achddou et al. [2]. The same authors introduce a novel technique called sequential sampling, that allows for fewer samples per iteration in cases where the sufficient decrease condition is satisfied or violated with a certain margin.
## 5. Handling constraints

In this section, we consider direct-search methods designed for constrained optimization problems of the form

$$ \min_{x \in \mathbb{R}^n} f(x) \quad \text{s.t.} \quad x \in \Omega, \qquad (25) $$

where the feasible set $\Omega$ is a subset of $\mathbb{R}^n$ that encodes constraints on the problem variables. As in the rest of the paper, information about the objective function is only available through function evaluations. As for the constraint set, access to derivative information is common for simple constraints, but is not necessary to design a direct-search scheme. In fact, constraints in derivative-free information are subject to a careful classification, according to the information they provide regarding feasibility or lack thereof [71]. In the context of this survey, we classify a constraint as *relaxable* or *unrelaxable* depending on whether the objective function and all other constraints can be queried at points that do not satisfy the first constraint. We also distinguish *quantifiable* and *non-quantifiable* constraints: the former corresponds to constraints that produce a numerical value, while the latter can be viewed as a binary certificate of constraint satisfaction.

The first part of this section is concerned with feasible direct-search methods, that aim at preserving feasibility of the iterates throughout the optimization process. Section 5.1 surveys theoretical guarantees that have been obtained in the presence of linear constraints, manifold constraints or non-relaxable constraints. We then move to infeasible direct-search methods, for which feasibility can only be guaranteed asymptotically, in Section 5.2.

### 5.1. Feasible methods

In this section, we consider feasible direct-search algorithms, that are particularly suited for unrelaxable constraints. Algorithm 6 describes a basic variation on Algorithm 1 that ensures the objective to only be evaluated at feasible points.

In Section 5.1.1, we discuss the case of linear constraints, that has been tackled by the main classes of direct-search algorithms. Sections 5.1.2 and 5.1.3 survey other feasible iterate strategies, based on Riemannian manifolds and extreme barrier, respectively.

<caption>**Algorithm 6: Feasible direct-search scheme.**</caption>

1: **Inputs:** Feasible point $x_0 \in \Omega$, initial stepsize $\alpha_0 > 0$, stepsize update parameters $0 < \theta < 1 \le \gamma$, forcing function $\rho: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$.
2: **for** $k = 0, 1, 2, \dots$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Select a set $D_k$ of poll directions.
4: &nbsp;&nbsp;&nbsp;&nbsp;**if** $x_k + s \in \Omega$ and $f(x_k + s) < f(x_k) - \rho(\alpha_k)$ hold for some $s \in \{\alpha_k \mid d: d \in D_k\}$ **then**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Declare the iteration as successful, set $x_{k+1} = x_k + s$ and $\alpha_{k+1} = \gamma\alpha_k$.
6: &nbsp;&nbsp;&nbsp;&nbsp;**else**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Declare the iteration as unsuccessful, set $x_{k+1} = x_k$ and $\alpha_{k+1} = \theta\alpha_k$.
8: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
9: **end for**

#### 5.1.1. Linear constraints

Linear constraints, and in particular bound constraints, are the simplest kind of constraints in optimization. In direct-search techniques, linear constraints can be easily handled by a careful choice of polling directions [65]. We present below the main tools and results for this purpose in the context of directional direct-search schemes, then discuss related results for mesh-based and line-search techniques.

Suppose that the feasible set $\Omega$ of Problem (25) is polyhedral and described by linear inequalities, i.e.,

$$ \Omega = \{x \in \mathbb{R}^n \mid Ax \le b\} \qquad (26) $$

for $A \in \mathbb{R}^{n \times n}$ and $b \in \mathbb{R}^l$. Without loss of generality, we assume that the rows of $A$ are unit vectors. Convergence results for Algorithm 6 in that setting aim at proving that the quantity $\chi(x_k)$ goes to zero, where

$$ \forall x \in \Omega, \chi(x) = \max_{\substack{x+d \in \Omega \\ \|d\| \le 1}} \langle d, [-\nabla f(x)] \rangle. \qquad (27) $$

The measure $\chi(x)$ is always nonnegative, and equals zero if and only if $x$ is a first-order KKT stationary point of Problem [65].

The analysis of direct-search methods in that setting relies heavily on the notion of approximate tangent and normal cones, that involve both a feasible point and a tentative stepsize.

**Definition 5.1.** Given $x \in \Omega$ and $\alpha > 0$, the approximate normal cone $N_{\Omega}(x, \alpha)$ is defined by

$$ N_{\Omega}(x, \alpha) = \text{pspan}(\{a_i \mid b_i - \langle a_i, x \rangle \le \alpha\}), $$

where pspan($\mathcal{A}$) denotes the set of nonnegative linear combinations of vectors in $\mathcal{A}$. The approximate normal cone is then defined as $T_{\Omega}(x, \alpha) = N_{\Omega}(x, \alpha)^{\circ}$.
It follows from Definition 5.1 that a displacement of length $\alpha$ along any direction in the approximate tangent cone preserves feasibility, while the normal cone is generated by normals corresponding to constraints that can be violated with such a displacement. A natural strategy for defining polling directions in Algorithm 6 then consists of using generators of the approximate tangent cone at the current iterate [66, Proposition 2.1].

Indeed, the cosine measure defined in Section 2 can be generalized to the linearly constrained setting [58,65] as follows.

**Definition 5.2.** Given $x \in \Omega$ and $\alpha > 0$, let $\pi_{T_{\Omega}(x, \alpha)}[-\cdot]$ denote the orthogonal projection onto the approximate tangent cone $T_{\Omega}(x, \alpha)$. Given $x \in \Omega$, $\alpha > 0$, and a direction set $\mathcal{D}$, the cosine measure generalized to linear constraints is defined by

$$ \mathrm{cm}_{T_{\Omega}(x, \alpha)}(\mathcal{D}) = \inf_{\mathbf{v} \in \mathbb{R}^n} \max_{d \in \mathcal{D}} \frac{\langle d, \mathbf{v} \rangle}{\|d\| \|\mathbf{v}\|} \qquad (28) $$

When $\mathcal{D}$ contains generators of $T(x, \alpha)$, it can be shown that this quantity is positive [66], and thus such directions guarantee feasible descent, i.e. preserve feasibility and yield descent from $x$ with a stepsize $\alpha$ as long as $\pi_{T_{\Omega}(x, \alpha)}[-\nabla f(x)] > 0$. Although the latter quantity is not a good measure of stationarity, it bears a natural connection with the measure $\chi$, described below.

**Lemma 5.1.** ([66, Proposition B.2]) Suppose that $\Omega$ is defined via the polyhedral description (26), and that there exists $B_g > 0$ such that $\|\nabla f(x)\| \le B_g$ for any $x \in \Omega$. Then, for any $x \in \Omega$ and any $\alpha > 0$, there exists $\eta > 0$ such that

$$ \chi(x) \le \pi_{T_{\Omega}(x, \alpha)}[-\nabla f(x)] + \frac{\alpha}{\eta} B_g. \qquad (29) $$

Note that $\eta$ can be computed a priori by considering all possible normal cones at $x$ [58].

Lemma 5.1 is instrumental to deriving global convergence guarantees for Algorithm 6. Assuming that generators for the approximate tangent cone are used, Kolda et al. [65] established a global convergence result of the form $\liminf_{k \to \infty} \chi(x_k) = 0$. A complexity bound was later established by Gratton et al. [58], and is summarized below.

**Theorem 5.1.** ([58, Theorem 4.1]) Suppose that $\Omega$ is defined via the polyhedral description (26), and that there exists $B_g > 0$ such that $\|\nabla f(x)\| \le B_g$ for any $x \in \Omega$. Suppose further that $f$ satisfies Assumptions 2.2 and 2.1. Finally, suppose that Algorithm 6 is applied using $\rho(\alpha) = c \alpha^2$ with $c > 0$ and choosing $\mathcal{D}_k$ as a set of unit norm generators for $T(x_k, a_k)$ at every iteration with $|\mathcal{D}_k| \le m$ and $\mathrm{cm}_{T_{\Omega}(x_k, a_k)}(\mathcal{D}_k) \ge \tau > 0$. Then, for any $\epsilon > 0$, the method reaches an iterate such that $\chi(x_k) \le \epsilon$ in at most

$$ \mathcal{O}(\tau^{-2}\epsilon^{-2}) \text{ iterations and } \mathcal{O}(m\tau^{-2}\epsilon^{-2}) \text{ function evaluations.} \qquad (30) $$

Although the bounds can be shown to hold for finite values of $m$ and positive values of $\tau$, precise values for these quantities are only available in the special case of bound constraints and linear equality constraints. When only bound constraints are enforced, using coordinate vectors and their negatives as poll directions yields $\tau = 1/\sqrt{n}$ and $m = 2n$. When the constraint set is defined using $1 \le n' \le n$ linearly independent equality constraints, the tangent cone is a linear subspace of $\mathbb{R}^n$. Using orthonormal basis vectors for that subspace and their negatives as poll directions yields $\tau = 1/\sqrt{n'}$ and $m = 2n'$ [58]. We note that a generalization of Theorem 5.1 through probabilistic analysis was described in Gratton et al. [58], showing that high-probability complexity results could also be derived akin to those of Theorem 3.4. Although the random polling sets to be used must be restricted to preserve feasibility, it is not clear that the complexity improves compared to the deterministic setting, except in the context of linear equality constraints. Nevertheless, by decomposing all tangent cones as the sum of pointed cones and subspaces, random directions can be drawn in the former, leading to a possible improvement of the number of directions used at every iteration.

We note that linear constraints have also been investigated in the context of mesh-based algorithms, with more emphasis on handling linear inequality constraints and bounds to obtain convergence results [1,85]. Although those result apply in both a smooth and a nonsmooth setting, no complexity guarantees are available. Furthermore, in a line-search setting, we note that the feasible direction framework based on approximate tangent cones was used in the context of bound constraints in Lucidi and Sciandrone [78] to establish convergence guarantees. The idea was then extended to general algebraic unrelaxable smooth constraints through the use of projections, leading again to global convergence guarantees [80]. Finally, Cristofari and Rinaldi [34] defined feasible direct-search methods for problems with structured feasible set (i.e., $\Omega = \mathrm{conv}(A)$ with $A = \{v_1, \dots, v_r\}$, $v_i \in \mathbb{R}^n$) that use directions pointing towards/away from the vertices of the feasible set, and are somehow related to zeroth-order Frank-Wolfe methods (see, e.g., [22] and references therein for further details on the Frank-Wolfe method).

### 5.1.2. Riemannian direct-search schemes

In this section, we discuss an emerging trend in feasible methods related to Riemannian optimization, a well-established subfield of nonlinear optimization that is concerned with finding the best value of a function over a smooth manifold [25]. Interestingly, the idea of extending direct-search methods to the Riemannian setting started early in the development of numerical Riemannian optimization. The first direct-search proposal for Riemannian optimization problems is due to Dreisigmeyer [44], and consists of an extension of mesh-based algorithm to Riemannian manifolds. The special case of reductive homogeneous spaces, that includes several
matrix manifolds, is discussed in a follow-up work, where both deterministic and probabilistic directional methods are proposed [45]. Those references focus on properly defining the algorithms rather than on establishing convergence guarantees.

More recently, Kungurtsev et al. [69] describe basic direct-search schemes of directional and line-search types for minimizing smooth and nonsmooth objectives on Riemannian manifolds. Global convergence properties are established, that match those of the unconstrained setting. More precisely, counterparts to Theorem 3.2 and 3.5 are established for directional variants on smooth and Lipschitz continuous objectives, respectively, while counterparts to Theorem 3.11 and 3.12 are obtained for variants using line-search extrapolation. We note, however, that complexity results have yet to be derived for Riemannian direct-search schemes.

### 5.1.3. Extreme barrier

The extreme barrier approach is arguably the most general way to handle constraints in direct-search methods so as to enforce feasibility. It consists of setting the objective equal to +∞ outside the feasible set. The first proof of convergence (to Clarke stationary points) of an extreme barrier approach was done for mesh adaptive direct-search algorithms [7]. Consequently, most of the theory in Section 3.2 extends to the case of a constrained set handled through extreme barrier. We present below the main ingredients to derive convergence results in the context of Algorithm 2 below, but note that similar arguments can be applied to any direct-search method with a dense set of poll directions. Akin to Section 5.1.1, we start by defining suitable cones that characterize stationarity.

**Definition 5.3.** The hypertangent cone $T_{\Omega}^{H}(x)$ to $\Omega$ at $x \in \Omega$ is the set of vectors $d \in \mathbb{R}^n$ such that there exist $\epsilon > 0$, $y \in \Omega \cap B_{\epsilon}(x)$ and $v \in B_{\epsilon}(d)$ satisfying $y + tv \in \Omega$ for any $t \in (0, \epsilon)$.

**Definition 5.4.** The Clarke tangent cone $T_{\Omega}(x)$ to $\Omega$ at $x \in \Omega$ is the set of vectors $d \in \mathbb{R}^n$ such that there exist sequences $\{t_k\}$, $\{x_k\}$, $\{d_k\}$ satisfying $t_k \to 0^+$, $x_k \to x$, $d_k \to d$ and $x_k + t_k d_k \in \Omega$ for every $k$.

The set $\Omega$ is called regular at $x$ if $T_{\Omega}(x)$ is the closure of $T_{\Omega}^{H}(x)$.

The use of hypertangent cones allows for generalizing Lemma 2.2 to the extreme barrier setting.

**Lemma 5.2.** Suppose that $f$ is Lipschitz continuous, and consider that Algorithm 2 is applied to (25) using an extreme barrier. Let $\mathcal{K} \subset \mathbb{N}$ be an index subsequence such that $\{x_k\}_{k \in \mathcal{K}}$ is a refining subsequence with refined point $\hat{x}$. Then, for any refining direction $d \in T_{\Omega}^{H}(\hat{x})$,

$$f^{\circ}(\hat{x}; d) \geq 0. \qquad (31)$$

Similarly, using the Clarke tangent cone together with the result of Lemma 5.2 leads to the following generalization of Theorem 3.5.

**Theorem 5.2.** Under the assumptions of Lemma 5.2, let $\mathcal{K}$ be an index set corresponding to unsuccessful iterations for Algorithm 2 applied with extreme barrier. Suppose further that $\lim_{k \in \mathcal{K}} \alpha_k = 0$, $x_k \to \hat{x}$, and that for every direction $d \in T_{\Omega}^{H}(\hat{x})$ there exists $\mathcal{L} \subset \mathcal{K}$ such that $\lim_{k \in \mathcal{L}} \frac{d_k}{\|d_k\|} = \frac{d}{\|d\|}$ with $d_k \in \mathbb{D}_k$. Finally, suppose that the set $\Omega$ is regular at $\hat{x}$. Then the point $\hat{x}$ is Clarke stationary for Problem (25).

The extreme barrier approach thus provides a straightforward way of extending a direct-search method to constraint sets while preserving convergence guarantees. On the other hand, when the constraints are quantifiable, using a progressive barrier prevents from getting any information regarding constraint violation. Infeasible approaches have been developed with that issue in mind, and we review them in the next section.

## 5.2. Infeasible methods

This section is concerned with direct-search schemes that allow for evaluating the objective outside of the feasible set. Such methods are termed infeasible, and apply to problems with quantitative relaxable constraints of the form

$$\min_{x \in \mathbb{R}^n} f(x) \quad \text{s.t.} \quad g(x) \le 0, \qquad (32)$$

where $g: \mathbb{R}^n \to \mathbb{R}^l$. Thanks to the quantitative nature of the constraints, it is possible to define an infeasibility measure $h(\cdot)$ that quantifies the distance to infeasibility. Typical choices for such measures are

$$h(x) = \sum_{i=1}^{l} \max\{0, g_i(x)\}, \qquad (33)$$

which is typically nonsmooth, and

$$h(x) = \sum_{i=1}^{l} \max\{0, g_i(x)^2\} \qquad (34)$$

which is a smooth penalty function for smooth functions $g_i$.
The theory for such progressive barrier functions was developed in Audet and Dennis [8]. A two-phase approach for solving Problem (32) using h consists of applying a suitable optimization method to the problem

$$ \min_{x \in \mathbb{R}^n} h(x) \qquad (35) $$

in order to find a feasible point. Once such a point is obtained, the original Problem (32) can be tackled using the extreme barrier approach of Section 5.1.3.

When the derivatives of g are not available, one may consider using Algorithm 1 or any of the methods discussed in Section 3 in the first phase. To this end, the sufficient decrease condition (2) must be modified so as to accept any feasible point that is encountered by the algorithm. However, note that the local nature of direct-search methods only guarantees convergence to a stationary point of the barrier function h. Still, provided the iterates converge towards the feasible set, a feasible point can be found in finitely many iterations. In that case, convergence results follow from that of the extreme barrier approach, and can be obtained for both progressive barrier functions (33) and (34).

Since the two-phase approach ultimately reduces to an extreme barrier approach, it suffers from the same shortcoming in that information at infeasible points eventually stops being used. The next three sections describe algorithms that exploit infeasibility information. Note that these variants are inspired by derivative-based constrained optimization techniques and, as such, possess the same advantages and drawbacks that their derivative-based counterparts. Thus, rather than formally stating their convergence results, we provide a summary of the key guarantees achievable by each algorithmic strategy, drawing from the main definitions and vocabulary of Section 5.1.

### 5.2.1. Filter-based methods

A filter approach consists of a biobjective view of the constrained optimization Problem (32), where the two objectives are the original function f and the progressive barrier h. In a filter-based algorithm, a set of “incumbent” points $S_k$ is maintained throughout the algorithm, and the corresponding values of f and h are stored. At every iteration, $S_k$ is updated by adding new points and discarding points dominated by new additions, where x is said to be dominated by y (written $x < y$) if

$$ [f(x) < f(y) \text{ and } h(x) \le h(y)] \text{ or } [f(x) \le f(y) \text{ and } h(x) < h(y)]. $$

With this definition, the set $S_k$ may consist of both feasible and infeasible points, the latter corresponding to better objective values.

The first direct-search method employing a filter is a mesh-based algorithm proposed by Audet and Dennis [6], but that method was not proven to converge towards stationary points for the original problem. The mesh-based algorithm of Dennis et al. [38] possess such a guarantee, provided both the constraint functions and the objective are differentiable (with approximate gradients available for the constraints). The progressive barrier of Audet and Dennis [8,11] combines the mesh adaptive direct-search paradigm with a filter based on a progressive barrier function. Algorithm 7 describes an analogous direct-search method based on Algorithm 1. Note that this choice is made for simplicity, and that one could have presented a variant on Algorithm 2 for the same purpose.

<caption>**Algorithm 7: Progressive barrier direct-search method.**</caption>

1: **Inputs:** $\emptyset \neq S_0 \subset \mathbb{R}^n$, $\alpha_0 > 0$, $0 < \theta < 1 \le \gamma$, $\rho: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$, $h_0 > 0$.
2: **for** $k = 0, 1, 2, \dots$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Select a set $D_k$ of poll directions. Select a feasible point $x_k^F \in S_k$ if such a point exists, and an infeasible point $x_k^I \in S_k$ if such a point exists.
4: &nbsp;&nbsp;&nbsp;&nbsp;**if** $y + s < \rho(\alpha_k)$ **y** holds for either $y = x_k^F$ or $y = x_k^I$ and $s \in \{\alpha_k d: d \in D_k\}$ **then**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Declare the iteration as successful, add $y + s$ to $S_k$, set $h_{k+1} = h_k$, and $\alpha_{k+1} = \gamma\alpha_k$.
6: &nbsp;&nbsp;&nbsp;&nbsp;**else**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $x_k^I$ was selected and there exists $s \in \{\alpha_k d: d \in D_k\}$ such that $h(x_k^I + s) < h_k - \rho(\alpha_k)$ and $f(x_k^I + s) > f(x_k^I)$ **then**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Declare the iteration as improving, add $x_k^I + s$ to $S_k$, and choose $h_{k+1} \le h_k - \rho(\alpha_k)$.
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $\alpha_{k+1} = \alpha_k$.
9: **else**
10: **end if**
11: **end if**
12: **end if**
13: Delete all the dominated points from $S_{k+1}$.
14: **end for**

At every iteration, Algorithm seeks improvement with respect to either a feasible point (if one has been found) or an infeasible point. A successful iteration corresponds to finding a point that dominates one of the two points. Note that Algorithm 7 uses special dominance relations of the form $x <_{\rho} y$ for some $\rho > 0$, that correspond to $x < y$ and $f(x) < f(y) - \rho$ or $h(x) < h(y) - \rho$. If one dominant point is found, this point is added to the set $S_k$. Otherwise, one checks whether an infeasible point with improved barrier value can be found, even if it increases the objective value, and such a point is added to $S_k$ if it exists. As a result, the method tries to improve both objectives or puts the emphasis on constraint violation.

In a mesh-based setting, Audet and Dennis [8] show convergence guarantees for refining sequences of the algorithm. If such a sequence consists of feasible points, convergence to a (Clarke) stationary point can be established under density of the refining directions. If the sequence consists of infeasible points, only convergence to a stationary point for Problem (35) is guaranteed.
### 5.2.2. Restoration-based methods

Restoration-based approaches rely on alternating steps that improve the objective with so-called restoration steps, which focus on improving the measure of infeasibility $h$. The approach of Bueno et al. [27] handles differentiable constraints, and alternates a gradient-based restoration phase with direct-search steps on a problem with linearized versions of the constraints. In this smooth setting, convergence to a KKT point for the original Problem (32) can be established.

The method of Gratton and Vicente [56] considers two conditions to accept a point at iteration $k$, namely

$$h(\mathbf{x}_k) - h(\mathbf{x}_k + \alpha_k \mathbf{d}_k) > \rho(\alpha_k) \quad \text{and} \quad h(\mathbf{x}_k) > C\rho(\alpha_k) \qquad (36)$$

for some $C > 1$, and

$$M(\mathbf{x}_k, \mu) - M(\mathbf{x}_k + \alpha_k \mathbf{d}_k, \mu) > \rho(\alpha_k) \qquad (37)$$

for some $\mu > 0$, where $M(\mathbf{x}, \mu) = f(\mathbf{x}) + \mu h(\mathbf{x})$ is referred to as a merit function. Within the restoration phase, only (36) is considered to accept a new point. Convergence results for this algorithm are established using the nonsmooth barrier (33). For nonsmooth objective, stationarity of refining points is established with a similar distinction than in the filter-based setting of Audet and Dennis [8]. In addition, results are specialized to the case of smooth objectives, thereby avoiding to require density of the refining directions.

### 5.2.3. Penalty approaches

In penalty based approaches, a direct-search scheme is applied to an auxiliary objective, that combines a measure of infeasibility with the objective function. The augmented Lagrangian technique proposed by Lewis and Torczon [72] handles nonlinear relaxable equality constraints as well as bounds on the objective. This technique was later extended to handle linear constraints in addition to bound constraints [67]. In both cases, a directional method is used to solve the augmented Lagrangian subproblems. The underlying theoretical guarantees rely on using the stepsize as a stopping criterion for solving the subproblems.

In another line of work, Liuzzi et al. [75] proposes an exact penalty approach to reformulate problem constrained with nonlinear, relaxable inequalities. The resulting penalized objective is defined on an open set, and tackled with a line-search method akin to those described in Section 3.3. Liuzzi et al. [76] considers a sequential penalty approach for problems with smooth objective, inequality constraints and bound constraints. In this approach, the inequality constraints other than bounds are relaxed and a line-search method is applied to minimize $f(\mathbf{x}) + \frac{1}{\epsilon}h(\mathbf{x})$ under bound constraints, where $h$ is the smooth infeasibility measure (34) and $\epsilon > 0$ is decreased based on a condition involving an upper bound on the stepizes. The case of nonsmooth relaxable constraints combined with bounds was addressed by Fasano et al. [50]. Under constraint qualification and for small enough $\epsilon$, it was then shown that all the stationary points corresponding to minimizing $f(\mathbf{x}) + \frac{1}{\epsilon}h(\mathbf{x})$ with $h$ defined by (33) are also (Clarke) stationary for the original problem.

## 6. Multiobjective direct-search methods

In this section, we survey the main developments in extending direct-search schemes to Multiobjective optimization (MOO) problems. The theory of multiobjective continuous optimization has undergone significant development since the 2000s, allowing for further development in the context of direct-search methods. Moreover, multiobjective problems are ubiquitous in engineering problems, a primary source of applications for direct-search methods.

Directional methods were extended to the multiobjective setting by Custódio et al. [36], along with convergence to a stationary point in the sense of Pareto. More recently, complexity bounds were derived for a directional direct-search scheme [37]. Mesh-based algorithms were considered by Audet et al. [10,12], the bi-objective version being successfully implemented in the NOMAD software [15]. Finally, line-search techniques are amenable to multiobjective extension. Liuzzi et al. [77] propose an analysis based on new optimality conditions. We consider MOO problems of the form

$$\min_{\mathbf{x} \in \mathbb{R}^n} F(\mathbf{x}) \equiv (f_1(\mathbf{x}), \dots, f_l(\mathbf{x}))^\top, \qquad (38)$$

where $f_i : \mathbb{R}^n \to \mathbb{R}$, with $i = 1, \dots, l$. In this context the concept of Pareto dominance is used to compare two points.

**Definition 6.1 (Pareto dominance).** Given two points $\mathbf{x}, \mathbf{y} \in \mathbb{R}^n$, we say that $\mathbf{x}$ dominates $\mathbf{y}$, when $F(\mathbf{x}) \le F(\mathbf{y})$, i.e., $f_i(\mathbf{x}) \le f_i(\mathbf{y})$ for all $i = 1, \dots, l$ and $F(\mathbf{x}) \ne F(\mathbf{y})$.

A point $\hat{\mathbf{x}}$ is then a local Pareto minimizer if there exists no point $\mathbf{y} \in B(\hat{\mathbf{x}}, \delta)$ s.t. $F(\mathbf{y}) \le F(\hat{\mathbf{x}})$, for some $\delta > 0$. In both the nonsmooth and the smooth settings, one can define necessary optimality conditions as follows.

**Definition 6.2.** Let $F$ be Lipschitz continuous near $\hat{\mathbf{x}}$. We say that $\hat{\mathbf{x}}$ is a Pareto-Clarke stationary point if

$$\forall \mathbf{d} \in \mathbb{R}^n, \exists i_d \in [1: l] : f_{i_d}^{\circ}(\hat{\mathbf{x}}; \mathbf{d}) \ge 0. \qquad (39)$$

**Definition 6.3.** Let $F$ be continuously differentiable at $\hat{\mathbf{x}}$. We say that $\hat{\mathbf{x}}$ is Pareto stationary, if

$$\forall \mathbf{d} \in \mathbb{R}^n, \exists i_d \in [1: l] : \langle \nabla f_{i_d}(\hat{\mathbf{x}}), \mathbf{d} \rangle \ge 0. \qquad (40)$$
It follows that the quantity

$$ \mu(\mathbf{x}) := - \min_{\|\mathbf{d}\| \le 1} \max_{i \in \{1:1\}} \langle \nabla f_i(\mathbf{x}), \mathbf{d} \rangle, \quad (41) $$

can be used to characterize Pareto stationarity. In our context, a finite set of polling directions is however used at each iteration. An approximation of $\mu(\mathbf{x})$ is hence often used in the analysis to take into account this aspect [37]:

$$ \mu_{\mathbb{D}_k}(\mathbf{x}) := - \min_{\mathbf{d} \in \mathbb{D}_k, \|\mathbf{d}\| \le 1} \max_{i \in \{1:1\}} \langle \nabla f_i(\mathbf{x}), \mathbf{d} \rangle. \quad (42) $$

We here focus on the Direct Multisearch (DMS) algorithm proposed in [36]. This approach, inspired by the search/poll paradigm of direct-search methods of directional type, does not aggregate any of the objective functions but rather handles a list $L_k$ of non-dominated points (from which the new iterates are chosen). The aim of DMS is to generate as many points in the Pareto front as possible from the polling procedure itself, while keeping the whole framework general enough to accommodate other strategies from the literature. DMS hence generalizes to multiobjective optimization (MOO) all direct-search methods of directional type. We report a simplified version of DMS in Algorithm 8. It uses the strict partial order induced by the cone $\mathbb{R}_+^m$. Following [37], we indicate with $D(L) \subset \mathbb{R}^m$ the image of the set points dominated by a list of points $L$ and $D(L; a)$ the image related to those set of points with an $\ell_\infty$-norm distance to $D(L)$ no larger than $a$. An iteration is successful in Algorithm 8 when there are modifications in the list of nondominated points, that is a new point was accepted, such that $F(\mathbf{x}) \notin D(L, \rho(\alpha))$. It is crucial to emphasize that, when dealing with MOO problems, the cone of descent directions related to all components of the objective $F$ can become arbitrarily narrow. We hence always need to consider density of the search directions at a given limit point [36,37] to guarantee convergence.

**Algorithm 8: Basic Direct Multisearch method.**

1: **Inputs:** Starting point/incumbent solution $\mathbf{x}_0 \in \mathbb{R}^n$, initial stepsize $\alpha_0 > 0$, $L_0 = \{(\mathbf{x}_0, \alpha_0)\}$, stepsize update parameters $0 < \theta < 1 \le \gamma$, forcing function $\rho: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$.

2: **for** $k = 0, 1, 2, \dots$ **do**

3: Select a set $\mathbb{D}_k$ of poll directions and evaluate $F(\mathbf{x})$ at the points $\{\mathbf{x}_k + \alpha_k \mathbf{d} : \mathbf{d} \in \mathbb{D}_k\}$

4: Compute $L_{k+1}$ by removing dominated points from $L_k \cup \{(\mathbf{x}_k + \alpha_k \mathbf{d}, \alpha_k) : \mathbf{d} \in \mathbb{D}_k\}$

5: **if** $L_{k+1} \neq L_k$ **then**

6: Declare the iteration as successful, replace new pairs with $(\mathbf{x}_k + \alpha_k \mathbf{d}, \gamma \alpha_k)$ in $L_{k+1}$.

7: **else**

8: Declare the iteration as unsuccessful, replace $(\mathbf{x}_k, \alpha_k)$ with $(\mathbf{x}_k, \theta \alpha_k)$ in $L_{k+1}$.

9: **end if**

10: **end for**

Convergence results can be established under the following assumptions, that generalize those made in single-objective optimization.

**Assumption 6.1.** For all $i \in [1: l]$, the function $f_i$ is continuously differentiable with Lipschitz continuous gradient with constant $M_i$. Set $M_{\max} = \max_{i \in [1: l]} M_i$.

**Assumption 6.2.** For all $i \in [1: l]$, the function $f_i$ is lower and upper bounded in $\{\mathbf{x} \in \mathbb{R}^n : F(\mathbf{x}) \notin D(\{\mathbf{x}_0\})\}$, with lower bound $f_i^{\min}$ and upper bound $f_i^{\max}$. Let

$$ F^{\min} := \min\{f_1^{\min}, \dots, f_l^{\min}\} \quad \text{and} \quad F^{\max} := \max\{f_1^{\max}, \dots, f_l^{\max}\}. $$

**Assumption 6.3.** The set $\{\mathbf{x} \in \mathbb{R}^n : F(\mathbf{x}) \notin D(\{\mathbf{x}_0\})\}$ is compact.

We further assume, without any loss of generality that all the positive spanning sets considered have normalized directions. Similarly to the analysis carried out in the single-objective case, we now report a result that connects at each iteration $k$ the approximate stationarity measure $\mu_{\mathbb{D}_k}(\mathbf{x}_k)$ and the stepsize $\alpha_k$, whose proof is given in [37]:

**Proposition 6.1.** Let Assumption 6.1 hold. Consider that Algorithm 8 is applied to the minimization of $F$, with $\mathbb{D}_k$ positive spanning set. Suppose that the $k$th iteration of Algorithm 8 is unsuccessful. Then,

$$ \mu_{\mathbb{D}_k}(\mathbf{x}_k) \le \alpha_k \frac{M_{\max}}{2} + \frac{\rho(\alpha_k)}{\alpha_k}. \quad (43) $$

The following assumption, which connects the exact and approximate stationarity measure, is also needed in the analysis.

**Assumption 6.4.** There exists $C_1 > 0$ such that

$$ |\mu_{\mathbb{D}_k}(\mathbf{x}_k) - \mu(\mathbf{x}_k)| \le C_1 \mu_{\mathbb{D}_k}(\mathbf{x}_k), \quad \forall k \ge 0. $$
A direct application of Proposition 6.1 under this assumption yields the following convergence result.

**Theorem 6.1.** Suppose that Algorithm 8 is applied to F under Assumptions 6.1, 6.2 and 6.4. Suppose also that the sequence ${D_k}$ is such that $D_k$ is a positive spanning set for all $k \ge 0$. Let $\mathcal{K}$ be a index set of unsuccessful iterations such that $\lim_{k \in \mathcal{K}} \alpha_k \to 0$. Then,

$$ \lim_{k \in \mathcal{K}} \mu(\mathbf{x}_k) = 0. \qquad (44) $$

As noticed in [37], Assumption 6.3 and the use of sufficient decrease to accept new nondominated points guarantee the existence of a refining subsequence ${\mathbf{x}_k}_\mathcal{K}$ such that ${\alpha_k}_\mathcal{K}$ converges to zero (see Definition 2.11). We would also like to highlight that Assumption 6.4 basically requires nonnegativity of the approximate measure $\mu_{D_k}(\mathbf{x}_k)$ at each $k$. This can only be obtained if a rich enough set $D_k$ of poll directions is available.

We now state the main convergence result for the basic DMS algorithmic scheme we reported above. Similarly to what we have seen in the single-objective case, we consider refining subsequences and assume density in the unit sphere for the related refining directions. The proof is directly derived from [36, Theorem 4.9].

**Theorem 6.2.** Suppose that Algorithm 8 is applied to a Lipschitz continuous function F. Let $\mathcal{K}$ be a sequence of unsuccessful steps such that ${\mathbf{x}_k}_\mathcal{K}$ is a refining subsequence converging to $\hat{\mathbf{x}}$. If the corresponding set of refining directions is dense in the unit sphere, then $\hat{\mathbf{x}}$ is a Pareto-Clarke stationary point. If in addition F is continuously differentiable, then $\hat{\mathbf{x}}$ is a Pareto stationary point.

It is important to note that the density of the refining directions is also needed in the smooth case. As mentioned before, this is obviously due to the fact that the cone of descent directions for the components of the function F can become arbitrarily narrow.

In addition to guaranteeing global convergence, the sufficient decrease condition also provides DMS with complexity bounds that quantify how many iterations or function evaluations are needed in the worst case scenario to reach an approximate optimality criterion, such as $\mu(\mathbf{x}_k) \le e$. This is once again comparable to what we observed in the single-objective case, though we note that the analysis is somewhat more involved. Indeed, it is based on linked sequences, i.e., sequences of pairs $(\mathbf{x}_{j_k}, \alpha_{j_k})$, such that for any $k = 1, 2, \dots$ the pair $(\mathbf{x}_{j_k}, \alpha_{j_k}) \in L_k$ is generated at iteration $k-1$ of Algorithm 8 by the pair $(\mathbf{x}_{j_{k-1}}, \alpha_{j_{k-1}}) \in \tilde{L}_{k-1}$.

**Theorem 6.3 ([37]).** Let the assumptions of Theorem 6.1 hold. Suppose further that $\rho(\alpha) = c \alpha^2$ for some $c > 0$, and that $|D_k| \le m$ for every $k \in \mathbb{N}$. Then, for any $e > 0$, Algorithm 8 reaches an iterate $\mathbf{x}_k$ such that $\mu(\mathbf{x}_k) \le e$ in at most

$$ \mathcal{O}(|L(e)|e^{-2m}) \text{ iterations and } \mathcal{O}(m|L(e)|e^{-2m}) \text{ function evaluations,} \qquad (45) $$

where $L(e)$ is the set of linked sequences between the first unsuccessful iteration and the iteration before the one that satisfies the criticality condition.

The reported bounds do not match the $\mathcal{O}(e^{-2})$ complexity bound that the gradient method for MOO problems has (see [51] for further details). The main reason behind this difference is that an iteration is declared successful in DMS when a sufficient decrease is guaranteed for at least one component of F, while in the gradient method the iterate is updated when a suitable decrease is obtained with respect to all components. The dependence on the term $|L(e)|$ is due instead to the fact that DMS tries to approximate the whole Pareto front using the list $L_k$, while the gradient method in [51] only finds one Pareto stationary point. As a final remark, we would like to highlight that for the case of (strongly) convex optimization problems, i.e., where all components of F are (strongly) convex, a complexity analysis should ideally expand upon the findings presented in [41], which focus on the single-objective scenario. However, it still remains as an open question whether similar complexity bounds to those derived in [41] could be established for the multiobjective case.

The handling of constraints in a MOO problem was already addressed in [10,12,36], where an extreme barrier approach was used for this purpose. An exact penalty approach was analyzed in [77] for handling general constraints in a linesearch-based algorithmic framework. To the best of the authors' knowledge, no existing work provides a detailed account of a Direct-Search algorithm for stochastic MOO problems.

## 7. Conclusion

Since their formal introduction more than 60 years ago, direct-search algorithms have been endowed with a rich convergence theory. In addition to possessing convergence and complexity guarantees in the smooth, unconstrained setting, direct-search methods can also be designed to cope with nonsmoothness, noisy function values, constraints of various types as well as multiple objectives. In this survey, we provided key results as well as pointers to the literature that we hope have convinced the readers of the strong theoretical foundations upon which modern direct-search algorithms are built.

Our presentation has highlighted several lines of research that will in our opinion shape the future of direct-search techniques. Recent advances in analyzing randomized direct-search techniques pave the way for scaling up those algorithms, possibly through combination with other tools from high-dimensional statistics. Modern constrained optimization techniques such as manifold theory provide principled ways of handling common constraints, and may lead to novel analyzes departing from the use of standard nonlinear
programming tools. Tackling noise in function values is likely to keep growing in importance to address problems involving massive amounts of data samples such as hyperparameter tuning, while much remains to be done in obtaining a comprehensive, probabilistic analysis of direct-search algorithms in that setting. Overall, we believe that these three research directions are likely to form the basis of next-generation direct-search methods.

## CRediT authorship contribution statement

K.J. Dzahini: Writing – review & editing, Writing – original draft, Supervision, Conceptualization. F. Rinaldi: Writing – review & editing, Writing – original draft, Supervision, Conceptualization. C.W. Royer: Writing – review & editing, Writing – original draft, Supervision, Conceptualization. D. Zeffiro: Writing – original draft, Supervision, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## References

[1] M.A. Abramson, O.A. Brezhneva, J.E. Dennis Jr., R.L. Pingel, Pattern search in the presence of degenerate linear constraints, Optim. Methods Softw. 23 (3) (2008) 297-319, https://doi.org/10.1080/10556780701521670.

[2] J. Achddou, O. Cappé, A. Garivier, Regret analysis of the stochastic direct search method for blind resource allocation, https://arxiv.org/abs/2210.05222, 2022.

[3] S. Alarie, C. Audet, P.-Y. Bouchet, S. Le Digabel, Optimisation of stochastic blackboxes with adaptive precision, SIAM J. Optim. 31 (4) (2021) 3127-3156, https://doi.org/10.1137/20M1318894.

[4] C. Audet, A survey on direct search methods for blackbox optimization and their applications, in: P. Pardalos, T. Rassias (Eds.), Mathematics Without Boundaries: Surveys in Interdisciplinary Research, Springer, New York, NY, 2014, pp. 31-56, chapter 2.

[5] C. Audet, J.E. Dennis Jr., Analysis of generalized pattern searches, SIAM J. Optim. 13 (3) (2002) 889-903, https://doi.org/10.1137/S1052623400378742.

[6] C. Audet, J.E. Dennis Jr., A pattern search filter method for nonlinear programming without derivatives, SIAM J. Optim. 14 (4) (2004) 980-1010, https://doi.org/10.1137/S105262340138983X.

[7] C. Audet, J.E. Dennis Jr., Mesh adaptive direct search algorithms for constrained optimization, SIAM J. Optim. 17 (1) (2006) 188-217, https://doi.org/10.1137/040603371.

[8] C. Audet, J.E. Dennis Jr., A progressive barrier for derivative-free nonlinear programming, SIAM J. Optim. 20 (1) (2009) 445-472, https://doi.org/10.1137/070692662.

[9] C. Audet, W. Hare, Derivative-Free and Blackbox Optimization, Springer Series in Operations Research and Financial Engineering., Springer, Cham, Switzerland, 2017.

[10] C. Audet, G. Savard, W. Zghal, Multiobjective optimization through a series of single-objective formulations, SIAM J. Optim. 19 (1) (2008) 188-210, https://doi.org/10.1137/060677513.

[11] C. Audet, J.E. Dennis Jr., S. Le Digabel, Globalization strategies for mesh adaptive direct search, Comput. Optim. Appl. 46 (2) (2010) 193-215, https://doi.org/10.1007/s10589-009-9266-1.

[12] C. Audet, G. Savard, W. Zghal, A mesh adaptive direct search algorithm for multiobjective optimization, Eur. J. Oper. Res. 204 (3) (2010) 545-556, https://doi.org/10.1016/j.ejor.2009.11.010.

[13] C. Audet, K.J. Dzahini, M. Kokkolaras, S. Le Digabel, Stochastic mesh adaptive direct search for blackbox optimization using probabilistic estimates, Comput. Optim. Appl. 79 (1) (2021) 1-34, https://doi.org/10.1007/s10589-020-00249-0.

[14] C. Audet, A. Batailly, S. Kojtych, Escaping unknown discontinuous regions in blackbox optimization, SIAM J. Optim. 32 (3) (2022) 1843-1870, https://doi.org/10.1137/21M1420915.

[15] C. Audet, S. Le Digabel, V. Rochon Montplaisir, C. Tribes, Algorithm 1027: NOMAD version 4: nonlinear optimization with the MADS algorithm, ACM Trans. Math. Softw. 48 (3) (2022) 35, https://doi.org/10.1145/3544489.

[16] T. Bäck, Evolutionary Algorithms in Theory and Practice: Evolution Strategies, Evolutionary Programming, Genetic Algorithms, Oxford University Press, New York, NY, USA, 1996.

[17] A.S. Berahas, R.H. Byrd, J. Nocedal, Derivative-free optimization of noisy functions via quasi-Newton methods, SIAM J. Optim. 29 (2) (2019) 965-993, https://doi.org/10.1137/18M1177718.

[18] A.S. Berahas, L. Cao, K. Scheinberg, Global convergence rate analysis of a generic line search algorithm with noise, SIAM J. Optim. 31 (2) (2021) 1489-1518, https://doi.org/10.1137/19M1291832.

[19] E.H. Bergou, E. Gorbunov, P. Richtářík, Stochastic three points method for unconstrained smooth minimization, SIAM J. Optim. 30 (4) (2020) 2726-2749, https://doi.org/10.1137/19M1244378.

[20] J. Blanchet, C. Cartis, M. Menickelly, K. Scheinberg, Convergence rate analysis of a stochastic trust-region method via submartingales, INFORMS J. Optim. 1 (2) (2019) 92-119, https://doi.org/10.1287/ijoo.2019.0016.

[21] R. Bollapragada, S.M. Wild, Adaptive sampling quasi-Newton methods for zeroth-order stochastic optimization, Math. Program. Comput. 15 (2) (2023) 327-364, https://doi.org/10.1007/s12532-023-00233-9.

[22] I.M. Bomze, F. Rinaldi, D. Zeffiro, Frank-Wolfe and friends: a journey into projection-free first-order optimization methods, 4OR 19 (3) (2021) 313-345, https://doi.org/10.1007/s10288-021-00493-y.

[23] A.J. Booker, J.E. Dennis Jr., P.D. Frank, D.B. Serafini, V. Torczon, M.W. Trosset, A rigorous framework for optimization of expensive functions by surrogates, Struct. Multidiscip. Optim. 17 (1) (1999) 1-13, https://doi.org/10.1007/BF01197708.

[24] P.-Y. Bouchet, Théorie de l'optimisation sans dérivées dans le cas discontinu, PhD thesis, Polytechnique Montréal, December 2023.

[25] N. Boumal, An Introduction to Optimization on Smooth Manifolds, Cambridge University Press, Cambridge, United Kingdom, 2023.

[26] A. Brilli, M. Kimiaiei, G. Liuzzi, S. Lucidi, Worst case complexity bounds for linesearch-type derivative-free algorithms, https://arxiv.org/abs/2302.05274, 2023.

[27] L.F. Bueno, A. Friedlander, J.M. Martinez, F.N.C. Sobral, Inexact restoration method for derivative-free optimization with smooth constraints, SIAM J. Optim. 23 (2) (2013) 1189-1213, https://doi.org/10.1137/110856253.

[28] C. Cartis, N.I.M. Gould, P.L. Toint, Evaluation Complexity of Algorithms for Nonconvex Optimization: Theory, Computation and Perspectives, MOS-SIAM Series on Optimization, vol. MO30, SIAM, 2022.

[29] K.H. Chang, Stochastic Nelder-Mead simplex method – a new globally convergent direct search method for simulation optimization, Eur. J. Oper. Res. 220 (3) (2012) 684-694, https://doi.org/10.1016/j.ejor.2012.02.028.
[30] R. Chen, M. Menickelly, K. Scheinberg, Stochastic optimization using a trust-region method and random models, Math. Program. 169 (2) (2018) 447–487,
https://doi.org/10.1007/s10107-017-1141-8.

[31] F.H. Clarke, Optimization and Nonsmooth Analysis, Classics in Applied Mathematics, vol. 5, John Wiley and Sons, New York, ISBN 9780898712568, 1983,
Reissued in 1990 by SIAM Publications, Philadelphia.

[32] A.R. Conn, K. Scheinberg, L.N. Vicente, Introduction to Derivative-Free Optimization, SIAM, Philadelphia, ISBN 978-0-898716-68-9, 2009.

[33] I.D. Coope, C.J. Price, On the convergence of grid-based methods for unconstrained optimization, SIAM J. Optim. 11 (4) (2001) 859–869, https://doi.org/10.1137/S1052623499354989.

[34] A. Cristofari, F. Rinaldi, A derivative-free method for structured optimization problems, SIAM J. Optim. 31 (2) (2021) 1079–1107, https://doi.org/10.1137/20M1337417.

[35] A.L. Custódio, L.N. Vicente, Analysis of direct searches for discontinuous functions, Math. Program. 133 (2012) 299–325, https://doi.org/10.1007/s10107-010-0429-8.

[36] A.L. Custódio, J.F. Madeira, A.I.F. Vaz, L.N. Vicente, Direct multiresearch for multiobjective optimization, SIAM J. Optim. 21 (3) (2011) 1109–1140, https://doi.org/10.1137/10079731X.

[37] A.L. Custódio, Y. Diouane, R. Garmanjani, E. Riccietti, Worst-case complexity bounds of directional direct-search methods for multiobjective optimization, J. Optim. Theory Appl. 188 (2021) 73–93, https://doi.org/10.1007/s10957-020-01781-z.

[38] J.E. Dennis Jr., C.J. Price, I.D. Coope, Direct search methods for nonlinearly constrained optimization using filters and frames, Optim. Eng. 5 (2) (2004) 123–144,
https://doi.org/10.1023/B:OPTE0000033371.04406.e0.

[39] M.A. Diniz-Ehrhardt, J.M. Martínez, M. Raydán, A derivative-free nonmonotone line-search technique for unconstrained optimization, J. Comput. Appl. Math. 219 (2) (2008) 383–397, https://doi.org/10.1016/j.cam.2007.07.017.

[40] Y. Diouane, V. Kungurtsev, F. Rinaldi, D. Zeffiro, Inexact direct-search methods for bilevel optimization problems, https://arxiv.org/abs/2307.09924, 2023.

[41] M. Dodangeh, L.N. Vicente, Worst case complexity of direct search under convexity, Math. Program. 155 (1–2) (2016) 307–332, https://doi.org/10.1007/s10107-014-0847-0.

[42] M. Dodangeh, L.N. Vicente, Z. Zhang, On the optimal order of worst case complexity of direct search, Optim. Lett. 10 (2016) 699–708, https://doi.org/10.1007/s11590-015-0908-1.

[43] E.D. Dolan, R.M. Lewis, V. Torczon, On the local convergence of pattern search, SIAM J. Optim. 14 (2) (2003) 567–583, https://doi.org/10.1137/S1052623400374495.

[44] D.W. Dreisigmeyer, Equality constraints, Riemannian manifolds and direct search methods, Technical Report LA-UR-06-7406, Los Alamos National Laboratory,
2006.

[45] D.W. Dreisigmeyer, Direct search methods on reductive homogeneous spaces, J. Optim. Theory Appl. 176 (3) (2018) 585–604, https://doi.org/10.1007/s10957-018-1225-5.

[46] K.J. Dzahini, Expected complexity analysis of stochastic direct-search, Comput. Optim. Appl. 81 (1) (2022) 179–200, https://doi.org/10.1007/s10589-021-00329-9.

[47] K.J. Dzahini, S.M. Wild, Direct search for stochastic optimization in random subspaces with zeroth-, first-, and second-order convergence and expected complexity, arXiv, arXiv:2403.13320, 2024, https://arxiv.org/abs/2403.13320.

[48] K.J. Dzahini, S.M. Wild, A class of sparse Johnson–Lindenstrauss transforms and analysis of their extreme singular values, SIAM J. Matrix Anal. Appl. 46 (1) (2025) 416–438, https://doi.org/10.1137/23M1605661.

[49] K.J. Dzahini, M. Kokkolaras, S. Le Digabel, Constrained stochastic blackbox optimization using a progressive barrier and probabilistic estimates, Math. Program. 198 (1) (2023) 675–732, https://doi.org/10.1007/s10107-022-01787-7.

[50] G. Fasano, G. Liuzzi, S. Lucidi, F. Rinaldi, A linesearch-based derivative-free approach for nonsmooth constrained optimization, SIAM J. Optim. 24 (3) (2014) 959–992, https://doi.org/10.1137/130940037.

[51] J. Fliege, A.I.F. Vaz, L.N. Vicente, Complexity of gradient descent for multiobjective optimization, Optim. Methods Softw. 34 (5) (2019) 949–959, https://doi.org/10.1080/10556788.2018.1510928.

[52] U.M. García-Palomares, J.F. Rodríguez, New sequential and parallel derivative-free algorithms for unconstrained minimization, SIAM J. Optim. 13 (1) (2002) 79–96, https://doi.org/10.1137/S1052623400370606.

[53] R. Garmanjani, L.N. Vicente, Smoothing and worst-case complexity for direct-search methods in nonsmooth optimization, IMA J. Numer. Anal. 33 (3) (2013) 1008–1028, https://doi.org/10.1093/imanum/drs027.

[54] M.G. Gasparo, A. Papini, A. Pasquali, Nonmonotone algorithms for pattern search methods, Numer. Algorithms 28 (2001) 171–186, https://doi.org/10.1023/A:1014046817188.

[55] A.A. Goldstein, Optimization for Lipschitz continuous functions, Math. Program. 13 (1) (1977) 14–22, https://doi.org/10.1007/BF01584320.

[56] S. Gratton, L.N. Vicente, A merit function approach for direct search, SIAM J. Optim. 24 (4) (2014) 1980–1998, https://doi.org/10.1137/130917661.

[57] S. Gratton, C.W. Royer, L.N. Vicente, Z. Zhang, Direct search based on probabilistic descent, SIAM J. Optim. 25 (3) (2015) 1515–1541, https://doi.org/10.1137/140961602.

[58] S. Gratton, C.W. Royer, L.N. Vicente, Z. Zhang, Direct search based on probabilistic feasible descent for bound and linearly constrained problems, Comput. Optim. Appl. 72 (3) (2019) 525–559, https://doi.org/10.1007/s10589-019-00062-4.

[59] L. Grippo, F. Rinaldi, A class of derivative-free nonmonotone optimization algorithms employing coordinate rotations and gradient approximations, Comput. Optim. Appl. 60 (2015) 1–33, https://doi.org/10.1007/s10589-014-9665-9.

[60] L. Grippo, F. Lampariello, S. Lucidi, Global convergence and stabilization of unconstrained minimization methods, J. Optim. Theory Appl. 56 (1988) 385–406,
https://doi.org/10.1007/BF00939550.

[61] R. Hooke, T.A. Jeeves, "Direct search" solution of numerical and statistical problems, J. Assoc. Comput. Mach. 8 (2) (1961) 212–229, https://doi.org/10.1145/321062.321069.

[62] B. Jin, K. Scheinberg, M. Xie, High probability complexity bounds for line search based on stochastic oracles, NeurIPS 34 (2021) 9193–9203.

[63] B. Jin, K. Scheinberg, M. Xie, Sample complexity analysis for adaptive optimization algorithms with stochastic oracles, https://arxiv.org/abs/2303.06838, 2023.

[64] C.T. Kelley, Implicit Filtering. Software Environment and Tools, Society for Industrial and Applied Mathematics, Philadelphia, 2011.

[65] T.G. Kolda, R.M. Lewis, V. Torczon, Optimization by direct search: new perspectives on some classical and modern methods, SIAM Rev. 45 (3) (2003) 385–482,
https://doi.org/10.1137/S003614450242889.

[66] T.G. Kolda, R.M. Lewis, V. Torczon, Stationarity results for generating set search for linearly constrained optimization, SIAM J. Optim. 17 (4) (2006) 943–968,
https://doi.org/10.1137/S1052623403433638.

[67] T.G. Kolda, R.M. Lewis, V. Torczon, A generating set direct search augmented Lagrangian algorithm for optimization with a combination of general and linear constraints, Technical report, Technical Report SAND2006-5315, Sandia National Laboratories, 2006.

[68] J. Konečný, P. Richtářík, Simple complexity analysis of simplified direct search, https://arxiv.org/abs/1410.0390, 2014.

[69] V. Kungurtsev, F. Rinaldi, D. Zeffiro, Retraction based direct search methods for derivative free Riemannian optimization, https://arxiv.org/abs/2202.11052, 2022.

[70] J. Larson, M. Menickelly, S.M. Wild, Derivative-free optimization methods, Acta Numer. 28 (2019) 287–404, https://doi.org/10.1017/s0962492919000060.

[71] S. Le Digabel, S.M. Wild, A taxonomy of constraints in black-box simulation-based optimization, Optim. Eng. (2023) 1–19, https://doi.org/10.1007/s11081-023-09839-3.
[72] R.M. Lewis, V. Torczon, A globally convergent augmented Lagrangian pattern search algorithm for optimization with general constraints and simple bounds, SIAM J. Optim. 12 (4) (2002) 1075-1089, https://doi.org/10.1137/S105262349833972.

[73] R.M. Lewis, V. Torczon, M.W. Trosset, Direct search methods: then and now, J. Comput. Appl. Math. 124 (1-2) (2000) 191-207, https://doi.org/10.1016/S0377-0427(00)00423-4.

[74] T. Lin, Z. Zheng, M. Jordan, Gradient-free methods for deterministic and stochastic nonsmooth nonconvex optimization, Adv. Neural Inf. Process. Syst. 35 (2022) 26160-26175.

[75] G. Liuzzi, S. Lucidi, A derivative-free algorithm for inequality constrained nonlinear programming via smoothing of an $l_{\infty}$ penalty function, SIAM J. Optim. 20 (1) (2009) 1-29, https://doi.org/10.1137/070711451.

[76] G. Liuzzi, S. Lucidi, M. Sciandrone, Sequential penalty derivative-free methods for nonlinear constrained optimization, SIAM J. Optim. 20 (5) (2010) 2614-2635, https://doi.org/10.1137/090750639.

[77] G. Liuzzi, S. Lucidi, F. Rinaldi, A derivative-free approach to constrained multiobjective nonsmooth optimization, SIAM J. Optim. 26 (4) (2016) 2744-2774, https://doi.org/10.1137/15M1037810.

[78] S. Lucidi, M. Sciandrone, A derivative-free algorithm for bound constrained optimization, Comput. Optim. Appl. 21 (2) (2002) 119-142, https://doi.org/10.1023/A:1013735414984.

[79] S. Lucidi, M. Sciandrone, On the global convergence of derivative-free methods for unconstrained optimization, SIAM J. Optim. 13 (1) (2002) 97-116, https://doi.org/10.1137/S1052623497330392.

[80] S. Lucidi, M. Sciandrone, P. Tseng, Objective-derivative-free methods for constrained optimization, Math. Program. 92 (2002) 37-59, https://doi.org/10.1007/s101070100266.

[81] J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (4) (1965) 308-313, https://doi.org/10.1093/comjnl/7.4.308.

[82] Y. Nesterov, V. Spokoiny, Random gradient-free minimization of convex functions, Found. Comput. Math. 17 (2017) 527-566, https://doi.org/10.1007/s10208-015-9296-2.

[83] C. Paquette, K. Scheinberg, A stochastic line search method with expected complexity analysis, SIAM J. Optim. 30 (1) (2020) 349-376, https://doi.org/10.1137/18M1216250.

[84] B. Pauwels, Optimisation sans dérivées sous incertitudes appliquées à des simulateurs coûteux, PhD thesis, Université de Toulouse, March 2016, http://thesepsups.ups-tlse.fr/3103/1/2016TOU30035.pdf.

[85] C.J. Price, I.D. Coope, Frames and grids in unconstrained and linearly constrained optimization: a nonsmooth approach, SIAM J. Optim. 14 (2) (2003) 415-438, https://doi.org/10.1137/S1052623402407084.

[86] F. Rinaldi, L.N. Vicente, D. Zeffiro, Stochastic trust-region and direct-search methods: a weak tail bound condition and reduced sample sizing, https://arxiv.org/abs/2202.11074, 2023.

[87] L. Roberts, C.W. Royer, Direct search based on probabilistic descent in reduced spaces, SIAM J. Optim. 33 (4) (2023) 3057-3082, https://doi.org/10.1137/22M1488569.

[88] V. Torczon, On the convergence of pattern search algorithms, SIAM J. Optim. 7 (1) (1997) 1-25, https://doi.org/10.1137/S1052623493250780.

[89] L.N. Vicente, Worst case complexity of direct search, EURO J. Comput. Optim. 1 (1) (2013) 143-153, https://doi.org/10.1007/s13675-012-0003-7.

[90] M.H. Wright, Direct Search Methods: Once Scorned, Now Respectable, Pitman Research Notes in Mathematics Series, vol. 344, 1996, pp. 191-208, https://api.semanticscholar.org/CorpusID:64988980.