{'fraction_non_alphanumeric': 0.04553995799712612, 'fraction_numerical': 0.01134077594782801, 'mean_word_length': 4.318128379967082, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 13, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SMGO: A Set Membership Approach to Data-Driven Global Optimization

Many science and engineering applications feature non-convex optimization problems where the objective function can not be handled analytically, i.e. it is a black box. Examples include design optimization via experiments, or via costly finite elements simulations. To solve these problems, global optimization routines are used. These iterative techniques must trade-off exploitation close to the current best point with exploration of unseen regions of the search space. In this respect, a new global optimization strategy based on a Set Membership (SM) framework is proposed. Assuming Lipschitz continuity of the cost function, the approach employs SM concepts to decide whether to switch from an exploitation mode to an exploration one, and vice-versa. The resulting algorithm, named SMGO (Set Membership Global Optimization) is presented. Theoretical properties regarding convergence and computational complexity are derived, and implementation aspects are discussed. Finally, the SMGO performance is evaluated on a set of benchmark non-convex problems and compared with those of other global optimization approaches.


Introduction
In many science and engineering fields, such as mechanical design, fluid-dynamics, electromagnetics, multi-physics simulations, control systems tuning, and chemical experiments, black-box optimization problems arise. Black-box functions are named as such due to the fact that an explicit mathematical model is unavailable, or too complicated to be handled analytically. In these cases, the optimization strategy can rely only on function values obtained through empirical tests. Moreover, these problems may present several local minima, and the time and resources required to carry out a single function evaluation are rather large, so that the solution method shall make the most efficient use of the available trials. To highlight this fact, we refer to the process of performing a new test as long function evaluation.
Seeking a minimizer in this framework is referred to as black-box optimization [1], or derivative-free optimization [2]. This is in contrast for example to gradient-based and Newton-type methods where the (first, or even further) derivatives are also assumed to be available. If the black-box function is assumed differentiable, one could resort to gradient estimation techniques to still employ these methods, which however require a rather large number of long function evaluations to estimate a local quantity (i.e., the gradient at the currently evaluated point) [3] and are designed to converge to a local optimum by always improving from the initial guess, not to explore the decision space searching for a global solution.
Due to time/resource limitations, in global optimization there are two conflicting aspects that must be considered when choosing the next test point for a long function evaluation [3]. The first, exploitation, pertains to improving the result by testing more points in the vicinity of the current best known one. The second, exploration, aims to gain more information about the cost function in other regions of the search space, seeking other, possibly better, minima. Exploitation and exploration strategies are at the core of most global optimization approaches, including stochastic global search, response surface, and Lipschitz-based methods.
Stochastic global search algorithms, such as random search methods [4], particle swarm algorithms [5], and its offshoot variants [6,7], use a population of search points per iteration. According to the corresponding function values, heuristics are applied to calculate the search point locations for the next iteration (or "generation"). These methods are widespread but require numerous long function evaluations per generation, which can be impractical. Moreover, aspects such as guaranteed convergence and optimality gap may be difficult to analyze. Response surface methods, on the other hand, generate an approximation of the black-box function via, e.g., Gaussian process regression [1], radial basis functions [8,9], or neural networks [10]. These approximations can then be used for the exploitation or exploration routines, whichever is appropriate. Albeit the number of long function evaluation is lower in this case, optimality gap is again not easily assessed, and there is an additional exponential computational burden to derive and refine the approximated objective function.
Another category of black-box optimizers are socalled Lipschitz-based algorithms. They rely on the assumption that the long function is Lipschitz continuous. An early method, proposed independently by Piyavskii [11] and Shubert [12], uses a known Lipschitz constant to draw the lower bounds of the black-box function given the existing samples, starting with the bounds (corners) of the search space. In this context, the Piyavskii-Shubert ("sawtooth") method chooses the next point by searching the location of the minimum lower bound. However, the method is proposed for one-dimensional function optimization only; furthermore, Lipschitz constants are rarely known in practical cases, and need to be estimated in the process.
A later proposal named DIRECT [13,14,15] is a modification to Piyavskii-Shubert method, removing the need for a known Lipschitz constant and starting samples from the corners. Instead, it invariably performs a first sample at the center of the search space, which is afterwards subdivided into three intervals. The Lipschitz constant estimate is used to draw the lower bounds, which are used to select potentially optimal intervals (hyperrectan-gles). Such hyperrectangles are iteratively subdivided, and their corresponding centers are sampled by batch, and potential intervals for sampling are selected again. This cycle proceeds until the maximum number of iterations is reached, with each iteration performing multiple long function evaluations. An offshoot method named DISIMPL [16,17] employs simplices instead of hyperrectangles, thus extending the usage from rectangular search spaces only to more general polytopic ones. On one hand, DIRECT and DISIMPL are completely deterministic and reproducible (assuming no noise on function evaluation), however, they only admit a fixed starting point, and do not propose how to handle existing sampled points at the beginning of the algorithms. Furthermore, the batch-based sampling means that the next set of trial points is not chosen until the current batch of sampling is finished, which can mean overlooking possible precious information added by individual samples as they come in.
LIPO and AdaLIPO [18] use a random binary variable to decide between exploitation and exploration; the only difference being that the former assumes a fixed Lipschitz constant, while the latter estimates it from existing data. For exploitation, they choose a random point within the set of potential optimizing points, i.e. those whose lower bound on the black-box function is below the current best sample. For exploration, it chooses a random point from the entire search space. These algorithms, owing to the simple implementation, are fast in terms of choosing the next sampling point, and can easily scale to high dimensional problems without significant hit on computational burden. However, such a completely randomized exploitation/exploration sampling may be inefficient in finding new local minima, because it does not fully utilize the shape of lower bounds. Furthermore, the results of such an algorithm are unrepeatable because of its randomized design.
In the described context, we propose a new Lipschitz-based algorithm that exploits Set Membership (SM) nonlinear function approximation theory [19]. So far, SM theory proved effective for nonlinear system identification [20], filter design [21], and controller design [22], among others. Here, we employ SM theory for the first time in global optimization. Our paper delivers the following contributions: (SMGO) algorithm for Lipschitz-continuous black-box functions is introduced. It employs the concepts of SM lower-and upper bounds in exploitation and exploration routines. Unlike the above-mentioned Lipschitz-based algorithms, we address the exploration problem systematically, by casting it as minimization of the uncertainty as inferred using SM techniques. Furthermore, our proposed algorithm is completely repeatable, and allows one to consider custom starting points and initial samples when available; • Theoretical properties of SMGO are derived, including convergence, optimality gap, and computational complexity; • Practical aspects are discussed, and ways to improve the computational efficiency are described. In particular, an iterative implementation is described, which exploits results from previous iterations to alleviate the computations at the current one; • The SMGO algorithm is compared with other representative optimizers with different test functions, showing that the performance of SMGO is very competitive w.r.t. the state of the art, especially on functions with many global minima.
A preliminary version of this work appeared in [23]. With respect to that contribution, the algorithm presented here is more efficient, moreover the theoretical results and implementation analysis are new, and more benchmark results are presented. The MATLAB code for SMGO, as well as the other scripts used to generate the results, are available on GitHub under the following URL: https://github.com/lorenzosabugjr/smgo. This paper is organized in seven sections. Section 2 gives the general problem statement and assumptions. Section 3 describes the proposed SMGO algorithm. The convergence properties and optimality gap calculations are discussed in Section 4, with implementation notes in Section 5. Section 6 compares the performance of the proposed algorithm with other global optimization techniques on representative test functions, and Section 7 concludes this paper and provides insights on future directions.

Problem Statement
Consider a cost function z = f o (x), f o : X → R, where x ∈ X is the vector of decision variables, X ⊂ R D is a compact and convex polytope ("search set"). No analytical form of f o is assumed available. The a priori knowledge about f o is given by the following assumption: Being a Lipschitz continuous function on a compact set, f o presents a global minimum z * : Let us denote with the corresponding set of global minimizers. We further assume that it is possible to acquire information about f o by sampling (long function evaluation): Assumption 2. Given a vector of decision variables x ∈ X , it is possible to sample the cost function without noise: The set of collected data is denoted as: where n ∈ N is the number of data points and z (i) = f o (x (i) ). For simplicity, with a slight abuse of notation, we will also write x (i) ∈ X n when the pair (x (i) , z (i) ) belongs to the data set (3).
We denote with (x * n , z * n ) the best pair in X n , where (x * n , z * n ) = arg min If the result of (4) is not unique, a lexicographic criterion is used to sort the minimizers, and the first one is picked. We further denote with δ n the optimality gap: We assume that a starting number n 0 ≥ 1 of data points is also available, for example collected by the user in a first testing campaign, forming the set X n0 ⊂ X .
We can now state the problem addressed in this paper.
Problem 1. Design an algorithm that, under Assumptions 1-2, generates a sequence of points {x (n0+1) , x (n0+2) , . . .}, x (i) ∈ X , such that: Moreover, for a finite sequence length N , provide a method to compute a bound on the optimality gap δ N .

Set Membership Global Optimization (SMGO): Algorithm
SMGO addresses Problem 1 with a sequential procedure, as common in global optimization [1,13,18]. At each iteration n ≥ n 0 , the next test point x (n+1) is chosen with a strategy that exploits the prior knowledge on function f o , given by Assumption 1, and the current data set X n . In particular, a valid (i.e., consistent with data) estimate of the Lipschitz constant γ o is derived at n 0 and updated at each n > n 0 , in order to estimate upper and lower bounds on f o . SMGO then leverages the latter to compute x (n+1) , choosing between an exploitation mode and an exploration one.

Preliminaries: Lipschitz constant estimation, cost function bounds, and search mesh
At the start of each iteration (n > n 0 ), one long function evaluation z (n) = f o (x (n) ) is carried out, where x (n) is the test point selected at the previous iteration. This resulting new data point (x (n) , z (n) ) is added to the set of collected samples: At the first iteration, n = n 0 , an initial estimate of the Lipschitz constant, γ n0 , is computed as: If n 0 = 1, the Lipschitz constant estimate can be initialized as an arbitrarily small number. When n > n 0 , the minimum feasible Lipschitz constant estimate γ n is updated as follows: (8) Such an estimate is guaranteed to be compatible with the available data and prior information on f o , given Assumption 2 [20]. Moreover, γ n is monotonically increasing with n and convergence to the true Lipschitz constant γ o is guaranteed if the test points densely cover the search set [22]. This is the case for the SMGO algorithm as shown in Section 4.
Given the data set X n and Lipschitz constant estimate γ n we can compute the following lower and upper bounds on f o [20]: When µ = 1, the functions z n , z n represent, respectively, the minimum and maximum values that the objective function can take for a given vector of decision variables x ∈ X , given the information collected up to iteration n.
In SMGO, an (user-defined) overestimation parameter µ > 1 is used in (9)-(10) to compensate for the fact that the γ n estimated from the currently available data points is an underestimator of the true Lipschitz constant γ o .
The set is referred to as the feasible functions set [20] at iteration n. The functions z n (x) and z n (x) are intersections of hypercones (referred simply as cones in the remainder), each one with vertex at a tested point. A representation of these bounds is depicted in Fig. 1 for a function f o with scalar input argument. It should be emphasized that z n and z n are calculated from the data in X n , without any need for additional evaluations of f o . The uncertainty about the value of f o at the point x ∈ X is computed as: The SMGO algorithm proceeds in two steps. First, an exploitation routine searches for a local minimizer of z n . Then, an exploration step seeks a maximizer of λ n . Exact solutions for such subproblems can be computed by utilizing Hyperbolic Voronoi Diagrams (HVDs) [20,24], however at rather high computational cost. To increase computational efficiency, we thus introduce a mesh given by segments l i,j connecting any two points Such a mesh allows us to derive analytic solutions to the computation of the next test point x (n+1) , according to the two strategies presented next.

Exploitation (Mode θ)
Consider the best data point (x * n , z * n ) (4) at iteration n. Mode θ searches a point x n θ * with the highest predicted improvement w.r.t. x * n , according to the lower bound z n (9). Specifically, the exploitation routine considers the set of segments l n * ,i between x * n and any other sampled point x (i) ∈ X n to search for the candidate point x n θ * . Here, index n * is the data-point index of the pair (x * n , z * n ) in the data set X n .
We first evaluate the lower bounds on f o (x) : x ∈ l n * ,i considering only the hypercones defined by the data points (x * n , z * n ) and (x (i) , z (i) ), i.e.: and The intersection of these two cones on the segment l n * ,i is the smallest feasible value of f o (x) : x ∈ l n * ,i considering the information given by the two points (x * n , z * n ), (x (i) , z (i) ) and the Lipschitz constant estimate γ n . Such an intersection, denoted with x (i) θ , can be derived analytically on the basis of geometrical considerations (see Fig. 2). In fact, denoting it follows that By evaluating (16) for all possible segments, we obtain the set X where the constraint (17b) ensures that the point selected is such that the lower bound z n considering the information provided by all the sampled points coincides with the lower bound provided by the current best point (x * n , z * n ) alone. This usually happens in the vicinity of x * n , in particular at the border of the hyperbolic Voronoi cell pertaining to x * n , see [20]. The feasibility of problem (17) is always guaranteed, as shown in Lemma 1 in Section 4.
To decide whether x n θ * shall be eventually taken as the next test point x (n+1) for the long function, or the exploration routine shall be called instead, the following condition is evaluated: where αγ n is referred to as the expected improvement threshold and α ∈ [0, 1) is the SMGO exploitation parameter, adjusted by the user. If (18) is fulfilled (as in the example of Fig. 2), it means that the potential improvement provided by x n θ * is larger than the threshold. In this case, we set Otherwise, SMGO switches to exploration mode, described in the following subsection.

Exploration (Mode ψ)
In this mode, the candidate point x n ψ * ∈ X is chosen as the one with largest λ n value among selected candidates. Sampling x n ψ * will yield λ n+1 (x n ψ * ) = 0 (due to noiseless function evaluation, see Assumption 2), and the uncertainty in its vicinity will correspondingly decrease.
Similar to Mode θ, the search for x n ψ * considers the segments among the samples in X n , in order to limit the computational complexity. In this case, we consider as candidates all the midpoints between any pair (x (i) , x (j) ) ∈ X n . The reason is that midpoints generally feature larger uncertainty than other points on the segments, which are closer to one of the two extremes. Hence, we form the set of midpoints: and, according to the described rationale, the exploration routine picks the next test point as: Also in this case, if (20) does not have a unique solution, a lexicographic criterion is added to select the maximizer.

Algorithm summary
A summarized flow of the SMGO method is given as pseudo-code in Algorithm 1.

Algorithm 1: SMGO Algorithm
Input: Long function fo, initial data points X n 0 , maximum number of long function evaluations N , exploitation parameter α, overestimation parameter µ 1 while n0 < n ≤ N do // Long function evaluation, data update 2 Evaluate the long function fo at x (n) and measure output z (n) . Update the data set Update current best sample (x * n , z * n ) and Lipschitz constant γ n from X n // Exploitation (Mode θ) 4 Evaluate lower bounds on segments around the current best sample x * n to choose the exploitation point x Assign test point for next iteration Go to next iteration n ← n + 1 11 Final optimal point and value: get the best sample (x * N , z * N ) from the set X N

On the choice of α and µ
The SMGO parameter α affects only Mode θ: a smaller α value leads to higher tendency for exploitation. Regarding µ, using a value close to one is recommended, because a higher Lipschitz constant overestimator leads to x , which would not fully utilize the Set Membership approach in exploitation. At the same time, excessively large values of µ increase the conservativeness of λ n . Moreover, in the computation of . . , n gain a higher relative importance with respect to the corresponding values of z (i) (see (9)-(10)), so that the exploration mode will tend to select the midpoints that are farthest away from the available data-points.

Theoretical Properties
We now analyze the convergence and optimality properties of the SMGO algorithm. First, we prove that problem (17) is always feasible, thus showing that the exploitation routine always returns a valid candidate next point x n θ * . Lemma 1. At any iteration n, problem (17) admits at least one feasible point.
Proof. Take the pair (x (i) , z (i) ) such that and consider the corresponding candidate exploitation point x (i) θ given by (16). By construction we have that and we want now to prove that this point also coincides with z n (x (i) θ ), i.e. that: Assume now, for the purpose of contradiction, that: thus invalidating (22). Then, denoting with a = 1 2 (1 − si µγ n ) (see (16)), we would have: which contradicts equation (21). Then, (23) is false, while (22) holds true, implying that point x (i) θ satisfies constraint (17b), hence proving the result.

Convergence of SMGO algorithm
We next show that after finite iterations, ansuboptimal point is obtained, as required by Problem 1. For simplicity, we consider the following technical assumption, which in practice can be replaced by a vertex estimation mechanism presented in Section 5.1. Let us denote with V ∈ N the number of vertices of the polytope X , and with c (v) the v-th vertex.
We start by proving four lemmas that are instrumental to show SMGO convergence.
In the following, for a point x and positive scalar r, denote with B(x, r) the closed hyper-ball: Proof. At iteration n, denote with G n the subset of the search space where a candidate test point is accepted to perform an exploitation, i.e., G n x ∈ X : z n (x) < z * n − αγ n .
When a candidate point inside this set is chosen as x (n+1) in Mode θ, two cases arise: 1. The new sample does not change the current best point x * n . In this case, z (n+1) ≥ z * n .
Denote as B(x (n+1) , r θ ), the hyper-ball centered at x (n+1) with finite radius It follows that ∀x ∈ B(x (n+1) , r θ ) Therefore all the points inside the hyper-ball are not eligible for exploitation. Hence, we have: i.e., G diminishes by a finite amount, which is valid even when γ n updates, because γ n is bounded by γ o . 2. The new sample replaces the current best point x * n . In this case z (n+1) < z * n and therefore G n+1 ⊂ G n due to the new threshold. Moreover, the hyper-ball B(x (n+1) , r θ ), with r θ = α µ around the new sample is also removed from G.
Hence, G n = ∅ after finite iterations, and Mode θ will fail, proving the lemma.
Then, there exists at least one sample x (j) ∈ X n ∩ Q, i = j.
Proof. Let us assume, for the purpose of contradiction, that the lemma does not hold. That is, there is no sample in the open half-space Q. Due to Assumption 3, this would also apply to all vertices c (v) , v = 1, . . . , V of X . Therefore, we would have: Two cases could then arise: 1. Inequality in (24) applies strictly for at least one vertex c (v) . In this case, any such vertex c (v) andx would reside on opposite sides of the hyper-plane However, consider any w ∈ X . Since the latter is a convex polytope, we can write Hence, it would apply that which would mean that all w ∈ X reside on the opposite side of P with respect tox . This can only happen ifx / ∈ X , falling in a contradiction with the assumptions. (24) applies for all c (v) , v = 1, . . . , V . This would mean that all vertices (and consequently, all w ∈ X ) belong to the hyper-plane P containing x (i) . In turn, this would imply that eitherx / ∈ X (if x − x (i) > 0), orx = x (i) , both of which fall in contradiction with the assumptions.

Equality in
Since both cases 1. and 2. lead to contradiction, the claim of the lemma must be true, thus completing the proof. Lemma 4. Consider any unsampled pointx and any radius r such that X n ∩ B(x, r) = ∅. Then, given Assumption 2 and a Lipschitz constant estimate γ n obtained with (7)- (8), it applies that Proof. Consider the pointx and denote with (x (a) , z (a) ), (x (b) , z (b) ) ∈ X n the data pairs that determine the upper and lower bounds z n (x), z n (x), respectively, i.e: Then, considering that X n ∩ B(x, r) = ∅, that (µ − 1) > 0, and taking into account (7)-(8) we have: which proves the result.
Note that in Lemma 4 we are not assuming that the estimate γ n is larger than γ o : the result follows by how the estimate is computed, combined with the use of µ > 1.
Lemma 5. Let Assumption 3 hold, and assume that Mode ψ is undertaken infinitely often as n → +∞. Then, for anyx ∈ X and any σ > 0, ∃n σ < ∞ such that min Proof. Consider any pointx ∈ X . Ifx ∈ X n for some n < ∞, we have min x (i) ∈X n x (i) −x = 0 and the claim is trivially proven. Consider then the casê x / ∈ X n , ∀n ∈ [1, +∞). Pick the nearest samplē x n = min and consider the open half-space: Now, we select samplex n x n = arg min Such a sample is always guaranteed to exist, in force of Lemma 3. We note that due to our selection criterion forx n , the set does not have any sample in its interior, otherwise any such sample should have been selected asx n in (28 Furthermore, we consider the set of points not belonging to H, whose uncertainty is larger than that ofm R n (m) = x / ∈ H : λ n (x) > λ n (m) .
We then sort all candidate midpoints in the set for some x (i) , x (j) ∈ X n (see (19)) and M = n(n − 1) 2 is the total number of segments at iteration n. In Mode ψ, the first-ranking candidate m n 1 is taken for sampling. Its uncertainty is necessarily larger than λ n (m) (otherwisem would have been ranked higher), i.e. m n 1 ∈ R n (m). Assume, for the purpose of contradiction, that points in H are never sampled.
At iteration n + 1, the midpoint m n 1 / ∈ H is thus sampled and its uncertainty λ n+1 (m n 1 ) becomes zero. Since this point does not belong to H, the value of h in (29) is still valid, while the uncertainty bound (30) may either be the same, or increase in case γ n+1 > γ n (see (8)). Moreover, for any point x ∈ X we have as shown by direct application of (9)- (11). Consider now the ball B(m n 1 , r λ ), where On the basis of (30), (31) and (32) we have that: This means that the volume of set R n (m) diminishes by a finite quantity at each iteration, since the value of h does not depend on n; and becomes null after finite iterations, since Mode ψ is assumed to be undertaken infinitely often. Denote withn − 1 < ∞ the iteration at which this happens. Then, no midpoint outside the set H would have uncertainty larger than λ n (m), so that a midpoint will be ranked first in M n−1 and sampled at iterationn, thus falling in contradiction. We thus demonstrated that, in finite iterations, a point inside H is sampled. Moreover, denoting such a point as x (n) , by construction we have Then, two cases may occur: either x (n) is closer tô x thanx n , or not. In the former case, repeating the process from (27) we'll have thatx n = x (n) andx n =x n . In the latter case,x n =x n and x n = x (n) . In both cases, the radius of the halfhyper-ball H shrinks because of (34). By induction we thus have: lim n→∞ x n −x n = 0.
Next, still consideringx n as the nearest sample tox (27), we are going to show that there exists a radius a > 0 such that To this end, consider any point x ∈ Q and assume, for the sake of contradiction, that (36) never holds. This would mean that the nearest point tô x on the segment connectingx n to x is always x n , no matter how close the two points are. In turn, this would imply that: i.e., which would imply that x / ∈ Q, thus falling in a contradiction. Hence, the directional derivative (37) must be negative for all x ∈ Q, meaning that there exists a scalar a > 0 such that in the neighborhood Q ∩ B(x n , a) the property (36) holds.
Combining (35) and (36), we obtain that, after a finite number k of iterations, the following condition holds: This implies that lim n→∞ x n −x = 0, and, for any σ > 0, we finally have ∃n σ : min thus proving the lemma.

Theorem 1. Let Assumptions 1 and 3 hold. Then,
Proof. Consider a pointx ∈ X * (see (2)) and take σ = γ o . For any n, denotē x n = arg min In virtue of Lemma 2, the exploration mode will be called infinitely often. Now, by applying Lemma 5, we have that ∃n σ < ∞ : Then, in virtue of Assumption 1 we have: Thus proving the result with n = n σ .

Optimality gap
Considering a run of the SMGO algorithm with N iterations, we now analyze the gap between the best sampled value and the actual global minimum, denoted as the optimality gap. In calculating such a difference, the SM-based guarantees z N (x) ≤ f o (x) ≤ z N (x) can be utilized in the case of a known Lipschitz constant γ o . Hence, the upper bound of the optimality gap δ N can be considered as the difference between z N * and the minimum lower bound z N defined in (9) generated by data set X N , i.e.
The calculation of min x∈X z N (x) in (39) entails calculation of the hyperbolic Voronoi cells and of the lower bounds at their vertices. Denoting with V N the set of such vertices generated from X N , (39) can be expressed as The acquisition of HVD vertices using the samples can be approached as described in [24], which is referred for more information. The calculation of HVD cells and vertices entail computations with exponential complexity w.r.t. D [24]. However, an interested user can obtain an upper bound on the optimality gap by carrying out (40).

Computational Aspects
This section presents an analysis of the computational aspects of the proposed algorithm and possible improvements. The first aspect discussed is on ensuring the coverage of X , followed by a discussion on the computational complexity.

On the coverage of the search set
The selection method for candidate points described in Section 3 allows SMGO to only explore the volume inside the convex hull of the sampled points in X n0 . Hence, in order to ensure exploration throughout all the search set X , all its vertices c (v) should also be considered in the candidate points generation, as stated in Assumption 3. One approach is to first evaluate f o on all c (v) before actually starting with exploitation/exploration. However, in the case of a hyperrectangle, this entails an exponential number of long function calls w.r.t. dimensionality D.
To address this issue, a simple approach is proposed in order to ensure coverage in X without explicitly calling f o to evaluate the corners. Moreover, with this approach Assumption 3 is not needed anymore for the theoretical results to hold. At each iteration n, the midpoints between any sample and all the vertices of X are included as candidate exploration samples, even if the vertices do not belong to X n . Moreover, for the purpose of computing the uncertainty of each midpoint, the cost function value f o (c (v) ) at the v-th vertex c (v) is estimated as: i.e., equal to the cost function value at the nearest sampled point. The concept of this strategy is shown in Fig. 3. Thanks to Assumption 1 and Lemma 5, with an increasing number of sample points in X the estimate (41) approximates more and more accurately the true function value at the vertex.

On the computational burden
There are three main routines within an iteration of SMGO: Lipschitz constant update, exploitation by Mode θ routine, and eventually, if the expected improvement threshold is not achieved, exploration by Mode ψ. The worst-case SMGO complexity analysis assumes that each iteration eventually proceeds to Mode ψ.
The first routine compares the incoming n-th sample with n − 1 existing ones, resulting in O(n) complexity.
The Modes θ and ψ can be designed as iterative routines as well to significantly improve the computational burden, as discussed in the following.
The exploitation Mode θ comprises the selection of n − 1 candidate points x (i) θ for each segment l n * ,i from x * n to all the other sampled points according to (15)- (16), and the calculation of z n for each x (i) θ using (9). Both are O(n) operations, leading to a compounded complexity of O(n 2 ). However, when γ n and x * n do not change from iteration n − 1 to n, a cache can store x (i) θ locations and corresponding z n . From (9), it is convenient to In practice, information are saved 1 about the hypercone that actually generated z n for each x (i) θ . As long as x * n stays the same, previously existing x (i) θ locations do not need to be computed again. Furthermore, calculating z n for an existing exploitation point x θ − x (n) for an incoming sample (x (n) , z (n) ).
Note that when z n (x (i) θ ) changes due to (42), . As a result, computing z n for all existing candidate points is O(n), i.e. O(1) for each. Now, consider the situation when γ n is updated, assuming that x * n andk do not change. This implies rescaling ∆ (i) such that∆ 1 It can be understood that z (k) is the value at the tip, while∆ (i) is the height of the hypercone. and updating z n with the new sample (x (n) , z (n) ) using (42), still resulting in O(1) per candidate point. This implies a Mode θ aggregate complexity of O(n).
Finally, one new candidate point x (n) θ is evaluated using (9), which is O(n). Hence, as long as x * n does not change, the introduction of a cache reduces Mode θ to O(n). However, when x * n changes, i.e. x * n = x (n) , computing all x (i) θ locations must be repeated from the start, leading to O(n 2 ) worst-case complexity for Mode θ as stated before.
As defined, Mode ψ generates n(n − 1) 2 candidate points (midpoints among all existing samples) at each iteration. Furthermore, the calculation of λ n for each candidate point is O(n), resulting in O(n 3 ) for each exploration iteration. Note that behavior in Mode ψ does not depend on the current best point, because n new segments are added to the search set at iteration n. A cache solution can be introduced also in this mode, updating λ n -related values iteratively.
For each existing x (i,j) ψ , the values to be saved are z (k) and∆ (i,j) (for z n ), and z (ǩ) and∆ (i,j) (for z n ), wherek andǩ are defined aŝ The iterative updates are performed using eq. (42) and (43). This then results in O(1) for each existing midpoint, and O(n) for each new midpoint generated by new segments from incoming x (n) to the n − 1 existing points. This results in O(n 2 ) for Mode ψ, and in turn, a O(n 2 ) worst-case complexity for SMGO. Remark 1. Note that when γ n is updated (which always implies a growth of its value), the proposed iterative implementation might result in more conservative bounds z n and z n than if repeatedly recalculated at every iteration, i.e. lower z n for Mode θ, and larger λ n in Mode ψ. This occurs because with higher γ n -values the HVD layout increasingly approaches that of a non-hyperbolic Voronoi tessellation, i.e. the cones providing the tightest bounds at any x ψ . These effects are not accounted for by the iterative update procedure, which assumes that the bound-generating cones remain the same. However, the theoretical properties of SMGO are still valid with the iterative implementation. In fact, all the arguments for Lemma 2 are valid also with iteratively-computed lower bounds. Regarding the exploration mode, inequality (26) in Lemma 4 is consistent with the more conservative iteratively-computed γ n . Furthermore, when sampling a new point m n 1 , for which λ n+1 (m n 1 ) is exactly computed, (33) still holds. Hence, all arguments laid out in Lemma 5 still hold.

Performance Test Results
In this section, the proposed SMGO algorithm is evaluated on well-known benchmarks, usually employed in black-box optimization, comparing its performance with representative Lipschitz-based methods: DIRECT and AdaLIPO. In addition, the Bayesian optimization approach (BayesianO) is also considered in the tests.
The results are discussed in terms of iteration-based optimization performance and computational times. All scripts used for this section are available at https://github.com/lorenzosabugjr/smgo.

Test parameters
Seven test functions of varying structure and dimensions are chosen for this test, covering a variety of characteristics relevant to black-box optimization benchmarks. Their search bounds, optimal values and relevant features are summarized in Table 1. A reasonable assumption in real-life applications is that the black-box function has no available optimal value. Hence, each algorithm is given a budget of N = 500 long function evaluations for each optimization run. The optimal results after the evaluation budget are then measured and compared among the different algorithms, in what is referred as a fixed-cost/budget comparison [25]. All computations were performed in MATLAB 2020b on a system with AMD Ryzen 9 3900X (3.80 GHz) and 32 GB RAM.
To assess the sensitivity of each algorithm to different initial information and randomized decisions, 100 independent runs (trials) were performed. For each trial, the same randomly-generated starting point is given to AdaLIPO, SMGO, and BayesianO.
On the other hand, DIRECT considers a fixed set of search points, not allowing to randomize the starting samples.
SMGO is implemented as described in Algorithm 1 with parameters α = 0.001 and µ = 1.025. AdaLIPO is implemented as described in [18] with parameter p = 0.1. DIRECT is implemented as described in [14], not requiring any tuning parameter. Finally, the algorithm available in the Matlab Statistics and Machine Learning Toolbox is employed for the BayesianO implementation, using the provided default parameters.

Iteration-based performance
The evolution of the best sample z * n along the iterations (n ∈ [50, 500]) for each algorithm is shown in Figs. 4 and 5 for the Rosenbrock and Deb's functions, respectively. The graphs show the distribution of the best sampled cost for the 100 independent runs. Note that the results of the DI-RECT algorithm do not show dispersion due to the fact that it is a deterministic search method. DI-RECT follows a batch sampling strategy, hence the graphs show for each n value, the best results for the previously-sampled batches.
The Rosenbrock function is a standard test case for gradient-based algorithms, with a unique global minimum (unimodal) for D = 2; however it has 2 minima for D = 4 ∼ 30 [26]. The location of x * is at (1, 1, . . . , 1), hence the classical search bounds −10 ≤ x i ≤ 10 would be advantageous for DIRECT, which invariably samples the center of the search hyper-box. For this reason, the search bounds are changed to move x * away from the center.
As seen in Fig. 4, DIRECT has only improved marginally the initial sampled cost after the 500 iterations. Also AdaLIPO achieved marginal improvements after 500 iterations. This might be associated to the completely randomized exploration and exploitation decisions. On the other hand, both SMGO and BayesianO achieved large improvements during the first 100 iterations, while slight enhancements are observed during the final 400 iterations. SMGO exhibits a wider distribution among different runs and, on average, achieves a better result after 500 iterations.
In the case of Deb's function #1, as shown in Fig. 5, the presence of numerous global minima causes a bad performance for DIRECT, which barely improved its initial best cost. On the other   hand, AdaLIPO, BayesianO, and SMGO have resulted in similar results, achieving very similar average costs after the iteration limit, while BayesianO shows more concentrated values among runs. In this case the randomized mechanism of AdaLIPO allows it to discover low-cost regions to improve the best sample. Furthermore, Bayesian's usage of a response surface, and SMGO's geometry-based candidate points derivation have become useful in looking for low-cost regions in the search space. A summary (mean best values) of the algorithms' outcomes after 500 iterations are shown in Table 2. Note that the results tabulated for DIRECT were achieved at iterations which are a bit larger than the pre-defined maximum N = 500, owing to its batch-based samples generation. In most test cases, BayesianO achieved the best average results, while SMGO and DIRECT attained the best average performance in 2 cases each. On the other hand, AdaLIPO achieved fair results, even if it ranked second-place or worse for most of the tested functions.
Furthermore, for each test function (with the same starting points across all algorithms), the number of trials (out of 100) resulting in best results for each respective algorithm is shown in Table 3. As with the average performance, BayesianO achieved the most best trials in 9 out of the 13 tested functions. Furthermore, in Deb's #1 (5D), Schwefel (5D and 10D), and Salomon (5D and 10D), it resulted in best results in almost all of the given trials. On the other hand, SMGO had the most best results in Rosenbrock and Deb's #2 (10D), and had the second largest count of best trials in Deb's #1 (10D) and Deb's #2 (5D). AdaLIPO registered the second count of best trials in Deb's #2 (10D), and had no best trials for 4 functions. Lastly, DIRECT had no best trials for 8 out of 13 test functions, but took almost all the best results in the Brown function (5D and 10D).
We have also compared the outcomes of SMGO with the competitor methods using the Wilcoxon and Kruskal-Wallis non-parametric statistical tests, pairing the results by the shared random initial point. Both test were performed with 5% significance level and with the null hypothesis that SMGO results are statistically similar to those from the respective competitor method. The two tests agreed that SMGO results were statistically similar to AdaLIPO for Styblinski-Tang (5D), Deb's #2 (5D, 10D), Schwefel (10D), and Brown (10D). Furthermore, similarities were found between SMGO and BayesianO for Deb's #2 (10D). Given these results, SMGO is found to have competitive results compared with the considered global optimization algorithms on the test functions employed in the analysis.

Computational time
The computational time required to complete an iteration for each considered algorithms is analyzed for the case of the 5D Deb's function #1 optimization. Fig. 6 (in log scale) shows the time required to complete an iteration as n increases. The mechanism of DIRECT calculates recursive hyperboxes to be sampled by batches, hence the computational time per iteration is not well defined, therefore results are not shown in the plot. For AdaLIPO an almost flat line is observed, with 10 to 100 µs per iteration, implying a computational time independent of the iteration number. SMGO has resulted in a polynomial complexity w.r.t. iterations, which is expected from the computational analysis presented in Section 5. Finally, BayesianO is the most demanding algorithm, requiring computational times per iteration between 2 and 3 orders of magnitude larger than SMGO in this test case.
The computational times observed for all the considered test functions have shown trends similar to those on the 5D Debs function #1. Table 4 summarizes the results, showing the total computational time required to complete an optimization run. The AdaLIPO algorithm, based on a simple mechanism of randomized generation of exploitation and exploration points, always shows very fast computational times, not affected by the dimensionality of the decision variable, scarcely affected by n, and taking 0.01-0.02 s per optimization run for all the test functions.
From the tests, BayesianO requires the largest computational times, which achieved around 4-5 s per iteration at n = 500 for most of the 10 dimensions functions. On the other hand, SMGO exhibits an intermediate per-iteration time around 100 µs to 10 ms for 5 dimensions functions and tens to hundreds of ms for 10 dimensions functions. Furthermore, for SMGO the total optimization run time when applied to 10 dimensions functions is around 10 times those employed for 5 dimensions functions. However, this is significantly faster than BayesianO, ranging from around 15 times shorter time for 10D functions, to 70 times shorter time for 5D functions.
The observed computational times, combined with the competitive optimization results presented Test function 9.78 E -1 9.61 E -1 1.65 E -1    in Figs. 4-5, and Table 2, indicate that the proposed SMGO algorithm shows a good trade-off between optimization performance and computational speed, compared to the other considered methods.

Conclusions and Further Work
In this work a sequential algorithm for global optimization of black-box functions has been proposed. The SMGO algorithm assumes a Lipschitz- continuous cost function and is based on a nonlinear Set Membership function approximation. The selection of the test points where the black-box function is evaluated is carried out by solving simplified geometry-based problems on the guaranteed lower bounds or uncertainty intervals of the function, provided by a Set Membership model. It is shown that the SMGO algorithm finds a suboptimal solution with desired tolerance after finite iterations and converges asymptotically to the global minimum. Furthermore, the computational complexity of the algorithm is discussed and a cache-based solution is introduced to improve the computational burden. The proposed method is evaluated against other representative global optimization methods on several test functions, comparing the quality of the best solutions found after a fixed number of calls to the black-box function. The results show the competitiveness of the SMGO algorithm compared to state of the art methods found in literature. Ongoing research aims to extend the algorithm to include constraints and its evaluation on experimental applications.