2.1.1 Simplex methods

Simplex methods (not to be confused with Dantzig’s simplex method for linear programming) move and manipulate a collection of $n+1$ affinely independent points (i.e. the vertices of a simplex in $\mathbb{R}^{n}$ ) when solving (DET). The method of Spendley, Hext and Himsworth (Reference Spendley, Hext and Himsworth1962) involves either taking the point in the simplex with the largest function value and reflecting it through the hyperplane defined by the remaining $n$ points or moving the $n$ worst points toward the best vertex of the simplex. In this manner, the geometry of all simplices remains the same as that of the starting simplex. (That is, all simplices are similar in the geometric sense.)

Nelder and Mead (Reference Nelder and Mead1965) extend the possible simplex operations, as shown in Figure 2.1, by allowing the ‘expansion’ and ‘contraction’ operations in addition to the ‘reflection’ and ‘shrink’ operations of Spendley et al. (Reference Spendley, Hext and Himsworth1962). These operations enable the Nelder–Mead simplex method to distort the simplex in order to account for possible curvature present in the objective function.

Nelder and Mead (Reference Nelder and Mead1965) propose stopping further function evaluations when the standard error of the function values at the simplex vertices is small. Others, Woods (Reference Woods1985) for example, propose stopping when the size of the simplex’s longest side incident to the best simplex vertex is small.

Figure 2.1.

Primary Nelder–Mead simplex operations: original simplex, reflection, expansion, inner contraction, and shrink.

Nelder–Mead is an incredibly popular method, in no small part due to its inclusion in Numerical Recipes (Press, Teukolsky, Vetterling and Flannery Reference Press, Teukolsky, Vetterling and Flannery2007), which has been cited over $125\,000$ times and no doubt used many times more. The method (as implemented by Lagarias, Poonen and Wright Reference Lagarias, Poonen and Wright2012) is also the algorithm underlying fminsearch in MATLAB. Benchmarking studies highlight Nelder–Mead performance in practice (Moré and Wild Reference Moré and Wild2009, Rios and Sahinidis Reference Rios and Sahinidis2013).

The method’s popularity from its inception was not diminished by the lack of theoretical results proving its ability to identify stationary points. Woods (Reference Woods1985) presents a non-convex, two-dimensional function where Nelder–Mead converges to a non-stationary point (where the function’s Hessian is singular). Furthermore, McKinnon (Reference McKinnon1998) presents a class of thrice-continuously differentiable, strictly convex functions on $\mathbb{R}^{2}$ where the Nelder–Mead simplex fails to converge to the lone stationary point. The only operation that Nelder–Mead performs on this relatively routine function is repeated ‘inner contraction’ of the initial simplex.

Researchers have continued to develop convergence results for modified or limited versions of Nelder–Mead. Kelley (Reference Kelley1999a ) addresses Nelder–Mead’s theoretical deficiencies by restarting the method when the objective decrease on consecutive iterations is not larger than a multiple of the simplex gradient norm. Such restarts do not ensure that Nelder–Mead will converge: Kelley (Reference Kelley1999a ) shows an example of such behaviour. Price, Coope and Byatt (Reference Price, Coope and Byatt2002) embed Nelder–Mead in a different (convergent) algorithm using positive spanning sets. Nazareth and Tseng (Reference Nazareth and Tseng2002) propose a clever, though perhaps superfluous, variant that connects Nelder–Mead to golden-section search.

Lagarias, Reeds, Wright and Wright (Reference Lagarias, Reeds, Wright and Wright1998) show that Nelder–Mead (with appropriately chosen reflection and expansion coefficients) converges to the global minimizer of strictly convex functions when $n=1$ . Gao and Han (Reference Gao and Han2012) show that the contraction and expansion steps of Nelder–Mead satisfy a descent condition on uniformly convex functions. Lagarias, Poonen and Wright (Reference Lagarias, Poonen and Wright2012) show that a restricted version of the Nelder–Mead method – one that does not allow an expansion step – can converge to minimizers of any twice-continuously differentiable function with a positive-definite Hessian and bounded level sets. (Note that the class of functions from McKinnon (Reference McKinnon1998) have singular Hessians at only one point – their minimizers – and not at the point to which the simplex vertices are converging.)

The simplex method of Rykov (Reference Rykov1980) includes ideas from model-based methods. Rykov varies the number of reflected vertices from iteration to iteration, following one of three rules that depend on the function value at the simplex centroid $\boldsymbol{x}_{c}$ . Rykov considers both evaluating $f$ at the centroid and approximating $f$ at the centroid using the values of $f$ at the vertex. The non-reflected vertices are also moved in parallel with the reflected subset of vertices. In general, the number of reflected vertices is chosen so that $\boldsymbol{x}_{c}$ moves in a direction closest to $-\unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{c})$ . This, along with a test of sufficient decrease in $f$ , ensures convergence of the modified simplex method to a minimizer of convex, continuously differentiable functions with bounded level sets and Lipschitz-bounded gradients. (The sufficient-decrease condition is also shown to be efficient for the classical Nelder–Mead algorithm.)

Tseng (Reference Tseng1999) proposes a modified simplex method that keeps the $b_{k}$ best simplex vertices on a given iteration $k$ and uses them to reflect the remaining vertices. Their method prescribes that ‘the rays emanating from the reflected vertices toward the $b_{k}$ best vertices should contain, in their convex hull, the rays emanating from a weighted centroid of the $b_{k}$ best vertices toward the to-be-reflected vertices’. Their method also includes a fortified descent condition that is stronger than common sufficient-decrease conditions. If $f$ is continuously differentiable and bounded below and $b_{k}$ is fixed for all iterations, Tseng (Reference Tseng1999) prove that every cluster point of the sequence of candidate points generated by their method is a stationary point.

Bűrmen, Puhan and Tuma (Reference Bűrmen, Puhan and Tuma2006) propose a convergent version of a simplex method that does not require a sufficient descent condition to be satisfied. Instead, they ensure that evaluated points lie on a grid of points, and they show that this grid will be refined as the method proceeds.

2.1.2 Directional direct-search methods

Broadly speaking, each iteration of a directional direct-search (DDS) method generates a finite set of points near the current point $\boldsymbol{x}_{k}$ ; these poll points are generated by taking $\boldsymbol{x}_{k}$ and adding terms of the form $\unicode[STIX]{x1D6FC}_{k}\boldsymbol{d}$ , where $\unicode[STIX]{x1D6FC}_{k}$ is a positive step size and $\boldsymbol{d}$ is an element from a finite set of directions $\boldsymbol{D}_{k}$ . Kolda et al. (Reference Kolda, Lewis and Torczon2003) propose the term generating set search methods to encapsulate this class of methods.Footnote ³ The objective function $f$ is then evaluated at all or some of the poll points, and $\boldsymbol{x}_{k+1}$ is selected to be some poll point that produces a (sufficient) decrease in the objective and the step size is possibly increased. If no poll point provides a sufficient decrease, $\boldsymbol{x}_{k+1}$ is set to $\boldsymbol{x}_{k}$ and the step size is decreased. In either case, the set of directions $\boldsymbol{D}_{k}$ can (but need not) be modified to obtain $\boldsymbol{D}_{k+1}$ .

A general DDS method is provided in Algorithm 2, which includes a search step where $f$ is evaluated at any finite set of points $\boldsymbol{Y}_{k}$ , including $\boldsymbol{Y}_{k}=\emptyset$ . The search step allows one to (potentially) improve the performance of Algorithm 2. For example, points could be randomly sampled during the search step from the domain in the hope of finding a better local minimum, or a person running the algorithm may have problem-specific knowledge that can generate candidate points given the observed history of evaluated points and their function values. While the search step allows for this insertion of such heuristics, rigorous convergence results are driven by the more disciplined poll step. When testing for objective decrease in Algorithm 1, one can stop evaluating points in $\boldsymbol{P}$ (line 6) as soon as the first point is identified where there is (sufficient) decrease in $f$ . In this case, the polling (or search) step is considered opportunistic.

DDS methods are largely distinguished by how they generate the set of poll directions $\boldsymbol{D}_{k}$ at line 7of Algorithm 2. Perhaps the first approach is coordinate search, in which the poll directions are defined as $\boldsymbol{D}_{k}=\{\pm \boldsymbol{e}_{i}:i=1,2,\ldots ,n\}$ , where $\boldsymbol{e}_{i}$ denotes the $i$ th elementary basis vector (i.e. column $i$ of the identity matrix in $n$ dimensions). The first known description of coordinate search appears in the work of Fermi and Metropolis (Reference Fermi and Metropolis1952) where the smallest positive integer $l$ is sought such that $f(\boldsymbol{x}_{k}+l\unicode[STIX]{x1D6FC}\boldsymbol{e}_{1}/2)>f(\boldsymbol{x}_{k}+(l-1)\unicode[STIX]{x1D6FC}\boldsymbol{e}_{1}/2)$ . If an increase in $f$ is observed at $\boldsymbol{e}_{1}/2$ then $-\boldsymbol{e}_{1}/2$ is considered. After such an integer $l$ is identified for the first coordinate direction, $\boldsymbol{x}_{k}$ is updated to $\boldsymbol{x}_{k}\pm l\boldsymbol{e}_{1}/2$ and the second coordinate direction is considered. If $\boldsymbol{x}_{k}$ is unchanged after cycling through all coordinate directions, then the method is repeated but with $\pm \boldsymbol{e}_{i}/2$ replaced with $\pm \boldsymbol{e}_{i}/16$ , terminating when no improvement is observed for this smaller $\unicode[STIX]{x1D6FC}$ . In terms of Algorithm 2the search set $\boldsymbol{Y}_{k}=\emptyset$ at line 4, and the descent test at line 4of Algorithm 1merely tests for simple decrease, that is, $f(\boldsymbol{p}_{i})-f(\boldsymbol{x})<0$ . Other versions of acceptability in line 4of Algorithm 1are employed by methods discussed later.

Proofs that DDS methods converge first appeared in the works of Céa (Reference Céa1971) and Yu (Reference Yu1979), although both require the sequence of step-size parameters to be non-increasing. Lewis, Torczon and Trosset (Reference Lewis, Torczon and Trosset2000) attribute the first global convergence proof for coordinate search to Polak (Reference Polak1971, p. 43). In turn, Polak cites the ‘method of local variation’ of Banichuk, Petrov and Chernous’ko (Reference Banichuk, Petrov and Chernous’ko1966); although Banichuk et al. (Reference Banichuk, Petrov and Chernous’ko1966) do develop parts of a convergence proof, they state in Remark 1 that ‘the question of the strict formulation of the general sufficient conditions for convergence of the algorithm to a minimum remains open’.

Typical convergence results for DDS require that the set $\boldsymbol{D}_{k}$ is a positive spanning set (PSS) for the domain $\boldsymbol{\unicode[STIX]{x1D6FA}}$ ; that is, any point $\boldsymbol{x}\in \boldsymbol{\unicode[STIX]{x1D6FA}}$ can be written as

$$\begin{eqnarray}\displaystyle \boldsymbol{x}=\mathop{\sum }_{i=1}^{|\boldsymbol{D}_{k}|}\unicode[STIX]{x1D706}_{i}\boldsymbol{d}_{i}, & & \displaystyle \nonumber\end{eqnarray}$$

where $\boldsymbol{d}_{i}\in \boldsymbol{D}_{k}$ and $\unicode[STIX]{x1D706}_{i}\geq 0$ for all $i$ . Some of the first discussions of properties of positive spanning sets were presented by Davis (Reference Davis1954) and McKinney (Reference McKinney1962), but recent treatments have also appeared in Regis (Reference Regis2016). In addition to requiring positive spanning sets during the poll step, earlier DDS convergence results depended on $f$ being continuously differentiable. When $f$ is non-smooth, no descent direction is guaranteed for these early DDS methods, even when the step size is arbitrarily small. See, for example, the modification of the Dennis–Woods (Dennis, Jr and Woods Reference Dennis, Woods and Wouk1987) function by Kolda et al. (Reference Kolda, Lewis and Torczon2003, Figure 6.2) and a discussion of why coordinate-search methods (for example) will not move when started at a point of non-differentiability; moreover, when started at differentiable points, coordinate-search methods tend to converge to a point that is not (Clarke) stationary.

The pattern-search method of Torczon (Reference Torczon1991) revived interest in direct-search methods. The method therein contains ideas from both DDS and simplex methods. Given a simplex defined by $\boldsymbol{x}_{k},\boldsymbol{y}_{1},\ldots ,\boldsymbol{y}_{n}$ (where $\boldsymbol{x}_{k}$ is the simplex vertex with smallest function value), the polling directions are given by $\boldsymbol{D}_{k}=\{\boldsymbol{y}_{i}-\boldsymbol{x}_{k}:i=1,\ldots ,n\}$ . If a decrease is observed at the best poll point in $\boldsymbol{x}_{k}+\boldsymbol{D}_{k}$ , the simplex is set to either $\boldsymbol{x}_{k}\bigcup \boldsymbol{x}_{k}+\boldsymbol{D}_{k}$ or some expansion thereof. If no improvement is found during the poll step, the simplex is contracted. Torczon (Reference Torczon1991) shows that if $f$ is continuous on the level set of $\boldsymbol{x}_{0}$ and this level set is compact, then a subsequence of $\{\boldsymbol{x}_{k}\}$ converges to a stationary point of $f$ , a point where $f$ is non-differentiable, or a point where $f$ is not continuously differentiable.

A generalization of pattern-search methods is the class of generalized pattern-search (GPS) methods. Early GPS methods did not allow for a search step; the search-poll paradigm was introduced by Audet (Reference Audet2004). GPS methods are characterized by fixing a positive spanning set $\boldsymbol{D}$ and selecting $\boldsymbol{D}_{k}\subseteq \boldsymbol{D}$ during the poll step at line 7on each iteration of Algorithm 2. Torczon (Reference Torczon1997) assumes that the test for decrease in line 4in Algorithm 1is simple decrease, that is, that $f(\boldsymbol{p}_{i})<f(\boldsymbol{x})$ . Early analysis of GPS methods using simple decrease required the step size $\unicode[STIX]{x1D6FC}_{k}$ to remain rational (Audet and Dennis, Jr Reference Audet and Dennis2002, Torczon Reference Torczon1997). Audet (Reference Audet2004) shows that such an assumption is necessary by constructing small-dimensional examples where GPS methods do not converge if $\unicode[STIX]{x1D6FC}_{k}$ is irrational. Works below show that if a sufficient (instead of simple) decrease is ensured, $\unicode[STIX]{x1D6FC}_{k}$ can take irrational values.

A refinement of the analysis of GPS methods was made by Dolan, Lewis and Torczon (Reference Dolan, Lewis and Torczon2003), which shows that when $\unicode[STIX]{x1D6FB}f$ is Lipschitz-continuous, the step-size parameter $\unicode[STIX]{x1D6FC}_{k}$ scales linearly with $\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{k})\Vert$ . Therefore $\unicode[STIX]{x1D6FC}_{k}$ can be considered a reliable measure of first-order stationarity and justifies the traditional approach of stopping a GPS method when $\unicode[STIX]{x1D6FC}_{k}$ is small. Second-order convergence analyses of GPS methods have also been considered. Abramson (Reference Abramson2005) shows that, when applied to a twice-continuously differentiable $f$ , a GPS method that infinitely often has $\boldsymbol{D}_{k}$ include a fixed orthonormal basis and its negative will have a limit point satisfying a ‘pseudo-second-order’ stationarity condition. Building off the use of curvature information in Frimannslund and Steihaug (Reference Frimannslund and Steihaug2007), Abramson, Frimannslund and Steihaug (Reference Abramson, Frimannslund and Steihaug2013) show that a modification of the GPS framework that constructs approximate Hessians of $f$ will converge to points that are second-order stationary provided that certain conditions on the Hessian approximation hold (and a fixed orthonormal basis and its negative are in $\boldsymbol{D}_{k}$ infinitely often).

