Meta-Optimization of Dimension Adaptive Parameter Schema for Nelder–Mead Algorithm in High-Dimensional Problems

Abstract

Although proposed more than half a century ago, the Nelder–Mead simplex search algorithm is still widely used. Four numeric constants define the operations and behavior of the algorithm. The algorithm with the original constant values performs fine on most low-dimensional, but poorly on high-dimensional, problems. Therefore, to improve its behavior in high dimensions, several adaptive schemas setting the constants according to the problem dimension were proposed in the past. In this work, we present a novel adaptive schema obtained by a meta-optimization procedure. We describe a schema candidate with eight parameters subject to meta-optimization and define an objective function evaluating the candidate’s performance. The schema is optimized on up to 100-dimensional problems using the Parallel Simulated Annealing with Differential Evolution global method. The obtained global minimum represents the proposed schema. We compare the performance of the optimized schema with the existing adaptive schemas. The data profiles on the Gao–Han modified quadratic, Moré–Garbow–Hilstrom, and CUTEr (Constrained and Unconstrained Testing Environment, revisited) benchmark problem sets show that the obtained schema outperforms the existing adaptive schemas in terms of accuracy and convergence speed.

1. Introduction

Today, we can find optimization algorithms in almost every field of science and technology. A number of optimization algorithms exist, invented to fulfill various requirements regarding convergence rate, precision, robustness, and more. One of them (the Nelder–Mead simplex algorithm [1]), although more than half a century old, is still extensively used for solving a wide range of continuous optimization problems. The algorithm’s popularity is due to its simplicity and reasonably good performance observed in practical optimization cases.

Despite its popularity, the algorithm has proven convergence issues. McKinnon [2] presented a family of two-dimensional functions that cause the Nelder–Mead Algorithm (NMA) to converge to a non-stationary point. Galántai [3] provided a sufficient condition for repeated inside contractions in two dimensions, causing non-convergence. A relatively poor theoretical background on algorithm convergence is available for up to two-dimensional problems. Lagarias et al. [4] proved various limited convergence results for two-dimensional strictly convex objective functions. Further, Lagarias et al. [5] proved convergence for a restricted NMA, a version of the NMA without expansion steps, on two-dimensional strictly convex $C^2$ functions with bounded level sets.

Many modifications of the original NMA were proposed to avoid the algorithm’s known deficiencies and provide convergence. Kelley [6] proposed an oriented restart in case of detected stagnation. Tseng’s [7] and Nazareth and Tseng’s [8] versions of the algorithm guarantee convergence with a sufficient descent approach. Price et al. [9] used simplex reshaping to achieve convergence on $C^1$ functions, again satisfying sufficient descent conditions. Bűrmen et al. [10], on the other hand, introduced a grid-restrained version of the NMA, thus making the algorithm a pattern search method. The approach was generalized to a successive approximation of the objective function by Bűrmen and Tuma [11].

However, convergence analysis of the unmodified, original NMA stays a hard mathematical problem. No theoretical background is available above two dimensions. Torczon [12] discovered that the NMA fails because search direction and downhill gradient become orthogonal when the problem dimension is large enough. Wright [13] reported that several scholars observed how the NMA deteriorates with dimensionality but without any explanation. Further, Han and Neumann [14] showed that the NMA makes less and less progress per iteration with an increasing problem dimension. Gao and Han [15] suggested that poor performance in high dimensions is due to an increasing fraction of reflection steps.

Researchers addressed the poor performance of the original NMA in high dimensions in two ways. First, they proposed various algorithm modifications to improve the convergence rate in high-dimensional parameter spaces. Fajfar et al. [16] used genetic programming to evolve a direct search procedure using reflection, expansion, and contraction steps. Musafer’s [17] modification adjusts simplex size and direction by performing different NMA steps on various axial combinations. Fajfar et al. [18] proposed random centroid perturbation for improving the search direction, to name a few.

The second approach deals with NMA parameters and does not modify the algorithm itself in any way. Gao and Han [15] proposed the first schema of dimension-dependent NMA parameters to maintain the algorithm’s descent property in high dimensions. Kumar and Suri [19] suggested another schema obtained from parameter sensitivity analysis on five test functions. Mehta [20], on the other hand, observed that two schemas based on Chebyshev spaced points outperform Gao–Han’s and Kumar–Suri’s schemas.

This paper presents the global minimum of the meta-optimization of the NMA adaptive parameter schema. We used the Parallel Simulated Annealing with Differential Evolution (PSADE) robust global optimization method [21] running on a cluster of personal computers as the meta-optimization method. The subjects of meta-optimization are eight coefficients whose values define an individual adaptive parameter Schema Candidate (SC). The schema’s mathematical formulation is set in advance and does not evolve during the procedure. We run the NMA using the SC on several test functions and evaluate the SC in each iteration of the meta-optimization. The global minimum represents the best adaptive parameter schema corresponding to the predefined mathematical formulation of the schema and used objective function. We compare the performance of the NMA using the optimized schema with the NMA using the existing schemas on modified quadratic, i.e., Gao–Han (GH) [15], Moré–Garbow–Hilstrom (MGH) [22], and Constrained and Unconstrained Testing Environment, revisited (CUTEr) [23], sets of benchmark problems. Since the proposed adaptive parameter schema results from a meta-optimization procedure, we do not provide a mathematical background explaining the schema’s performance. However, the proposed schema outperforms all the other schemas and thus, to the best of our knowledge, currently represents the best dimension adaptive parameter schema for the NMA.

2. Adaptive Parameter Schemas for NMA

Each mathematical symbol used in Section 2 and Section 3 is explained at first use. However, for clarification, all the symbols are also listed in the Abbreviations section at the end of the paper.

The original NMA [1,4] is well known, therefore we provide only a brief introduction. The NMA is an unconstrained minimization algorithm operating on an objective function f of $n_p$ variables handling $(n_p+1)$ points ${\mathbf{P}}_i$, i.e., simplex vertices, with objective function values $f\left({\mathbf{P}}_i\right)$, $i=0,1\dots n_p$. In one iteration of the NMA, the following steps are performed:

• Order and relabel simplex vertices to satisfy $f\left({\mathbf{P}}_0\right)\le f\left({\mathbf{P}}_1\right)\le \dots \le f\left({\mathbf{P}}_{n_p}\right)$. Calculate the centroid point $\overline{\mathbf{P}}=\frac{1}{n_p}{\sum }_{i=0}^{n_p-1}{\mathbf{P}}_i$ excluding the highest objective function value point.
• Reflect ${\mathbf{P}}_{n_p}$ over $\overline{\mathbf{P}}$ to obtain the reflected point ${\mathbf{P}}_r=\overline{\mathbf{P}}+\alpha (\overline{\mathbf{P}}-{\mathbf{P}}_{n_p})$, $\alpha >0$.
• If $f\left({\mathbf{P}}_r\right)<f\left({\mathbf{P}}_0\right)$, expand ${\mathbf{P}}_r$ to obtain the expanded point
  ${\mathbf{P}}_e=\overline{\mathbf{P}}+\frac{\beta }{\alpha }({\mathbf{P}}_r-\overline{\mathbf{P}})=\overline{\mathbf{P}}+\beta (\overline{\mathbf{P}}-{\mathbf{P}}_{n_p})$, $\beta >\alpha$.
  If $f\left({\mathbf{P}}_e\right)<f\left({\mathbf{P}}_r\right)$, replace ${\mathbf{P}}_{n_p}$ with ${\mathbf{P}}_e$ and end the iteration.
• If $f\left({\mathbf{P}}_r\right)<f\left({\mathbf{P}}_{n_p-1}\right)$, replace ${\mathbf{P}}_{n_p}$ with ${\mathbf{P}}_r$ and end the iteration.
• If $f\left({\mathbf{P}}_r\right)<f\left({\mathbf{P}}_{n_p}\right)$, contract ${\mathbf{P}}_r$ towards $\overline{\mathbf{P}}$ to obtain the contracted point
  ${\mathbf{P}}_{rc}=\overline{\mathbf{P}}+\frac{\gamma }{\alpha }({\mathbf{P}}_r-\overline{\mathbf{P}})=\overline{\mathbf{P}}+\gamma (\overline{\mathbf{P}}-{\mathbf{P}}_{n_p})$, $\gamma <\alpha$.
  If $f\left({\mathbf{P}}_{rc}\right)<f\left({\mathbf{P}}_{n_p}\right)$, replace ${\mathbf{P}}_{n_p}$ with ${\mathbf{P}}_{rc}$ and end the iteration.
• If $f\left({\mathbf{P}}_r\right)\ge f\left({\mathbf{P}}_{n_p}\right)$, contract ${\mathbf{P}}_{n_p}$ towards $\overline{\mathbf{P}}$ to obtain the contracted point
  ${\mathbf{P}}_{nc}=\overline{\mathbf{P}}+\gamma ({\mathbf{P}}_{n_p}-\overline{\mathbf{P}})=\overline{\mathbf{P}}-\gamma (\overline{\mathbf{P}}-{\mathbf{P}}_{n_p})$.
  If $f\left({\mathbf{P}}_{nc}\right)<f\left({\mathbf{P}}_{n_p}\right)$, replace ${\mathbf{P}}_{n_p}$ with ${\mathbf{P}}_{nc}$ and end the iteration.
• Shrink the entire simplex towards ${\mathbf{P}}_0$, ${\mathbf{P}}_i:={\mathbf{P}}_0+\delta ({\mathbf{P}}_i-{\mathbf{P}}_0)$, $\delta <1$, $i=1,2,\dots n_p$.

NMA iterations are repeated until convergence is achieved. Algorithm behavior depends on $\alpha$ (reflection), $\beta$ (expansion), $\gamma$ (contraction), and $\delta$ (shrink) parameter values. The NMA default values are

$$\alpha =1,\beta =2,\gamma =\frac{1}{2},\delta =\frac{1}{2}.$$

(1)

Adaptive parameter schemas define NMA parameter values as functions of the number of variables $n_p$. The existing schemas (Gao–Han Schema (GHS) [15], Kumar–Suri Schema (KSS) [19], Chebyshev Crude Schema (CCS) [20], and Chebyshev Refined Schema (CRS) [20]) considered in this paper are

$$\begin{array}{ccc}\hfill \mathrm{GHS}:& \hfill & \alpha =1,\beta =1+\frac{2}{n_p},\gamma =\frac{3}{4}-\frac{1}{2n_p},\delta =1-\frac{1}{n_p}\hfill \\ \hfill \mathrm{KSS}:& \hfill & \alpha =1+\frac{3}{5n_p},\beta =\frac{6}{5},\gamma =\frac{19}{20}-\frac{3}{n_p}-\frac{3}{n_p^2},\delta =1-\frac{1}{n_p}\hfill \\ \hfill \mathrm{CCS}:& \hfill & \alpha =1+cos\frac{(n_p-1-n_p\%2)\pi }{2n_p},\beta =1+cos\frac{(n_p-3-n_p\%2)\pi }{2n_p},\hfill \\ \hfill & \hfill & \gamma =1+cos\frac{(n_p+3+n_p\%2)\pi }{2n_p},\delta =1+cos\frac{(n_p+1+n_p\%2)\pi }{2n_p}\hfill \\ \hfill \mathrm{CRS}:& \hfill & \alpha =1+cos\frac{(n_c-1)\pi }{2n_c},\beta =1+cos\frac{(n_c-3)\pi }{2n_c},\hfill \\ \hfill & \hfill & \gamma =1+cos\frac{(n_c+5)\pi }{2n_c},\delta =1+cos\frac{(n_c+3)\pi }{2n_c},n_c=2(9+\lfloor \frac{n_p-1}{5}\rfloor )\hfill \end{array}$$

(2)

where % denotes modulo operation, and $\lfloor \rfloor$ the floor function.

In general, the initial simplex vertices ${\mathbf{P}}_i$ are random. In this paper, however, to assure repeatability of the results, the initial simplex is generated from the starting point ${\mathbf{x}}_0$ using Pfeffer’s method [15]. The first vertex is starting point ${\mathbf{P}}_0={\mathbf{x}}_0$. The remaining vertices are generated by varying the $i\mathrm{th}$ component ${\mathbf{P}}_i={\mathbf{x}}_0+{ϵ}_i{\mathbf{e}}_i$. ${\mathbf{e}}_i$ is the $i\mathrm{th}$ component unit vector, and ${ϵ}_i$ is given by