In general, first-order convergence results (there exists a limit point $\boldsymbol{x}_{\ast }$ of $\{\boldsymbol{x}_{k}\}$ generated by a GPS method such that $\unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{\ast })=\mathbf{0}$ ) for GPS methods can be demonstrated when $f$ is continuously differentiable. For general Lipschitz-continuous (but non-smooth) functions $f$ , however, one can only demonstrate that on a particular subsequence ${\mathcal{K}}$ , satisfying $\{{\boldsymbol{x}_{k}\}}_{k\in {\mathcal{K}}}\rightarrow \boldsymbol{x}_{\ast }$ , for each $\boldsymbol{d}$ that appears infinitely many times in $\{{\boldsymbol{D}_{k}\}}_{k\in {\mathcal{K}}}$ , it holds that $f^{\prime }(\boldsymbol{x}_{\ast };\boldsymbol{d})\geq 0$ ; that is, the directional derivative at $\boldsymbol{x}_{\ast }$ in the direction $\boldsymbol{d}$ is non-negative.

The flexibility of GPS methods inspired various extensions. Abramson, Audet and Dennis, Jr (Reference Abramson, Audet and Dennis2004) consider adapting GPS to utilize derivative information when it is available in order to reduce the number of points evaluated during the poll step. Abramson, Audet, Dennis, Jr and Le Digabel (Reference Abramson, Audet, Dennis and Le Digabel2009b ) and Frimannslund and Steihaug (Reference Frimannslund and Steihaug2011) re-use previous function evaluations in order to determine the next set of directions. Custódio and Vicente (Reference Custódio and Vicente2007) consider re-using previous function evaluations to compute simplex gradients; they show that the information obtained from simplex gradients can be used to reorder the poll points $\boldsymbol{P}$ in Algorithm 1. A similar use of simplex gradients in the non-smooth setting is considered by Custódio, Dennis, Jr and Vicente (Reference Custódio, Dennis and Vicente2008). Hough, Kolda and Torczon (Reference Hough, Kolda and Torczon2001) discuss modifications to Algorithm 2that allow for increased efficiency when concurrent, asynchronous evaluations of $f$ are possible; an implementation of the method of Hough et al. (Reference Hough, Kolda and Torczon2001) is presented by Gray and Kolda (Reference Gray and Kolda2006).

The early analysis of Torczon (Reference Torczon1991, Section 7) of pattern-search methods when $f$ is non-smooth carries over to GPS methods as well; such methods may converge to a non-stationary point. This motivated a further generalization of GPS methods, mesh adaptive direct search (MADS) methods (Audet and Dennis, Jr Reference Audet and Dennis2006, Abramson and Audet Reference Abramson and Audet2006). Inspired by Coope and Price (Reference Coope and Price2000), MADS methods augment GPS methods by incorporating a mesh parametrized by a mesh parameter $\unicode[STIX]{x1D6FD}_{k}^{m}>0$ . In the $k$ th iteration, given the fixed PSS $\boldsymbol{D}$ and the mesh parameter $\unicode[STIX]{x1D6FD}_{k}^{m}$ , the MADS mesh around the current point $\boldsymbol{x}_{k}$ is

$$\begin{eqnarray}\displaystyle {\mathcal{M}}_{k}=\mathop{\bigcup }_{\boldsymbol{x}\in \boldsymbol{S}_{k}}\left\{\boldsymbol{x}+\unicode[STIX]{x1D6FD}_{k}^{m}\mathop{\sum }_{j=1}^{|\boldsymbol{D}|}\unicode[STIX]{x1D706}_{j}\boldsymbol{d}_{j}:\boldsymbol{d}_{j}\in \boldsymbol{D},\unicode[STIX]{x1D706}_{j}\in \mathbb{N}\bigcup \{0\}\right\}, & & \displaystyle \nonumber\end{eqnarray}$$

where $\boldsymbol{S}_{k}$ is the set of points at which $f$ has been evaluated prior to the $k$ th iteration of the method.

MADS methods additionally define a frame

$$\begin{eqnarray}\displaystyle {\mathcal{F}}_{k}=\{\boldsymbol{x}_{k}+\unicode[STIX]{x1D6FD}_{k}^{m}\boldsymbol{d}^{f}:\boldsymbol{d}^{f}\in \boldsymbol{D}_{k}^{f}\}, & & \displaystyle \nonumber\end{eqnarray}$$

where $\boldsymbol{D}_{k}^{f}$ is a finite set of directions, each of which is expressible as

$$\begin{eqnarray}\displaystyle \boldsymbol{d}^{f}=\mathop{\sum }_{j=1}^{|\boldsymbol{D}|}\unicode[STIX]{x1D706}_{j}\boldsymbol{d}_{j}, & & \displaystyle \nonumber\end{eqnarray}$$

with each $\unicode[STIX]{x1D706}_{j}\in \mathbb{N}\bigcup \{0\}$ and $\boldsymbol{d}_{j}\in \boldsymbol{D}$ . Additionally, MADS methods define a frame parameter $\unicode[STIX]{x1D6FD}_{k}^{f}$ and require that each $\boldsymbol{d}^{f}\in \boldsymbol{D}_{k}^{f}$ satisfies $\unicode[STIX]{x1D6FD}_{k}^{m}\Vert \boldsymbol{d}^{f}\Vert \leq \unicode[STIX]{x1D6FD}_{k}^{f}\max \{\Vert \boldsymbol{d}\Vert :\boldsymbol{d}\in \boldsymbol{D}\}$ . Observe that in each iteration, ${\mathcal{F}}_{k}\subsetneq {\mathcal{M}}_{k}$ . Note that the mesh is never explicitly constructed nor stored over the domain. Rather, points are evaluated only at what would be nodes of some implicitly defined mesh via the frame.

In the poll step of Algorithm 2, the set of poll directions $\boldsymbol{D}_{k}$ is chosen as $\{\boldsymbol{y}-\boldsymbol{x}_{k}:\boldsymbol{y}\in {\mathcal{F}}_{k}\}$ . The role of the step-size parameter $\unicode[STIX]{x1D6FC}_{k}$ in Algorithm 2is completely replaced by the behaviour of $\unicode[STIX]{x1D6FD}_{k}^{f},\unicode[STIX]{x1D6FD}_{k}^{m}$ . If there is no improvement at a candidate solution during the poll step, $\unicode[STIX]{x1D6FD}_{k}^{m}$ is decreased, resulting in a finer mesh; likewise $\unicode[STIX]{x1D6FD}_{k}^{f}$ is decreased, resulting in a finer local mesh around $\boldsymbol{x}_{k}$ . MADS intentionally allows the parameters $\unicode[STIX]{x1D6FD}_{k}^{m}$ and $\unicode[STIX]{x1D6FD}_{k}^{f}$ to be decreased at different rates; roughly speaking, by driving $\unicode[STIX]{x1D6FD}_{k}^{m}$ to zero faster than $\unicode[STIX]{x1D6FD}_{k}^{f}$ is driven to zero, and by choosing the sequence $\{\boldsymbol{D}_{k}^{f}\}$ to satisfy certain conditions, the directions in ${\mathcal{F}}_{k}$ become asymptotically dense around limit points of $\boldsymbol{x}_{k}$ . That is, it is possible to decrease $\unicode[STIX]{x1D6FD}_{k}^{m},\unicode[STIX]{x1D6FD}_{k}^{f}$ at rates such that poll directions will be arbitrarily close to any direction. This ensures that the Clarke directional derivative is non-negative in all directions around any limit point of the sequence of $\boldsymbol{x}_{k}$ generated by MADS; that is,

(2.1) $$\begin{eqnarray}f_{C}^{\prime }(\boldsymbol{x}_{\ast };\boldsymbol{d})\geq 0\quad \text{for all directions }\boldsymbol{d},\end{eqnarray}$$

with an analogous result also holding for constrained problems, with (2.1) reduced to all feasible directions $\boldsymbol{d}$ . (DDS methods for constrained optimization will be discussed in Section 7.) This powerful result highlights the ability of directional direct-search methods to address non-differentiable functions $f$ .

MADS does not prescribe any one approach for adjusting $\unicode[STIX]{x1D6FD}_{k}^{m},\unicode[STIX]{x1D6FD}_{k}^{f}$ so that the poll directions are dense, but Audet and Dennis, Jr (Reference Audet and Dennis2006) demonstrate an approach where randomized directions are completed to be a PSS and $\unicode[STIX]{x1D6FD}_{k}^{f}$ either is $n\sqrt{\unicode[STIX]{x1D6FD}_{k}^{m}}$ or $\sqrt{\unicode[STIX]{x1D6FD}_{k}^{m}}$ results in a asymptotically dense poll directions for any convergent subsequence of $\{\boldsymbol{x}_{k}\}$ . MADS does not require a sufficient-decrease condition.

Recent advances to MADS-based algorithms have focused on reducing the number of function evaluations required in practice by adaptively reducing the number of poll points queried; see, for example, Audet, Ianni, Le Digabel and Tribes (Reference Audet, Ianni, Le Digabel and Tribes2014) and Alarie et al. (Reference Alarie, Amaioua, Audet, Le Digabel and Leclaire2018). Smoothing-based extensions to noisy deterministic problems include Audet, Ihaddadene, Le Digabel and Tribes (Reference Audet, Ihaddadene, Le Digabel and Tribes2018b ). Vicente and Custódio (Reference Vicente and Custódio2012) show that MADS methods converge to local minima even for a limited class of discontinuous functions that satisfy some assumptions concerning the behaviour of the disconnected regions of the epigraph at limit points.

Worst-case complexity analysis.

Throughout this survey, when discussing classes of methods, we will refer to their worst-case complexity (WCC). Generally speaking, WCC refers to an upper bound on the number of function evaluations $N_{\unicode[STIX]{x1D716}}$ required to attain an $\unicode[STIX]{x1D716}$ -accurate solution to a problem drawn from a problem class. Correspondingly, the definition of $\unicode[STIX]{x1D716}$ -accurate varies between different problem classes. For instance, and of particular immediate importance, if an objective function is assumed Lipschitz-continuously differentiable (which we denote by $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ ), then an appropriate notion of first-order $\unicode[STIX]{x1D716}$ -accuracy is

(2.2) $$\begin{eqnarray}\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{k})\Vert \leq \unicode[STIX]{x1D716}.\end{eqnarray}$$

That is, the WCC of a method applied to the class ${\mathcal{L}}{\mathcal{C}}^{1}$ is characterized by $N_{\unicode[STIX]{x1D716}}$ , an upper bound on the number of function evaluations the method requires before (2.2) is satisfied for any $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ . Similarly, we can define a notion of second-order $\unicode[STIX]{x1D716}$ -accuracy as

(2.3) $$\begin{eqnarray}\max \{\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{k})\Vert ,-\unicode[STIX]{x1D706}_{k}\}\leq \unicode[STIX]{x1D716},\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{k}$ denotes the minimum eigenvalue of $\unicode[STIX]{x1D6FB}^{2}f(\boldsymbol{x}_{k})$ .

Note that WCCs can only be derived for methods for which convergence results have been established. Indeed, in the problem class ${\mathcal{L}}{\mathcal{C}}^{1}$ , first-order convergence results canonically have the form

(2.4) $$\begin{eqnarray}\lim _{k\rightarrow \infty }\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{k})\Vert =0.\end{eqnarray}$$

The convergence in (2.4) automatically implies the weaker lim-inf-type result

(2.5) $$\begin{eqnarray}\liminf _{k\rightarrow \infty }\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{k})\Vert =0,\end{eqnarray}$$

from which it is clear that for any $\unicode[STIX]{x1D716}>0$ , there must exist finite $N_{\unicode[STIX]{x1D716}}$ so that (2.2) holds. In fact, in many works, demonstrating a result of the form (2.5) is a stepping stone to proving a result of the form (2.4). Likewise, demonstrating a second-order WCC of the form (2.3) depends on showing

(2.6) $$\begin{eqnarray}\lim _{k\rightarrow \infty }\max \{\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{k})\Vert ,-\unicode[STIX]{x1D706}_{k}\}=0,\end{eqnarray}$$

which guarantees the weaker lim-inf-type result

(2.7) $$\begin{eqnarray}\liminf _{k\rightarrow \infty }\max \{\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}_{k})\Vert ,-\unicode[STIX]{x1D706}_{k}\}=0.\end{eqnarray}$$

Proofs of convergence for DDS methods applied to functions $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ often rely on a (sub)sequence of positive spanning sets $\{\boldsymbol{D}_{k}\}$ satisfying

(2.8) $$\begin{eqnarray}\text{cm}(\boldsymbol{D}_{k})=\min _{\boldsymbol{v}\in \mathbb{R}^{n}\setminus \{\mathbf{0}\}}\max _{\boldsymbol{d}\in \boldsymbol{D}_{k}}\displaystyle \frac{\boldsymbol{d}^{\top }\boldsymbol{v}}{\Vert \boldsymbol{d}\Vert \Vert \boldsymbol{v}\Vert }\geq \unicode[STIX]{x1D705}>0,\end{eqnarray}$$

where $\text{cm}(\cdot )$ is the cosine measure of a set. Under Assumption (2.8), Vicente (Reference Vicente2013) obtains a WCC of type (2.2) for a method in the Algorithm 2framework. In that work, it is assumed that $\boldsymbol{Y}_{k}=\emptyset$ at every search step. Moreover, sufficient decrease is tested at line 4of Algorithm 1; in particular, Vicente (Reference Vicente2013) checks in this line whether $f(\boldsymbol{p}_{i})<f(\boldsymbol{x})-c\unicode[STIX]{x1D6FC}_{k}^{2}$ for some $c>0$ , where $\unicode[STIX]{x1D6FC}_{k}$ is the current step size in Algorithm 2. Under these assumptions, Vicente (Reference Vicente2013) demonstrates a WCC in $O(\unicode[STIX]{x1D716}^{-2})$ . Throughout this survey, we will refer to Table A.1 for more details concerning specific WCCs. In general, though, we will often summarize WCCs in terms of their $\unicode[STIX]{x1D716}$ -dependence, as this provides an asymptotic characterization of a method’s complexity in terms of the accuracy to which one wishes to solve a problem.

When $f\in {\mathcal{L}}{\mathcal{C}}^{2}$ , work by Gratton et al. (Reference Gratton, Royer and Vicente2016) essentially augments the DDS method analysed by Vicente (Reference Vicente2013), but forms an approximate Hessian via central differences from function evaluations obtained (for free) by using a particular choice of $\boldsymbol{D}_{k}$ . Gratton et al. (Reference Gratton, Royer and Vicente2016) then demonstrate that this augmentation of Algorithm 2has a subsequence that converges to a second-order stationary point. That is, they prove a convergence result of the form (2.7) and demonstrate a WCC result of type (2.3) in $O(\unicode[STIX]{x1D716}^{-3})$ (see Table A.1).

We are unaware of WCC results for MADS methods; this situation may be unsurprising since MADS methods are motivated by non-smooth problems, which depend on the generation of a countably infinite number of poll directions. However, WCC results are not necessarily impossible to obtain in structured non-smooth cases, which we discuss in Section 5. We will discuss a special case where smoothing functions of a non-smooth function are assumed to be available in Section 5.3.2.

2.2.1 Quality of smooth model approximation

A natural first indicator of the quality of a model used for optimization is the degree to which the model locally approximates the function $f$ and its derivatives. To say anything about the quality of such approximation, one must make an assumption about the smoothness of both the model and function. For the moment, we leave this assumption implicit, but it will be formalized in subsequent sections.

A function $m:\mathbb{R}^{n}\rightarrow \mathbb{R}$ is said to be a $\boldsymbol{\unicode[STIX]{x1D705}}$ -fully linear model of $f$ on ${\mathcal{B}}(\boldsymbol{x};\unicode[STIX]{x1D6E5})=\{\boldsymbol{y}~:~\Vert \boldsymbol{x}-\boldsymbol{y}\Vert \leq \unicode[STIX]{x1D6E5}\}$ if

(2.9a) $$\begin{eqnarray}\displaystyle |f(\boldsymbol{x}+\boldsymbol{s})-m(\boldsymbol{x}+\boldsymbol{s})| & \leq & \displaystyle \unicode[STIX]{x1D705}_{\text{ef}}\unicode[STIX]{x1D6E5}^{2},\quad \text{for all}~\boldsymbol{s}\in {\mathcal{B}}(\mathbf{0};\unicode[STIX]{x1D6E5}),\end{eqnarray}$$

(2.9b) $$\begin{eqnarray}\displaystyle \Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}+\boldsymbol{s})-\unicode[STIX]{x1D6FB}m(\boldsymbol{x}+\boldsymbol{s})\Vert & \leq & \displaystyle \unicode[STIX]{x1D705}_{\text{eg}}\unicode[STIX]{x1D6E5},\quad \text{for all}~\boldsymbol{s}\in {\mathcal{B}}(\mathbf{0};\unicode[STIX]{x1D6E5}),\end{eqnarray}$$

$\boldsymbol{\unicode[STIX]{x1D705}}=(\unicode[STIX]{x1D705}_{\text{ef}},\unicode[STIX]{x1D705}_{\text{eg}})$

$\boldsymbol{\unicode[STIX]{x1D705}}=(\unicode[STIX]{x1D705}_{\text{ef}},\unicode[STIX]{x1D705}_{\text{eg}},\unicode[STIX]{x1D705}_{\text{eH}})$

$m$

$\boldsymbol{\unicode[STIX]{x1D705}}$

$f$

${\mathcal{B}}(\boldsymbol{x};\unicode[STIX]{x1D6E5})$

(2.10a) $$\begin{eqnarray}\displaystyle |f(\boldsymbol{x}+\boldsymbol{s})-m(\boldsymbol{x}+\boldsymbol{s})| & \leq & \displaystyle \unicode[STIX]{x1D705}_{\text{ef}}\unicode[STIX]{x1D6E5}^{3},\quad \text{for all}~\boldsymbol{s}\in {\mathcal{B}}(\mathbf{0};\unicode[STIX]{x1D6E5}),\end{eqnarray}$$

(2.10b) $$\begin{eqnarray}\displaystyle \Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x}+\boldsymbol{s})-\unicode[STIX]{x1D6FB}m(\boldsymbol{x}+\boldsymbol{s})\Vert & \leq & \displaystyle \unicode[STIX]{x1D705}_{\text{eg}}\unicode[STIX]{x1D6E5}^{2},\quad \text{for all}~\boldsymbol{s}\in {\mathcal{B}}(\mathbf{0};\unicode[STIX]{x1D6E5}),\end{eqnarray}$$

(2.10c) $$\begin{eqnarray}\displaystyle \Vert \unicode[STIX]{x1D6FB}^{2}f(\boldsymbol{x}+\boldsymbol{s})-\unicode[STIX]{x1D6FB}^{2}m(\boldsymbol{x}+\boldsymbol{s})\Vert & \leq & \displaystyle \unicode[STIX]{x1D705}_{\text{eH}}\unicode[STIX]{x1D6E5},\quad \text{for all}~\boldsymbol{s}\in {\mathcal{B}}(\mathbf{0};\unicode[STIX]{x1D6E5}).\end{eqnarray}$$

Models satisfying (2.9) or (2.10) are called Taylor-like models. To understand why, consider the second-order Taylor model

(2.11) $$\begin{eqnarray}m(\boldsymbol{x}+\boldsymbol{s})=f(\boldsymbol{x})+\unicode[STIX]{x1D6FB}f(\boldsymbol{x})^{\text{T}}\boldsymbol{s}+\displaystyle \frac{1}{2}\boldsymbol{s}^{\text{T}}\unicode[STIX]{x1D6FB}^{2}f(\boldsymbol{x})\boldsymbol{s}.\end{eqnarray}$$

This model is a $\boldsymbol{\unicode[STIX]{x1D705}}$ -fully quadratic model of $f$ , with

$$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D705}_{\text{ef}},\unicode[STIX]{x1D705}_{\text{eg}},\unicode[STIX]{x1D705}_{\text{eH}})=(L_{\text{H}}/6,L_{\text{H}}/2,L_{\text{H}}), & & \displaystyle \nonumber\end{eqnarray}$$

on any ${\mathcal{B}}(\boldsymbol{x};\unicode[STIX]{x1D6E5})$ , where $f$ has a Lipschitz-continuous second derivative with Lipschitz constant $L_{\text{H}}$ .

As illustrated in the next section, one also can guarantee that models that do not employ derivative information satisfy these approximation bounds in (2.9) or (2.10). This approximation quality is used by derivative-free algorithms to ensure that a sufficient reduction predicted by the model $m$ yields an attainable reduction in the function $f$ as $\unicode[STIX]{x1D6E5}$ becomes smaller.

2.2.2 Polynomial models

Polynomial models are the most commonly used models for derivative-free local optimization. We let ${\mathcal{P}}^{d,n}$ denote the space of polynomials of $n$ variables of degree $d$ and $\boldsymbol{\unicode[STIX]{x1D719}}~:~\mathbb{R}^{n}\rightarrow \mathbb{R}^{\text{dim}({\mathcal{P}}^{d,n})}$ define a basis for this space. For example, quadratic models can be obtained by using the monomial basis

(2.12) $$\begin{eqnarray}\boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{x})=[1,x_{1},\ldots ,x_{n},x_{1}^{2},\ldots x_{n}^{2},x_{1}x_{2},\ldots ,x_{n-1}x_{n}]^{\text{T}},\end{eqnarray}$$

for which $\text{dim}({\mathcal{P}}^{2,n})=(n+1)(n+2)/2$ ; linear models can be obtained by using the first $\text{dim}({\mathcal{P}}^{1,n})=n+1$ components of (2.12); quadratic models with diagonal Hessians, which are considered by Powell (Reference Powell2003), can be obtained by using the first $2n+1$ components of (2.12).

Any polynomial model $m\in {\mathcal{P}}^{d,n}$ is defined by $\boldsymbol{\unicode[STIX]{x1D719}}$ and coefficients $\boldsymbol{a}\in \mathbb{R}^{\text{dim}({\mathcal{P}}^{d,n})}$ through

(2.13)

Given a set of $p$ points $\boldsymbol{Y}=\{\boldsymbol{y}_{1},\ldots ,\boldsymbol{y}_{p}\}$ , a model that interpolates $f$ on $\boldsymbol{Y}$ is defined by the solution $\boldsymbol{a}$ to

(2.14) $$\begin{eqnarray}\boldsymbol{\unicode[STIX]{x1D6F7}}(\boldsymbol{Y})\boldsymbol{a}=\left[\begin{array}{@{}ccc@{}}\boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{y}_{1}) & \cdots \, & \boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{y}_{p})\end{array}\right]^{\text{T}}\boldsymbol{a}=\left[\begin{array}{@{}c@{}}f(\boldsymbol{y}_{1})\\ \vdots \\ f(\boldsymbol{y}_{p})\end{array}\right]\!.\end{eqnarray}$$

The existence, uniqueness and conditioning of a solution to (2.14) depend on the location of the sample points $\boldsymbol{Y}$ through the matrix $\boldsymbol{\unicode[STIX]{x1D6F7}}(\boldsymbol{Y})$ . We note that when $n>1$ , $|\boldsymbol{Y}|=\text{dim}({\mathcal{P}}^{d,n})$ is insufficient for guaranteeing that $\boldsymbol{\unicode[STIX]{x1D6F7}}(\boldsymbol{Y})$ is non-singular (Wendland Reference Wendland2005). Instead, additional conditions, effectively on the geometry of the sample points $\boldsymbol{Y}$ , must be satisfied.

Simplex gradients and linear interpolation models.

The geometry conditions needed to uniquely define a linear model are relatively straightforward: the sample points $\boldsymbol{Y}$ must be affinely independent; that is, the columns of

(2.15) $$\begin{eqnarray}\boldsymbol{Y}_{-1}=\left[\begin{array}{@{}ccc@{}}\boldsymbol{y}_{2}-\boldsymbol{y}_{1} & \cdots \, & \boldsymbol{y}_{n+1}-\boldsymbol{y}_{1}\end{array}\right]\end{eqnarray}$$

must be linearly independent. Such sample points define what is referred to as a simplex gradient $\boldsymbol{g}$ through $\boldsymbol{g}=[a_{2},\ldots ,a_{n+1}]^{\text{T}}$ , when the monomial basis $\unicode[STIX]{x1D719}$ is used in (2.14).

Simplex gradients can be viewed as a generalization of first-order finite-difference estimates (e.g. the forward differences based on evaluations at the points $\{\boldsymbol{y}_{1},\boldsymbol{y}_{1}+\unicode[STIX]{x1D6E5}\boldsymbol{e}_{1},\ldots ,\boldsymbol{y}_{1}+\unicode[STIX]{x1D6E5}\boldsymbol{e}_{n}\}$ ); their use in optimization algorithms dates at least back to the work of Spendley et al. (Reference Spendley, Hext and Himsworth1962) that inspired Nelder and Mead (Reference Nelder and Mead1965). Other example usage includes pattern search (Custódio and Vicente Reference Custódio and Vicente2007, Custódio et al. Reference Custódio, Dennis and Vicente2008) and noisy optimization (Kelley Reference Kelley1999b , Bortz and Kelley Reference Bortz, Kelley and Borggaard1998); the study of simplex gradients continues with recent works such as those of Regis (Reference Regis2015) and Coope and Tappenden (Reference Coope and Tappenden2019).