$${ϵ}_i=\left\{\begin{array}{ccc}0.05{\mathbf{x}}_0{\mathbf{e}}_i\hfill & \hfill & {\mathbf{x}}_0{\mathbf{e}}_i\ne 0\hfill \\ 0.00025\hfill & \hfill & {\mathbf{x}}_0{\mathbf{e}}_i=0\hfill \end{array}\right.,i=1,2\dots n_p.$$

(3)

The starting point ${\mathbf{x}}_0$ is ${[1,1\dots 1]}^{\mathrm{T}}$ for GH benchmarks [15], and as given in [22] for MGH benchmarks.

A Nelder–Mead run terminates when the simplex becomes too flat or shrinks below a certain size. In this paper, we use tolerances $To{l}_{\mathrm{f}}$ for simplex flatness and $To{l}_{\mathrm{X}}$ for simplex size. A Nelder–Mead run stops when both criteria (4) are met:

$$\underset{i=1}{\overset{n_p}{max}}|f\left({\mathbf{P}}_i\right)-f\left({\mathbf{P}}_0\right)|<To{l}_{\mathrm{f}}\underset{j=0}{\overset{n_p-1}{max}}\underset{i=1}{\overset{n_p}{max}}|P_{ij}-P_{0j}|<To{l}_{\mathrm{X}},$$

(4)

where ${\mathbf{P}}_i={[P_{i0},P_{i1},\dots P_{i(n_p-1)}]}^T$ is the $i\mathrm{th}$ vertex of the simplex.

3. Optimization of the Adaptive Parameter Schema

The default (1) and adaptive parameter schema functions (2) shown in Figure 1 have in general similar behavior. By choosing appropriate values $c_{0p}$ and $c_{1p}$, an individual parameter c from a particular schema could be closely fitted with function $c=c_{0p}+\frac{c_{1p}}{n}$.

At this point, a question arises: can a better adaptive parameter schema of the formulation (5) be obtained by choosing the right values for $c_{0\alpha }$, $c_{1\alpha }$, $c_{0\beta }$, $c_{1\beta }$, etc.? Do such values exist, and what are they? A meta-optimization procedure could provide some answers.

$$\alpha =c_{0\alpha }+\frac{c_{1\alpha }}{n_p},\beta =c_{0\beta }+\frac{c_{1\beta }}{n_p},\gamma =c_{0\gamma }+\frac{c_{1\gamma }}{n_p},\delta =c_{0\delta }+\frac{c_{1\delta }}{n_p}$$

(5)

Let us first define the meta-optimization procedure. As in any optimization, we need optimization parameters, an objective function, and an optimization method. Optimization parameters ($c_{0\alpha }\dots c_{1\delta }$) follow from the mathematical formulation describing an SC (5). Therefore, we have an eight-dimensional meta-optimization parameter space.

The meta-optimization objective function measures the weighted difference in data profiles [24] between the best reference schema and an SC. For a mathematical formulation of the objective function, some definitions are needed. A data profile function of a schema s over a set of benchmark problems $\mathcal{P}$

$$d_{s\mathcal{P}}\left(\kappa \right)=\frac{|p\in \mathcal{P}:\frac{t_{ps}}{n_p+1}\le \kappa |}{\left|\mathcal{P}\right|}$$

(6)

defines the fraction of problems in a set $\mathcal{P}$ that schema s solves in $\kappa$ simplex gradient estimates. s is a schema from the set of schemas $\mathcal{S}$ ($s\in \mathcal{S}=\{$NMA, GHS, KSS, CCS, CRS, SC}). $|·|$ denotes cardinality of a set. $t_{ps}$ is the number of objective function evaluations needed by schema s to achieve convergence on problem p, and $n_p$ is the problem dimension. Since one simplex gradient estimate corresponds to $n_p+1$ objective function evaluations, fraction $\frac{t_{ps}}{n_p+1}$ is the number of simplex gradient estimates required for convergence. Convergence is achieved when

$$f\left(\mathbf{x}\right)\le f_{\mathrm{L}}+\tau (f\left({\mathbf{x}}_{\mathbf{0}}\right)-f_{\mathrm{L}}),$$

(7)

where $\mathbf{x}$ is an evaluated point in $n_p$-dimensional parameter space, and $f_{\mathrm{L}}$ is the lowest objective function value reached by any of the schemas $s\in \mathcal{S}$ in ${\kappa }_{\max }$ simplex gradient estimates. Tolerance $\tau$ specifies accuracy level. Moré and Wild [24] use tolerance values ${10}^{-1}$, ${10}^{-3}$, ${10}^{-5}$, and ${10}^{-7}$. We set convergence condition tolerance $\tau$ to ${10}^{-7}$ as Mehta did in his work [20]. If a particular schema s fails to satisfy condition (7) for problem p, then $t_{ps}$ is set to infinity.

The final meta-optimization objective function $h($SC) evaluates a particular SC defined by eight meta-optimization parameters $c_{0\alpha }\dots c_{1\delta }$ with

$$h\left(\mathrm{SC}\right)=\sum _{\mathcal{P}\in \mathcal{X}}\sum _{\kappa =1}^{{\kappa }_{\max }}\left((\underset{r\in \mathcal{R}}{max}{d}_{r\mathcal{P}}\left(\kappa \right)-d_{\mathrm{SC}\mathcal{P}}\left(\kappa \right))\times \left\{\begin{array}{cc}{w}_{{\mathcal{P}}_{+}}\hfill & {\displaystyle \underset{r\in \mathcal{R}}{max}{d}_{r\mathcal{P}}\left(\kappa \right)>d_{\mathrm{SC}\mathcal{P}}\left(\kappa \right)}\hfill \\ w_{{\mathcal{P}}_{-}}\hfill & {\displaystyle \underset{r\in \mathcal{R}}{max}{d}_{r\mathcal{P}}\left(\kappa \right)\le d_{\mathrm{SC}\mathcal{P}}\left(\kappa \right)}\hfill \end{array}\right.\right).$$

(8)

The GH and MGH benchmark problems are treated separately, $\mathcal{X}\in \left\{\mathrm{GH},\mathrm{MGH}\right\}$. For every number of simplex gradient estimates $\kappa$, the data profile of the SC is compared with the best of the reference profiles, $r\in \mathcal{R}=\mathcal{S}-\{$SC}. The difference is weighted with $w_{{\mathcal{P}}_{+}}$ when the best reference is better, and $w_{{\mathcal{P}}_{-}}$ otherwise. Weights $w_{{\mathcal{P}}_{+}}$ and $w_{{\mathcal{P}}_{-}}$ define a trade-off between under- and over-achieving the optimization goal (which in turn is the best schema’s performance). Usually, under-achieving is penalized more than over-achieving is rewarded. In our case, the weight values were set to $w_{{\mathrm{GH}}_{+}}=w_{{\mathrm{MGH}}_{+}}=10$, and $w_{{\mathrm{GH}}_{-}}=w_{{\mathrm{MGH}}_{-}}=1$.

Finally, we have to choose an optimization method to perform our meta-optimization procedure. We are searching for a global minimum in eight-dimensional parameter space. Therefore, a global optimization method with proven convergence can do the job. We chose PSADE [21], a parallel version of [25] since it is available in the PyOPUS Python package [26]. Among others, the PyOPUS package includes optimization algorithms (original NMA included), parallel processing support, and benchmark problems (GH, MGH, and CUTEr problems included), all the ingredients needed in our meta-optimization procedure. It can be found in the Python Package Index (PyPI) software repository. The PSADE method exhibited good performance on global benchmark functions as well as on real optimization problems [27,28,29]. Further, PSADE is an asynchronous global method achieving speedups up to the number of slave computational cores when run in parallel. We ran PSADE on a cluster of 25 personal computers. Instead of PSADE, one of a plethora of newer global methods, e.g., [30,31,32,33], could be used. However, besides faster convergence, we do not expect significantly better results.

A meta-optimization search for a better NMA parameter schema requires significant computational power. An individual SC, represented by eight meta-optimization parameter values $c_{0\alpha }\dots c_{1\delta }$, has to be evaluated against reference parameter schemas in every meta-optimization iteration. One SC evaluation requires as many Nelder–Mead runs as there are problems included in the objective function evaluation. In general, a single Nelder–Mead run stops when the termination criteria are achieved, e.g., when the simplex becomes flat or shrinks below the tolerance. However, additional improvement is possible if the procedure runs further. When driven beyond tolerances, a non-convergent SC may become convergent, although rather slow. Gao and Han [15], Kumar and Suri [19], and Mehta [20] all set an absolute limit to the number of objective function evaluations to ${10}^6$, i.e., 9900–90,909 simplex gradient estimates. However, in their results, 1000–2000 estimates are needed on average for the NMA using an adaptive parameter schema to converge. They set termination tolerances $To{l}_{\mathrm{f}}$, $To{l}_{\mathrm{X}}$ in range ${10}^{-10}$–${10}^{-4}$. Therefore, after some preliminary tests, we set ${\kappa }_{\max }$ to 5000. More would be better. And—after some preliminary experimental optimization runs—we established that around ${10}^6$ meta-optimization iterations are required to achieve convergence in an unconstrained eight-dimensional parameter space using the robust PSADE global optimization method [21]. Again, the more, the better. Thus, the total number of required objective function evaluations can be estimated as the number of simplex gradient estimates per benchmark problem times the number of simplex vertices summed over all benchmark problems times the number of meta-optimization iterations. For GH and MGH benchmarks altogether, this is over ${10}^{13}$ evaluations.

Therefore, brute force is not very promising. Instead, we tuned the meta-optimization parameters by conducting a series of shorter meta-optimization runs. We varied parameter space constraints, NMA termination tolerances $To{l}_{\mathrm{f}}$ and $To{l}_{\mathrm{X}}$, the gradient estimates limit ${\kappa }_{\max }$ per Nelder–Mead run, meta-optimization iteration limit, and we also experimented with the objective function definition at the beginning. By observing the results, we gradually eliminated parts that were not essential, e.g., setting NMA termination tolerances $To{l}_{\mathrm{f}}$ and $To{l}_{\mathrm{X}}$ too low can significantly extend the meta-optimization procedure without producing a significantly better result. The final meta-optimization parameter values are as follows: the parameter space is a discrete eight-dimensional box with $0.01$ grid and constraints set to ${[c_{0\alpha },c_{1\alpha },c_{0\beta },c_{1\beta },c_{0\gamma },c_{1\gamma },c_{0\delta },c_{1\delta }]}^{\mathrm{T}}\in [[0.80,1.20],$$[0.20,0.60],$$[0.85,1.25],$$[0.35,0.75],$$[0.65,1.05],$$[-0.50,-0.10],$$[0.05,0.45],$${[-0.40,0.00]]}^{\mathrm{T}}$. NMA termination tolerances are set to $To{l}_{\mathrm{f}}=To{l}_{\mathrm{X}}={10}^{-4}$. The simplex gradient estimates limit was set to ${\kappa }_{\max }=5000$, enough to catch degenerated SCs. The final objective function formulation is given in (8). Meta-optimization iteration limit was set to ${10}^6$.

The final meta-optimized parameter values, i.e., the global minimum of (8), are

$$\begin{array}{cccc}{c}_{0\alpha }=1.02,\hfill & c_{1\alpha }=0.31,\hfill & c_{0\beta }=1.06,\hfill & c_{1\beta }=0.53,\hfill \\ c_{0\gamma }=0.82,\hfill & c_{1\gamma }=-0.27,\hfill & c_{0\delta }=0.28,\hfill & c_{1\delta }=-0.19.\hfill \end{array}$$

(9)

This is the global minimum when (5) is used as the adaptive parameter schema. We chose (5) as the adaptive parameter schema because of its similarity to the existing schemas.

4. Evaluation of the Optimized Schema with Discussion

In this section, we present and discuss the properties of the optimized schema (9). We analyze its performance and compare it with the other schemas, (1) and (2). Since the schema is a result of the meta-optimization procedure, we do not deal with the issue of why the schema performs well. A mathematical evaluation of an arbitrary schema is given with the objective function definition (8). However, an analytical solution to the defined meta-optimization problem exceeds the scope of this paper.

The optimized parameter schema (9) is compared with the original NMA (1) and the existing adaptive parameter schemas (2), GHS, KSS, CCS, and CRS, in Figure 2. We can observe that the optimized schema follows the same pattern as the other schemas. However, the optimized schema approaches its high-dimensional value faster which is the consequence of smaller $|c_{1p}/c_{0p}|$ ratios. Otherwise, the optimized schema curves do not significantly deviate from the others. Shrink parameter $\delta$ is an exception. While $\delta$ in all the other adaptive schemas approaches 1, the optimized schema’s high-dimensional $\delta$ value is $c_{0\delta }=0.28$ which is far lower, even than $0.5$ of the original NMA. On the other hand, the shrink parameter turns out to be insignificant for GH and MGH benchmarks. The fraction of shrink steps is at most 0.5% for Penalty I and Penalty II problems. However, the fraction of shrink steps in the optimized schema is 0.0% for all GH and MGH benchmarks.

To evaluate the obtained schema further, we compare it in terms of accuracy and speed. Note that the meta-optimization objective function (8) in combination with NMA termination tolerances rewards speed. Therefore, accuracy was not a subject of meta-optimization.

Table 1. Accuracy of the original NMA, the existing adaptive schemas (GHS, KSS, CCS, CRS), and the optimized schema on Gao–Han (GH) modified quadratic benchmark problems.

	${\mathit{n}}_{\mathit{p}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{NMA}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{GHS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{KSS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{CCS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{CRS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{opt}.}$
$ϵ=$ 0.0	10	3.5 × 10−323	0.0	0.0	3.5 × 10−323	0.0	0.0
$\sigma =$ 0.0	20	2 × 10−322	0.0	0.0	0.0	0.0	10−323
	30	1.14 × 10−11	0.0	0.0	0.0	0.0	0.0
	40	2.03 × 10−4	0.0	0.0	0.0	0.0	0.0
	50	5.54 × 10−4	0.0	0.0	0.0	0.0	0.0
	60	1.38 × 10−5	0.0	0.0	5 × 10−324	0.0	5 × 10−324
	70	5.76 × 10−5	0.0	0.0	0.0	0.0	0.0
	80	4.87 × 10−6	5 × 10−323	0.0	0.0	0.0	0.0
	90	2.75 × 10−6	1.4 × 10−322	0.0	5 × 10−324	0.0	0.0
	100	3.19 × 10−6	6 × 10−323	0.0	2 × 10−323	0.0	0.0
$ϵ=$ 0.05	10	0.0	0.0	0.0	0.0	0.0	0.0
$\sigma =$ 0.0	20	6.23 × 10−322	0.0	0.0	0.0	5 × 10−324	0.0
	30	5.31 × 10−3	0.0	0.0	0.0	0.0	5 × 10−324
	40	1.32 × 10−2	0.0	0.0	0.0	0.0	0.0
	50	1.62 × 10−1	0.0	0.0	0.0	0.0	0.0
	60	12.7	0.0	0.0	0.0	0.0	0.0
	70	8.24	2 × 10−323	0.0	0.0	0.0	0.0
	80	32.2	1.24 × 10−322	0.0	5 × 10−324	0.0	0.0
	90	3.77	5.4 × 10−323	5 × 10−324	5 × 10−324	0.0	0.0
	100	278	6 × 10−323	0.0	10−323	0.0	0.0
$ϵ=$ 0.0	10	0.0	0.0	0.0	0.0	0.0	0.0
$\sigma =$ 0.0001	20	2.05 × 10−3	0.0	0.0	0.0	0.0	0.0
	30	1.91 × 10−5	0.0	0.0	0.0	0.0	0.0
	40	16.7	0.0	0.0	0.0	0.0	0.0
	50	2.63	0.0	0.0	0.0	0.0	0.0
	60	10.9	0.0	0.0	10−323	0.0	0.0
	70	276	5 × 10−324	0.0	0.0	0.0	0.0
	80	292	8 × 10−323	0.0	0.0	0.0	0.0
	90	11.7	6 × 10−323	5 × 10−324	5 × 10−324	0.0	0.0
	100	48.7	7.4 × 10−323	0.0	10−323	0.0	0.0
$ϵ=$ 0.05	10	1.5 × 10−323	0.0	0.0	0.0	0.0	0.0
$\sigma =$ 0.0001	20	1.93 × 10−4	0.0	0.0	0.0	5 × 10−324	5 × 10−324
	30	1.12 × 10−2	0.0	0.0	0.0	0.0	5 × 10−324
	40	7.31 × 10−1	0.0	0.0	0.0	0.0	0.0
	50	37.2	0.0	0.0	0.0	0.0	5 × 10−324
	60	179	0.0	0.0	0.0	0.0	0.0
	70	18.7	3 × 10−323	0.0	0.0	0.0	0.0
	80	16.4	7.4 × 10−323	0.0	5 × 10−324	0.0	0.0
	90	1480	1.3 × 10−322	1.5 × 10−323	0.0	0.0	0.0
	100	3802	5 × 10−323	0.0	10−323	0.0	0.0
accurate		7/40	40/40	40/40	40/40	40/40	40/40

shows accuracy results for a GH modified quadratic function (10) up to problem dimension $n_p=100$. Parameter $ϵ\ge 0$ defines condition number of matrix $\mathbf{D}$, and $\sigma \ge 0$ specifies deviation from quadratic form. The minimum value of (10) is zero (${min}_{\mathbf{x}\in {\mathbb{R}}^{n_p}}f\left(\mathbf{x}\right)=0$) for any $ϵ$, $\sigma$, or $n_p$. An individual Nelder–Mead run was stopped after ${\kappa }_{\max }$ = 25,000 simplex gradient estimates, that is, after $N_{\mathrm{eval}}=25,000(n_p+1)$ objective function evaluations. We applied no tolerance-based termination criteria, i.e., $To{l}_{\mathrm{f}}=To{l}_{\mathrm{X}}=0$. The table shows the lowest achieved objective function values. A schema is considered to be accurate if the achieved minimum value is correct to at least six decimal places, i.e., when $f\left(\mathbf{x}\right)<$ 5 × 10${}^{-7}$. The accurate schema values are shown in bold. If a schema does not converge according to condition (7) for $\tau ={10}^{-7}$, the value is shown in italics.

$$\begin{array}{c}f\left(\mathbf{x}\right)={\mathbf{x}}^{\mathrm{T}}\mathbf{D}\mathbf{x}+\sigma {\left({\mathbf{x}}^{\mathrm{T}}\mathbf{B}\mathbf{x}\right)}^2\\ \mathbf{D}=\mathrm{diag}[(1+ϵ),{(1+ϵ)}^2,\dots {(1+ϵ)}^{n_p}],\mathbf{B}={\mathbf{U}}^{\mathrm{T}}\mathbf{U},\mathbf{U}=\left[\begin{array}{ccc}1& \cdots & 1\\ & \ddots & ⋮\\ & & 1\end{array}\right]\end{array}$$

(10)

Given 25,000 simplex gradient estimates, the optimized schema is accurate and convergent for all 40 GH benchmark problems, as all the existing schemas also are. As expected, the original NMA encounters accuracy and convergence problems in high dimensions.

Accuracy results for MGH benchmark problems [22] are shown in

Table 2. Accuracy of the original NMA, the existing adaptive schemas (GHS, KSS, CCS, CRS), and the optimized schema on Moré–Garbow–Hilstrom (MGH) benchmark problems.

Function	${\mathit{n}}_{\mathit{p}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{NMA}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{GHS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{KSS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{CCS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{CRS}}$	$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{opt}.}$
Extended	12	2.91 × 10−28	2.35 × 10−29	1.73 × 10−28	4.64 × 10−29	5.18 × 10−29	6.88 × 10−29
Rosenbrock	18	20.0	6.97 × 10−29	6.51 × 10−29	1.44 × 10−28	5.66 × 10−28	1.35 × 10−28
	24	12.5	1.72 × 10−28	2.58 × 10−28	1.86 × 10−28	3.72 × 10−28	6.21 × 10−28
	30	34.5	4.09 × 10−28	6.70 × 10−28	7.28 × 10−28	3.70 × 10−27	2.76 × 10−28
	36	49.1	8.72 × 10−28	6.64 × 10−28	6.81 × 10−28	4.77 × 10−28	1.11 × 10−27
Extended	12	8.34 × 10−55	3.33 × 10−57	7.09 × 10−59	3.09 × 10−57	1.07 × 10−59	1.18 × 10−58
Powell	24	1.33 × 10−9	1.83 × 10−54	3.45 × 10−56	5.37 × 10−56	1.67 × 10−53	3.39 × 10−53
singular	40	1.69 × 10${}^{-6}$	1.06 × 10−50	2.34 × 10−52	1.46 × 10−52	5.22 × 10−52	2.33 × 10−53
	60	4.16 × 10−4	9.71 × 10${}^{-6}$	3.43 × 10−50	2.88 × 10−52	1.28 × 10−37	1.16 × 10−46
Penalty I	10	7.57 × 10−5	7.09 × 10−5	7.09 × 10−5	7.60 × 10${}^{-5}$	7.09 × 10−5	7.09 × 10−5
Penalty II	10	2.98 × 10−4	2.94 × 10−4	2.94 × 10−4	2.98 × 10−4	2.95 × 10−4	2.94 × 10−4
Variably	12	4.77	3.72 × 10−30	1.47 × 10−29	3.64 × 10−29	2.30 × 10−29	1.78 × 10−29
dimensioned	18	4.22	8.96 × 10−30	2.06 × 10−29	1.52 × 10−29	4.74 × 10−29	4.25 × 10−29
	24	11.5	8.22 × 10−29	8.37 × 10−29	7.52 × 10−29	9.23 × 10−29	2.27 × 10−28
	30	40.5	8.08 × 10−29	1.08 × 10−28	1.06 × 10−28	3.38 × 10−28	4.49 × 10−28
	36	60.1	4.21 × 10−28	1.46 × 10−28	8.82 × 10−29	8.35 × 10−28	7.60 × 10−28
Trigonometric	10	2.80 × 10${}^{-5}$	2.80 × 10−5	2.80 × 10−5	2.80 × 10−5	2.80 × 10−5	2.80 × 10−5
	20	1.35 × 10−6	1.35 × 10−6	6.03 × 10−6	6.86 × 10−6	1.35 × 10−6	1.35 × 10−6
	30	2.20 × 10−5	9.90 × 10−7	9.90 × 10−7	5.65 × 10−6	9.90 × 10−7	5.98 × 10−7
	40	1.41 × 10−5	1.55 × 10−6	3.95 × 10−6	1.68 × 10−7	5.58 × 10−7	1.55 × 10−6
	50	2.52 × 10−5	2.24 × 10−7	3.41 × 10−6	9.23 × 10−7	1.11 × 10−6	2.24 × 10−7
	60	3.87 × 10−5	8.68 × 10−7	7.57 × 10−7	7.57 × 10−7	1.27 × 10−6	4.62 × 10−7
Discrete	10	6.85 × 10−32	3.03 × 10−33	8.36 × 10−32	2.20 × 10−31	3.07 × 10−32	1.59 × 10−32
boundary	20	4.69 × 10−30	7.24 × 10−32	2.51 × 10−32	2.39 × 10−32	1.05 × 10−31	3.92 × 10−32
value	30	9.87 × 10−6	1.10 × 10−31	1.43 × 10−31	1.19 × 10−31	2.02 × 10−31	9.29 × 10−32
	40	6.46 × 10−6	4.58 × 10−31	3.55 × 10−31	1.37 × 10−31	4.76 × 10−31	3.45 × 10−31
	50	5.72 × 10−6	6.02 × 10−31	5.70 × 10−31	2.84 × 10−31	5.35 × 10−31	4.35 × 10−31
	60	3.19 × 10−6	2.46 × 10−30	8.09 × 10−31	6.64 × 10−31	1.11 × 10−30	7.39 × 10−31
Discrete	10	1.91 × 10−31	4.24 × 10−33	1.44 × 10−32	2.27 × 10−31	2.56 × 10−32	3.08 × 10−33
integral	20	7.69 × 10−30	4.62 × 10−32	2.90 × 10−32	6.27 × 10−32	3.40 × 10−32	2.37 × 10−32
equation	30	7.11 × 10−4	2.22 × 10−31	2.50 × 10−31	3.45 × 10−31	8.55 × 10−32	2.50 × 10−31
	40	3.63 × 10−4	3.82 × 10−31	3.07 × 10−31	3.21 × 10−31	4.25 × 10−31	3.04 × 10−31
	50	3.05 × 10−3	8.51 × 10−31	1.34 × 10−30	5.95 × 10−31	1.47 × 10−30	7.34 × 10−31
	60	4.46 × 10−4	2.24 × 10−30	1.30 × 10−30	1.59 × 10−30	9.74 × 10−31	6.12 × 10−31
Broyden	10	3.99 × 10−30	3.12 × 10−30	7.31 × 10−30	3.28 × 10−29	2.92 × 10−30	2.92 × 10−30
tridiagonal	20	3.20 × 10−26	1.63 × 10−29	3.34 × 10−29	6.15 × 10−29	2.45 × 10−29	3.15 × 10−29
	30	4.70 × 10−26	1.68 × 10−28	1.19 × 10−28	9.66 × 10−29	6.50 × 10−29	8.18 × 10−29
	40	9.11 × 10−14	2.24 × 10−28	6.73 × 10−28	3.70 × 10−28	2.45 × 10−28	2.24 × 10−28
	50	2.67 × 10−13	5.82 × 10−28	6.98 × 10−28	4.61 × 10−28	6.86 × 10−28	4.54 × 10−28
	60	3.78 × 10−11	9.12 × 10−28	1.69 × 10−27	1.01 × 10−27	1.34 × 10−27	1.10 × 10−27
Broyden	10	4.18 × 10−28	4.61 × 10−30	4.81 × 10−30	6.82 × 10−29	7.43 × 10−30	2.36 × 10−30
banded	20	1.85 × 10−26	2.63 × 10−29	6.04 × 10−29	7.47 × 10−29	1.60 × 10−28	4.77 × 10−29
	30	12.2	2.25 × 10−28	1.34 × 10−28	2.08 × 10−28	3.15 × 10−28	2.89 × 10−28
	40	2.02 × 10${}^{-6}$	3.48 × 10−28	7.41 × 10−28	9.32 × 10−28	1.34 × 10−28	3.25 × 10−28
	50	9.33 × 10${}^{-5}$	6.08 × 10−28	1.38 × 10−27	6.78 × 10−28	5.98 × 10−28	1.04 × 10−27
	60	5.13 × 10${}^{-6}$	2.93 × 10−27	3.44 × 10−27	4.09 × 10−27	7.98 × 10−28	7.59 × 10−28
accurate		15/46	40/46	40/46	39/46	39/46	42/46

. The same rules (25,000 simplex gradient estimates, $To{l}_{\mathrm{f}}=To{l}_{\mathrm{X}}=0$, etc.) apply. Minima of the used MGH functions are all zero with the exception of Penalty I and Penalty II functions. For $n_p=10$, the corresponding minima are 7.0876515…×10${}^{-5}$ and 0.00029366054…, respectively. Thus, the six digit accuracy criteria are given by $f_{\mathrm{Penalty}\mathrm{I}}\left(\mathbf{x}\right)<$ 7.087655 ×10${}^{-5}$, and $f_{\mathrm{Penalty}\mathrm{II}}\left(\mathbf{x}\right)<$ 0.0002936615.

Again, with 25,000 simplex gradient estimates available, the optimized schema is accurate and convergent for all MGH benchmarks, except for four trigonometric functions, where the result is approaching the minimum, although not reaching it. We can observe similar behavior for the other schemas as well. Furthermore, for trigonometric function benchmarks, the schemas achieved relatively small objective function value reduction, which is reflected in relation $\tau \frac{f\left({\mathbf{x}}_0\right)}{f_{\mathrm{L}}}\ll 1$. Consequently, $f_{\mathrm{L}}\gg \tau (f\left({\mathbf{x}}_{\mathbf{0}}\right)-f_{\mathrm{L}})$, and convergence condition (7) degenerates into $f\left(\mathbf{x}\right)\lesssim f_{\mathrm{L}}$. Strictly following the convergence condition, only schemas reaching the lowest value $f_{\mathrm{L}}$ are considered convergent. Nevertheless, the optimized schema managed to produce the lowest objective function value in all trigonometric function benchmarks, except one.

Convergence speed of the obtained schema (9) is compared to the other schemas (1) and (2) using data profiles (6) and (7) [24]. Figure 3 shows data profiles for GH and MGH benchmark sets separately and combined. The profiles are calculated at $\tau ={10}^{-7}$ with ${\kappa }_{\max }$ = 25,000 simplex gradient estimates and without tolerance-based algorithm termination ($To{l}_{\mathrm{f}}=To{l}_{\mathrm{X}}=0$). The graphs are shown for up to 15,000 simplex gradient estimates.

Graphs in Figure 3 reveal a slight advantage in terms of speed for the optimized schema over the existing adaptive parameter schemas. The schemas perform similarly when only GH benchmarks are considered. Although the optimized schema can be declared as the fastest (solves the highest percentage of problems at almost any given $\kappa$), the remaining schemas quickly follow. All the adaptive schemas solve 100% of GH problems after approximately 2500 simplex gradient estimates. The CCS is the first to achieve this goal at ~2100 simplex gradient estimates. The original NMA, on the other hand, manages to be competitive with the adaptive schemas for the first 10% of problems. With the lower-dimensional problems solved, it starts to lag, finally solving less than 30% of GH problems.

The advantage of the optimized schema becomes notable when observing the MGH benchmark set. At ~1600 simplex gradient estimates, the optimized schema solves more than 90% of the problems while the remaining adaptive schemas solve up to 83%, and the original NMA merely 50% of the problems. It is clearly the fastest schema and has the highest final percentage of solved problems (98%). Other adaptive schemas solve up to 93%, and the original NMA solves 59% of the problems.

The optimized schema remains the fastest when all GH and MGH benchmark problems are considered. Its advantage over the existing adaptive schemas is 6% at 1000, 5% at 2000, and 3% at 3000 simplex gradient estimates. The original NMA manages to keep up for less than 20% of the problems at ~150 simplex gradient estimates. The final percentage of solved problems is 99% for the optimized schema, up to 97% for the remaining adaptive schemas, and 44% for the original NMA.

In Figure 4, the same measurement of convergence speed is repeated with tolerance based algorithm termination set to $To{l}_{\mathrm{f}}=To{l}_{\mathrm{X}}={10}^{-4}$. The profiles are once more calculated for $\tau ={10}^{-7}$. The maximum number of simplex gradient estimates ${\kappa }_{\max }$ = 25,000 does not play any role in this experiment because none of the individual Nelder–Mead runs ever uses the entire budget of available simplex gradient estimates. The algorithm is always terminated earlier by the tolerance-based criterion. The graphs are shown for up to 8000 simplex gradient estimates. No further progress is made beyond that point.

With tolerance-based termination enabled, the optimized schema performs even better than the other schemas. When we consider only GH benchmarks, the schemas again perform similarly, although, in general, the optimized schema is still slightly faster. The CCS is the first that solves 100% of GH problems in ~2100 estimates. Other adaptive schemas quickly follow. The only notable alteration can be observed for the original NMA which now performs worse and solves less than 20% of GH problems. The original NMA manages to converge in an additional ~10% of GH problems when it is allowed to run beyond the tolerance-based stopping criterion.

The advantage of the optimized schema becomes apparent in MGH data profiles. The KSS, CCS, and CRS keep pace for up to ~370 simplex gradient estimates where ~70% of the MGH problems are solved. The GHS starts to lag at ~240 estimates with ~40% of the problems solved. It catches up at 3000 estimates, finally solving 70% of the MGH problems. KSS, the best of the existing adaptive schemas, ends at 71%. The optimized schema is clearly better with 80% of the problems solved in 730 estimates, ending with 82% at 1660 estimates. The original NMA again performs worse compared to the run without tolerance-based termination, ending with 23% of solved problems.

For all-inclusive GH and MGH benchmark data profiles, the optimized schema starts to stand out from the rest of the existing adaptive schemas at ~60% of solved problems, achieved within ~300 simplex gradient estimates. The optimized schema reaches its final result, i.e., 90% of solved problems, in 2400 estimates. Its advantage over the best existing adaptive schema, i.e., KSS, is 6%. The KSS comes to 84% in 7020 estimates. As expected, the original NMA ends at a modest 21% of solved problems within 1200 estimates.

Figure 5 shows our last convergence speed measurement on GH, MGH, and CUTEr benchmark problems, 169 problems in total. The data profiles are shown for the CUTEr benchmark set, and GH, MGH, and CUTEr benchmark sets combined. They are calculated at convergence condition tolerance $\tau ={10}^{-7}$ with ${\kappa }_{\max }$ = 25,000 simplex gradient estimates limit. Cases without the tolerance-based algorithm termination ($To{l}_{\mathrm{f}}=To{l}_{\mathrm{X}}=0$), and with it ($To{l}_{\mathrm{f}}=To{l}_{\mathrm{X}}={10}^{-4}$), are shown. When the tolerance-based criterion is applied, the Nelder–Mead runs are always terminated before the limit of ${\kappa }_{\max }$ simplex gradient estimates is reached.

Data profiles in Figure 5 confirm our previous observations. Although not meta-optimized on CUTEr benchmark problems, the optimized schema solves the highest percentage of the problems in all shown cases at any given $\kappa$. It makes a difference of 4 to 6% compared to the first follower at ~400 simplex gradient estimates when tolerance-based algorithm termination is applied, and at ~1100 estimates when it is not. The optimized schema reaches or comes close to its final result in ~4000 estimates. It solves 91 to 97% of problems, GHS 84 to 93%, CRS 84 to 93%, KSS 86 to 91%, CCS 81 to 88%, and the original NMA manages 13 to 36%.

Besides the 40 GH and 46 MGH benchmarks, the following problems are included in data profiles in Figure 5: BrownAlmostLinear with dimensions $n_p=\{20,30,40,50,70,100\}$ from MGH set, HilbertQuadratic with dimensions $n_p=\{10,30,60,90\}$, OrenPower [34] with dimensions $n_p=\{30,50,60,70,80,90,100\}$, and ARWHEAD_100, DQDRTIC_50, DQDRTIC_100, SPARSINE_50, SPARSINE_100, CHNROSNB_25, CHNROSNB_50, SCOSINE_10, LIARWHD_100, FLETCHBV_100, DIXON3DQ_100, OSCIGRAD_25, OSCIGRAD_100, NONCVXUN_10, NONCVXUN_100, PENALTY1_50, PENALTY1_100, SINQUAD_50, SINQUAD_100, FLETCBV3_100, PENALTY2_100, TOINTGSS_50, TOINTGSS_100, ARGLINC_50, EXTROSNB_100, COSINE_100, TRIDIA_50, TRIDIA_100, NONDQUAR_100, QUARTC_25, QUARTC_100, FREUROTH_100, WATSON_31, ERRINROS_25, ERRINROS_50, NONDIA_20, NONDIA_30, NONDIA_50, NONDIA_90, NONDIA_100, MANCINO_20, MANCINO_30, DQRTIC_100, ENGVAL1_50, ENGVAL1_100, HILBERTA_10, FLETCBV2_100, TQUARTIC_10, EDENSCH_36, ARGLINA_50, ARGLINA_100, BOX_100, POWELLSG_36, POWELLSG_40, POWELLSG_60, POWELLSG_80, POWELLSG_100, POWER_75, POWER_100, HILBERTB_50, ARGLINB_50, MOREBV_50, BDQRTIC_100, SCURLY10_100, VAREIGVL_49, VAREIGVL_99 from CUTEr benchmark problem set.

Speed of convergence of a particular adaptive parameter schema is mirrored in the simplex’s best value descent during an individual Nelder–Mead run. The descent rate can be expressed by $cos\theta$, where $\theta$ is an angle between the search direction $\mathbf{d}$ and the gradient of the objective function $\nabla f\left(\mathbf{x}\right)$:

$$cos\theta =\frac{{\mathbf{d}}^{\mathrm{T}}\nabla f\left(\mathbf{x}\right)}{\left|\mathbf{d}\right||\nabla f(\mathbf{x}\left)\right|}.$$

(11)

Search direction $\mathbf{d}$ is locally descending when $cos\theta <0$. The fastest descent is achieved at $cos\theta =-1$. According to the NMA definition, the search direction is $\mathbf{d}=c({\mathbf{P}}_{n_p}$−$\overline{\mathbf{P}})$, where c is the reflection ($\alpha$), expansion ($\beta$), or contraction ($\gamma$) NMA parameter. The descent rate in a non-shrinking NMA iteration is therefore calculated as

$$cos\theta =\frac{{({\mathbf{P}}_{n_p}-\overline{\mathbf{P}})}^{\mathrm{T}}\nabla f\left({\mathbf{P}}_{n_p}\right)}{|{\mathbf{P}}_{n_p}-\overline{\mathbf{P}}\left|\right|\nabla f\left({\mathbf{P}}_{n_p}\right)|}.$$

(12)

The simplex’s best value descents and corresponding descent rates for $n_p=$ 100-dimensional GH benchmark problems are shown in Figure 6. The figure shows that all the existing schemas as well as the optimized adaptive schema manage to maintain some level of descent during the entire Nelder–Mead run. A higher descent rate, i.e., more negative $cos\theta$, ensures faster objective function descent and fulfillment of the termination criteria. The optimized schema is the fastest or near fastest in all shown cases except for the $ϵ=0.0$, $\sigma =0.0001$ case. This is partly reflected in the poorer descent rate of the optimized schema for that particular case.

The tolerance boundary intersections shown in Figure 6 are $t_{ps}$ values from (6). Data profiles in Figure 3, Figure 4 and Figure 5 summarize tolerance boundary intersections over the entire benchmark problem set.

The absence of a sufficient descent rate is fatal for the original NMA. The original NMA manages some slow descent only in the $ϵ=0.0$, $\sigma =0.0$ case. In all remaining cases, $cos\theta$ approaches $0^{\circ }$. Search direction $\mathbf{d}$ becomes orthogonal to the negative gradient which was first observed by Torczon [12]. As a consequence, the original NMA stops descending and does not achieve the convergence boundary.

In [4], the authors prove that NMA does not perform shrink iteration when the objective function is strictly convex. Furthermore, for a uniformly convex objective function, the descent rate is provided by expansion and contraction iterations [15], although the effect diminishes with problem dimension $n_p$. In other words, to maintain a sufficient descent rate, an adequate share of expansion and contraction iterations is needed. Note that a uniformly convex function is also strictly convex. Since the modified quadratic function (10) is uniformly convex, the above applies to the GH benchmark set. The share of non-reflection iterations, i.e., expansion and contraction iterations combined, is shown in Figure 7. In general, it declines with the problem dimension for all schemas and ($ϵ$,$\sigma$) pairs. Nevertheless, all the adaptive parameter schemas manage to keep the share above 5%, which provides a sufficient descent rate. The CCS stands out with its lowest non-reflection share above 35%, yet, such a high share is not reflected in better performance. The non-reflection share alone, therefore, does not guarantee high convergence speed. The lowest non-reflection share of the optimized schema is 12% at $n_p=100$ in the $ϵ=0.0$, $\sigma =0.0001$ case. The original NMA’s share is, on the other hand, never greater than 26% (which in turn is achieved for lower-dimensional problems). With problem dimension increase, it quickly drops as low as to 0.56% in the worst case, which confirms the poor performance and convergence problems of the original NMA schema (1).

5. Conclusions

Adaptive parameter schemas address poor NMA performance on high-dimensional problems. We used a meta-optimization procedure to find a novel adaptive parameter schema presented in this paper. Although the meta-optimization problem seems simple, brute force optimization is not feasible due to the immense computing power required. To set up the problem, we defined a mathematical formulation of the adaptive parameter schema and an objective function evaluating a schema’s performance. We tuned the meta-optimization parameters in a series of shorter meta-optimization runs. The final settings constrain the meta-optimization parameter space, define a single NMA run termination criteria to evaluate an SC’s performance, limit the number of NMA iterations to catch non-convergent SCs, and limit the number of meta-optimization iterations. We used PSADE, a robust global parallel asynchronous method.

The performance of the proposed adaptive parameter schema is discussed and compared with the existing schemas. The share of non-reflection iterations and the descent rate do not show any significant deviation of the proposed schema from the existing ones. However, data profiles on GH modified quadratic, MGH, and CUTEr benchmark problems show that the proposed schema outperforms the existing ones in both accuracy and convergence speed. We performed the evaluation with and without tolerance-based termination of the NMA.

The proposed schema is a result of a meta-optimization procedure. We evaluate its performance but, on the other hand, provide no mathematical explanation for why the schema performs so well. The proposed schema is the global minimum determined by the schema’s mathematical formulation and meta-optimization objective function definition.