Provided that (2.15) is non-singular, it is straightforward to show that linear interpolation models are $\boldsymbol{\unicode[STIX]{x1D705}}$ -fully linear model of $f$ in a neighbourhood of $\boldsymbol{y}_{1}$ . In particular, if $\boldsymbol{Y}\subset {\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ and $f$ has an $L_{\text{g}}$ -Lipschitz-continuous first derivative on an open domain containing ${\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ , then (2.9) holds on ${\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ with

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D705}_{\text{eg}}=L_{\text{g}}(1+\sqrt{n}\unicode[STIX]{x1D6E5}\Vert \boldsymbol{Y}_{-1}^{-1}\Vert /2)\quad \text{and}\quad \unicode[STIX]{x1D705}_{\text{e f}}=L_{\text{g}}/2+\unicode[STIX]{x1D705}_{\text{eg}}.\end{eqnarray}$$

The expressions in (2.16) also provide a recipe for obtaining a model with a potentially tighter error bound over ${\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ : modify $\boldsymbol{Y}\subset {\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ to decrease $\Vert \boldsymbol{Y}_{-1}^{-1}\Vert$ . We note that when $\boldsymbol{Y}_{-1}$ contains orthonormal directions scaled by $\unicode[STIX]{x1D6E5}$ , one recovers $\unicode[STIX]{x1D705}_{\text{eg}}=L_{\text{g}}(1+\sqrt{n}/2)$ and $\unicode[STIX]{x1D705}_{\text{ef}}=L_{\text{g}}(3+\sqrt{n})/2$ , which is the least value one can obtain from (2.16) given the restriction that $\boldsymbol{Y}\subset {\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ . Hence, by performing LU or QR factorization with pivoting, one can obtain directions (which are then scaled by $\unicode[STIX]{x1D6E5}$ ) in order to improve the conditioning of $\boldsymbol{Y}_{-1}^{-1}$ and hence the approximation bound. Such an approach is performed by Conn, Scheinberg and Vicente (Reference Conn, Scheinberg and Vicente2008a ) for linear models and by Wild and Shoemaker (Reference Wild and Shoemaker2011) for fully linear radial basis function models.

The geometric conditions on $\boldsymbol{Y}$ , induced by the approximation bounds in (2.9) or (2.10), can be viewed as playing a similar role to the geometric conditions (e.g. positive spanning) imposed on $\boldsymbol{D}$ in directional direct-search methods. Naturally, the choice of basis function used for any model affects the quantitative measure of that model’s quality.

Note that many practical methods employ interpolation sets contained within a constant multiple of the trust-region radius (i.e.  $\boldsymbol{Y}\subset {\mathcal{B}}(\boldsymbol{y}_{1};c_{1}\unicode[STIX]{x1D6E5})$ for a constant $c_{1}\in [1,\infty )$ ).

Quadratic interpolation models.

Quadratic interpolation models have been used for derivative-free optimization for at least fifty years (Winfield Reference Winfield1969, Reference Winfield1973) and were employed by a series of methods that revitalized interest in model-based methods; see, for example, Conn and Toint (Reference Conn, Toint, Di Pillo and Giannessi1996), Conn, Scheinberg and Toint (Reference Conn, Scheinberg and Toint1997b ), Conn, Scheinberg and Toint (Reference Conn, Scheinberg, Toint, Iserles and Buhmann1997a ) and Powell (Reference Powell and Yuan1998b , Reference Powell2002).

Of course, the quality of an interpolation model (quadratic or otherwise) in a region of interest is determined by the position of the underlying points being interpolated. For example, if a model $m$ interpolates a function $f$ at points far away from a certain region of interest, the model value may differ greatly from the value of $f$ in that region. $\unicode[STIX]{x1D6EC}$ -poisedness is a concept to measure how well a set of points is dispersed through a region of interest, and ultimately how well a model will estimate the function in that region.

The most commonly used metric for quantifying how well points are positioned in a region of interest is based on Lagrange polynomials. Given a set of $p$ points $\boldsymbol{Y}=\{\boldsymbol{y}_{1},\ldots ,\boldsymbol{y}_{p}\}$ , a basis of Lagrange polynomials satisfies

(2.17) $$\begin{eqnarray}\ell _{j}(\boldsymbol{y}_{i})=\left\{\begin{array}{@{}ll@{}}1\quad & \text{if }i=j,\\ 0\quad & \text{if }i\neq j.\end{array}\right.\end{eqnarray}$$

We now define $\unicode[STIX]{x1D6EC}$ -poisedness. A set of points $\boldsymbol{Y}$ is said to be $\unicode[STIX]{x1D6EC}$ -poised on a set $\boldsymbol{B}$ if $\boldsymbol{Y}$ is linearly independent and the Lagrange polynomials $\{\ell _{1},\ldots ,\ell _{p}\}$ associated with $\boldsymbol{Y}$ satisfy

(2.18) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}\geq \max _{1\leq i\leq p}\max _{\boldsymbol{x}\in B}|\ell _{i}(\boldsymbol{x})|.\end{eqnarray}$$

(For an equivalent definition of $\unicode[STIX]{x1D6EC}$ -poisedness, see Conn et al. (Reference Conn, Scheinberg and Vicente2009b , Definition 3.6).) Note that the definition of $\unicode[STIX]{x1D6EC}$ -poisedness is independent of the function being modelled. Also, the points $\boldsymbol{Y}$ need not necessarily be elements of the set $\boldsymbol{B}$ . Also, note that if a model is poised on a set $\boldsymbol{B}$ , it is poised on any subset of $\boldsymbol{B}$ . One is usually interested in the least value of $\unicode[STIX]{x1D6EC}$ so that (2.18) holds.

Powell’s unconstrained optimization by quadratic approximation method (UOBYQA) follows such an approach in maximizing the Lagrange polynomials. In Powell (Reference Powell and Yuan1998b ), Powell (Reference Powell2001) and Powell (Reference Powell2002), significant care is given to the linear algebra expense associated with this maximization and the associated change of basis as the methods change their interpolation sets. For example, in Powell (Reference Powell and Yuan1998b ), particular sparsity in the Hessian approximation is employed with the aim of capturing curvature while keeping linear algebraic expenses low.

Maintaining, and the question of to what extent it is necessary to maintain, this geometry for quadratic models has been intensely studied; see, for example, Fasano, Morales and Nocedal (Reference Fasano, Morales and Nocedal2009), Marazzi and Nocedal (Reference Marazzi and Nocedal2002), D’Ambrosio, Nannicini and Sartor (Reference D’Ambrosio, Nannicini and Sartor2017) and Scheinberg and Toint (Reference Scheinberg and Toint2010).

Underdetermined quadratic interpolation models.

A fact not to be overlooked in the context of derivative-free optimization is that employing an interpolation set $\boldsymbol{Y}$ requires availability of the $|\boldsymbol{Y}|$ function values $\{f(\boldsymbol{y}_{i}):\boldsymbol{y}_{i}\in \boldsymbol{Y}\}$ . When the function $f$ is computationally expensive to evaluate, the $(n+1)(n+2)/2$ points required by fully quadratic models can be a burden, potentially with little benefit, to obtain repeatedly in an optimization algorithm.

Beginning with Powell (Reference Powell2003), Powell investigated quadratic models constructed from fewer than $(n+1)(n+2)/2$ points. The most successful of these strategies was detailed in Powell (Reference Powell2004a ) and Powell (Reference Powell2004b ) and resolved the $(n+1)(n+2)/2-|\boldsymbol{Y}|$ remaining degrees of freedom by solving problems of the form

(2.19)

to obtain a model $m$ about a point of interest $\check{\boldsymbol{x}}$ . Solutions to (2.19) are models with a Hessian closest in Frobenius norm to a specified $\boldsymbol{H}=\boldsymbol{H}^{\text{T}}$ among all models that interpolate $f$ on $\boldsymbol{Y}$ . A popular implementation of this strategy is the NEWUOA solver (Powell Reference Powell, Di Pillo and Roma2006).

By using the basis

(2.20) $$\begin{eqnarray}\displaystyle \begin{array}{@{}lll@{}}\boldsymbol{\unicode[STIX]{x1D719}}(\check{\boldsymbol{x}}+\boldsymbol{x}) & = & \displaystyle [\boldsymbol{\unicode[STIX]{x1D719}}_{\text{fg}}(\check{\boldsymbol{x}}+\boldsymbol{x})^{\text{T}}|\,\boldsymbol{\unicode[STIX]{x1D719}}_{\text{H}}(\check{\boldsymbol{x}}+\boldsymbol{x})^{\text{T}}]^{\text{T}}\\ & = & \displaystyle \biggl[1,\,x_{1},\,\ldots ,\,x_{n}\bigg|\,\displaystyle \frac{1}{2}x_{1}^{2},\,\ldots ,\,\displaystyle \frac{1}{2}x_{n}^{2},\,\displaystyle \frac{1}{\sqrt{2}}x_{1}x_{2},\,\ldots ,\displaystyle \frac{1}{\sqrt{2}}x_{n-1}x_{n}\biggr]^{\text{T}},\end{array} & & \displaystyle\end{eqnarray}$$

the problem (2.19) is equivalent to the problem

(2.21)

Existence and uniqueness of solutions to (2.21) again depend on the positioning of the points in $\boldsymbol{Y}$ . Notably, a necessary condition for there to be a unique minimizer of the seminorm is that at least $n+1$ of the points in $\boldsymbol{Y}$ be affinely independent. Lagrange polynomials can be defined for this case; Conn, Scheinberg and Vicente (Reference Conn, Scheinberg and Vicente2008b ) establish conditions for $\unicode[STIX]{x1D6EC}$ -poisedness (and hence a fully linear, or better, approximation quality) of such models.

Powell (Reference Powell2004c , Reference Powell2007, Reference Powell2008) develops efficient solution methodologies for (2.21) when $\boldsymbol{H}$ and $m$ are constructed from interpolation sets that differ by at most one point, and employ these updates in NEWUOA and subsequent solvers. Wild (Reference Wild2008b ) and Custódio et al. (Reference Custódio, Rocha and Vicente2009) use $\boldsymbol{H}=\mathbf{0}$ in order to obtain tighter fully linear error bounds of models resulting from (2.21). A strategy of using even fewer interpolation points (including those in a proper subspace of $\mathbb{R}^{n}$ ) is developed by Powell (Reference Powell2013) and Zhang (Reference Zhang2014). In Section 5.2, we summarize approaches that exploit knowledge of sparsity of the derivatives of $f$ in building quadratic models that interpolate fewer than $(n+1)(n+2)/2$ points.

Figure 2.2.

(a) Minimum-norm-Hessian model through five points in ${\mathcal{B}}(\boldsymbol{x}_{k};\unicode[STIX]{x1D6E5}_{k})$ and its minimizer. (b) Absolute value of a sixth Lagrange polynomial for the five points. (c) Minimum-norm-Hessian model through five points in ${\mathcal{B}}(\boldsymbol{x}_{k+1};\unicode[STIX]{x1D6E5}_{k+1})$ and its minimizer. (d) Absolute value of a sixth Lagrange polynomial for the five points.

Figure 2.2 shows quadratic models in two dimensions that interpolate $(n+1)(n+2)/2-1=5$ points as well as the associated magnitude of the remaining Lagrange polynomial (note that this polynomial vanishes at the five interpolated points).

Regression models.

Just as one can establish approximation bounds and geometry conditions when $\boldsymbol{Y}$ is linearly independent, the same can be done for overdetermined regression models (Conn et al. Reference Conn, Scheinberg and Vicente2008b , Reference Conn, Scheinberg and Vicente2009b ). This can be accomplished by extending the definition of Lagrange polynomials from (2.17) to the regression case. That is, given a basis $\boldsymbol{\unicode[STIX]{x1D719}}~:~\mathbb{R}^{n}\rightarrow \mathbb{R}^{\text{dim}({\mathcal{P}}^{d,n})}$ and points $\boldsymbol{Y}=\{\boldsymbol{y}_{1},\ldots ,\boldsymbol{y}_{p}\}\subset \mathbb{R}^{n}$ with $p>\text{dim}({\mathcal{P}}^{d,n})$ , the set of polynomials satisfies

(2.22) $$\begin{eqnarray}\ell _{j}(\boldsymbol{y}_{i})\overset{\text{l.s.}}{=}\left\{\begin{array}{@{}ll@{}}1\quad & \text{if }i=j,\\ 0\quad & \text{if }i\neq j,\end{array}\right.\end{eqnarray}$$

where $\overset{\text{l.s.}}{=}$ denotes the least-squares solution. The regression model can be recovered finding the least-squares solution (now overdetermined) system from (2.14), and the definition of $\unicode[STIX]{x1D6EC}$ -poisedness (in the regression sense) is equivalent to (2.18). Ultimately, given a linear regression model through a set of $\unicode[STIX]{x1D6EC}$ -poised points $\boldsymbol{Y}\subset {\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ , and if $f$ has an $L_{\text{g}}$ -Lipschitz-continuous first derivative on an open domain containing ${\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ , then (2.9) holds on ${\mathcal{B}}(\boldsymbol{y}_{1};\unicode[STIX]{x1D6E5})$ with

(2.23) $$\begin{eqnarray}\unicode[STIX]{x1D705}_{\text{eg}}=\displaystyle \frac{5}{2}\sqrt{p}L_{\text{g}}\unicode[STIX]{x1D6EC}\quad \text{and}\quad \unicode[STIX]{x1D705}_{\text{ef}}=\displaystyle \frac{1}{2}L_{\text{g}}+\unicode[STIX]{x1D705}_{\text{eg}}.\end{eqnarray}$$

Conn et al. (Reference Conn, Scheinberg and Vicente2008b ) note the fact that the extension of Lagrange polynomials does not apply to the 1-norm or infinity-norm case. Billups, Larson and Graf (Reference Billups, Larson and Graf2013) show that the definition of Lagrange polynomials can be extended to the weighted regression case. Verdério, Karas, Pedroso and Scheinberg (Reference Verdério, Karas, Pedroso and Scheinberg2017) show that (2.9) can also be recovered for support vector regression models.

Efficiently minimizing the model (regardless of type) over a trust region is integral to the usefulness of such models within an optimization algorithm. In fact, this necessity is a primary reason for the use of low-degree polynomial models by the majority of derivative-free trust-region methods. For quadratic models, the resulting subproblem remains one of the most difficult non-convex optimization problems solvable in polynomial time, as illustrated by Moré and Sorensen (Reference Moré and Sorensen1983). As exemplified by Powell (Reference Powell and Griffiths1997), the implementation of subproblem solvers is a key concern in methods seeking to perform as few algebraic operations between function evaluations as possible.

2.2.3 Radial basis function interpolation models

An additional way to model non-linearity with potentially less restrictive geometric conditions is by using radial basis functions (RBFs). Such models take the form

(2.24) $$\begin{eqnarray}m(\boldsymbol{x})=\mathop{\sum }_{i=1}^{|\boldsymbol{Y}|}b_{i}\unicode[STIX]{x1D713}(\Vert \boldsymbol{x}-\boldsymbol{y}_{i}\Vert )+\boldsymbol{a}^{\text{T}}\boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{x}),\end{eqnarray}$$

where $\unicode[STIX]{x1D713}~:~\mathbb{R}_{+}\rightarrow \mathbb{R}$ is a conditionally positive-definite univariate function and $\boldsymbol{a}^{\text{T}}\boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{x})$ represents a (typically low-order) polynomial as before; see, for example, Buhmann (Reference Buhmann2000). Given a sample set $\boldsymbol{Y}$ , RBF model coefficients $(\boldsymbol{a},\boldsymbol{b})$ can be obtained by solving the augmented interpolation equations

(2.25) $$\begin{eqnarray}\left[\begin{array}{@{}cccc@{}}\unicode[STIX]{x1D713}(\Vert \boldsymbol{y}_{1}-\boldsymbol{y}_{1}\Vert ) & \cdots \, & \unicode[STIX]{x1D713}(\Vert \boldsymbol{y}_{1}-\boldsymbol{y}_{|\boldsymbol{Y}|}\Vert ) & \boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{y}_{1})^{\text{T}}\\ \vdots & & \vdots & \vdots \\ \unicode[STIX]{x1D713}(\Vert \boldsymbol{y}_{|\boldsymbol{Y}|}-\boldsymbol{y}_{1}\Vert ) & \cdots \, & \unicode[STIX]{x1D713}(\Vert \boldsymbol{y}_{|\boldsymbol{Y}|}-\boldsymbol{y}_{|\boldsymbol{Y}|}\Vert ) & \boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{y}_{|\boldsymbol{Y}|})^{\text{T}}\\ \boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{y}_{1}) & \cdots \, & \boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{y}_{|\boldsymbol{Y}|}) & \mathbf{0}\end{array}\right]\left[\begin{array}{@{}c@{}}\boldsymbol{b}\\ \boldsymbol{a}\end{array}\right]=\left[\begin{array}{@{}c@{}}f(\boldsymbol{y}_{1})\\ \vdots \\ f(\boldsymbol{y}_{|\boldsymbol{Y}|})\\ \mathbf{0}\end{array}\right]\!.\end{eqnarray}$$

That RBFs are conditionally positive-definite ensures that (2.25) is non-singular provided that the degree $d$ of the polynomial $\boldsymbol{\unicode[STIX]{x1D719}}$ is sufficiently large and that $\boldsymbol{Y}$ is poised for degree- $d$ polynomial interpolation. For example, cubic ( $\unicode[STIX]{x1D713}(r)=r^{3}$ ) RBFs require a linear polynomial; multiquadric ( $\unicode[STIX]{x1D713}(r)=-(\unicode[STIX]{x1D6FE}^{2}+r^{2})^{1/2}$ ) RBFs require a constant polynomial; and inverse multiquadric ( $\unicode[STIX]{x1D713}(r)=(\unicode[STIX]{x1D6FE}^{2}+r^{2})^{-1/2}$ ) and Gaussian ( $\unicode[STIX]{x1D713}(r)=\exp (-\unicode[STIX]{x1D6FE}^{-2}r^{2})$ ) RBFs do not require a polynomial. Consequently, RBFs have relatively unrestrictive geometric requirements on the interpolation points $\boldsymbol{Y}$ while allowing for modelling a wide range of non-linear behaviour.

This feature is typically exploited in global optimization (see e.g. Björkman and Holmström Reference Björkman and Holmström2000, Gutmann Reference Gutmann2001 and Regis and Shoemaker Reference Regis and Shoemaker2007), whereby an RBF surrogate model is employed to globally approximate $f$ . However, works such as Oeuvray and Bierlaire (Reference Oeuvray and Bierlaire2009), Oeuvray (Reference Oeuvray2005), Wild (Reference Wild2008a ) and Wild and Shoemaker (Reference Wild and Shoemaker2013) establish and use local approximation properties of these models. This approach is typically performed by relying on a linear polynomial $\boldsymbol{a}^{\text{T}}\boldsymbol{\unicode[STIX]{x1D719}}(\boldsymbol{x})$ , which can be used to establish that the RBF model in (2.24) can be a fully linear local approximation of smooth  $f$ .

2.2.4 Trust-region methods

Having discussed issues of model construction, we are now ready to present a general statement of a model-based trust-region method in Algorithm 3.

A distinguishing characteristic of derivative-free model-based trust-region methods is how they manage $\boldsymbol{Y}_{k}$ , the set of points used to construct the model $m_{k}$ . Some methods ensure that $\boldsymbol{Y}_{k}$ contains a scaled stencil of points around $\boldsymbol{x}_{k}$ ; such an approach can be attractive since the objective at such points can be evaluated in parallel. A fixed stencil can also ensure that all models sufficiently approximate the objective. Other methods construct $\boldsymbol{Y}$ by using previously evaluated points near $\boldsymbol{x}_{k}$ , for example, those points within ${\mathcal{B}}(\boldsymbol{x}_{k};c_{1}\unicode[STIX]{x1D6E5}_{k})$ for some constant $c_{1}\in [1,\infty )$ . Depending on the set of previously evaluated points, such methods may need to add points to $\boldsymbol{Y}_{k}$ that most improve the model quality. Determining which additional points to add to $\boldsymbol{Y}_{k}$ can be computationally expensive, but the method should be willing to do so in the hope of needing fewer evaluations of the objective function at new points in $\boldsymbol{Y}_{k}$ . Most methods do not ensure that models are valid on every iteration but rather make a single step toward improving the model. Such an approach can ensure a high-quality model in a finite number of improvement steps. (Exceptional methods that ensure model quality before $\boldsymbol{s}_{k}$ is calculated are the methods of Powell and manifold sampling of Khan, Larson and Wild (Reference Khan, Larson and Wild2018).) The ORBIT method (Wild, Regis and Shoemaker Reference Wild, Regis and Shoemaker2008) places a limit on the size of $\boldsymbol{Y}_{k}$ (e.g. in order to limit the amount of linear algebra or to prevent overfitting). In the end, such restrictions on $\boldsymbol{Y}_{k}$ may determine whether $m_{k}$ is an interpolation or regression model.

Derivative-free trust-region methods share many similarities with traditional trust-region methods, for example, the use of a $\unicode[STIX]{x1D70C}$ -test to determine whether a step is taken or rejected. As in a traditional trust-region method, the $\unicode[STIX]{x1D70C}$ -test measures the ratio of actual decrease observed in the objective versus the decrease predicted by the model.

On the other hand, the management of the trust-region radius parameter $\unicode[STIX]{x1D6E5}_{k}$ in Algorithm 3differs remarkably from traditional trust-region methods. Derivative-free variants require an additional test of model quality, the failure of which results in shrinking $\unicode[STIX]{x1D6E5}_{k}$ . When derivatives are available, Taylor’s theorem ensures model accuracy for small $\unicode[STIX]{x1D6E5}_{k}$ . In the derivative-free case, such a condition must be explicitly checked in order to ensure that $\unicode[STIX]{x1D6E5}_{k}$ does not go to zero merely because the model is poor, hence the inclusion of tests of model quality. As a direct result of these considerations, $\unicode[STIX]{x1D6E5}_{k}\rightarrow 0$ as Algorithm 3converges; this is generally not the case in traditional trust-region methods.

As in derivative-based trust-region methods, the solution to the trust-region subproblem in line 11of Algorithm 3must satisfy a Cauchy decrease condition. Given the model $m_{k}$ used in Algorithm 3,we define the optimal step length in the direction $-\unicode[STIX]{x1D6FB}m_{k}(\boldsymbol{x}_{k})$ by

and the corresponding Cauchy step

$$\begin{eqnarray}\displaystyle \boldsymbol{s}_{k}^{C}=-t_{k}^{C}\unicode[STIX]{x1D6FB}m_{k}(\boldsymbol{x}_{k}). & & \displaystyle \nonumber\end{eqnarray}$$

It is straightforward to show (see e.g. Conn, Scheinberg and Vicente Reference Conn, Scheinberg and Vicente2009b , Theorem 10.1) that

(2.26) $$\begin{eqnarray}m_{k}(\boldsymbol{x}_{k})-m_{k}(\boldsymbol{x}_{k}+\boldsymbol{s}_{k}^{C})\geq \displaystyle \frac{1}{2}\Vert \unicode[STIX]{x1D6FB}m_{k}(\boldsymbol{x}_{k})\Vert \min \left\{\displaystyle \frac{\Vert \unicode[STIX]{x1D6FB}m_{k}(\boldsymbol{x}_{k})\Vert }{\Vert \unicode[STIX]{x1D6FB}^{2}m_{k}(\boldsymbol{x}_{k})\Vert },\unicode[STIX]{x1D6E5}_{k}\right\}.\end{eqnarray}$$

That is, (2.26) states that, provided that both $\unicode[STIX]{x1D6E5}_{k}\approx \Vert \unicode[STIX]{x1D6FB}m_{k}(\boldsymbol{x}_{k})\Vert$ and a uniform bound exists on the norm of the model Hessian, the model decrease attained by the Cauchy step $\boldsymbol{s}_{k}^{C}$ is of the order of $\unicode[STIX]{x1D6E5}_{k}^{2}$ . In order to prove convergence, it is desirable to ensure that each step $\boldsymbol{s}_{k}$ generated in line 11of Algorithm 3decreases the model $m_{k}$ by no less than $\boldsymbol{s}_{k}^{C}$ does, or at least some fixed positive fraction of the decrease achieved by $\boldsymbol{s}_{k}^{C}$ . Because successful iterations ensure that the actual decrease attained in an iteration is at least a constant fraction of the model decrease, the sequence of decreases of Algorithm 3are square-summable, provided that $\unicode[STIX]{x1D6E5}_{k}\rightarrow 0$ . (This is indeed the case for derivative-free trust-region methods.) Hence, in most theoretical treatments of these methods, it is commonly stated as an assumption that the subproblem solution $\boldsymbol{s}_{k}$ obtained in line 11of Algorithm 3satisfies

(2.27) $$\begin{eqnarray}m_{k}(\boldsymbol{x}_{k})-m_{k}(\boldsymbol{x}_{k}+\boldsymbol{s}_{k})\geq \unicode[STIX]{x1D705}_{\text{fcd}}(m_{k}(\boldsymbol{x}_{k})-m_{k}(\boldsymbol{x}_{k}+\boldsymbol{s}_{k}^{C})),\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{\text{fcd}}\in (0,1]$ is the fraction of the Cauchy decrease. In practice, when $m_{k}$ is a quadratic model, subproblem solvers have been well studied and often come with guarantees concerning the satisfaction of (2.27) (Conn, Gould and Toint Reference Conn, Gould and Toint2000). Wild et al. (Reference Wild, Regis and Shoemaker2008, Figure 4.3) demonstrate the satisfaction of an assumption like (2.27) when the model $m_{k}$ is a radial basis function.

Under reasonable smoothness assumptions, most importantly $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ , algorithms in the Algorithm 3framework have been shown to be first-order convergent (i.e. (2.4)) and second-order convergent (i.e. (2.6)), with the (arguably) most well-known proof given by Conn, Scheinberg and Vicente (Reference Conn, Scheinberg and Vicente2009a ). In more recent work, Garmanjani, Jùdice and Vicente (Reference Garmanjani, Jùdice and Vicente2016) provide a WCC bound of the form (2.2) for Algorithm 3, recovering essentially the same upper bound on the number of function evaluations required by DDS methods found in Vicente (Reference Vicente2013), that is, a WCC bound in $O(\unicode[STIX]{x1D716}^{-2})$ (see Table A.1). When $f\in {\mathcal{L}}{\mathcal{C}}^{2}$ , Gratton, Royer and Vicente (Reference Gratton, Royer and Vicente2019a ) demonstrate a second-order WCC bound of the form (2.3) in $O(\unicode[STIX]{x1D716}^{-3})$ ; in order to achieve this result, fully quadratic models $m_{k}$ are required. In Section 3.3, a similar result is achieved by using randomized variants that do not require a fully quadratic model in every iteration.

Early analysis of Powell’s UOBYQA method shows that, with minor modifications, the algorithm can converge superlinearly in neighbourhoods of strict convexity (Han and Liu Reference Han and Liu2004). A key distinction between Powell’s methods and other model-based trust-region methods is the use of separate neighbourhoods for model quality and trust-region steps, with each of these neighbourhoods changing dynamically. Convergence of such methods is addressed by Powell (Reference Powell2010, Reference Powell2012).

The literature on derivative-free trust-region methods is extensive. We mention in passing several additional classes of trust-region methods that have not fallen neatly into our discussion thus far. Wedge methods (Marazzi and Nocedal Reference Marazzi and Nocedal2002) explicitly enforce geometric properties ( $\unicode[STIX]{x1D6EC}$ -poisedness) of the sample set between iterations by adding additional constraints to the trust-region subproblem. Alexandrov, Dennis, Jr, Lewis and Torczon (Reference Alexandrov, Dennis, Lewis and Torczon1998) consider a trust-region method utilizing a hierarchy of model approximations. In particular, if derivatives can be obtained but are expensive, then the method of Alexandrov et al. (Reference Alexandrov, Dennis, Lewis and Torczon1998) uses a model that interpolates not only zeroth-order information but also first-order (gradient) information. For problems with deterministic noise, Elster and Neumaier (Reference Elster and Neumaier1995) propose a method that projects the solutions of a trust-region subproblem onto a dynamically refined grid, encouraging better practical behaviour. Similarly, for problems with deterministic noise, Maggiar, Wächter, Dolinskaya and Staum (Reference Maggiar, Wächter, Dolinskaya and Staum2018) propose a model-based trust-region method that implicitly convolves the objective function with a Gaussian kernel, again yielding better practical behaviour.

2.3.1 Finite differences

As noted in Section 1.1.2, many of the earliest derivative-free methods employed finite-difference-based estimates of derivatives. The most popular first-order directional derivative estimates include the forward/reverse difference

(2.28) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{\text{f}}(f;\boldsymbol{x};\boldsymbol{d};h)=\displaystyle \frac{f(\boldsymbol{x}+h\boldsymbol{d})-f(\boldsymbol{x})}{h}\end{eqnarray}$$

and central difference

(2.29) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{\text{c}}(f;\boldsymbol{x};\boldsymbol{d};h)=\displaystyle \frac{f(\boldsymbol{x}+h\boldsymbol{d})-f(\boldsymbol{x}-h\boldsymbol{d})}{2h},\end{eqnarray}$$

where $h\neq 0$ is the difference parameter and the non-trivial $\boldsymbol{d}\in \mathbb{R}^{n}$ defines the direction. Several recent methods, including the methods described in Sections 2.3.2, 2.3.3 and 3.1.2, use such estimates and employ difference parameters or directions that dynamically change.

As an example of a potentially dynamic choice of difference parameter, we consider the usual case of roundoff errors. We denote by $f_{\infty }^{\prime }(\boldsymbol{x};\boldsymbol{d})$ the directional derivative at $\boldsymbol{x}$ of the infinite-precision (i.e. based on real arithmetic) objective function $f_{\infty }$ in the unit direction $\boldsymbol{d}$ (i.e.  $\Vert \boldsymbol{d}\Vert =1$ ). We then have the following error for forward or reverse finite-difference estimates based on the function $f$ available through computation:

(2.30) $$\begin{eqnarray}|\unicode[STIX]{x1D6FF}_{\text{f}}(f;\boldsymbol{x};\boldsymbol{d};h)-f_{\infty }^{\prime }(\boldsymbol{x};\boldsymbol{d})|\leq \displaystyle \frac{1}{2}L_{\text{g}}(\boldsymbol{x})|h|+2\displaystyle \frac{\unicode[STIX]{x1D716}_{\infty }(\boldsymbol{x})}{|h|},\end{eqnarray}$$

provided that $|f_{\infty }^{\prime \prime }(\cdot ;\boldsymbol{d})|\leq L_{\text{g}}(\boldsymbol{x})$ and $|f_{\infty }(\cdot )-f(\cdot )|\leq \unicode[STIX]{x1D716}_{\infty }(\boldsymbol{x})$ on the interval $[\boldsymbol{x},\boldsymbol{x}+h\boldsymbol{d}]$ . In Gill, Murray and Wright (Reference Gill, Murray and Wright1981) and Gill, Murray, Saunders and Wright (Reference Gill, Murray, Saunders and Wright1983), the recommended difference parameter is $h=2\sqrt{\unicode[STIX]{x1D716}_{\infty }(\boldsymbol{x})/L_{\text{g }}(\boldsymbol{x})}$ , which yields the minimum value $2\sqrt{\unicode[STIX]{x1D716}_{\infty }(\boldsymbol{x})L_{\text{g }}(\boldsymbol{x})}$ of the upper bound in (2.30); when $\unicode[STIX]{x1D716}_{\infty }$ is a bound on the roundoff error and $L_{\text{g}}$ is of order one, then the familiar $h\in O(\sqrt{\unicode[STIX]{x1D716}_{\infty }})$ is obtained.

Similarly, if one models the error between $f_{\infty }$ and $f$ as a stationary stochastic process (through the ansatz denoted by $f_{\boldsymbol{\unicode[STIX]{x1D709}}}$ ) with variance $\unicode[STIX]{x1D716}_{\text{f}}(\boldsymbol{x})^{2}$ , minimizing the upper bound on the mean-squared error,

(2.31) $$\begin{eqnarray}\mathbb{E}_{\boldsymbol{\unicode[STIX]{x1D709}}}[(\unicode[STIX]{x1D6FF}_{\text{f}}(f_{\boldsymbol{\unicode[STIX]{x1D709}}};\boldsymbol{x};\boldsymbol{d};h)-f_{\infty }^{\prime }(\boldsymbol{x};\boldsymbol{d}))^{2}]\leq \displaystyle \frac{1}{4}L_{\text{g}}(\boldsymbol{x})^{2}h^{2}+2\displaystyle \frac{\unicode[STIX]{x1D716}_{\text{f}}(\boldsymbol{x})^{2}}{h^{2}},\end{eqnarray}$$

yields the choice $h=(\sqrt{8}\unicode[STIX]{x1D716}_{\text{f}}(\boldsymbol{x})/L_{\text{g}}(\boldsymbol{x}))^{1/2}$ with an associated root-mean-squared error of $(\sqrt{2}\unicode[STIX]{x1D716}_{\infty }(\boldsymbol{x})L_{\text{g}}(\boldsymbol{x}))^{1/2}$ ; see, for example, Moré and Wild (Reference Moré and Wild2012, Reference Moré and Wild2014). A rough procedure for computing $\unicode[STIX]{x1D716}_{\text{f}}$ is provided in Moré and Wild (Reference Moré and Wild2011) and used in recent methods such as that of Berahas, Byrd and Nocedal (Reference Berahas, Byrd and Nocedal2019).

In both cases (2.30) and (2.31), the first-order error is $c\sqrt{\unicode[STIX]{x1D716}(\boldsymbol{x})L_{\text{g }}(\boldsymbol{x})}$ (for a constant $c\leq 2$ ), which can be used to guide the decision on whether the derivatives estimates are of sufficient accuracy.

2.3.2 Implicit filtering

Implicit filtering is a hybrid of a grid-search algorithm (evaluating all points on a lattice) and a Newton-like local optimization method. The gradient (and possible Hessian) estimates for local optimization are approximated by the central differences $\{\unicode[STIX]{x1D6FF}_{\text{c}}(f;\boldsymbol{x}_{k};\boldsymbol{e}_{i};\unicode[STIX]{x1D6E5}_{k}):i=1,\ldots ,n\}$ . The difference parameter $\unicode[STIX]{x1D6E5}_{k}$ decreases when implicit filtering encounters a stencil failure at $\boldsymbol{x}_{k}$ , that is,

(2.32) $$\begin{eqnarray}f(\boldsymbol{x}_{k})\leq f(\boldsymbol{x}_{k}\pm \unicode[STIX]{x1D6E5}_{k}\boldsymbol{e}_{i}),\end{eqnarray}$$

where $\boldsymbol{e}_{i}$ is the $i$ th elementary basis vector. This is similar to direct-search methods, but notice that implicit filtering is not polling opportunistically: all polling points are evaluated on each iteration.

The basic version of implicit filtering from Kelley (Reference Kelley2011) is outlined in Algorithm 4. Note that most implementations of implicit filtering require a bound-constrained domain.

Considerable effort has been devoted to extensions of Algorithm 4when $f$ is ‘noisy’. Gilmore and Kelley (Reference Gilmore and Kelley1995) show that implicit filtering converges to local minima of (DET) when the objective $f$ is the sum of a smooth function $f_{s}$ and a high-frequency, low-amplitude function $f_{n}$ , with $f_{n}\rightarrow 0$ quickly in a neighbourhood of all minimizers of $f_{s}$ . Under similar assumptions, Choi and Kelley (Reference Choi and Kelley2000) show that Algorithm 4converges superlinearly if the step sizes $\unicode[STIX]{x1D6E5}_{k}$ are defined as a power of the norm of the previous iteration’s gradient approximation.

2.3.3 Adaptive regularized methods

Cartis, Gould and Toint (Reference Cartis, Gould and Toint2012) perform an analysis of adaptive regularized cubic (ARC) methods and propose a derivative-free method, ARC-DFO. ARC-DFO is an extension of ARC whereby gradients are replaced with central finite differences of the form (2.29), with the difference parameter monotonically decreasing within a single iteration of the method. ARC-DFO is an intrinsically model-based method akin to Algorithm 3, but the objective within each subproblem regularizes third-order behaviour of the model. Thus, like a trust-region method, ARC-DFO employs trial steps and model gradients. During the main loop of ARC-DFO, if the difference parameter exceeds a constant factor of the minimum of the trial step norm or the model gradient norm, then the difference parameter is shrunk by a constant factor, and the iteration restarts to obtain a new trial step. This mechanism is structurally similar to a derivative-free trust-region method’s checks on model quality. Cartis et al. (Reference Cartis, Gould and Toint2012) show that ARC-DFO demonstrates a WCC result of type (2.2) in $O(\unicode[STIX]{x1D716}^{-3/2})$ , the same asymptotic result (in terms of $\unicode[STIX]{x1D716}$ -dependence) that the authors demonstrate for derivative-based variants of ARC methods. In terms of dependence on $\unicode[STIX]{x1D716}$ , this result is a strict improvement over the WCC results of the same type demonstrated for DDS and trust-region methods, although this result is proved under the stronger assumption that $f\in {\mathcal{L}}{\mathcal{C}}^{2}$ .

In a different approach, Hare and Lucet (Reference Hare and Lucet2013) show convergence of a derivative-free method that penalizes large steps via a proximal regularizer, thereby removing the necessity for a trust region. Lazar and Jarre (Reference Lazar and Jarre2016) regularize their line-search with a term seeking to minimize a weighted change of the model’s third derivatives.

2.3.4 Line-search-based methods

Several line-search-based methods for derivative-free optimization have been developed. Grippo, Lampariello and Lucidi (Reference Grippo, Lampariello and Lucidi1988) and De Leone, Gaudioso and Grippo (Reference De Leone, Gaudioso and Grippo1984) (two of the few papers appearing in the 1980s concerning derivative-free optimization) both analyse conditions on the step sizes used in a derivative-free line-search algorithm, and provide methods for constructing such steps. Lucidi and Sciandrone (Reference Lucidi and Sciandrone2002b ) present methods that combine pattern-search and line-search approaches in a convergent framework. The VXQR method of Neumaier, Fendl, Schilly and Leitner (Reference Neumaier, Fendl, Schilly and Leitner2011) performs a line search on a direction computed from a QR factorization of previously evaluated points. Neumaier et al. (Reference Neumaier, Fendl, Schilly and Leitner2011) apply VXQR to problems with $n=1000$ , a large problem dimension among the methods considered here.

Consideration has also been given to non-monotone line-search-based derivative-free methods. Since gradients are not available in derivative-free optimization, the search direction in a line-search method may not be a descent direction. Non-monotone methods allow one to still employ such directions in a globally convergent framework. Grippo and Sciandrone (Reference Grippo and Sciandrone2007) extend line-search strategies based on coordinate search and the method of Barzilai and Borwein (Reference Barzilai and Borwein1988) to develop a globally convergent non-monotone derivative-free method. Grippo and Rinaldi (Reference Grippo and Rinaldi2014) extend such non-monotone strategies to broader classes of algorithms that employ simplex gradients, hence further unifying direct-search and model-based methods. Another non-monotone line-search method is proposed by Diniz-Ehrhardt, Martínez and Raydan (Reference Diniz-Ehrhardt, Martínez and Raydan2008), who encapsulate early examples of randomized DDS methods (Section 3.2).

2.3.5 Methods for non-smooth optimization

In Section 2.1.2, we discuss how MADS handles non-differentiable objective functions by densely sampling directions on a mesh, thereby ensuring that all Clarke directional derivatives are non-negative (i.e. (2.1)). Another early analysis of a DDS method on a class of non-smooth objectives was performed by García-Palomares and Rodríguez (Reference García-Palomares and Rodríguez2002).

Gradient sampling methods are a developing class of algorithms for general non-smooth non-convex optimization; see the recent survey by Burke et al. (Reference Burke, Curtis, Lewis, Overton and Simões2018). These methods attempt to estimate the $\unicode[STIX]{x1D716}$ -subdifferential at a point $\boldsymbol{x}$ by evaluating a random sample of gradients in the neighbourhood of $\boldsymbol{x}$ and constructing the convex hull of these gradients. In a derivative-free setting, the approximation of these gradients is not as immediately obvious in the presence of non-smoothness, but there exist gradient-sampling methods that use finite-difference estimates with specific smoothing techniques (Kiwiel Reference Kiwiel2010).

In another distinct line of research, Bagirov, Karasözen and Sezer (Reference Bagirov, Karasözen and Sezer2007) analyse a derivative-free variant of subgradient descent, where subgradients are approximated via so-called discrete gradients. In Section 5.3, we will further discuss methods for minimizing composite non-smooth objective functions of the form $f=h\circ F$ , where $h$ is non-smooth but a closed-form expression is known and $F$ is assumed smooth. These methods are characterized by their exploitation of the knowledge of $h$ , making them less general than the methods for non-smooth optimization discussed so far.

3.1.1 Pure random search

Pure random search is a natural method to start with for randomized derivative-free optimization. Pure random search is popular for multiple reasons; in particular, it is easy to implement (with few or no user-defined tolerances), and (if the points generated are independent of one another) it exhibits perfect scaling in terms of evaluating $f$ at many points simultaneously.

A pure random-search method is given in Algorithm 5, where points are generated randomly from $\boldsymbol{\unicode[STIX]{x1D6FA}}$ . For example, if $\boldsymbol{\unicode[STIX]{x1D6FA}}=\{\boldsymbol{x}:\boldsymbol{c}(\boldsymbol{x})\leq \mathbf{0},\boldsymbol{l}\leq \boldsymbol{x}\leq \boldsymbol{u}\}$ , line 3of Algorithm 5may involve drawing points uniformly at random from $[\boldsymbol{l},\boldsymbol{u}]$ and checking whether they satisfy $\boldsymbol{c}(\boldsymbol{x})\leq \mathbf{0}$ .

If the procedure for generating points in line 3of Algorithm 5is independent of the function values observed, then the entire set of points used within pure random search can be generated beforehand: we intentionally omit the index $k$ in the statement of Algorithm 5.

Nevertheless, for the sake of analysis, it is useful to consider an ordering of the sequence of random points generated by Algorithm 5.With such a sequence $\{\boldsymbol{x}_{k}\}$ , one can analyse the best points after $N$ evaluations,

If $f_{\ast }$ is the global minimum value, then

$$\begin{eqnarray}\displaystyle \mathbb{P}[f(\hat{\boldsymbol{x}}_{N})\leq f_{\ast }+\unicode[STIX]{x1D716}]=1-\mathop{\prod }_{k=1}^{N}(1-\mathbb{P}[\boldsymbol{x}_{k}\in \boldsymbol{{\mathcal{L}}}_{f_{\ast }+\unicode[STIX]{x1D716}}(f)]), & & \displaystyle \nonumber\end{eqnarray}$$

where $\unicode[STIX]{x1D716}\geq 0$ and $\boldsymbol{{\mathcal{L}}}_{\unicode[STIX]{x1D6FC}}(f)=\{\boldsymbol{x}:f(\boldsymbol{x})\leq \unicode[STIX]{x1D6FC}\}$ . Provided that the procedure used to generate points at line 3of Algorithm 5satisfies

$$\begin{eqnarray}\displaystyle \lim _{N\rightarrow \infty }\mathop{\prod }_{k=1}^{N}(1-\mathbb{P}[\boldsymbol{x}_{k}\in \boldsymbol{{\mathcal{L}}}_{f_{\ast }+\unicode[STIX]{x1D716}}(f)])=0 & & \displaystyle \nonumber\end{eqnarray}$$

for all $\unicode[STIX]{x1D716}>0$ , then $f(\hat{\boldsymbol{x}}_{k})$ converges in probability to $f_{\ast }$ . For example, if each $\boldsymbol{x}_{k}$ is drawn independently and uniformly over $\boldsymbol{\unicode[STIX]{x1D6FA}}$ , then one can calculate the number of evaluations required to ensure that the $\hat{\boldsymbol{x}}_{k}$ returned by Algorithm 5satisfies $\hat{\boldsymbol{x}}_{k}\in \boldsymbol{{\mathcal{L}}}_{f_{\ast }+\unicode[STIX]{x1D716}}$ with probability $p\in (0,1)$ , that is,

$$\begin{eqnarray}\displaystyle N\geq \displaystyle \frac{\log (p)}{\log \biggl(1-\displaystyle \frac{\unicode[STIX]{x1D707}(\boldsymbol{{\mathcal{L}}}_{f_{\ast }+\unicode[STIX]{x1D716}}\bigcap \boldsymbol{\unicode[STIX]{x1D6FA}})}{\unicode[STIX]{x1D707}(\boldsymbol{\unicode[STIX]{x1D6FA}})}\biggr)}, & & \displaystyle \nonumber\end{eqnarray}$$

provided $\unicode[STIX]{x1D707}(\boldsymbol{{\mathcal{L}}}_{f(\boldsymbol{x}_{\ast })+\unicode[STIX]{x1D716}}\cap \boldsymbol{\unicode[STIX]{x1D6FA}})>0$ and $\boldsymbol{\unicode[STIX]{x1D6FA}}$ is measurable.

Random-search methods typically make few assumptions about $f$ ; see Zhigljavsky (Reference Zhigljavsky1991) for further discussion about the convergence of pure random search. Naturally, a method that assumes only that $f$ is measurable on $\boldsymbol{\unicode[STIX]{x1D6FA}}$ is likely to produce function values that converge more slowly to $f_{\ast }$ when applied to an $f\in {\mathcal{C}}^{0}$ than does a method that exploits the continuity of $f$ . Heuristic modifications of random search have sought to improve empirical performance on certain classes of problems, while still maintaining random search’s global optimization property; see, for example, the work of Zabinsky and Smith (Reference Zabinsky and Smith1992) and Patel, Smith and Zabinsky (Reference Patel, Smith and Zabinsky1989).

3.1.2 Nesterov random search

We refer to the method discussed in this section as Nesterov random search because of the seminal article by Nesterov and Spokoiny (Reference Nesterov and Spokoiny2017), but the idea driving this method is much older. A similar method, for instance, is discussed in Polyak (Reference Polyak1987, Chapter 3.4).

The method of Nesterov random search is largely motivated by Gaussian smoothing. In particular, given a covariance (i.e. symmetric positive-definite) matrix $\boldsymbol{B}$ and a smoothing parameter $\unicode[STIX]{x1D707}>0$ , consider the Gaussian smoothed function

$$\begin{eqnarray}\displaystyle f_{\unicode[STIX]{x1D707}}(\boldsymbol{x})=\sqrt{\displaystyle \frac{\text{det}(\boldsymbol{B})}{(2\unicode[STIX]{x1D70B})^{n}}}\int _{\mathbb{R}^{n}}f(\boldsymbol{x}+\unicode[STIX]{x1D707}\boldsymbol{u})\exp \biggl(-\displaystyle \frac{1}{2}\boldsymbol{u}^{\text{T}}\boldsymbol{B}\boldsymbol{u}\biggr)\text{d}\boldsymbol{u}. & & \displaystyle \nonumber\end{eqnarray}$$

This smoothing has many desirable properties; for instance, if $f$ is Lipschitz-continuous with constant $L_{\text{f}}$ , then $f_{\unicode[STIX]{x1D707}}$ is Lipschitz-continuous with a constant no worse than $L_{\text{f}}$ for all $\unicode[STIX]{x1D707}>0$ . Likewise, if $f$ has Lipschitz-continuous gradients with constant $L_{\text{g}}$ , then $f_{\unicode[STIX]{x1D707}}$ has Lipschitz-continuous gradients with a constant no worse than $L_{\text{g}}$ for all $\unicode[STIX]{x1D707}>0$ . If $f$ is convex, then $f_{\unicode[STIX]{x1D707}}$ is convex.

One can show that

$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}f_{\unicode[STIX]{x1D707}}(\boldsymbol{x})=\displaystyle \frac{1}{\unicode[STIX]{x1D707}}\sqrt{\displaystyle \frac{\text{det}(\boldsymbol{B})}{(2\unicode[STIX]{x1D70B})^{n}}}\int _{\mathbb{R}^{n}}(f(\boldsymbol{x}+\unicode[STIX]{x1D707}\boldsymbol{u})-f(\boldsymbol{x}))\exp \biggl(-\displaystyle \frac{1}{2}\boldsymbol{u}^{\text{T}}\boldsymbol{B}\boldsymbol{u}\biggr)\boldsymbol{B}\boldsymbol{u}\text{d}\boldsymbol{u}. & & \displaystyle \nonumber\end{eqnarray}$$

In other words, $\unicode[STIX]{x1D6FB}f_{\unicode[STIX]{x1D707}}(\boldsymbol{x})$ , which can be understood as an approximation of $\unicode[STIX]{x1D6FB}f(\boldsymbol{x})$ in the smooth case, can be computed via an expectation over $\boldsymbol{u}\in \mathbb{R}^{n}$ weighted by the finite difference $f(\boldsymbol{x}+\unicode[STIX]{x1D707}\boldsymbol{u})-f(\boldsymbol{x})$ and inversely weighted by a radial distance from $\boldsymbol{x}$ . With this interpretation in mind, Nesterov and Spokoiny (Reference Nesterov and Spokoiny2017) propose a collection of random gradient-free oracles, where one first generates a Gaussian random vector $\boldsymbol{u}\in {\mathcal{N}}(\mathbf{0},\boldsymbol{B}^{-1})$ and then uses one of

(3.1) $$\begin{eqnarray}\begin{array}{@{}lll@{}}\boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u}) & = & \unicode[STIX]{x1D6FF}_{\text{f}}(f;\boldsymbol{x};\boldsymbol{u};\unicode[STIX]{x1D707})\boldsymbol{B}\boldsymbol{u},\quad \text{or}\\ \hat{\boldsymbol{g}}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u}) & = & \unicode[STIX]{x1D6FF}_{\text{c}}(f;\boldsymbol{x};\boldsymbol{u};\unicode[STIX]{x1D707})\boldsymbol{B}\boldsymbol{u},\end{array}\end{eqnarray}$$

for a difference parameter $\unicode[STIX]{x1D707}>0$ . Nesterov and Spokoiny also propose a third oracle, $\boldsymbol{g}_{0}(\boldsymbol{x};\boldsymbol{u})=f^{\prime }(\boldsymbol{x};\boldsymbol{u})\boldsymbol{B}\boldsymbol{u}$ , intended for the optimization of non-smooth functions; this oracle assumes the ability to compute directional derivatives $f^{\prime }(\boldsymbol{x};\boldsymbol{u})$ . For this reason, Nesterov and Spokoiny refer to all oracles as gradient-free instead of derivative-free. Given the scope of this survey, we focus on the derivative-free oracles $\boldsymbol{g}_{\unicode[STIX]{x1D707}}$ and $\hat{\boldsymbol{g}}_{\unicode[STIX]{x1D707}}$ displayed in (3.1).

With an oracle $\boldsymbol{g}$ chosen as either oracle in (3.1), Nesterov random-search methods are straightforward to define, and we do so in Algorithm 6. In Algorithm 6, $\text{proj}(\cdot ;\boldsymbol{\unicode[STIX]{x1D6FA}})$ denotes projection onto a domain $\boldsymbol{\unicode[STIX]{x1D6FA}}$ .

A particularly striking result proved in Nesterov and Spokoiny (Reference Nesterov and Spokoiny2017) was perhaps the first WCC result for an algorithm (Algorithm 6) in the case where $f\in {\mathcal{L}}{\mathcal{C}}^{0}$ – that is, $f$ may be both non-smooth and non-convex. Because of the randomized nature of iteratively sampling from a Gaussian distribution in Algorithm 6, complexity results are given as expectations. That is, letting $\boldsymbol{U}_{k}=\{\boldsymbol{u}_{0},\boldsymbol{u}_{1},\ldots ,\boldsymbol{u}_{k}\}$ denote the random variables associated with the first $k$ iterations of Algorithm 6, complexity results are stated in terms of expectations with respect to the filtration defined by these variables. A WCC is given as an upper bound on the number of $f$ evaluations needed to attain the approximate ( $\unicode[STIX]{x1D716}>0$ ) optimality condition

(3.2) $$\begin{eqnarray}\mathbb{E}_{\boldsymbol{U}_{k-1}}[\Vert \unicode[STIX]{x1D6FB}f_{\check{\unicode[STIX]{x1D707}}}(\hat{\boldsymbol{x}}_{k})\Vert ]\leq \unicode[STIX]{x1D716},\end{eqnarray}$$

where $\hat{\boldsymbol{x}}_{k}=\text{arg}\,\text{min}_{j=0,1,\ldots ,k-1}f(\boldsymbol{x}_{j})$ . By fixing a particular choice of $\check{\unicode[STIX]{x1D707}}$ (dependent on $\unicode[STIX]{x1D716}$ , $n$ and Lipschitz constants), Nesterov and Spokoiny (Reference Nesterov and Spokoiny2017) demonstrate that the number of $f$ evaluations needed to attain (3.2) is in $O(\unicode[STIX]{x1D716}^{-3})$ ; see Table A.1. For $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ (but still non-convex), Nesterov and Spokoiny (Reference Nesterov and Spokoiny2017) prove a WCC result of type (3.2) in $O(\unicode[STIX]{x1D716}^{-2})$ for the same method. WCC results of Algorithm 6under a variety of stronger assumptions on the convexity and differentiability of $f$ are also shown in Table A.1 and discussed in Section 4. We further note that some randomized methods of the form Algorithm 6have also been developed for (STOCH), which we discuss in Section 6.

We remark on an undesirable feature of the convergence analysis for variants of Algorithm 6: the analysis of these methods supposes that the sequence $\{\unicode[STIX]{x1D6FC}_{k}\}$ is chosen as a constant that depends on parameters, including $L_{\text{f}}$ , that may not be available to the method. Similar assumptions concerning the preselection of $\{\unicode[STIX]{x1D6FC}_{k}\}$ also appear in the convex cases discussed in Section 4, and we highlight these dependencies in Table A.1.

4.2.1 One-point bandit feedback

In one-point bandit feedback, a method is assumed to have access to an oracle $\tilde{f}$ that returns unbiased estimates of $f$ . In particular, given two points $\boldsymbol{x}_{1}$ , $\boldsymbol{x}_{2}\in \boldsymbol{\unicode[STIX]{x1D6FA}}$ , two separate calls to the oracle will return $\tilde{f}(\boldsymbol{x}_{1};\boldsymbol{\unicode[STIX]{x1D709}}_{1})$ and $\tilde{f}(\boldsymbol{x}_{2};\boldsymbol{\unicode[STIX]{x1D709}}_{2})$ ; methods do not have control over the selection of the random variables $\boldsymbol{\unicode[STIX]{x1D709}}_{1}$ and $\boldsymbol{\unicode[STIX]{x1D709}}_{2}$ . Many one-point bandit methods do not fall neatly into the frameworks discussed in this survey (Agarwal et al. Reference Agarwal, Foster, Hsu, Kakade, Rakhlin and Shawe-Taylor2011, Belloni, Liang, Narayanan and Rakhlin Reference Belloni, Liang, Narayanan and Rakhlin2015, Bubeck, Lee and Eldan Reference Bubeck, Lee and Eldan2017). See Table 4.1 for a summary of best known WCC results of type (4.9) for one-point bandit feedback methods.

One example of a method using one-point bandit feedback, whose development falls naturally into our discussion thus far, is given by Flaxman, Kalai and McMahan (Reference Flaxman, Kalai and McMahan2005); they analyse a method resembling Algorithm 6, but the gradient-free oracle is chosen as

(4.10) $$\begin{eqnarray}\boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})=\displaystyle \frac{\tilde{f}(\boldsymbol{x}+\unicode[STIX]{x1D707}\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})}{\unicode[STIX]{x1D707}}\boldsymbol{u},\end{eqnarray}$$

where $\boldsymbol{u}$ is drawn uniformly from the unit $n$ -dimensional sphere. That is, given a realization $\boldsymbol{\unicode[STIX]{x1D709}}\in \boldsymbol{\unicode[STIX]{x1D6EF}}$ , a stochastic gradient estimator based on a stochastic evaluation at the point $\boldsymbol{x}+\unicode[STIX]{x1D707}\boldsymbol{u}$ is computed via (4.10). For this method, and for general convex $f$ , Flaxman et al. (Reference Flaxman, Kalai and McMahan2005) demonstrate a WCC bound of type (4.9) in $O(n^{2}\unicode[STIX]{x1D716}^{-4})$ for general convex $f$ . For strongly convex $f$ , Flaxman et al. (Reference Flaxman, Kalai and McMahan2005) demonstrate a WCC bound of type (4.9) in $O(n^{2}\unicode[STIX]{x1D716}^{-3})$ . Various extensions of these results are given in Saha and Tewari (Reference Saha and Tewari2011), Dekel, Eldan and Koren (Reference Dekel, Eldan, Koren and Cortes2015) and Gasnikov et al. (Reference Gasnikov, Krymova, Lagunovskaya, Usmanova and Fedorenko2017). To the best of our knowledge, the best known upper bound on the WCC for smooth, strongly convex problems is in $O(n^{2}\unicode[STIX]{x1D716}^{-2})$ and the best known upper bound on the WCC for smooth, convex problems is in $O(n\unicode[STIX]{x1D716}^{-3})$ (Gasnikov et al. Reference Gasnikov, Krymova, Lagunovskaya, Usmanova and Fedorenko2017).

Table 4.1.

Best known WCCs for $N_{\unicode[STIX]{x1D716}}$ , the number of evaluations required to bound (4.9), for one-point bandit feedback. $N_{\unicode[STIX]{x1D716}}$ is given only in terms of $n$ and $\unicode[STIX]{x1D716}$ . See text for the definition of $\unicode[STIX]{x1D6FD}$ -smooth; here $\unicode[STIX]{x1D6FD}>2$ . Method types include random search (RS), mirror descent (MD) and ellipsoidal.

For the solution of (STOCH) when $\boldsymbol{\unicode[STIX]{x1D6FA}}=\mathbb{R}^{n}$ , Bach and Perchet (Reference Bach, Perchet and Feldman2016) analyse a method resembling Algorithm 6wherein the gradient-free oracle is chosen as

(4.11) $$\begin{eqnarray}\boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}}_{+},\boldsymbol{\unicode[STIX]{x1D709}}_{-})=\displaystyle \frac{\tilde{f}(\boldsymbol{x}+\unicode[STIX]{x1D707}\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}}_{+})-\tilde{f}(\boldsymbol{x}-\unicode[STIX]{x1D707}\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}}_{-})}{\unicode[STIX]{x1D707}}\boldsymbol{u},\end{eqnarray}$$

where $\boldsymbol{u}$ is again drawn uniformly from the unit $n$ -dimensional sphere. We remark that $\boldsymbol{\unicode[STIX]{x1D709}}_{+},\boldsymbol{\unicode[STIX]{x1D709}}_{-}\in \boldsymbol{\unicode[STIX]{x1D6EF}}$ in (4.11) are different realizations; this distinction will become particularly relevant in Section 4.2.2. For the solution of (STOCH) when $\boldsymbol{\unicode[STIX]{x1D6FA}}\subsetneq \mathbb{R}^{n}$ , Bach and Perchet (Reference Bach, Perchet and Feldman2016) also consider a gradient-free oracle like (4.10). Bach and Perchet (Reference Bach, Perchet and Feldman2016) are particularly interested in how the smoothness of $f$ can be exploited to improve complexity bounds in the one-point bandit feedback paradigm; they define a parameter $\unicode[STIX]{x1D6FD}$ and say that a function $f$ is $\unicode[STIX]{x1D6FD}$ -smooth provided $f$ is almost everywhere $(\unicode[STIX]{x1D6FD}-1)$ -times Lipschitz-continuously differentiable (a strictly weaker condition than assuming $f\in {\mathcal{L}}{\mathcal{C}}^{\unicode[STIX]{x1D6FD}-1}$ ). When $\unicode[STIX]{x1D6FD}=2$ , Bach and Perchet (Reference Bach, Perchet and Feldman2016) recover similar results to those seen in Table 4.1 for when $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ . When $\unicode[STIX]{x1D6FD}>2$ , however, in both the constrained and unconstrained case, Bach and Perchet (Reference Bach, Perchet and Feldman2016) prove a WCC result of type (4.9) in $O(n^{2}\unicode[STIX]{x1D716}^{2\unicode[STIX]{x1D6FD}/(1-\unicode[STIX]{x1D6FD})})$ when $f$ is convex and $\unicode[STIX]{x1D6FD}$ -smooth. Further assuming that $f$ is $c$ -strongly convex, Bach and Perchet (Reference Bach, Perchet and Feldman2016) prove a WCC result of type (4.9) in $O(n^{2}\unicode[STIX]{x1D716}^{(\unicode[STIX]{x1D6FD}+1)/(1-\unicode[STIX]{x1D6FD})}/c)$ . Notice that asymptotically, as $\unicode[STIX]{x1D6FD}\rightarrow \infty$ , this bound is in $O(n^{2}\unicode[STIX]{x1D716}^{-1}/c)$ . This asymptotic result is particularly important because it attains the lower bound on optimization error demonstrated by Shamir (Reference Shamir, Shalev-Shwartz and Steinwart2013) for strongly convex $\infty$ -smooth (quadratic) $f$ .

We also note that Bubeck et al. (Reference Bubeck, Lee and Eldan2017) conjecture that a particular kernel method can achieve a WCC of type (4.9) in $O(n^{3}\unicode[STIX]{x1D716}^{-2})$ for general convex functions. In light of well-known results in deterministic convex optimization, the WCCs summarized in Table 4.1 may be surprising. In particular, for any $c$ -strongly convex function $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ , the best known WCC results are in $O(\unicode[STIX]{x1D716}^{-2})$ . We place particular emphasis on this result because it illustrates a gap between derivative-free and derivative-based optimization that is not just a factor of $n$ . In this particular one-point bandit feedback setting, there do not seem to exist methods that achieve the optimalFootnote ⁹ $O(\unicode[STIX]{x1D716}^{-1})$ convergence rate attainable by gradient-based methods for smooth strongly convex stochastic optimization. Hu, Prashanth, György and Szepesvári (Reference Hu, Prashanth, György and Szepesvári2016) partially address this issue concerning one-point bandit feedback, which they refer to as ‘uncontrolled noise’. These observations motivated the study of two-point (multi-point) bandit feedback, which we will discuss in the next section, Section 4.2.2.

We further remark that every WCC in Table 4.1 has a polynomial dependence on the dimension $n$ , raising natural questions about the applicability of these methods in high-dimensional settings. Wang, Du, Balakrishnan and Singh (Reference Wang, Du, Balakrishnan, Singh, Storkey and Perez-Cruz2018) consider a mirror descent method employing a special gradient-free oracle computed via a compressed sensing technique. They prove a WCC of type (4.9) in $O(\log (d)^{3/2}n_{z}^{2}\unicode[STIX]{x1D716}^{-3})$ , under additional assumptions on derivative sparsity, most importantly, that for all $\boldsymbol{x}\in \boldsymbol{\unicode[STIX]{x1D6FA}}$ , $\Vert \unicode[STIX]{x1D6FB}f(\boldsymbol{x})\Vert _{0}\leq n_{z}$ for some $n_{z}$ . Thus, provided $n_{z}\ll n$ , the polynomial dependence on $n$ becomes a logarithmic dependence on $n$ , at the expense of a WCC with a strictly worse dependence on $\unicode[STIX]{x1D716}$ than $\unicode[STIX]{x1D716}^{-2}$ . In Section 5.2, we discuss methods that similarly exploit known sparsity of objective function derivatives.

In an extension of Algorithm 6, Chen (Reference Chen2015, Chapter 3) dynamically updates the difference parameter by exploiting knowledge of the changing variance of $\boldsymbol{\unicode[STIX]{x1D709}}$ .

4.2.2 Two-point (multi-point) bandit feedback

We now focus on the stochastic paradigm of two-point (or multi-point) bandit feedback. In this setting of bandit feedback, we do not encounter the same gaps in WCC results between derivative-free and derivative-based optimization that exist for one-point bandit feedback.

The majority of methods analysed in the two-point bandit feedback setting are essentially random-search methods of the form Algorithm 6. The gradient-free oracles from (3.1) in the two-point setting takes one of the two forms

(4.12) $$\begin{eqnarray}\begin{array}{@{}l@{}}\boldsymbol{g}_{\unicode[STIX]{x1D707}_{k}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})=\unicode[STIX]{x1D6FF}_{\text{f}}(\tilde{f}(\cdot ,\boldsymbol{\unicode[STIX]{x1D709}});\boldsymbol{x};\boldsymbol{u};\unicode[STIX]{x1D707}_{k})\boldsymbol{B}\boldsymbol{u},\quad \text{or}\\ \hat{\boldsymbol{g}}_{\unicode[STIX]{x1D707}_{k}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})=\unicode[STIX]{x1D6FF}_{\text{c}}(\tilde{f}(\cdot ,\boldsymbol{\unicode[STIX]{x1D709}});\boldsymbol{x};\boldsymbol{u};\unicode[STIX]{x1D707}_{k})\boldsymbol{B}\boldsymbol{u}.\end{array}\end{eqnarray}$$

The key observation in (4.12) is that $\boldsymbol{\unicode[STIX]{x1D709}}$ denotes a single realization used in the computation of both function values in the definitions of the oracles (see (2.28) and (2.29)). This assumption of ‘controllable realizations’ separates two-point bandit feedback from the more pessimistic one-point bandit feedback. This property of being able to recall a single realization $\boldsymbol{\unicode[STIX]{x1D709}}$ for two (or more) evaluations of $\tilde{f}$ is precisely why this setting of bandit feedback is called ‘two-point’ (or ‘multi-point’). This property will also be exploited in Section 6.1.

The early work of Agarwal, Dekel and Xiao (Reference Agarwal, Dekel, Xiao, Kalai and Mohri2010) directly addresses the discussed gap in WCC results and demonstrates that a random-search method resembling Algorithm 6, but applied in the two-point (or multi-point) setting as opposed to the one-point setting, attains a WCC of type (4.9) in $O(\unicode[STIX]{x1D716}^{-1})$ ; this is a complexity result matching the optimal rate (in terms of $\unicode[STIX]{x1D716}$ -dependence) shown by Agarwal et al. (Reference Agarwal, Wainwright, Bartlett, Ravikumar and Bengio2009). See Table 4.2 for a summary of best known WCC results of type (4.9) for two-point bandit feedback methods.

Nesterov and Spokoiny (Reference Nesterov and Spokoiny2017) provide a WCC result for Algorithm 6using the stochastic gradient-free oracles (4.12), but strictly better WCCs have since been established. We also note that, in contrast with the gradient-free oracles of (3.1) in Section 3.1.2, the difference parameter $\unicode[STIX]{x1D707}$ is written as $\unicode[STIX]{x1D707}_{k}$ in (4.12), indicating that a selection for $\unicode[STIX]{x1D707}_{k}$ must be made in the $k$ th iteration. In the works that we discuss here, $\unicode[STIX]{x1D707}_{k}$ in (4.12) is either chosen as a constant sufficiently small or else $\unicode[STIX]{x1D707}_{k}\rightarrow 0$ at a rate typically of the order of $1/k$ . We also remark that many of the results discussed in this section trivially hold for deterministic convex problems, and can be seen as an extension of results concerning methods of the form of Algorithm 6.

Provided $f\in {\mathcal{L}}{\mathcal{C}}^{0}$ , the best known WCCs of type (4.9) (with variations in constants) can be found in Duchi, Jordan, Wainwright and Wibisono (Reference Duchi, Jordan, Wainwright and Wibisono2015), Gasnikov et al. (Reference Gasnikov, Krymova, Lagunovskaya, Usmanova and Fedorenko2017) and Shamir (Reference Shamir2017). These works consider variants of mirror-descent methods with approximate gradients given by estimators of the form (4.12); see Algorithm 7for a description of a basic mirror-descent method. Algorithm 7depends on the concept of a Bregman divergence $D_{\unicode[STIX]{x1D713}}$ , used in line 6of Algorithm 7to define a proximal-point subproblem. The Bregman divergence is defined by a function $\unicode[STIX]{x1D713}~:~\boldsymbol{\unicode[STIX]{x1D6FA}}\rightarrow \mathbb{R}$ , which is assumed to be $1$ -strongly convex with respect to the norm $\Vert \cdot \Vert _{p}$ . To summarize many findings in this area (Duchi et al. Reference Duchi, Jordan, Wainwright and Wibisono2015, Gasnikov et al. Reference Gasnikov, Krymova, Lagunovskaya, Usmanova and Fedorenko2017, Shamir Reference Shamir2017), if $\Vert \cdot \Vert _{q}$ is the dual norm to $\Vert \cdot \Vert _{p}$ (i.e.  $p^{-1}+q^{-1}=1$ ) where $p\in \{1,2\}$ and $R_{p}$ denotes the radius of the feasible set in the $\Vert \cdot \Vert _{p}$ -norm, then a bound on WCC of type (4.9) in $O(n^{2/q}R_{p}^{2}\unicode[STIX]{x1D716}^{-2})$ can be established for a method like Algorithm 7in the two-point feedback setting where $f\in {\mathcal{L}}{\mathcal{C}}^{0}$ . These WCCs for methods like Algorithm 7are responsible for the popularity of mirror-descent methods in machine learning. For many machine learning problems, solutions are typically sparse, and so, in some sense, $R_{1}\leq R_{2}$ . Thus, using a function $\unicode[STIX]{x1D713}$ that is $1$ -strongly convex with respect to the $\Vert \cdot \Vert _{1}$ -norm (e.g. simply letting $\unicode[STIX]{x1D713}=\Vert \cdot \Vert _{1}$ ) may be preferable to $p=2$ , in both theory and practice.

Duchi et al. (Reference Duchi, Jordan, Wainwright and Wibisono2015) also provide an information-theoretic lower bound on convergence rates for any method in the (non-strongly) convex, $f\in {\mathcal{L}}{\mathcal{C}}^{0}$ , two-point feedback setting. This bound matches the best known WCCs up to constants, demonstrating that these results are tight. This lower bound is still of the order of $\unicode[STIX]{x1D716}^{-2}$ , matching the result of Agarwal et al. (Reference Agarwal, Wainwright, Bartlett, Ravikumar and Bengio2009) in $\unicode[STIX]{x1D716}$ -dependence in the case where $f\in {\mathcal{L}}{\mathcal{C}}^{0}$ . It is also remarkable that this result is only a factor of $\sqrt{n}$ worse than the bounds provided by Agarwal et al. (Reference Agarwal, Wainwright, Bartlett, Ravikumar and Bengio2009) for the derivative-based case, as opposed to the factor of $n$ that one may expect.

Table 4.2.

Best known WCCs of type (4.9) for two-point bandit feedback. $N_{\unicode[STIX]{x1D716}}$ is given only in terms of $n$ and $\unicode[STIX]{x1D716}$ . See text for the definition of $p,q$ . $R_{p}$ denotes the size of the feasible set in the $\Vert \cdot \Vert _{p}$ -norm. If $f$ is $c_{p}$ -strongly convex, then $f$ is strongly convex with respect to the $\Vert \cdot \Vert _{p}$ -norm with constant $c_{p}$ . $\unicode[STIX]{x1D70E}$ is the standard deviation on the gradient estimator $\boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})$ (i.e.  $\mathbb{E}_{\boldsymbol{\unicode[STIX]{x1D709}}}[\Vert \boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})-\unicode[STIX]{x1D6FB}f(\boldsymbol{x})\Vert ^{2}]\leq \unicode[STIX]{x1D70E}^{2}$ ). The Lipschitz constant of the gradient $L_{\text{g}}$ is defined by the $\Vert \cdot \Vert _{2}$ -norm. $\star$ denotes the additional assumption that $\mathbb{E}_{\boldsymbol{\unicode[STIX]{x1D709}}}[\Vert \boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})\Vert ]<\infty$ . Method types include random search (RS), mirror descent (MD) and accelerated mirror descent (AMD).

Additionally assuming $f(\boldsymbol{x})$ is strongly convex (but still assuming $f\in {\mathcal{L}}{\mathcal{C}}^{0}$ ), Agarwal et al. (Reference Agarwal, Dekel, Xiao, Kalai and Mohri2010) prove, for a method like Algorithm 6, a WCC in $O(n^{2}\unicode[STIX]{x1D716}^{-1})$ . Using a method like Algorithm 7, Gasnikov et al. (Reference Gasnikov, Krymova, Lagunovskaya, Usmanova and Fedorenko2017) improve the dependence on $n$ in this WCC to $O(n^{2/q}c_{p}^{-1}\unicode[STIX]{x1D716}^{-1})$ , provided $f$ is $c_{p}$ -strongly convex with respect to $\Vert \cdot \Vert _{p}$ .

We now address the case where $f$ is convex and $f\in {\mathcal{L}}{\mathcal{C}}^{1}$ . Given an assumption that $\Vert \boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})\Vert$ is uniformly bounded, Agarwal et al. (Reference Agarwal, Dekel, Xiao, Kalai and Mohri2010) demonstrate a WCC of type (4.9) in $O(n^{2}\unicode[STIX]{x1D716}^{-1})$ . Dropping this somewhat restrictive assumption on the gradient-free oracle and assuming instead that the oracle used in a method like Algorithm 6has bounded variance (i.e. the oracle satisfies $\mathbb{E}_{\boldsymbol{\unicode[STIX]{x1D709}}}[\Vert \boldsymbol{g}_{\unicode[STIX]{x1D707}}(\boldsymbol{x};\boldsymbol{u};\boldsymbol{\unicode[STIX]{x1D709}})-\unicode[STIX]{x1D6FB}f(\boldsymbol{x})\Vert ^{2}]\leq \unicode[STIX]{x1D70E}^{2}$ ), Ghadimi and Lan (Reference Ghadimi and Lan2013) prove a WCC of a type similar to (4.9) (we avoid a discussion of randomized stopping) in $O(\max \{nL_{\text{g}}\Vert \boldsymbol{x}_{0}-\boldsymbol{x}_{\ast }\Vert \unicode[STIX]{x1D716}^{-1},n\unicode[STIX]{x1D70E}^{2}L_{g}\Vert \boldsymbol{x}_{0}-\boldsymbol{x}_{\ast }\Vert \unicode[STIX]{x1D716}^{-2}\})$ . We mention that Nesterov and Spokoiny (Reference Nesterov and Spokoiny2017) hinted at a similar WCC result, but with a strictly worse dependence on $n$ , and different assumptions on $\boldsymbol{\unicode[STIX]{x1D709}}$ .

5.3.1 Convex $h$

When $h$ in (5.5) is convex (note that $f$ may still be non-convex due to non-convexity in $\boldsymbol{F}$ ), one thrust of research extends the techniques in derivative-based composite non-smooth optimization; see the works of Yuan (Reference Yuan1985) and Fletcher (Reference Fletcher1987, Chapter 14). For example, Yuan (Reference Yuan1985) use derivatives to construct convex first-order approximations of $f$ near $\boldsymbol{x}$ ,

(5.6) $$\begin{eqnarray}\ell (\boldsymbol{x}+\boldsymbol{s})=h(\boldsymbol{F}(\boldsymbol{x})+\unicode[STIX]{x1D6FB}\boldsymbol{F}(\boldsymbol{x})\boldsymbol{s}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FB}\boldsymbol{F}$ denotes the Jacobian of $\boldsymbol{F}$ ; see Hare (Reference Hare2017) for properties of such approximations in the derivative-free setting. By replacing $\unicode[STIX]{x1D6FB}\boldsymbol{F}$ in (5.6) with the matrix $\boldsymbol{M}(\boldsymbol{x}_{k})$ containing the gradients of a fully linear approximation to $\boldsymbol{F}$ at $\boldsymbol{x}_{k}$ , Grapiglia, Yuan and Yuan (Reference Grapiglia, Yuan and Yuan2016) and Garmanjani et al. (Reference Garmanjani, Jùdice and Vicente2016) independently analyse a model-based trust-region method similar to Algorithm 3from Section 2.2.4 that uses the non-smooth trust-region subproblem

(5.7)

Note that only $\boldsymbol{F}$ is assumed to be a black-box; these methods exploit the fact that $h$ is convex with a known form in order to appropriately solve (5.6). Both Grapiglia et al. (Reference Grapiglia, Yuan and Yuan2016) and Garmanjani et al. (Reference Garmanjani, Jùdice and Vicente2016) use the stationarity measure of Yuan (Reference Yuan1985) in their analysis,

(5.8) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}(\boldsymbol{x})=\ell (\boldsymbol{x})-\min _{\Vert \boldsymbol{s}\Vert \leq 1}\ell (\boldsymbol{x}+\boldsymbol{s}),\end{eqnarray}$$

for which it is known that $\unicode[STIX]{x1D6F9}(\boldsymbol{x}_{\ast })=0$ if and only if $\boldsymbol{x}_{\ast }$ is a critical point of $f$ in the sense that $\ell (\boldsymbol{x}_{\ast })\leq \ell (\boldsymbol{x}_{\ast }+\boldsymbol{s})$ for all $\boldsymbol{s}\in \mathbb{R}^{n}$ . Worst-case complexity results that bound the effort required to attain $\Vert \unicode[STIX]{x1D6F9}(\boldsymbol{x}_{k})\Vert \leq \unicode[STIX]{x1D716}$ are included in Table A.1.

Methods using the local approximation (5.6) require convexity of $h$ in order for $\unicode[STIX]{x1D6F9}$ in (5.8) to be interpreted as a stationarity measure. From a practical perspective, the form of $h$ directly affects the difficulty of solving the trust-region subproblem. Grapiglia et al. (Reference Grapiglia, Yuan and Yuan2016) demonstrate this approach on a collection of problems of the form

(5.9) $$\begin{eqnarray}f(\boldsymbol{x})=\max _{i=1,\ldots ,p}F_{i}(\boldsymbol{x}),\end{eqnarray}$$

where each $F_{i}$ is assumed smooth. Garmanjani et al. (Reference Garmanjani, Jùdice and Vicente2016) test their method on a collection of problems of the form

(5.10) $$\begin{eqnarray}f(\boldsymbol{x})=\Vert \boldsymbol{F}(\boldsymbol{x})\Vert _{1}=\mathop{\sum }_{i=1}^{p}|F_{i}(\boldsymbol{x})|.\end{eqnarray}$$

For objectives of the form (5.9) or (5.10), the associated trust-region subproblems (5.7) can be cast as linear programs when the $\infty$ -norm defines the trust region. (An early example of such an approach appears in Madsen (Reference Madsen, Balinski and Wolfe1975), where linear approximations to each $F_{i}$ in (5.9) are constructed.) Although more general convex $h$ could fit into this framework, one must be wary of the difficulty of the resulting subproblems.

Direct-search methods have also been adapted for composite non-smooth functions of specific forms. In these variants, knowledge of $h$ in (5.5) informs the selection of search directions in a manner similar to that described in Section 5.2. Ma and Zhang (Reference Ma and Zhang2009) and Bogani, Gasparo and Papini (Reference Bogani, Gasparo and Papini2009) consider the cases of $f$ of the form (5.9) and (5.10), respectively.

Objectives of the form (5.9) are also addressed by Hare, Planiden and Sagastizábal (Reference Hare, Planiden and Sagastizábal2019), who develop an algorithm that decomposes such problems into orthogonal subspaces associated with directions of non-smoothness and directions of smoothness. The resulting derivative-free $VU$ -algorithm employs model-based estimates of gradients to form and update this decomposition (Hare Reference Hare2014). Liuzzi, Lucidi and Sciandrone (Reference Liuzzi, Lucidi and Sciandrone2006) address finite minimax problems by converting the original problem into a smooth problem using an exponential penalty function. Their DDS method adjusts the penalty parameter via a rule that depends on the current step size in order to guarantee convergence to a Clarke stationary point.

Approximate gradient-sampling methods are developed and analysed by Hare and Macklem (Reference Hare and Macklem2013) and Hare and Nutini (Reference Hare and Nutini2013) for the finite minimax problem (5.9). These methods effectively exploit the subdifferential structure of $h(\boldsymbol{y})=\max _{i=1,\ldots ,p}y_{i}$ and employ derivative-free approximations of each $\unicode[STIX]{x1D6FB}F_{i}(\boldsymbol{x})$ . Larson, Menickelly and Wild (Reference Larson, Menickelly and Wild2016) propose a variant of gradient sampling, called manifold sampling, for objectives of the form (5.10). Unlike (approximate) gradient sampling, manifold sampling does not depend on a random sample of points to estimate the $\unicode[STIX]{x1D716}$ -subdifferential.

5.3.2 Non-convex $h$

When $h$ is non-convex, minimization of (5.5) is considerably more challenging than when $h$ is convex. Few methods exist that exploit the structure of non-convex $h$ . One of the many challenges is that the model in (5.6) may no longer be an underestimator of $h$ . Khan et al. (Reference Khan, Larson and Wild2018) propose a manifold sampling method for piecewise linear $h$ ; in contrast to the previously discussed methods, this method does not require that $h$ be convex. Other methods applicable for non-convex $h$ employ smoothing functions.

As mentioned in Section 2.1.2, the worst-case complexity of DDS methods applied to non-smooth (Lipschitz-continuous) objective functions is difficult to analyse. The reason that DDS methods generate an asymptotically dense set of polling directions is to ensure that no descent directions exist. An exception to this generality, however, is functions for which an appropriate smoothing function exists. Given a locally Lipschitz-continuous $f$ , we say that $f_{\unicode[STIX]{x1D707}}:\mathbb{R}^{n}\rightarrow \mathbb{R}$ is a smoothing function for $f$ provided that for any $\unicode[STIX]{x1D707}\in (0,\infty )$ , $f_{\unicode[STIX]{x1D707}}$ is continuously differentiable and that

$$\begin{eqnarray}\displaystyle \lim _{\boldsymbol{z}\rightarrow \boldsymbol{x},\unicode[STIX]{x1D707}\rightarrow 0^{+}}f_{\unicode[STIX]{x1D707}}(\boldsymbol{z})=f(\boldsymbol{x}), & & \displaystyle \nonumber\end{eqnarray}$$

for all $\boldsymbol{x}\in \mathbb{R}^{n}$ .

Thus, if a smoothing function $f_{\unicode[STIX]{x1D707}}$ exists for $f$ , it is natural to iteratively apply a method for smooth unconstrained optimization to obtain approximate solutions $\boldsymbol{x}_{k}$ to $\min _{\boldsymbol{x}}f_{\unicode[STIX]{x1D707}_{k}}(\boldsymbol{x})$ while decreasing $\unicode[STIX]{x1D707}_{k}$ . We roughly prescribe such a smoothing method in Algorithm 8.

Garmanjani and Vicente (Reference Garmanjani and Vicente2012) consider the DDS framework analysed by Vicente (Reference Vicente2013) as the method $\mathfrak{M}$ in Algorithm 8.They terminate $\mathfrak{M}$ when the step-size parameter $\unicode[STIX]{x1D6FC}$ of Algorithm 2is sufficiently small, where the notion of sufficiently small scales with $\unicode[STIX]{x1D707}_{k}$ in Algorithm 8.Garmanjani and Vicente (Reference Garmanjani and Vicente2012) prove a first-order stationarity result of the form

(5.11) $$\begin{eqnarray}\liminf _{k\rightarrow \infty }\Vert \unicode[STIX]{x1D6FB}f_{\unicode[STIX]{x1D707}_{k}}(\boldsymbol{x}_{k})\Vert =0.\end{eqnarray}$$

Under certain assumptions (for instance, that $h$ satisfies some regularity conditions at $\boldsymbol{F}(\boldsymbol{x}_{\ast })$ ) this first-order stationarity result is equivalent to $0\in \unicode[STIX]{x2202}f(\boldsymbol{x}_{\ast })$ ; that is, $\boldsymbol{x}_{\ast }$ is Clarke stationary.

Garmanjani and Vicente (Reference Garmanjani and Vicente2012) consider the decrease to be sufficient in line 4of Algorithm 1if $f(\boldsymbol{p}_{i})<f(\boldsymbol{x})-c_{1}\unicode[STIX]{x1D6FC}^{3/2}$ for some $c_{1}>0$ . If, furthermore, $\mathfrak{M}$ terminates in each iteration of Algorithm 8when $\unicode[STIX]{x1D6FC}<c_{2}\unicode[STIX]{x1D707}_{k}^{2}$ for some $c_{2}>0$ , then an upper bound on the number of function evaluations needed to obtain

(5.12) $$\begin{eqnarray}\Vert \unicode[STIX]{x1D6FB}f_{\unicode[STIX]{x1D707}_{\ast }}(\boldsymbol{x}_{k})\Vert \leq \unicode[STIX]{x1D716},\end{eqnarray}$$

for $\unicode[STIX]{x1D716}\in (0,1)$ and $\unicode[STIX]{x1D707}_{\ast }\in O(n^{-1/2}\unicode[STIX]{x1D716})$ , is in $O(\unicode[STIX]{x1D716}^{-3})$ ; see Table A.1. We note that while the sequence of smoothing parameters $\unicode[STIX]{x1D707}_{k}$ induces a type of limiting behaviour of the gradients (as seen in (5.12)) returned by the method $\mathfrak{M}$ used in Algorithm 8, this still does not necessarily recover elements of the Clarke subdifferential of $f$ . The smoothing functions $f_{\unicode[STIX]{x1D707}_{k}}$ must satisfy an additional gradient consistency property in order for Algorithm 8to produce a sequence of points $\boldsymbol{x}_{k}$ converging to Clarke stationary points (Rockafellar and Wets Reference Rockafellar and Wets2009, Theorem 9.67).

Garmanjani et al. (Reference Garmanjani, Jùdice and Vicente2016) consider the use of a model-based trust-region method $\mathfrak{M}$ in Algorithm 8. The authors demonstrate the first-order convergence result (5.11); they also prove the same WCC as is proved by Garmanjani and Vicente (Reference Garmanjani and Vicente2012).