# A Radial Basis Function-Based Optimization Algorithm with Regular Simplex Set Geometry in Ellipsoidal Trust-Regions

Authors: Tom Lefebvre; Frederik De Belie; Guillaume Crevecoeur
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, we investigate two ideas in the context of the interpolation-based optimization paradigm tailored to derivative-free black-box optimization problems. The proposed architecture maintains a radial basis function interpolation model of the actual objective that is managed according to a trust-region globalization scheme. We focus on two distinctive ideas. Firstly, we explore an original sampling strategy to adapt the interpolation set to the new trust region. Abetter-than-linear interpolation model is guaranteed by maintaining a well-poised supporting subset that pursues a near regular simplex geometry of n+1 points plus the trust-region center. This strategy improves the geometric distribution of the interpolation points whilst also optimally exploiting the existing interpolation set. On account of the associated minimal interpolation set size, the better-than-linear interpolation model will exhibit curvature, which is a necessary condition for the second idea. Therefore, we explore the generalization of the classic spherical to an ellipsoidal trust-region geometry by matching the contour ellipses with the inverse of the local problem hessian. This strategy is enabled by the certainty of a curved interpolation model and is introduced to accounts for the local output anisotropy of the objective function when generating new interpolation points. Instead of adapting the sampling strategy to an ellipsoid, we carry out the sampling in an affine transformed space. The combination of both methods is validated on a set of multivariate benchmark problems and compared with ORBIT.


---

## Body

## 1. Introduction

We consider the following optimization problem with the objective functionf:X⊂ℝn⟶ℝ. We further assume that f is continuously differentiable, smooth, and bounded from below and that ∇f is Lipschitz continuous on the subset X. We restrict the minimization to a convex subset X⊂ℝn, however, we assume that none of these constraints are active at the minimum so that (1) can be considered unconstrained in practice (the subset X is mainly introduced to emphasize that the search is limited to a range of practical interest. It may also be a means to delimit the search space to the definition domain of the function).(1)minx∈X⊂ℝnfx.Such problems are common in engineering, and consequently, a large family of search algorithms has been tailored to solve such problems. To motivate the development of derivative-free optimization methods, we premise additional assumptions on the computational tractability of the objective.As is the case in many engineering design applications, we assume that any evaluation of the objective is associated with the execution of a complex simulation or legacy code. Such evaluations are time-demanding and usually render local derivative information inaccessible. These assumptions rule out the use of the members of the family of gradient-based optimization algorithms that exploit the first- and second-order derivative information to search for points that satisfy the first- and second-order optimality conditions. Under these conditions, it is also intractable to approximate derivative information indirectly using, for example, a finite difference scheme, given the tight evaluation budget that is available. In short, these assumptions imply that the search algorithm has to rely exclusively on direct evaluations of the objective function to probe the local geometry of the problem. The denominator derivative-free optimization (DFO) refers to any search algorithm that uses only the direct evaluations of the objective function.


### 1.1. Overview of DFO Methods

In this family, there are two main classes to be recognized, which are as follows: stochastic search methods and deterministic search methods. Either can be further subdivided into population-based and evolutionary strategies and direct search methods and interpolation-based search methods, respectively. The distinction between stochastic or deterministic here refers to the property that, excluding therandom initialization of the algorithm, the execution and convergence history is the same for deterministic algorithms whilst being subject to variation for stochastic search algorithms.Generally speaking, stochastic search algorithms maintain a belief over the solution space representing possible candidate solutions. This belief is either represented with an actual population or with a parametric distribution, e.g., the multivariate normal distribution. Members of the latter class are denoted as evolutionary strategies (ESs). Stochastic search methods proceed with the evaluation of any new members of the population and use stochastic operations to generate a new population. In the past few decades, large varieties of population-based algorithms have been invented that are inspired by the swarm intelligence exhibited by animal populations. Examples are monarch butterfly optimization [1], earthworm optimization algorithm [2], elephant herding optimization [3], moth search algorithm [4], slime mould algorithm [5], Harris hawks optimization [6], artificial bee colony optimization [7], and grey wolf optimization [8], amongst others, which, in some way or another, are all related to the classic particle swarm optimization algorithm [9]. These algorithms have found applications in structural design optimization mostly and are particularly appealing to find global solutions at the cost of sampling the objective function abundantly. Instead of a population, the class of evolutionary search algorithms maintains a parametric distribution, which is used to spawn a sample set, which, in turn, is evaluated and assessed to infer the next parametric belief. The benefit of these evolutionary strategies is that they may be derived from much more principled arguments. Again, their main appeal is that they solve for a global solution and are particularly well-equipped to treat noisy or multi-global objective functions, where other algorithms, even population-based algorithms, are prone to get stuck in local minimal. Examples popular in the machine learning and reinforcement learning communities are the covariance-matrix adaption evolutionary strategy [10–13], natural evolutionary strategies [14–17], and information geometric optimization [18]. On account of their inherent stochasticity, the execution of stochastic search methods may be subject to excessive covariance, and their actual performance may differ for each individual run. Basically, they battle noise with noise.Deterministic search methods eliminate that problem. The subclass of direct search methods encompasses the historical Nelder-Mead simplex method [19] and modern pattern search methods [20], such as mesh adaptive direct search (MADS) methods and related variants [21–24]. These direct methods are characterized by the fact that they try to satisfy the zeroth-order optimality criterion by comparing the ordering of evaluated points but usually do not exploit any directional preferences implied by them. The advantage here is that they will not be misled by false positives or false negatives (e.g., a strong objective slope that suddenly flattens) but consequently can also not exploit, or at least to a lesser extent, true negatives and true positives (e.g., a strong objective slope that actually points into the direction of the minimum). The other subclass of direct DFO methods can be described as interpolation-based optimization (IBO) methods. The basic idea here is to approximate the objective function by fitting an interpolation model to the available objective function evaluations. The hypothesis is that the underlying objective function is subjected to tendencies that are picked up by the interpolation model, which, in turn, can be exploited to explore the potential regions of interest. Within this scope, there are basically two main branches that can be classified based on the fact whether they maintain a global or a local interpolation function. The global branch was founded when the efficient global optimization (EGO) algorithm was first published [25]. Members of this class maintain a so-called acquisition function that can be optimized to determine the next best point to be evaluated and added to the interpolation set. A large family of methods has been developed based on this idea. The other branch of strictly deterministic interpolation-based methods is based on the trust-region framework and will be discussed in the following subsection.In conclusion, we note that the principal idea of interpolation-based optimization, i.e., replacing the objective function with an interpolation-based proxy, has been subject to cross-fertilization with ideas from the population-based stochastic search family and are commonly referred to as surrogate-assisted stochastic search methods. As opposed to the acquisition function used in EGO and related methods, here, a so-called batch infill criterion is used that can generate a whole set of points instead of a single point at once. Examples are TLRBF [26], RBFBS [27], SATLBO [28], SAGWO [29], DSCPSO-EMM [30], amongst others. Further exposition is out of scope.


### 1.2. Trust-Region IBO and Contributions

This article focuses on the particular deterministic DFO subclass of trust-region-based IBO methods, or TRIBO for short [31–34]. As opposed to the global IBO methods surveyed above, TRIBOs target a compact neighborhood in each iteration. Such a strategy limits the required amount of sample points to obtain a reliable surrogate in the compact neighborhood considered during that iteration [35]. A standard gradient-based optimizer can then be used to solve the local subproblem confined to the compact neighborhood.This analysis allows one to identify three key components. A framework is required to maintain a sequence of compact neighborhoods. Secondly, a strategy should be present to populate the new compact neighborhood with interpolation points. In turn, this interpolation set can be used to determine a new local interpolation model that is sufficiently accurate in the compact neighborhood. In this work, we will make use of the radial basis functions (RBF) interpolation framework [36–38]. The first two components are discussed next.The trust-region globalization scheme offers a framework to generate a sequence of compact neighborhoods that will be referred to as trust-regions henceforth. When the surrogate satisfies Taylor-like error bounds on its function value and gradient for every iteration, convergence to first-order optimality can be guaranteed [39]. In the context of interpolation, such error bounds can be established, depending on the geometry of the interpolation set used [40]. The geometry of an interpolation set can be understood as the distribution of points with respect to a certain region. The quality of geometry is quantified by a set-dependent geometry metric. An arbitrary set of points is then said to be well-poised for interpolation within the specified region if the corresponding geometry measure is sufficiently large.To guarantee a valid interpolation model suited for optimization within a specified trust region, IBO methods maintain a subset of the available interpolation set for which the geometry constant is maximized. Consequently, the well-poisedness of the subset geometry is assured, and so is the model validity. In every iteration, this constant is maximized by removing points from the interpolation set and generating new ones according to the current size and location of the trust region whilst limiting the overall accumulation of objective evaluations. This delicate mechanism is referred to as interpolation set management.In [32, 41, 42], such a geometry measure was described using a basis of polynomials. Other authors have described this constant in the function of the QR pivots of the Vandermonde matrix determined by the interpolation points [31, 40]. The mentioned methods consider a subset of n+1 interpolation points to guarantee a valid model. In any of these cases, the interpolation set management pursues orthogonality between n interpolation points relative to the current trust-region center. The center then serves as the n+1th interpolation point.The contribution of this article is twofold. Firstly, in our algorithm, we pursue a subset ofn+2 instead of n+1 interpolation points, as seen in [31, 32]. An original sampling strategy is devised that is inspired by the volumetric extremity properties of the regular n-simplex. This strategy has two benefits. For one, we achieve an improved interpolation set distribution, which is heuristically justified for low-dimensional problems. Then, considering that we use the RBF framework to construct a local interpolation model, it is guaranteed that the surrogate will exhibit curvature.Secondly, we propose to leverage this curvature information. We, therefore, generalize the traditional spherical trust-regions to ellipsoidal trust regions. The motivation to use an ellipsoidal geometry is again twofold. An ellipsoidal trust-region will adapt to the geometry of the objective function and will, therefore, allow larger steps toward more promising regions compared to using spherical geometry, i.e., the trust-region boundaries will align with the contours of the objective function, and the objective will look like a spherical problem. This idea is already leveraged in the context of evolutionary computation [10, 43]. Secondly, our interpolation set management is input-space oriented, meaning it only takes into account the local distribution of the interpolation points and does not consider the associated output values. The introduction of an ellipsoidal trust-region, whose shape reflects the objective’s response behavior, allows to account for the output-space when executing the interpolation set management and will, therefore, benefit interpolation-based metamodeling.


### 1.3. Outline

To keep this article self-contained, we provide a brief review of the general trust-region framework for surrogate management in Section2. In Section 3, we discuss the radial basis function interpolation framework and provide details on the implementation. Section 4 motivates our choice for using n+2 interpolation points and introduces our interpolation set management. In Section 5, we discuss the transition toward ellipsoidal trust regions. In Section 6, we validate our methods empirically on a number of benchmark problems.


## 1.1. Overview of DFO Methods

In this family, there are two main classes to be recognized, which are as follows: stochastic search methods and deterministic search methods. Either can be further subdivided into population-based and evolutionary strategies and direct search methods and interpolation-based search methods, respectively. The distinction between stochastic or deterministic here refers to the property that, excluding therandom initialization of the algorithm, the execution and convergence history is the same for deterministic algorithms whilst being subject to variation for stochastic search algorithms.Generally speaking, stochastic search algorithms maintain a belief over the solution space representing possible candidate solutions. This belief is either represented with an actual population or with a parametric distribution, e.g., the multivariate normal distribution. Members of the latter class are denoted as evolutionary strategies (ESs). Stochastic search methods proceed with the evaluation of any new members of the population and use stochastic operations to generate a new population. In the past few decades, large varieties of population-based algorithms have been invented that are inspired by the swarm intelligence exhibited by animal populations. Examples are monarch butterfly optimization [1], earthworm optimization algorithm [2], elephant herding optimization [3], moth search algorithm [4], slime mould algorithm [5], Harris hawks optimization [6], artificial bee colony optimization [7], and grey wolf optimization [8], amongst others, which, in some way or another, are all related to the classic particle swarm optimization algorithm [9]. These algorithms have found applications in structural design optimization mostly and are particularly appealing to find global solutions at the cost of sampling the objective function abundantly. Instead of a population, the class of evolutionary search algorithms maintains a parametric distribution, which is used to spawn a sample set, which, in turn, is evaluated and assessed to infer the next parametric belief. The benefit of these evolutionary strategies is that they may be derived from much more principled arguments. Again, their main appeal is that they solve for a global solution and are particularly well-equipped to treat noisy or multi-global objective functions, where other algorithms, even population-based algorithms, are prone to get stuck in local minimal. Examples popular in the machine learning and reinforcement learning communities are the covariance-matrix adaption evolutionary strategy [10–13], natural evolutionary strategies [14–17], and information geometric optimization [18]. On account of their inherent stochasticity, the execution of stochastic search methods may be subject to excessive covariance, and their actual performance may differ for each individual run. Basically, they battle noise with noise.Deterministic search methods eliminate that problem. The subclass of direct search methods encompasses the historical Nelder-Mead simplex method [19] and modern pattern search methods [20], such as mesh adaptive direct search (MADS) methods and related variants [21–24]. These direct methods are characterized by the fact that they try to satisfy the zeroth-order optimality criterion by comparing the ordering of evaluated points but usually do not exploit any directional preferences implied by them. The advantage here is that they will not be misled by false positives or false negatives (e.g., a strong objective slope that suddenly flattens) but consequently can also not exploit, or at least to a lesser extent, true negatives and true positives (e.g., a strong objective slope that actually points into the direction of the minimum). The other subclass of direct DFO methods can be described as interpolation-based optimization (IBO) methods. The basic idea here is to approximate the objective function by fitting an interpolation model to the available objective function evaluations. The hypothesis is that the underlying objective function is subjected to tendencies that are picked up by the interpolation model, which, in turn, can be exploited to explore the potential regions of interest. Within this scope, there are basically two main branches that can be classified based on the fact whether they maintain a global or a local interpolation function. The global branch was founded when the efficient global optimization (EGO) algorithm was first published [25]. Members of this class maintain a so-called acquisition function that can be optimized to determine the next best point to be evaluated and added to the interpolation set. A large family of methods has been developed based on this idea. The other branch of strictly deterministic interpolation-based methods is based on the trust-region framework and will be discussed in the following subsection.In conclusion, we note that the principal idea of interpolation-based optimization, i.e., replacing the objective function with an interpolation-based proxy, has been subject to cross-fertilization with ideas from the population-based stochastic search family and are commonly referred to as surrogate-assisted stochastic search methods. As opposed to the acquisition function used in EGO and related methods, here, a so-called batch infill criterion is used that can generate a whole set of points instead of a single point at once. Examples are TLRBF [26], RBFBS [27], SATLBO [28], SAGWO [29], DSCPSO-EMM [30], amongst others. Further exposition is out of scope.


## 1.2. Trust-Region IBO and Contributions

This article focuses on the particular deterministic DFO subclass of trust-region-based IBO methods, or TRIBO for short [31–34]. As opposed to the global IBO methods surveyed above, TRIBOs target a compact neighborhood in each iteration. Such a strategy limits the required amount of sample points to obtain a reliable surrogate in the compact neighborhood considered during that iteration [35]. A standard gradient-based optimizer can then be used to solve the local subproblem confined to the compact neighborhood.This analysis allows one to identify three key components. A framework is required to maintain a sequence of compact neighborhoods. Secondly, a strategy should be present to populate the new compact neighborhood with interpolation points. In turn, this interpolation set can be used to determine a new local interpolation model that is sufficiently accurate in the compact neighborhood. In this work, we will make use of the radial basis functions (RBF) interpolation framework [36–38]. The first two components are discussed next.The trust-region globalization scheme offers a framework to generate a sequence of compact neighborhoods that will be referred to as trust-regions henceforth. When the surrogate satisfies Taylor-like error bounds on its function value and gradient for every iteration, convergence to first-order optimality can be guaranteed [39]. In the context of interpolation, such error bounds can be established, depending on the geometry of the interpolation set used [40]. The geometry of an interpolation set can be understood as the distribution of points with respect to a certain region. The quality of geometry is quantified by a set-dependent geometry metric. An arbitrary set of points is then said to be well-poised for interpolation within the specified region if the corresponding geometry measure is sufficiently large.To guarantee a valid interpolation model suited for optimization within a specified trust region, IBO methods maintain a subset of the available interpolation set for which the geometry constant is maximized. Consequently, the well-poisedness of the subset geometry is assured, and so is the model validity. In every iteration, this constant is maximized by removing points from the interpolation set and generating new ones according to the current size and location of the trust region whilst limiting the overall accumulation of objective evaluations. This delicate mechanism is referred to as interpolation set management.In [32, 41, 42], such a geometry measure was described using a basis of polynomials. Other authors have described this constant in the function of the QR pivots of the Vandermonde matrix determined by the interpolation points [31, 40]. The mentioned methods consider a subset of n+1 interpolation points to guarantee a valid model. In any of these cases, the interpolation set management pursues orthogonality between n interpolation points relative to the current trust-region center. The center then serves as the n+1th interpolation point.The contribution of this article is twofold. Firstly, in our algorithm, we pursue a subset ofn+2 instead of n+1 interpolation points, as seen in [31, 32]. An original sampling strategy is devised that is inspired by the volumetric extremity properties of the regular n-simplex. This strategy has two benefits. For one, we achieve an improved interpolation set distribution, which is heuristically justified for low-dimensional problems. Then, considering that we use the RBF framework to construct a local interpolation model, it is guaranteed that the surrogate will exhibit curvature.Secondly, we propose to leverage this curvature information. We, therefore, generalize the traditional spherical trust-regions to ellipsoidal trust regions. The motivation to use an ellipsoidal geometry is again twofold. An ellipsoidal trust-region will adapt to the geometry of the objective function and will, therefore, allow larger steps toward more promising regions compared to using spherical geometry, i.e., the trust-region boundaries will align with the contours of the objective function, and the objective will look like a spherical problem. This idea is already leveraged in the context of evolutionary computation [10, 43]. Secondly, our interpolation set management is input-space oriented, meaning it only takes into account the local distribution of the interpolation points and does not consider the associated output values. The introduction of an ellipsoidal trust-region, whose shape reflects the objective’s response behavior, allows to account for the output-space when executing the interpolation set management and will, therefore, benefit interpolation-based metamodeling.


## 1.3. Outline

To keep this article self-contained, we provide a brief review of the general trust-region framework for surrogate management in Section2. In Section 3, we discuss the radial basis function interpolation framework and provide details on the implementation. Section 4 motivates our choice for using n+2 interpolation points and introduces our interpolation set management. In Section 5, we discuss the transition toward ellipsoidal trust regions. In Section 6, we validate our methods empirically on a number of benchmark problems.


## 2. Trust-Region Methods

In this section, we give a brief review of the traditional trust-region globalization scheme and its applications in IBO. Therefore, we will introduce the so-called class of fully linear surrogate models that generalize the characteristics of the surrogate model used by the trust-region methodology from gradient-based models (such as second-order Taylor models) to more general model types (such as radial basis interpolation models).


### 2.1. Generic Trust-Region Method

The trust-region framework is a widely-used technique to attain the global convergence of gradient-based optimization algorithms. At every iterationk, a surrogate model, mk, is considered that approximates the actual objective within a neighborhood of iterate xk. Traditionally, this neighborhood, i.e., the trust-region, is modeled as a sphere.(2)ℬxk,Δk=x∈ℝn:x−xk≤Δk,where xk is the current iterate and Δk is referred to as the trust-region radius.A trial step,sk, is obtained by minimizing the surrogate model, mk, in the trust-region, ℬxk,Δk. This minimization problem is referred to as the trust-region subproblem. The resulting trial step, sk, gives rise to the trivial candidate iterate, xk+1=xk+sk.(3)sk=argmins∈ℬ0,Δkmkxk+s.In every iteration, the size and position of the trust-region are updated in such a way that the generated sequence of iterates converges. The configuration of the trust-region, which is determined by the couple,xk,Δk, is updated according to a fundamental scheme based on the ratio of actual to predicted decrease, ρk. The ratio constitutes a metric for the predictive quality of the surrogate with respect to the current trust-radius [35].(4)ρk=fxk−fxk+skmkxk−mkxk+sk.Based on the value ofρk, the trust region is updated according to the following general rules. If the prediction of the surrogate is good, the radius is increased. In case of a bad prediction, the radius is reduced. For other cases, the current radius is simply maintained.To regulate the expansion and contraction of the trust-region, one takes values0<Δm≤ΔM,0≤η1<η2<1 and 0<γ1<1<γ2, where Δm and ΔM represent the minimal and maximal radii, respectively. The value of Δk is then updated as follows:(5)Δk+1=maxγ1Δk,Δm,ρk≤η1Δk,η1<ρk≤η2minγ2Δk,ΔM,η2<ρk.To decide whether to accept or reject a trial step,sk, we adopt the rules from [35]. If the candidate iterate, xk+sk, is not accepted, then the trust-radius shrinks, given that ρk<0.(6)xk+1=xk+sk,fxk+sk<fxkxk,otherwise.The generic trust-region algorithm is presented in Algorithm 1.Algorithm 1: Generic trust-region algorithm.
(1)
choosex0∈X⊂ℝn, 0<Δm≤Δ0≤ΔM, 0≤η1<η2<1 and 0<γ1<1<γ2(2)
fork=0,1,… until convergence do(3)
obtain modelmk(4)
solve sub-problem (3) to obtain sk(5)
evaluatefxk+sk and compute ρk through (4)(6)
updateΔk+1 and xk+1 according to (5) and (6)(7)
end for(8)
returnxk


### 2.2. Trust-Regions and Fully Linear Models

In the gradient-based trust-region method, the surrogate,mk, is provided by a quadratic model of the following form:(7)mkx=fxk+gk⊤x−xk,+12x−xk⊤Hkx−xk,where gk and Hk are the evaluations of the first- and second-order derivatives of the objective function f in the iterate xk, respectively. According to Taylor’s theorem, a region exists within the vicinity of xk, where the assumption that model mk is a reliable approximation of f is true. A region of reliable approximation can be defined as a region where the surrogate modeling approximation error, reflected by ρk, is small enough such that the descend in the surrogate model, mk, restricted to this region corresponds to the descend of the objective function f. It is apparent that the strength of the trust-region method thus depends on the capacity of the updating rules to maintain a region, where the surrogate modeling approximation error is sufficiently bounded so that mk is suited for optimization purposes [44, 45].From that perspective, the trust-region method can be decomposed into two crucial mechanisms that interact to generate the sequence of iterates. A procedure must be provided that can generate a quadratic approximation modelmk of the goal function f for any given iterate xk. The updating rules will then maintain a series of regions, where the corresponding and consecutive models, mk, are suited to generate the descending sequence of iterates as the solution of subproblem (3). The latter mechanism is in that aspect generic and can be employed to maintain a series of trust-regions given any type of model. Therefore the mechanism extends from quadratic models to the interpolation based and thus derivative-free models that will be used here.In this work, we are interested in the class of so-called fully linear models. The definition of a fully linear model establishes Taylor-like bounds on arbitrary models, and therefore, it generalizes the notion of a reliable approximation. As a consequence, the trust-region framework can be generalized to include such fully linear models.We adopt the definition from [39, 40] with X^=x∈X:ℬx,ΔM⊆X.Definition 1.
For givenκf>0,κg>0, and for givenxk∈X^,Δk∈0,ΔM, definingℬk=ℬxk,Δk, surrogate modelmkis said to be fully linear onℬk, if for allx∈ℬk, (8)fx−mkx≤κfΔk2,∇fx−∇mkx≤κgΔk.
Let us assume two different procedures that can generate a fully linear surrogate modelmk on ℬxk,Δk or that are able to verify whether the surrogate model, mk, is fully linear on ℬxk,Δk and that for any given xk∈X^, Δk∈0,ΔM. If two such procedures are available, (1) can be modified to construct a trust-region algorithm with fully linear models that will converge to a first-order critical point of the function f [39]. Algorithm 2 presents such a trust-region algorithm using fully linear models.
The main modification is the introduction of the termination criteria,∇mkxkεg/2. This condition is a reformulation of the first-order optimality condition on the objective function, ∇fxkϵg. Assuming that mk is fully linear on ℬxk,Δkg, we have the following:(9)∇fxk≤∇mkxk+∇fxk−∇mkxk≤ϵg2+κgΔkg.
Algorithm 2 assures that upon termination, modelmk will be fully linear on ℬxk,ϵg/2κg. Substituting ϵg/2κg for Δkg in (9) reveals that this is equivalent to satisfying the first-order optimality criterion, ∇fxk≤ϵg. We remark that it is thus fundamental to provide procedures that can guarantee the tightest possible bound on κf and κg. In Section 4, such procedures are developed, considering that model mk is an interpolation model.Algorithm 2: Trust-region algorithm with fully linear property.
(1)
choosex0∈X⊂ℝn, 0<Δm≤Δ0≤ΔM, 0≤η1<η2<1 and 0<γ1<1<γ2(2)
fork=0,1,…do(3)
construct fully linear modelmk valid in ℬxk,Δk(4)
if∇mkxk≤ϵg/2then(5)
ifmk is fully linear in ℬxk,ϵg/2κgthen(6)
stop(7)
else(8)
setxk+1=xk and Δk+1=ϵg/2κg(9)
end if(10)
else(11)
solve sub-problem (3) to obtain sk(12)
evaluatefxk+sk and compute ρk through (4)(13)
updateΔk+1 and xk+1 according to (5) and (6)(14)
end if(15)
end for(16)
returnxk


### 2.3. Interpolation Models

When first-order derivative information,gk, is not directly available, it is possible to construct a surrogate model based on strategically generated data, i.e., the direct evaluations of the objective function. A model mk denotes an interpolation model if it satisfies the interpolation conditions in (10) for a given input set, Z=z1,…,zq, with interpolation points, zi∈ℝn, and a corresponding output set, ℱ=fz1,…,fzq.(10)mkzi=fzi,i=1,…,q.A quadratic model, as in (7), is capable of modeling the curvature of the underlying function f and allows, therefore, to proof global convergence to second-order critical points [39]. Correspondingly, several DFO methods have been developed that operate with quadratic interpolation models [34, 46, 47]. To determine a quadratic model, the number of interpolation points, q, must either equal or exceed q=1/2n+1n+2, thus requiring at least that number of objective function evaluations. This requirement contradicts with our intention to limit the number of objective function evaluations.A linear model can be constructed with as few asn+1 interpolation points. Algorithms that operate with minimal interpolation sets of q=n+1 are described by [31, 32, 48]. However, the drawback, in this case, is that a linear model is not capable of representing any curvature present in the objective.In this work, we impose the additional requirement that the minimal number of interpolation points must equalq=n+2 because it is then guaranteed that the resulting model can represent a curvature, given a proper interpolation framework (see Section 3). We will refer to such models as better-than-linear models, noting that these models classify as fully linear models but are constructed with at least one more interpolation point than the minimum for exact linear models q=n+1.


### 2.4. Interpolation-Based Optimization

Anticipating that a procedure is available that can generate an interpolation model,mk, that classifies as a fully linear model, an interpolation framework is trivially embedded in Algorithm 2.In Section3, we will discuss a multivariate interpolation method that can model curvature and is flexible in the sense that it does not impose strict requirements on the number of interpolation points. Within the context of this framework, it will be shown that an interpolation model can be classified as a fully linear model if the set, Z, satisfies conditions on the geometry of the points with respect to the trust-region [40]. The geometry can be understood as the scattering of interpolation points inside a specified region. In Section 4, the geometry of an interpolation set is quantified using a novel geometry metric.The procedures of verifying and generating fully linear interpolation models thus boil down to verifying whether the geometric conditions are satisfied with respect to having a fully-linear model for a given interpolation set and to complementing an existing or possibly empty set so that the condition is satisfied by construction. In the context of IBO, these two mechanisms are provided by a single procedure that is referred to as the interpolation set management or the sampling strategy. A sampling strategy for minimal interpolation sets ofq=n+2 is proposed in Section 4.


## 2.1. Generic Trust-Region Method

The trust-region framework is a widely-used technique to attain the global convergence of gradient-based optimization algorithms. At every iterationk, a surrogate model, mk, is considered that approximates the actual objective within a neighborhood of iterate xk. Traditionally, this neighborhood, i.e., the trust-region, is modeled as a sphere.(2)ℬxk,Δk=x∈ℝn:x−xk≤Δk,where xk is the current iterate and Δk is referred to as the trust-region radius.A trial step,sk, is obtained by minimizing the surrogate model, mk, in the trust-region, ℬxk,Δk. This minimization problem is referred to as the trust-region subproblem. The resulting trial step, sk, gives rise to the trivial candidate iterate, xk+1=xk+sk.(3)sk=argmins∈ℬ0,Δkmkxk+s.In every iteration, the size and position of the trust-region are updated in such a way that the generated sequence of iterates converges. The configuration of the trust-region, which is determined by the couple,xk,Δk, is updated according to a fundamental scheme based on the ratio of actual to predicted decrease, ρk. The ratio constitutes a metric for the predictive quality of the surrogate with respect to the current trust-radius [35].(4)ρk=fxk−fxk+skmkxk−mkxk+sk.Based on the value ofρk, the trust region is updated according to the following general rules. If the prediction of the surrogate is good, the radius is increased. In case of a bad prediction, the radius is reduced. For other cases, the current radius is simply maintained.To regulate the expansion and contraction of the trust-region, one takes values0<Δm≤ΔM,0≤η1<η2<1 and 0<γ1<1<γ2, where Δm and ΔM represent the minimal and maximal radii, respectively. The value of Δk is then updated as follows:(5)Δk+1=maxγ1Δk,Δm,ρk≤η1Δk,η1<ρk≤η2minγ2Δk,ΔM,η2<ρk.To decide whether to accept or reject a trial step,sk, we adopt the rules from [35]. If the candidate iterate, xk+sk, is not accepted, then the trust-radius shrinks, given that ρk<0.(6)xk+1=xk+sk,fxk+sk<fxkxk,otherwise.The generic trust-region algorithm is presented in Algorithm 1.Algorithm 1: Generic trust-region algorithm.
(1)
choosex0∈X⊂ℝn, 0<Δm≤Δ0≤ΔM, 0≤η1<η2<1 and 0<γ1<1<γ2(2)
fork=0,1,… until convergence do(3)
obtain modelmk(4)
solve sub-problem (3) to obtain sk(5)
evaluatefxk+sk and compute ρk through (4)(6)
updateΔk+1 and xk+1 according to (5) and (6)(7)
end for(8)
returnxk


## 2.2. Trust-Regions and Fully Linear Models

In the gradient-based trust-region method, the surrogate,mk, is provided by a quadratic model of the following form:(7)mkx=fxk+gk⊤x−xk,+12x−xk⊤Hkx−xk,where gk and Hk are the evaluations of the first- and second-order derivatives of the objective function f in the iterate xk, respectively. According to Taylor’s theorem, a region exists within the vicinity of xk, where the assumption that model mk is a reliable approximation of f is true. A region of reliable approximation can be defined as a region where the surrogate modeling approximation error, reflected by ρk, is small enough such that the descend in the surrogate model, mk, restricted to this region corresponds to the descend of the objective function f. It is apparent that the strength of the trust-region method thus depends on the capacity of the updating rules to maintain a region, where the surrogate modeling approximation error is sufficiently bounded so that mk is suited for optimization purposes [44, 45].From that perspective, the trust-region method can be decomposed into two crucial mechanisms that interact to generate the sequence of iterates. A procedure must be provided that can generate a quadratic approximation modelmk of the goal function f for any given iterate xk. The updating rules will then maintain a series of regions, where the corresponding and consecutive models, mk, are suited to generate the descending sequence of iterates as the solution of subproblem (3). The latter mechanism is in that aspect generic and can be employed to maintain a series of trust-regions given any type of model. Therefore the mechanism extends from quadratic models to the interpolation based and thus derivative-free models that will be used here.In this work, we are interested in the class of so-called fully linear models. The definition of a fully linear model establishes Taylor-like bounds on arbitrary models, and therefore, it generalizes the notion of a reliable approximation. As a consequence, the trust-region framework can be generalized to include such fully linear models.We adopt the definition from [39, 40] with X^=x∈X:ℬx,ΔM⊆X.Definition 1.
For givenκf>0,κg>0, and for givenxk∈X^,Δk∈0,ΔM, definingℬk=ℬxk,Δk, surrogate modelmkis said to be fully linear onℬk, if for allx∈ℬk, (8)fx−mkx≤κfΔk2,∇fx−∇mkx≤κgΔk.
Let us assume two different procedures that can generate a fully linear surrogate modelmk on ℬxk,Δk or that are able to verify whether the surrogate model, mk, is fully linear on ℬxk,Δk and that for any given xk∈X^, Δk∈0,ΔM. If two such procedures are available, (1) can be modified to construct a trust-region algorithm with fully linear models that will converge to a first-order critical point of the function f [39]. Algorithm 2 presents such a trust-region algorithm using fully linear models.
The main modification is the introduction of the termination criteria,∇mkxkεg/2. This condition is a reformulation of the first-order optimality condition on the objective function, ∇fxkϵg. Assuming that mk is fully linear on ℬxk,Δkg, we have the following:(9)∇fxk≤∇mkxk+∇fxk−∇mkxk≤ϵg2+κgΔkg.
Algorithm 2 assures that upon termination, modelmk will be fully linear on ℬxk,ϵg/2κg. Substituting ϵg/2κg for Δkg in (9) reveals that this is equivalent to satisfying the first-order optimality criterion, ∇fxk≤ϵg. We remark that it is thus fundamental to provide procedures that can guarantee the tightest possible bound on κf and κg. In Section 4, such procedures are developed, considering that model mk is an interpolation model.Algorithm 2: Trust-region algorithm with fully linear property.
(1)
choosex0∈X⊂ℝn, 0<Δm≤Δ0≤ΔM, 0≤η1<η2<1 and 0<γ1<1<γ2(2)
fork=0,1,…do(3)
construct fully linear modelmk valid in ℬxk,Δk(4)
if∇mkxk≤ϵg/2then(5)
ifmk is fully linear in ℬxk,ϵg/2κgthen(6)
stop(7)
else(8)
setxk+1=xk and Δk+1=ϵg/2κg(9)
end if(10)
else(11)
solve sub-problem (3) to obtain sk(12)
evaluatefxk+sk and compute ρk through (4)(13)
updateΔk+1 and xk+1 according to (5) and (6)(14)
end if(15)
end for(16)
returnxk


## 2.3. Interpolation Models

When first-order derivative information,gk, is not directly available, it is possible to construct a surrogate model based on strategically generated data, i.e., the direct evaluations of the objective function. A model mk denotes an interpolation model if it satisfies the interpolation conditions in (10) for a given input set, Z=z1,…,zq, with interpolation points, zi∈ℝn, and a corresponding output set, ℱ=fz1,…,fzq.(10)mkzi=fzi,i=1,…,q.A quadratic model, as in (7), is capable of modeling the curvature of the underlying function f and allows, therefore, to proof global convergence to second-order critical points [39]. Correspondingly, several DFO methods have been developed that operate with quadratic interpolation models [34, 46, 47]. To determine a quadratic model, the number of interpolation points, q, must either equal or exceed q=1/2n+1n+2, thus requiring at least that number of objective function evaluations. This requirement contradicts with our intention to limit the number of objective function evaluations.A linear model can be constructed with as few asn+1 interpolation points. Algorithms that operate with minimal interpolation sets of q=n+1 are described by [31, 32, 48]. However, the drawback, in this case, is that a linear model is not capable of representing any curvature present in the objective.In this work, we impose the additional requirement that the minimal number of interpolation points must equalq=n+2 because it is then guaranteed that the resulting model can represent a curvature, given a proper interpolation framework (see Section 3). We will refer to such models as better-than-linear models, noting that these models classify as fully linear models but are constructed with at least one more interpolation point than the minimum for exact linear models q=n+1.


## 2.4. Interpolation-Based Optimization

Anticipating that a procedure is available that can generate an interpolation model,mk, that classifies as a fully linear model, an interpolation framework is trivially embedded in Algorithm 2.In Section3, we will discuss a multivariate interpolation method that can model curvature and is flexible in the sense that it does not impose strict requirements on the number of interpolation points. Within the context of this framework, it will be shown that an interpolation model can be classified as a fully linear model if the set, Z, satisfies conditions on the geometry of the points with respect to the trust-region [40]. The geometry can be understood as the scattering of interpolation points inside a specified region. In Section 4, the geometry of an interpolation set is quantified using a novel geometry metric.The procedures of verifying and generating fully linear interpolation models thus boil down to verifying whether the geometric conditions are satisfied with respect to having a fully-linear model for a given interpolation set and to complementing an existing or possibly empty set so that the condition is satisfied by construction. In the context of IBO, these two mechanisms are provided by a single procedure that is referred to as the interpolation set management or the sampling strategy. A sampling strategy for minimal interpolation sets ofq=n+2 is proposed in Section 4.


## 3. Radial Basis Functions

Algorithm 2 allows to engage the trust-region globalization scheme whilst substituting interpolation models (10) for the surrogate model, mk. It requires a flexible, multivariate interpolation modeling framework that can capture curvature and does not impose strict conditions on the number of interpolation points since we will recycle interpolation points from previous iterations. The radial basis function (RBF) interpolation framework has readily shown its utility within IBO [31–33]. The traditional RBF framework decides arbitrarily on the shape of the modeling RBFs by fixing a shape parameter a priori. Alternatively, one can determine the shape parameter based on the interpolation data. In that case, the RBF framework coincides with universal kriging [36], where the shape parameter is referred to as a hyperparameter. The modeling assumptions of the framework and the determination of the model parameters and hyperparameters are discussed next.


### 3.1. Interpolation Framework

Consider a given set of interpolation pointsZ and corresponding objective function evaluations ℱ. The RBF multivariate interpolation framework then presumes a function model of the following form:(11)mx=∑zi∈Zαiϕx−zi+px,where ϕ:ℝ+⟶ℝ is an element of the class of RBFs and px is an element of an n-variate polynomial space Pdn of at most degree d, referred to as the polynomial tail. Let P^dn=π1x,…,πd^x be an n-variate basis for Pdn with d^being the size of this n-variate basis. Then, the polynomial tail can be represented by a linear combination of this basis, i.e., px=∑πj∈P^dnβjπjx, with βj being the expansion coefficient. Some popular choices for RBF are listed in Table 1. Notice that the function ϕx has isosurfaces coinciding with hyperspheres in ℝn, i.e., the functions radiate out isotropically.Table 1
Radial basis functions,γ>0, where d0 is the minimal polynomial degree to have CPD.
Typeϕrd0Cubicr31Thin plate spliner2logr1Linearmax0,1−γr0Multiquadraticr2+γ20Gaussianexp−γrβ−1By introducing vectors,ϕx and α, ∈ℝq and, πx and β, ∈ℝd^, which entries are given by ϕx−zi,αi,πjxandβj, respectively-, we can reformulate the RBF model as(12)mx=∑zi∈Zαiϕx−zi+∑πj∈P^dnβjπjx=ϕx⊤⋅α+πx⊤⋅β.


### 3.2. Determination of Model Parameters

The RBF function model (12) is uniquely defined by the parameters α and β. These parameters can be determined as follows:Let the Vandermonde matricesΦ∈ℝq×q and Π∈ℝq×d^ be defined as Φik=ϕzi−zk and Πij=πjzi. The set of q interpolation conditions in (10) can then be written as follows:(13)ΦΠ⋅αβ=f,where the entries of f are defined as the function evaluations contained in ℱ.These conditions do not render uniqueα and β. Therefore, d^ additional conditions are required to complement (13) so that the assembled system matrix is nonsingular. To this end, we require that the model mx is the smoothest model of the form described by (12) that complies with the interpolation condition (10). For RBFs, function smoothness can be defined in relation with the concept of conditional positive definiteness (CPD) [37]. RBFs have the specific property that to any RBF, ϕ, one can associate a value d0 such that if the maximum degree d of the polynomial tail satisfies d≥d0, the matrix Φ will be positive definite with respect to all vectors α∈nullΠ. Equivalently, we have that if d≥d0 and Π⊤α=0, it holds that −1d0+1α⊤Φα>0. In conclusion, when CPD is satisfied, i.e., α⊤Φα>0, it can be shown that the corresponding interpolator is the smoothest function of form (12) satisfying (10). For an elaborate proof and additional details on CPD, we refer to [37].Formally, the matrix representation of the interpolation condition can be complemented as such by the function smoothness condition to obtain the linear system of equations.(14)ΦΠΠ⊤0⋅αβ=f0.The model parameters are determined uniquely if the system matrix is nonsingular. It is easily verified that ifd≥d0, it is then sufficient that Π⊤ has a full column rank. If so, the interpolation set is said to be well-poised for interpolation. In Section 4, how the well-posedness corresponds with the geometrical properties of the points in Z is discussed. [4]Based on the previous observation, the explicit solution of the model parameters is as follows:(15)β=Π⊤Φ−1Π−1Π⊤Φ−1f,α=Φ−1f−Πβ,This result may also be obtained using the universal kriging (UK) framework. Here, a linear predictor of the form mx=ηx⊤f is considered. The objective function is modeled as the realization of a polynomial tail and a random process, fx=px+wx. The process’ stochasticity is quantified by a spatial correlation function that models the output correlation as a function of the spatial distance between the input values. The correlation is thus expressed in terms of RBF as corrwziwzj=ϕzi−zj. It can be shown, e.g., [36, 38], that the optimal linear unbiased predictor, η, is then determined by the following system of equations, where ξ can be interpreted as a Lagrangian multiplier.(16)ΦΠΠ⊤0⋅ηξ=ϕπ.For additional details, we refer to appendix A.


### 3.3. Hyperparameters

The value of the shape parameterγ (see Table 1) is arbitrary and affects the receding radial effect of RBF contribution in (12). Heuristics may be provided, however, they usually lack proper theoretical foundation. A suitable value is mostly determined by trial-and-error. Furthermore, we emphasize that each component of the interpolation points zi contributes equally to the distance metric, ⋅, i.e., the model presumes the isotropic dependence of the objective function on all of its variables. Such isotropic behavior can be achieved by normalizing the optimization space by performing a linear transformation of the input space. Nonetheless, such measure cannot account for any local deviations on the global trend. Since we are to construct local interpolation models whose validity is not to stretch beyond the limits of the trust-region boundary, we could redefine this linear transformation locally.To account for local anisotropy in the principal directions of the input space, one may define the distance metric as (17). The RBFs in (12) are then redefined as ϕx−ziγ. It is clear that in this case, the univariate hyperparameter γ has become redundant. In practical implementations of the UK framework, the determination of hyperparameters γi is automated using a maximum log likelihood estimate based on the observations, ℱ.(17)xγ2=∑i=1nγixi2.We refer to appendix A and [49] for elaborate details.


## 3.1. Interpolation Framework

Consider a given set of interpolation pointsZ and corresponding objective function evaluations ℱ. The RBF multivariate interpolation framework then presumes a function model of the following form:(11)mx=∑zi∈Zαiϕx−zi+px,where ϕ:ℝ+⟶ℝ is an element of the class of RBFs and px is an element of an n-variate polynomial space Pdn of at most degree d, referred to as the polynomial tail. Let P^dn=π1x,…,πd^x be an n-variate basis for Pdn with d^being the size of this n-variate basis. Then, the polynomial tail can be represented by a linear combination of this basis, i.e., px=∑πj∈P^dnβjπjx, with βj being the expansion coefficient. Some popular choices for RBF are listed in Table 1. Notice that the function ϕx has isosurfaces coinciding with hyperspheres in ℝn, i.e., the functions radiate out isotropically.Table 1
Radial basis functions,γ>0, where d0 is the minimal polynomial degree to have CPD.
Typeϕrd0Cubicr31Thin plate spliner2logr1Linearmax0,1−γr0Multiquadraticr2+γ20Gaussianexp−γrβ−1By introducing vectors,ϕx and α, ∈ℝq and, πx and β, ∈ℝd^, which entries are given by ϕx−zi,αi,πjxandβj, respectively-, we can reformulate the RBF model as(12)mx=∑zi∈Zαiϕx−zi+∑πj∈P^dnβjπjx=ϕx⊤⋅α+πx⊤⋅β.


## 3.2. Determination of Model Parameters

The RBF function model (12) is uniquely defined by the parameters α and β. These parameters can be determined as follows:Let the Vandermonde matricesΦ∈ℝq×q and Π∈ℝq×d^ be defined as Φik=ϕzi−zk and Πij=πjzi. The set of q interpolation conditions in (10) can then be written as follows:(13)ΦΠ⋅αβ=f,where the entries of f are defined as the function evaluations contained in ℱ.These conditions do not render uniqueα and β. Therefore, d^ additional conditions are required to complement (13) so that the assembled system matrix is nonsingular. To this end, we require that the model mx is the smoothest model of the form described by (12) that complies with the interpolation condition (10). For RBFs, function smoothness can be defined in relation with the concept of conditional positive definiteness (CPD) [37]. RBFs have the specific property that to any RBF, ϕ, one can associate a value d0 such that if the maximum degree d of the polynomial tail satisfies d≥d0, the matrix Φ will be positive definite with respect to all vectors α∈nullΠ. Equivalently, we have that if d≥d0 and Π⊤α=0, it holds that −1d0+1α⊤Φα>0. In conclusion, when CPD is satisfied, i.e., α⊤Φα>0, it can be shown that the corresponding interpolator is the smoothest function of form (12) satisfying (10). For an elaborate proof and additional details on CPD, we refer to [37].Formally, the matrix representation of the interpolation condition can be complemented as such by the function smoothness condition to obtain the linear system of equations.(14)ΦΠΠ⊤0⋅αβ=f0.The model parameters are determined uniquely if the system matrix is nonsingular. It is easily verified that ifd≥d0, it is then sufficient that Π⊤ has a full column rank. If so, the interpolation set is said to be well-poised for interpolation. In Section 4, how the well-posedness corresponds with the geometrical properties of the points in Z is discussed. [4]Based on the previous observation, the explicit solution of the model parameters is as follows:(15)β=Π⊤Φ−1Π−1Π⊤Φ−1f,α=Φ−1f−Πβ,This result may also be obtained using the universal kriging (UK) framework. Here, a linear predictor of the form mx=ηx⊤f is considered. The objective function is modeled as the realization of a polynomial tail and a random process, fx=px+wx. The process’ stochasticity is quantified by a spatial correlation function that models the output correlation as a function of the spatial distance between the input values. The correlation is thus expressed in terms of RBF as corrwziwzj=ϕzi−zj. It can be shown, e.g., [36, 38], that the optimal linear unbiased predictor, η, is then determined by the following system of equations, where ξ can be interpreted as a Lagrangian multiplier.(16)ΦΠΠ⊤0⋅ηξ=ϕπ.For additional details, we refer to appendix A.


## 3.3. Hyperparameters

The value of the shape parameterγ (see Table 1) is arbitrary and affects the receding radial effect of RBF contribution in (12). Heuristics may be provided, however, they usually lack proper theoretical foundation. A suitable value is mostly determined by trial-and-error. Furthermore, we emphasize that each component of the interpolation points zi contributes equally to the distance metric, ⋅, i.e., the model presumes the isotropic dependence of the objective function on all of its variables. Such isotropic behavior can be achieved by normalizing the optimization space by performing a linear transformation of the input space. Nonetheless, such measure cannot account for any local deviations on the global trend. Since we are to construct local interpolation models whose validity is not to stretch beyond the limits of the trust-region boundary, we could redefine this linear transformation locally.To account for local anisotropy in the principal directions of the input space, one may define the distance metric as (17). The RBFs in (12) are then redefined as ϕx−ziγ. It is clear that in this case, the univariate hyperparameter γ has become redundant. In practical implementations of the UK framework, the determination of hyperparameters γi is automated using a maximum log likelihood estimate based on the observations, ℱ.(17)xγ2=∑i=1nγixi2.We refer to appendix A and [49] for elaborate details.


## 4. Interpolation set Management

In Section2, we described a trust-region method with fully linear models. It was demonstrated how a sequence of fully linear models mk can be utilized to converge to a first-order critical point of the objective function f. When considering interpolation models satisfying (10), the question whether mk is a fully linear model boils down to a proper management of the interpolation set Z used to construct mk. This so-called sampling strategy should maintain a proper geometry of the interpolation set.In [40], the geometry of an arbitrary interpolation set, Z, was quantified by a single set-dependent constant directly related to the conditions in (8) to define a fully linear model. The lemma in [40], however, only applies to interpolation models for which the set size equals q=n+1 exactly. On this basis, the authors in [31] devised an IBO algorithm, employing RBF interpolation models, given at least n+1 interpolation points to construct a fully linear model. By applying a QR pivoting strategy directly related to the measure described by [40], their sampling strategy assures that a subset of n+1 interpolation points (amongst the entire available data set) exhibits the required geometry to guarantee a fully linear interpolation model. Possibly, additional samples need to be generated to have the verification of the full linearity of the interpolation model. We, henceforth, refer to these points to verify that mk is fully linear as the affine subset. Additional points from the entire data set are subsequently recycled to enhance the modeling capabilities of the RBF model as long as the well-poisedness of (14) is not corrupted.Having an interpolation modelmk representing the curvature requires an affine subset of q≥n+2 (see Section 3.2) that is adequately embedded in an interpolation-based trust-region algorithm. In other words, we aim to ensure the convergence of Algorithm 2 to second-order critical points whilst only adding a single point to the affine subset with respect to [31, 32]. Since an estimate of the curvature is now always available, the remark on the RBF hyperparameters in Section 3.3 can be generalized to arbitrary directions. The latter is investigated in Section 5.In Section4.1, we introduce a generalized version of the theorem in [40] for interpolation sets with size exceeding n+1. On that basis, we identify the optimal affine subset geometry that corresponds with n+2 interpolation points in 4.2. Finally, in 4.3, we propose a procedure capable to maintain an interpolation set with an affine subset of n+2 points so that the interpolation model, mk, is fully linear and exhibits curvature.


### 4.1. Better-Than-Linear Models

The uniqueness of the RBF interpolation model parameters implied well-poisedness, i.e., full column rank, of the matrixΠ⊤. Thus, the size of the interpolation set must equal or exceed the size of the polynomial basis. Given that we address a minimal interpolation set size n+2 and thus d^≤n+2, the polynomial tail will be linear or equivalently px∈P1n. A convenient choice for the basis P^1n is then clearly πx⊤=x⊤1. Since we will embed the interpolation framework in a trust-region framework, the current iterate xk is always available for interpolation since fxk is evaluated for calculating the predictive quality of the interpolation model (4). Henceforth, we will consider a trust-region centered at the origin, 0, by performing a shift of coordinates each time. Furthermore, we require that 0 is always an element of Z, i.e., Z=0,z1,…,zq−1.Introducing the interpolation set Vandermonde matrix,Z, matrix Π⊤ can be written as follows:(18)Π⊤=0z1⋯zq−111⋯1=0Z11⊤.As mentioned in [40], the conditions (8) were related to the interpolation set geometry, and more specifically, they were related to matrix Z. Here, a generalized lemma is presented, establishing similar Taylor-like conditions for better-than-linear models, i.e., when q≥n+2. A derivation is included in appendix B.Lemma 1.
Suppose thatmis an RBF model with linear tail and RBF contributionα⊤ϕthat satisfies the interpolation condition,mzi=fzi, for the interpolation setZ:=0,z1,…,zq−1⊂ℬ0,Δ, whereq≥n+2and satisfiesZ01Δ≤ΛZ. Furthermore, suppose thatfandmare continuously differentiable in the hypersphereℬ0,Δand that∇fand∇ϕ⊤αare Lipschitz continuous inℬ0,Δ, with Lipschitz constantsγfandγϕ, respectively, then the following two inequalities are satisfied for anyx∈ℬ0,Δ:(19)mx−fx≤52ΛZq−1+12γf+γϕΔ2∇mx−∇fx≤52ΛZq−1γf+γϕΔ.These conditions are similar to (8). The lemma provides, therefore, an admittedly conservative estimate of the constantsκf and κg, however, foremost it reveals the crucial role of the interpolation set geometry, which is reflected by the value of Z.It follows that our interpolation set management must ensure the boundedness ofΔk01Zk from below by the a priori fixed value ΛZ01 at each iteration.


### 4.2. Optimal Critical Set Geometry

The value ofΔ01Z is understood as a measure for the algebraic well-poisedness of a given set Z⊂ℬ0,Δ within ℬ0,Δ [40]. For any given number of interpolation points q, we can define the optimal critical set geometry that corresponds with the largest value of Z corresponding to the tightest bounds on κf and κg. Hence, assuring the set’s well-poisedness reduces to generating an affine subset that resembles the optimal set geometry as closely as possible whilst recycling as many points as possible. This strategy balances the requirement to minimize the total amount of function evaluations without corrupting the model quality. Figure 1 depicts optimal geometries for n=2 dimensions, maximizing the measure Z for affine subsets of size q=n+1, q=n+2, and q=2n+1, respectively. The optimal set geometry for q=n+1 complies with the normalized set Δ01Z∖0 being orthonormal with respect to the origin, or equivalently when Z coincides with a standard simplex. This set geometry is pursued by QR pivoting strategies as elaborated in [31] and Lagrange or Newton polynomials-based strategies as elaborated in [32, 40]. The optimal set geometry for q=n+2 corresponds with the regular simplex geometry (plus the origin) [19, 50] and that for q=2n+1 with that of an orthoplex (plus the origin) [51]. Orthoplex geometry has been recommended by Powell in [34]. This geometry also determines the fixed search directions in many pattern search optimization methods.Figure 1
Optimally poised interpolation sets for varyingq and n=2 next to a contour of ϑx;Z. Note that the center point, 0, is each time incorporated in the interpolation sets. (a) q=n+1. (b) q=n+2. (c) q=2n+1.
(a)(b)(c)In this context, it is relevant to inspect the spatial distribution of the 2D geometries in Figure1. The optimal set geometries for q=n+2 and q=2n+1 are unoriented: the center of mass of the corresponding sets coincide with the origin. It does not apply for q=n+1, where the set is biased. Note that the above is valid for arbitrary n. To quantify this more rigorously, we introduce the metric, ϑ∗Z, defined as follows:(20)ϑ∗Z=maxx∈ℬ0,Δϑx;Z=1q∑zi∈Zx−zi.Its value is equal to the maximum mean distance that any test point in the trust-region can be removed from the modelling setZ (including 0). The value gives an indication for the worst-case contribution that the radial basis term (that is directly related to the distance between the test point and the interpolation points) has compared to the linear tail. In Table 2, we compare ϑ∗Z values for different optimal set geometries. It can be observed that ϑ∗Z values can be improved by only adding a single point, referring to the shift from the standard to the regular simplex geometry. The contours of ϑx;Z are presented in Figure 1 that allow one to visualize the regions that are not represented well by the optimal set geometries. For n⟶∞, ϑ∗Z tends to 2 for all three given configurations, implying that all are as good, or as worse, for infinite dimensional problems.Table 2
ϑZ for different optimal set geometries contained within ℬ0,1.
qn+1n+22n+2n≥21+2n2+nn/n+13+2n2−n/n+21+2n2−nn+2n2+nn/2n+121.56521.251.245331.58211.29281.155241.58561.31651.214251.58441.33211.2513∞222The discussion based onϑ∗Z provides an additional argument to choose the regular simplex over the configuration in Figure 1(a), certainly for low-dimensional problems.


### 4.3. Regular Simplex Affine Subsets

According to the previous subsection, the optimal affine subset (without the origin) forq=n+2 is a regular simplex. 0 remains an element of the affine subset that affects the definition of a regular simplex affine subset based on the n+1 remaining interpolation points.Formally, the set of all regular simplices lying on the surface∂ℬ0,1 is given by (21). ln⋅ (see appendix C) is a scaling factor depending on the radius of the n-sphere (here 1).(21)ℛn0,1=s0,…,sn|si∈∂ℬ0,1,andsi−sj=ln1,i≠j.The geometry of the regular simplex can be related to the matrix norm,Z, by examining the hypervolume spanned by the convex hull of the affine subset, Z∖0. Consider, therefore, the hypervolume, VZ, of an arbitrary set, Z=z1,…,zn+1, defined in (22) with spanning vectors, z2−z1,…,zn+1−z1. As outlined in [50, 52], the regular simplex has volumetric extremity properties and maximizes the set well-poisedness metric, Z.(22)VZ=detz2−z1⋯zn+1−z1n!.Hence, the algebraic problem of maximizingZ coincides with the geometric problem of maximizing the volume of the convex hull of Z. Approaching the problem from this point of view helps to devise an intuitive interpolation set management, resulting in a critical subset that is (close to) an element of the set ℛnxk,Δk, whilst recycling the interpolation points was evaluated in previous iterations.Based on the hypervolumeVZ of the subset Z, a metric can be construed whether Z matches a regular simplex sufficiently. The principal idea is that VZ should be close to the theoretical maximal hypervolume of the regular simplex, V∗. In particular, we require the hypervolume that is not smaller than the maximal hypervolume multiplied with some predefined relative error factor 0≤μ1≤1. If VZ≥μ1V∗ is not satisfied, new interpolation points need to be generated. This procedure guarantees Δ01Z≥ΛZ01.Before we elaborate the particularities of our set management, we mention that each set of pointsℛ∈ℛmn0,1, is a lower order regular simplex itself, where m is the dimension of the regular subsimplex considered, and therefore, m≤n. ℛmn0,1 is defined in (23). This statement holds since any subset of a set ℛ∈ℛn0,1 simply inherits the property that the distance between every two points is the same.(23)ℛmn0,1=s0,…,sm,Vert∃sm+1,…,sn:s0,…,sn∈ℛn0,1.


#### 4.3.1. Algorithm Expand2volume

Based on the previous idea, Algorithm 3 constructs an affine subset in a systematic manner so thatmk classifies as a fully linear model for given ΛZ. The goal is to attain an affine subset with hypervolume VZ≥μ1V∗ whilst adding as few points as possible to the available set of interpolation points. The input is the set of all interpolations points available, scaled with respect to the trust radius, Δk, and unbiased with respect to the center, xk, i.e., Zk=zi=Δk01si−xk|si∈Sk∖xk. Sk is the set incorporating every point evaluated so far. The current origin is excluded such that we remain with only the candidate interpolation points to form the affine subset.The algorithm embarks on removing all points that are outside a certain periphery, parametrized by the variableθ1≥1. Note that we artificially enlarge the trust-region since we allow parameter θ1 to be greater than 1. We can always assure that the interpolation model is fully linear by redefining κf and κg. It ensures that in practice, the number of interpolation points actually accepted to be part of the affine subset is slightly increased, improving the overall convergence. Otherwise, the algorithm might end up expanding an empty set too regularly (certainly for larger n). If this operation yields a subset, Z1, that is still larger than the required affine subset size of n+1 (which, in practice, will almost never occur), all points that are furthest away from the trust-surface are additionally removed.If the size of the resulting subset,Z1, is smaller than n+1, it is expanded by calling Simexpand Algorithm 4. The procedure described by Simexpand exploits the idea that if the set Z1 were to be an element of ℛn−mn0,1, we can expand it with an element of ℛm−1n0,1 so that their union would then constitute an element from ℛnn0,1. The expansion element can be generated by rotating and translating a properly sized lower-order regular simplex.Our numerical implementation proceeds as follows: a set,X is determined to be an element of ℛm−10,1 embedded in the n-dimensional workspace. The direction e is determined as the orthogonal projection direction of the origin on the subspace spanned by the vertices of Z1. Subsequently, a rotation matrix Q is determined so that the basis spanned by the rotated element QX is orthogonal to the union of the basis spanned by Z1 and the direction e. In a final step, the rotated element QX is scaled according to rmn1 and translated along the direction e over a distance dmn1. The factors rmn. and dmn. are explained in appendix C and will only maximize the volume in the occasional coincidence, where Z1∈ℛn−mn0,1. Nonetheless, this same reasoning can be applied when Z1∉ℛn−mn0,1. In that case, an additional subproblem must be solved to identify the optimal value of σ to maximize the volume of their union. Here, σ is a scaling factor related to rmn. and dmn.. In the description of the algorithm, we make no distinction between a set X and its matrix column vector representation X.We also use the following operators: operator. constructs a matrix containing a set of column vectors that span the space determined by the vertices in Z=z0,…,zq, i.e., Z=z1−z0…zq−z0. The function regsimn,m outputs a set of vertices X=x0,…,xm that constitutes a regular m-simplex in an n-dimensional coordinate space. The matrix function QR. constructs orthonormal bases for span N and kernel K, determined by the columns of the input matrix, respectively. The mechanics of Simexpand are illustrated in Figure 2. Each subfigure depicts a set Z1 that differs in size from the previous set, the corresponding expansion set Z2, the two bases N, and NX and the vector e.Figure 2
Illustration of the expansion mechanism ofSimexpand for varying m and n=3. (a) m = 1, (b) m = 2, and (c) m = 3.
(a)(b)(c)If the unionZ1∪Z2 does not satisfy the criticality condition VZ1∪Z2≥μ1V∗, an iterative procedure is initiated, where points are removed from the set Z1 systematically. Then, the contracted set can be expanded by applying Simexpand again. To avoid that, expansion points are generated too close to a point that was removed in an earlier iteration. Information from the removed point is transferred to the contracted set Z1. Information about a point that gets removed is stored by transferring it to a memory set Z′. The mean of this memory set is used to complement the points in Z1 that remain before expanding the set with a new regular m+1-simplex. In this fashion, the added simplex is generated away from the mass center of the original set Z1.The algorithm ends after at mostn+1 iterations, i.e., when the entire original set Z1 is transferred to the memory set Z′. We remark that in this case, the mass center is expanded with an n−1-regular simplex that is, however, not necessarily an element of ℛn−1n0,1. As a consequence, the eventual affine subset needs to not be an element of ℛn0,1. As this does not ensure that the sufficiency condition is satisfied, the algorithm closes with a final verification of the hypervolume and possible expansion of the set Z2.Algorithm 3: Expand2volume.
Input:Z0=Z/0,θ1,μ1Output:Z2(1)
obtainZ1 with elements zi∈Z0 such that zi≤θ1(2)
defineZ2=∅(3)
if|Z1|>n+1 then(4)
orderzi such that 1−z1≤…≤1−zq(5)
obtainZ1=Z1/zn+2,…,zq(6)
else ifZ1<n+1 then(7)
getZ2=SimexpandZ1(8)
end if(9)
find initial memoryz′=zi∈Z1:zi≤zj,∀j(10)
define memory setZ′=z′(11)
whileVZ1∪Z2<μ1V∗ or |Z1|>0 do(12)
findzi∈Z1:zi−z′≤zj−z′,∀j(13)
contractZ1=Z1/zi(14)
redefine memory setZ′=Z′∪zi(15)
expandZ2=SimexpandZ1∪meanZ′(16)
end while(17)
ifVZ1∪Z2<μ1V∗ then(18)
Z2=Z2∪SimexpandZ2(19)
end if(20)
returnZ2Algorithm 4: Simexpand.
Input:Z1Output:Z2(1)
m=n−Z1(2)
X=regsimn,m−1(3)
N,∼=QRZ1, N0,K0=QRZ1∪0, NX,KX=QRX∪0(4)
∼,e=QRNK0(5)
Q=K0N0⋅NXKX⊤(6)
σ∗=argminσVσ⋅1m⊤⋅e+1−σ2⋅Q⋅X(7)
Z2=z1…zn−Z=σ∗⋅1m⊤⋅e+1−σ∗2⋅Q⋅X(8)
returnZ2


### 4.4. Remarks

With respect to algorithms 3 and 4 described above, we add the following: during the initialization ofExpand2-volume, all points that are furthest away from the trust-surface are removed. It is a heuristic. Alternatively, one might run through a more complex scheme, where every possible n+1 size set is verified with respect to their hypervolume, and the largest set is selected. These points could also be added to the memory set. However, it is a situation that will almost never occur in practice, partially because the optimization algorithm will guide the iterations to unexplored regions of the parameter space. It is also because the trust region will contract when the sequence of iterates starts to converge. Hence, by proper choice of μ1 and θ1, it can be avoided that the complete interpolation set of the previous iteration is within the current periphery.Furthermore, we remark that the alignment of the kernel associated with the space spanned by the unionZ1∪0 and the space spanned by the regular simplex X is underdetermined. Here, this additional degree of freedom remains unexploited. Alternatively, this degree of freedom could be employed by orienting the simplex that is addedso that the generated points are located furthest away from the available set of interpolation points or such that the added points correspond to a mean minimal predicted objective function value.


## 4.1. Better-Than-Linear Models

The uniqueness of the RBF interpolation model parameters implied well-poisedness, i.e., full column rank, of the matrixΠ⊤. Thus, the size of the interpolation set must equal or exceed the size of the polynomial basis. Given that we address a minimal interpolation set size n+2 and thus d^≤n+2, the polynomial tail will be linear or equivalently px∈P1n. A convenient choice for the basis P^1n is then clearly πx⊤=x⊤1. Since we will embed the interpolation framework in a trust-region framework, the current iterate xk is always available for interpolation since fxk is evaluated for calculating the predictive quality of the interpolation model (4). Henceforth, we will consider a trust-region centered at the origin, 0, by performing a shift of coordinates each time. Furthermore, we require that 0 is always an element of Z, i.e., Z=0,z1,…,zq−1.Introducing the interpolation set Vandermonde matrix,Z, matrix Π⊤ can be written as follows:(18)Π⊤=0z1⋯zq−111⋯1=0Z11⊤.As mentioned in [40], the conditions (8) were related to the interpolation set geometry, and more specifically, they were related to matrix Z. Here, a generalized lemma is presented, establishing similar Taylor-like conditions for better-than-linear models, i.e., when q≥n+2. A derivation is included in appendix B.Lemma 1.
Suppose thatmis an RBF model with linear tail and RBF contributionα⊤ϕthat satisfies the interpolation condition,mzi=fzi, for the interpolation setZ:=0,z1,…,zq−1⊂ℬ0,Δ, whereq≥n+2and satisfiesZ01Δ≤ΛZ. Furthermore, suppose thatfandmare continuously differentiable in the hypersphereℬ0,Δand that∇fand∇ϕ⊤αare Lipschitz continuous inℬ0,Δ, with Lipschitz constantsγfandγϕ, respectively, then the following two inequalities are satisfied for anyx∈ℬ0,Δ:(19)mx−fx≤52ΛZq−1+12γf+γϕΔ2∇mx−∇fx≤52ΛZq−1γf+γϕΔ.These conditions are similar to (8). The lemma provides, therefore, an admittedly conservative estimate of the constantsκf and κg, however, foremost it reveals the crucial role of the interpolation set geometry, which is reflected by the value of Z.It follows that our interpolation set management must ensure the boundedness ofΔk01Zk from below by the a priori fixed value ΛZ01 at each iteration.


## 4.2. Optimal Critical Set Geometry

The value ofΔ01Z is understood as a measure for the algebraic well-poisedness of a given set Z⊂ℬ0,Δ within ℬ0,Δ [40]. For any given number of interpolation points q, we can define the optimal critical set geometry that corresponds with the largest value of Z corresponding to the tightest bounds on κf and κg. Hence, assuring the set’s well-poisedness reduces to generating an affine subset that resembles the optimal set geometry as closely as possible whilst recycling as many points as possible. This strategy balances the requirement to minimize the total amount of function evaluations without corrupting the model quality. Figure 1 depicts optimal geometries for n=2 dimensions, maximizing the measure Z for affine subsets of size q=n+1, q=n+2, and q=2n+1, respectively. The optimal set geometry for q=n+1 complies with the normalized set Δ01Z∖0 being orthonormal with respect to the origin, or equivalently when Z coincides with a standard simplex. This set geometry is pursued by QR pivoting strategies as elaborated in [31] and Lagrange or Newton polynomials-based strategies as elaborated in [32, 40]. The optimal set geometry for q=n+2 corresponds with the regular simplex geometry (plus the origin) [19, 50] and that for q=2n+1 with that of an orthoplex (plus the origin) [51]. Orthoplex geometry has been recommended by Powell in [34]. This geometry also determines the fixed search directions in many pattern search optimization methods.Figure 1
Optimally poised interpolation sets for varyingq and n=2 next to a contour of ϑx;Z. Note that the center point, 0, is each time incorporated in the interpolation sets. (a) q=n+1. (b) q=n+2. (c) q=2n+1.
(a)(b)(c)In this context, it is relevant to inspect the spatial distribution of the 2D geometries in Figure1. The optimal set geometries for q=n+2 and q=2n+1 are unoriented: the center of mass of the corresponding sets coincide with the origin. It does not apply for q=n+1, where the set is biased. Note that the above is valid for arbitrary n. To quantify this more rigorously, we introduce the metric, ϑ∗Z, defined as follows:(20)ϑ∗Z=maxx∈ℬ0,Δϑx;Z=1q∑zi∈Zx−zi.Its value is equal to the maximum mean distance that any test point in the trust-region can be removed from the modelling setZ (including 0). The value gives an indication for the worst-case contribution that the radial basis term (that is directly related to the distance between the test point and the interpolation points) has compared to the linear tail. In Table 2, we compare ϑ∗Z values for different optimal set geometries. It can be observed that ϑ∗Z values can be improved by only adding a single point, referring to the shift from the standard to the regular simplex geometry. The contours of ϑx;Z are presented in Figure 1 that allow one to visualize the regions that are not represented well by the optimal set geometries. For n⟶∞, ϑ∗Z tends to 2 for all three given configurations, implying that all are as good, or as worse, for infinite dimensional problems.Table 2
ϑZ for different optimal set geometries contained within ℬ0,1.
qn+1n+22n+2n≥21+2n2+nn/n+13+2n2−n/n+21+2n2−nn+2n2+nn/2n+121.56521.251.245331.58211.29281.155241.58561.31651.214251.58441.33211.2513∞222The discussion based onϑ∗Z provides an additional argument to choose the regular simplex over the configuration in Figure 1(a), certainly for low-dimensional problems.


## 4.3. Regular Simplex Affine Subsets

According to the previous subsection, the optimal affine subset (without the origin) forq=n+2 is a regular simplex. 0 remains an element of the affine subset that affects the definition of a regular simplex affine subset based on the n+1 remaining interpolation points.Formally, the set of all regular simplices lying on the surface∂ℬ0,1 is given by (21). ln⋅ (see appendix C) is a scaling factor depending on the radius of the n-sphere (here 1).(21)ℛn0,1=s0,…,sn|si∈∂ℬ0,1,andsi−sj=ln1,i≠j.The geometry of the regular simplex can be related to the matrix norm,Z, by examining the hypervolume spanned by the convex hull of the affine subset, Z∖0. Consider, therefore, the hypervolume, VZ, of an arbitrary set, Z=z1,…,zn+1, defined in (22) with spanning vectors, z2−z1,…,zn+1−z1. As outlined in [50, 52], the regular simplex has volumetric extremity properties and maximizes the set well-poisedness metric, Z.(22)VZ=detz2−z1⋯zn+1−z1n!.Hence, the algebraic problem of maximizingZ coincides with the geometric problem of maximizing the volume of the convex hull of Z. Approaching the problem from this point of view helps to devise an intuitive interpolation set management, resulting in a critical subset that is (close to) an element of the set ℛnxk,Δk, whilst recycling the interpolation points was evaluated in previous iterations.Based on the hypervolumeVZ of the subset Z, a metric can be construed whether Z matches a regular simplex sufficiently. The principal idea is that VZ should be close to the theoretical maximal hypervolume of the regular simplex, V∗. In particular, we require the hypervolume that is not smaller than the maximal hypervolume multiplied with some predefined relative error factor 0≤μ1≤1. If VZ≥μ1V∗ is not satisfied, new interpolation points need to be generated. This procedure guarantees Δ01Z≥ΛZ01.Before we elaborate the particularities of our set management, we mention that each set of pointsℛ∈ℛmn0,1, is a lower order regular simplex itself, where m is the dimension of the regular subsimplex considered, and therefore, m≤n. ℛmn0,1 is defined in (23). This statement holds since any subset of a set ℛ∈ℛn0,1 simply inherits the property that the distance between every two points is the same.(23)ℛmn0,1=s0,…,sm,Vert∃sm+1,…,sn:s0,…,sn∈ℛn0,1.


### 4.3.1. Algorithm Expand2volume

Based on the previous idea, Algorithm 3 constructs an affine subset in a systematic manner so thatmk classifies as a fully linear model for given ΛZ. The goal is to attain an affine subset with hypervolume VZ≥μ1V∗ whilst adding as few points as possible to the available set of interpolation points. The input is the set of all interpolations points available, scaled with respect to the trust radius, Δk, and unbiased with respect to the center, xk, i.e., Zk=zi=Δk01si−xk|si∈Sk∖xk. Sk is the set incorporating every point evaluated so far. The current origin is excluded such that we remain with only the candidate interpolation points to form the affine subset.The algorithm embarks on removing all points that are outside a certain periphery, parametrized by the variableθ1≥1. Note that we artificially enlarge the trust-region since we allow parameter θ1 to be greater than 1. We can always assure that the interpolation model is fully linear by redefining κf and κg. It ensures that in practice, the number of interpolation points actually accepted to be part of the affine subset is slightly increased, improving the overall convergence. Otherwise, the algorithm might end up expanding an empty set too regularly (certainly for larger n). If this operation yields a subset, Z1, that is still larger than the required affine subset size of n+1 (which, in practice, will almost never occur), all points that are furthest away from the trust-surface are additionally removed.If the size of the resulting subset,Z1, is smaller than n+1, it is expanded by calling Simexpand Algorithm 4. The procedure described by Simexpand exploits the idea that if the set Z1 were to be an element of ℛn−mn0,1, we can expand it with an element of ℛm−1n0,1 so that their union would then constitute an element from ℛnn0,1. The expansion element can be generated by rotating and translating a properly sized lower-order regular simplex.Our numerical implementation proceeds as follows: a set,X is determined to be an element of ℛm−10,1 embedded in the n-dimensional workspace. The direction e is determined as the orthogonal projection direction of the origin on the subspace spanned by the vertices of Z1. Subsequently, a rotation matrix Q is determined so that the basis spanned by the rotated element QX is orthogonal to the union of the basis spanned by Z1 and the direction e. In a final step, the rotated element QX is scaled according to rmn1 and translated along the direction e over a distance dmn1. The factors rmn. and dmn. are explained in appendix C and will only maximize the volume in the occasional coincidence, where Z1∈ℛn−mn0,1. Nonetheless, this same reasoning can be applied when Z1∉ℛn−mn0,1. In that case, an additional subproblem must be solved to identify the optimal value of σ to maximize the volume of their union. Here, σ is a scaling factor related to rmn. and dmn.. In the description of the algorithm, we make no distinction between a set X and its matrix column vector representation X.We also use the following operators: operator. constructs a matrix containing a set of column vectors that span the space determined by the vertices in Z=z0,…,zq, i.e., Z=z1−z0…zq−z0. The function regsimn,m outputs a set of vertices X=x0,…,xm that constitutes a regular m-simplex in an n-dimensional coordinate space. The matrix function QR. constructs orthonormal bases for span N and kernel K, determined by the columns of the input matrix, respectively. The mechanics of Simexpand are illustrated in Figure 2. Each subfigure depicts a set Z1 that differs in size from the previous set, the corresponding expansion set Z2, the two bases N, and NX and the vector e.Figure 2
Illustration of the expansion mechanism ofSimexpand for varying m and n=3. (a) m = 1, (b) m = 2, and (c) m = 3.
(a)(b)(c)If the unionZ1∪Z2 does not satisfy the criticality condition VZ1∪Z2≥μ1V∗, an iterative procedure is initiated, where points are removed from the set Z1 systematically. Then, the contracted set can be expanded by applying Simexpand again. To avoid that, expansion points are generated too close to a point that was removed in an earlier iteration. Information from the removed point is transferred to the contracted set Z1. Information about a point that gets removed is stored by transferring it to a memory set Z′. The mean of this memory set is used to complement the points in Z1 that remain before expanding the set with a new regular m+1-simplex. In this fashion, the added simplex is generated away from the mass center of the original set Z1.The algorithm ends after at mostn+1 iterations, i.e., when the entire original set Z1 is transferred to the memory set Z′. We remark that in this case, the mass center is expanded with an n−1-regular simplex that is, however, not necessarily an element of ℛn−1n0,1. As a consequence, the eventual affine subset needs to not be an element of ℛn0,1. As this does not ensure that the sufficiency condition is satisfied, the algorithm closes with a final verification of the hypervolume and possible expansion of the set Z2.Algorithm 3: Expand2volume.
Input:Z0=Z/0,θ1,μ1Output:Z2(1)
obtainZ1 with elements zi∈Z0 such that zi≤θ1(2)
defineZ2=∅(3)
if|Z1|>n+1 then(4)
orderzi such that 1−z1≤…≤1−zq(5)
obtainZ1=Z1/zn+2,…,zq(6)
else ifZ1<n+1 then(7)
getZ2=SimexpandZ1(8)
end if(9)
find initial memoryz′=zi∈Z1:zi≤zj,∀j(10)
define memory setZ′=z′(11)
whileVZ1∪Z2<μ1V∗ or |Z1|>0 do(12)
findzi∈Z1:zi−z′≤zj−z′,∀j(13)
contractZ1=Z1/zi(14)
redefine memory setZ′=Z′∪zi(15)
expandZ2=SimexpandZ1∪meanZ′(16)
end while(17)
ifVZ1∪Z2<μ1V∗ then(18)
Z2=Z2∪SimexpandZ2(19)
end if(20)
returnZ2Algorithm 4: Simexpand.
Input:Z1Output:Z2(1)
m=n−Z1(2)
X=regsimn,m−1(3)
N,∼=QRZ1, N0,K0=QRZ1∪0, NX,KX=QRX∪0(4)
∼,e=QRNK0(5)
Q=K0N0⋅NXKX⊤(6)
σ∗=argminσVσ⋅1m⊤⋅e+1−σ2⋅Q⋅X(7)
Z2=z1…zn−Z=σ∗⋅1m⊤⋅e+1−σ∗2⋅Q⋅X(8)
returnZ2


## 4.3.1. Algorithm Expand2volume

Based on the previous idea, Algorithm 3 constructs an affine subset in a systematic manner so thatmk classifies as a fully linear model for given ΛZ. The goal is to attain an affine subset with hypervolume VZ≥μ1V∗ whilst adding as few points as possible to the available set of interpolation points. The input is the set of all interpolations points available, scaled with respect to the trust radius, Δk, and unbiased with respect to the center, xk, i.e., Zk=zi=Δk01si−xk|si∈Sk∖xk. Sk is the set incorporating every point evaluated so far. The current origin is excluded such that we remain with only the candidate interpolation points to form the affine subset.The algorithm embarks on removing all points that are outside a certain periphery, parametrized by the variableθ1≥1. Note that we artificially enlarge the trust-region since we allow parameter θ1 to be greater than 1. We can always assure that the interpolation model is fully linear by redefining κf and κg. It ensures that in practice, the number of interpolation points actually accepted to be part of the affine subset is slightly increased, improving the overall convergence. Otherwise, the algorithm might end up expanding an empty set too regularly (certainly for larger n). If this operation yields a subset, Z1, that is still larger than the required affine subset size of n+1 (which, in practice, will almost never occur), all points that are furthest away from the trust-surface are additionally removed.If the size of the resulting subset,Z1, is smaller than n+1, it is expanded by calling Simexpand Algorithm 4. The procedure described by Simexpand exploits the idea that if the set Z1 were to be an element of ℛn−mn0,1, we can expand it with an element of ℛm−1n0,1 so that their union would then constitute an element from ℛnn0,1. The expansion element can be generated by rotating and translating a properly sized lower-order regular simplex.Our numerical implementation proceeds as follows: a set,X is determined to be an element of ℛm−10,1 embedded in the n-dimensional workspace. The direction e is determined as the orthogonal projection direction of the origin on the subspace spanned by the vertices of Z1. Subsequently, a rotation matrix Q is determined so that the basis spanned by the rotated element QX is orthogonal to the union of the basis spanned by Z1 and the direction e. In a final step, the rotated element QX is scaled according to rmn1 and translated along the direction e over a distance dmn1. The factors rmn. and dmn. are explained in appendix C and will only maximize the volume in the occasional coincidence, where Z1∈ℛn−mn0,1. Nonetheless, this same reasoning can be applied when Z1∉ℛn−mn0,1. In that case, an additional subproblem must be solved to identify the optimal value of σ to maximize the volume of their union. Here, σ is a scaling factor related to rmn. and dmn.. In the description of the algorithm, we make no distinction between a set X and its matrix column vector representation X.We also use the following operators: operator. constructs a matrix containing a set of column vectors that span the space determined by the vertices in Z=z0,…,zq, i.e., Z=z1−z0…zq−z0. The function regsimn,m outputs a set of vertices X=x0,…,xm that constitutes a regular m-simplex in an n-dimensional coordinate space. The matrix function QR. constructs orthonormal bases for span N and kernel K, determined by the columns of the input matrix, respectively. The mechanics of Simexpand are illustrated in Figure 2. Each subfigure depicts a set Z1 that differs in size from the previous set, the corresponding expansion set Z2, the two bases N, and NX and the vector e.Figure 2
Illustration of the expansion mechanism ofSimexpand for varying m and n=3. (a) m = 1, (b) m = 2, and (c) m = 3.
(a)(b)(c)If the unionZ1∪Z2 does not satisfy the criticality condition VZ1∪Z2≥μ1V∗, an iterative procedure is initiated, where points are removed from the set Z1 systematically. Then, the contracted set can be expanded by applying Simexpand again. To avoid that, expansion points are generated too close to a point that was removed in an earlier iteration. Information from the removed point is transferred to the contracted set Z1. Information about a point that gets removed is stored by transferring it to a memory set Z′. The mean of this memory set is used to complement the points in Z1 that remain before expanding the set with a new regular m+1-simplex. In this fashion, the added simplex is generated away from the mass center of the original set Z1.The algorithm ends after at mostn+1 iterations, i.e., when the entire original set Z1 is transferred to the memory set Z′. We remark that in this case, the mass center is expanded with an n−1-regular simplex that is, however, not necessarily an element of ℛn−1n0,1. As a consequence, the eventual affine subset needs to not be an element of ℛn0,1. As this does not ensure that the sufficiency condition is satisfied, the algorithm closes with a final verification of the hypervolume and possible expansion of the set Z2.Algorithm 3: Expand2volume.
Input:Z0=Z/0,θ1,μ1Output:Z2(1)
obtainZ1 with elements zi∈Z0 such that zi≤θ1(2)
defineZ2=∅(3)
if|Z1|>n+1 then(4)
orderzi such that 1−z1≤…≤1−zq(5)
obtainZ1=Z1/zn+2,…,zq(6)
else ifZ1<n+1 then(7)
getZ2=SimexpandZ1(8)
end if(9)
find initial memoryz′=zi∈Z1:zi≤zj,∀j(10)
define memory setZ′=z′(11)
whileVZ1∪Z2<μ1V∗ or |Z1|>0 do(12)
findzi∈Z1:zi−z′≤zj−z′,∀j(13)
contractZ1=Z1/zi(14)
redefine memory setZ′=Z′∪zi(15)
expandZ2=SimexpandZ1∪meanZ′(16)
end while(17)
ifVZ1∪Z2<μ1V∗ then(18)
Z2=Z2∪SimexpandZ2(19)
end if(20)
returnZ2Algorithm 4: Simexpand.
Input:Z1Output:Z2(1)
m=n−Z1(2)
X=regsimn,m−1(3)
N,∼=QRZ1, N0,K0=QRZ1∪0, NX,KX=QRX∪0(4)
∼,e=QRNK0(5)
Q=K0N0⋅NXKX⊤(6)
σ∗=argminσVσ⋅1m⊤⋅e+1−σ2⋅Q⋅X(7)
Z2=z1…zn−Z=σ∗⋅1m⊤⋅e+1−σ∗2⋅Q⋅X(8)
returnZ2


## 4.4. Remarks

With respect to algorithms 3 and 4 described above, we add the following: during the initialization ofExpand2-volume, all points that are furthest away from the trust-surface are removed. It is a heuristic. Alternatively, one might run through a more complex scheme, where every possible n+1 size set is verified with respect to their hypervolume, and the largest set is selected. These points could also be added to the memory set. However, it is a situation that will almost never occur in practice, partially because the optimization algorithm will guide the iterations to unexplored regions of the parameter space. It is also because the trust region will contract when the sequence of iterates starts to converge. Hence, by proper choice of μ1 and θ1, it can be avoided that the complete interpolation set of the previous iteration is within the current periphery.Furthermore, we remark that the alignment of the kernel associated with the space spanned by the unionZ1∪0 and the space spanned by the regular simplex X is underdetermined. Here, this additional degree of freedom remains unexploited. Alternatively, this degree of freedom could be employed by orienting the simplex that is addedso that the generated points are located furthest away from the available set of interpolation points or such that the added points correspond to a mean minimal predicted objective function value.


## 5. Ellipsoidal Trust-Regions

The procedures described in the previous two sections can be employed to verify and construct a model,mk, that satisfies Definition 1 for a fully linear model in a given trust-region, ℬxk,Δk. Hence, these procedures can be used to embed an interpolation framework in the trust-region algorithm with fully linear models described by Algorithm 2, factually realizing the IBO algorithm.As for now, the geometry of the interpolation points is considered only with regard to the input space and with respect to the trust-region defined in it. Given a spherical trust-region, our input space-oriented sampling strategy will treat every direction equally in the function of a well-conditioned system matrixZ. However, from the perspective of optimization, it might be of interest to account for the geometry of the objective function as well.The gradient of the objective function will be more sensitive to certain search directions than to others as determined by the local curvature of the objective function. Currently, the set management, however, treats all directions equally (as a result of the spherical shape of the trust-region) and will try to homogenize the distance from the trust-region center and in between the interpolation points. Now, as a result of the local curvature, the corresponding function values may vary either faster or slower than that can be expected when solely considering the local objective gradient. Hence, the geometry of the objective function with respect to the arbitrary choice of the optimization space coordinates is reflected by its local curvature.Therefore, we propose to execute the set management in a transformed space by introducing an auxiliary optimization parameter,u. The goal of this transformation is to obtain a normalized curvature of the objective function in the transformed optimization space. We propose to use an adaptive affine mapping from the transformed space to the original space that changes with the trust-region iterations. Notice that the origin in the auxiliary space corresponds with the trust-region center in the original space.(24)xu=Akz+bk=Aku+xk.The trust-region will be stretched along those directions, where the local gradient is varying slowly and is shrunk in those directions, where the local gradient is varying rapidly. As a consequence, the objective function values evaluated along the trust-region boundary will increase or decrease in equal amounts with respect to the objective function value at the trust-region center. It will have a positive effect on the hyperparameter tuning described in 3.3 since the principal directions are now homogenized.Another consequence is that a spherical trust-region in the auxiliary spaceu will now have an ellipsoidal shape in the original optimization space illustrated by figures 3(a) and 3(b). As a result, a step of length 1 taken in the auxiliary space represented in the original space may exceed the original spherical trust-region radius. Two situations can be distinguished when inspecting the directional curvature. When having a large directional curvature on the one hand, the ellipsoidal geometry will have a stabilizing effect, considering that the set of admissible steps (i.e., the trust-region) becomes less aggressive in terms of large jumps in the objective function, which will also have an indirect effect on the modelling procedure. On the other hand, when having a small directional curvature, the ellipsoidal trust region will allow the optimization procedure to take large steps in the direction of a slowly varying gradient. In this manner, the original optimization space will be explored more rapidly. For instance, in figure 3(b), the trust-region in the auxiliary space encloses a larger part of the inner contour that contains the actual minimum. Similar ideas to exploit the curvature information have been used in the context of pattern search algorithms [?], however, to the extent of our knowledge, not to reshape the trust-region geometry in the context of IBO.Figure 3
Illustration of a spherical trust-region and its isovolumetric ellipsoidal counterpart for the two-dimensional Rosenbrock function. The curvature was normalized using the exact Hessian. (a) Spherical trust-region. (b) Ellipsoidal trust-region.
(a)(b)


### 5.1. Curvature Normalization

The goal of the coordinate transformation is to normalize the local curvature of the problem in the new coordinate space. This condition can be expressed mathematically as given. For reasons that will become apparent, we introduce an additional scaling factorλk2.(25)∇u2fxk=Ak⊤∇x2fxkAk=λk2In.Assume that we now have an estimateBk for the Hessian ∇x2fxk. We can formulate an expression for Ak by considering the singular value decomposition (SVD) of matrix Bk. The SVD decomposition of the positive definite matrix Bk is given by UkΛkUk⊤. Here, Uk is an orthonormal matrix and Λk is a diagonal matrix with the singular values σ1,…,σn on its main diagonal. It is then easily verified that,(26)Ak=λkUkΛk01/2.Finally, we will utilize the scaling factorλk to keep subproblem (3) as it manifests in the auxiliary coordinate space proportional to its manifestation in the original coordinate space. Mathematically, this proportionality can be realized by determining a measure that remains invariant under the transformation. Recall that as a result of the affine coordinate transformation, a spherical trust-region in the auxiliary optimization space will adopt an ellipsoidal geometry in the original optimization space. We propose the following measure of proportionality: the value of λk is determined such that the spherical and ellipsoidal trust-region representations remain isovolumetric.Consider that the volume of ann-dimensional hypersphere with radius Δk is equal to ηn⋅Δkn and that the volume of an n-dimensional ellipsoid with the given half principal axis lengths, a1,…,an, is equal to ηn⋅∏iai, where ηn is a mutual scaling factor. We also have that by the construction of Ak, half the principal axis lengths can be written in the function of the singular values as ai=λkσi01/2. An isovolumetric relation between the trust-region representations in both the original as well as the auxiliary coordinate space can be maintained by proper choice of the scaling factor λk, i.e.,(27)λk=Δk∏i=1nσi1/2n.The exact coordinate transformation of (24) can then be written explicitly as follows:(28)xu=Δk∏i=1nσi1/2nUkΛk01/2u+xk.


### 5.2. Exponential Smoothing

Several methods can be devised to provide an estimate of the Hessian,Bk. Since our sample strategy ensures that the local interpolation model is curved, we can calculate an approximation of the Hessian at the next iterate in the auxiliary coordinate space directly from the model, i.e., ∇u2mkuk+1.This estimate can be calculated back to the original coordinate space to provide an estimate for the actual Hessian,B^k+1=Ak∇u2mkuk+1Ak⊤. To get a smooth approximation, we propose to use exponential smoothing. Note that if μ2=1, this framework reduces to the original spherical trust-region framework.(29)Bk+1=μ2Bk+1−μ2B^k+1,B0=In.


## 5.1. Curvature Normalization

The goal of the coordinate transformation is to normalize the local curvature of the problem in the new coordinate space. This condition can be expressed mathematically as given. For reasons that will become apparent, we introduce an additional scaling factorλk2.(25)∇u2fxk=Ak⊤∇x2fxkAk=λk2In.Assume that we now have an estimateBk for the Hessian ∇x2fxk. We can formulate an expression for Ak by considering the singular value decomposition (SVD) of matrix Bk. The SVD decomposition of the positive definite matrix Bk is given by UkΛkUk⊤. Here, Uk is an orthonormal matrix and Λk is a diagonal matrix with the singular values σ1,…,σn on its main diagonal. It is then easily verified that,(26)Ak=λkUkΛk01/2.Finally, we will utilize the scaling factorλk to keep subproblem (3) as it manifests in the auxiliary coordinate space proportional to its manifestation in the original coordinate space. Mathematically, this proportionality can be realized by determining a measure that remains invariant under the transformation. Recall that as a result of the affine coordinate transformation, a spherical trust-region in the auxiliary optimization space will adopt an ellipsoidal geometry in the original optimization space. We propose the following measure of proportionality: the value of λk is determined such that the spherical and ellipsoidal trust-region representations remain isovolumetric.Consider that the volume of ann-dimensional hypersphere with radius Δk is equal to ηn⋅Δkn and that the volume of an n-dimensional ellipsoid with the given half principal axis lengths, a1,…,an, is equal to ηn⋅∏iai, where ηn is a mutual scaling factor. We also have that by the construction of Ak, half the principal axis lengths can be written in the function of the singular values as ai=λkσi01/2. An isovolumetric relation between the trust-region representations in both the original as well as the auxiliary coordinate space can be maintained by proper choice of the scaling factor λk, i.e.,(27)λk=Δk∏i=1nσi1/2n.The exact coordinate transformation of (24) can then be written explicitly as follows:(28)xu=Δk∏i=1nσi1/2nUkΛk01/2u+xk.


## 5.2. Exponential Smoothing

Several methods can be devised to provide an estimate of the Hessian,Bk. Since our sample strategy ensures that the local interpolation model is curved, we can calculate an approximation of the Hessian at the next iterate in the auxiliary coordinate space directly from the model, i.e., ∇u2mkuk+1.This estimate can be calculated back to the original coordinate space to provide an estimate for the actual Hessian,B^k+1=Ak∇u2mkuk+1Ak⊤. To get a smooth approximation, we propose to use exponential smoothing. Note that if μ2=1, this framework reduces to the original spherical trust-region framework.(29)Bk+1=μ2Bk+1−μ2B^k+1,B0=In.


## 6. Numerical Results

In this section, we present the numerical results of the presented derivative-free optimization algorithm applied on a set of unconstrained problems from the CUTEst test environment [53]. As indicated by the neighborhood representation metric, ϑZ, the method should be especially suited for low-dimensional problems and was developed from this point of view. We have applied the algorithm on a set of test functions with a maximum problem dimension n=8. All tested algorithms were initialized from the traditional starting point [54] and are run once provided that theirexecution is otherwise entirely deterministic. The following benchmark procedure has been adopted, and the results are presented in Table 3. The results as obtained with 4 different versions of the present algorithm have been compared with those obtained with its closest competitors in the context of derivative-free optimization, being ORBIT [31] and NMSMAX. Since ORBIT outperforms the NEWUOA algorithm (at least for small scale problems [31]), we did not compare our algorithm with NEWUOA [34]. The details of ORBIT are elaborated in the associated paper and are essentially similar to the details of the presented algorithm when the spherical trust-region framework is used and are apart from the novel regular simplex-based sampling strategy. We used the implementation made available at [55]. The NMSMAX algorithm is a frequently practiced implementation of the Nelder-Mead simplex method [19] and is made available in the matrix computation toolbox [56]. The NMSMAX algorithm was initialized with a regular simplex in accordance with the present algorithm and is included as a benchmark.Table 3
Comparative results between 4 versions of the algorithms ORBIT and NMSMAX.
NamenTotal number of function evaluationsAbsolute difference to known minimumV1V2V3V4ORBITNMSV1V2V3V4ORBITNMSBEALE228253242419635,75e-076,25e-093,10e-071,34e-075,16e-102,93e-07BRANIN22425323653561,10e-075,35e-077,60e-086,14e-071,10e-057,81e-07ROSENBROCK2373733398381002,51e-079,56e-089,40e-071,69e-071,80e-033,71e-07SCHWEFEL3622625263269699,35e-074,29e-071,44e-077,43e-089,17e-074,39e-07GULF349415916522610003893,29e-077,71e-071,84e-098,95e-072,66e+005,60e-07HART33524649693481299,98e-098,33e-077,14e-072,64e-083,91e-066,08e-07HELIX3534248484431158,29e-075,97e-073,08e-074,50e-074,61e-044,31e-07ROSENBROCK37768808810001763,81e-074,41e-076,66e-071,18e-075,20e-035,12e-07SCHWEFEL435546416777931,72e-077,98e-087,20e-071,81e-088,05e-077,07e-07BROWNDEN4118841121213982366,77e-079,34e-065,16e-063,44e-069,60e-069,77e-06HART4450424838621223,86e-072,78e-073,95e-072,94e-075,56e-078,85e-07POWELL4201921041486153666,77e-077,43e-078,85e-074,51e-076,01e-059,45e-07ROSENBROCK413611114314410003328,67e-079,90e-073,10e-077,74e-072,75e-028,78e-07SCHWEFEL44636162741121562,46e-079,04e-075,55e-073,18e-079,25e-079,75e-07WOOD416679851092552176,35e-075,16e-073,78e-072,13e-084,04e-047,57e-07BIGGS666029424352915610007,54e-041,52e-079,65e-079,48e-076,60e-035,68e-05HART66104116951071391795,63e-075,43e-075,23e-074,08e-078,21e-079,39e-07ROSENBROCK630225726443510008183,02e-076,78e-078,09e-077,34e-077,22e-029,90e-07SCHWEFEL461241161111691732218,37e-077,10e-072,38e-077,37e-079,47e-078,81e-07TRID64052062051943995416,30e-036,67e-072,77e-071,79e-076,90e-037,94e-07WATSON611783971192223052,11e-078,63e-071,56e-074,16e-079,83e-077,21e-07POWELL81681351391512254107,61e-072,88e-075,77e-073,95e-071,00e-069,70e-07ROSENBROCK8912648624769100010007,41e-077,25e-079,35e-078,96e-074,81e-012,75e-01TRID81000648585778100010001,10e-049,81e-078,49e-073,68e-073,53e-042,40e-03Each algorithm was given a budget of 1000 function evaluations. Since every test function has a known global minimum,f∗, the tolerance of having a convergence of an algorithm was 10−6, i.e., fxk−f∗<10−6. The internal stopping conditions of each of the algorithms were set to 0 so that an algorithm only stops when this benchmark condition was satisfied or when it exceeded its evaluation budget. These internal stopping conditions include the gradient condition (9) and the minimal simplex size for NMSMAX. Only the total number of function evaluations and the difference between the final objective function value and the known minimum are recorded. We did not consider the CPU time and performance and data profiles that are gaining currency in the optimization community [57, 58]. Since the focus of the paper is primarily on the novelty of the sampling strategy and the ellipsoidal trust-region framework, we reason that the adopted approach was best-suited to illustrate the elaborated methods.Each algorithm was applied using its standard settings to mimic the realistic practical usage. The only parameter that was varied was the initial trust-region radius, or in case of the NMSMAX algorithm, the size of the initial simplex. It corresponds with an arbitrary uniform scaling of the original optimization space, which is a standard practice when performing optimization routines. The standard parameter settings of the presented algorithm are given in Table4, those of each of the competing algorithms can be found in the respective references. The parameters present in Table 4 were chosen via trial-and-error on a limited benchmark set. Those corresponding to the ORBIT algorithm were adopted from that reference. For future development, a short description of the influence of each parameter is included in table as well.Table 4
Standard settings of the algorithm.
θ11.25Determines the set of points that is considered in algorithm 3. Forθ1=1, all points must be inside the trust-region. For slightly larger values, the nearby points are also considered. Large values will positively bias the volume so that even for poor set distributions, the volumetric condition will be satisfied.η10Standard trust-region parameter. Adopted from [35].η20.6Seeη1γ10.5Seeη1γ22Seeη1β110−2Lower-limit for the hyperparameters in the RBF framework. Determines how fast the RBFs decay. For normalized search space and reasonably smooth objective functions, it was determined to be a good limit by trial-and-error.β210−1Upper-limit for the hyperparameters in the RBF framework. Seeβ1.Δm0Minimal trust-region size. Determines the smallest focus the trust-region framework can have.ΔM3Δ0Maximum trust-region size. Determines the largest focus the trust-region framework can have. Value was adopted from [31].qM100Maximum number of points to be included in the RBF model. The value was chosen arbitrarily and worked well for relatively low-dimensional problems. When the number is too large, the influence of far-away points may distort the local geometry of the objective. Furthermore, the fitting and evaluating of the interpolation model becomes computationally challenging as a result of hyperparameter optimization.We have considered 4 different versions of the present algorithm. Details about the numerical implementation are given in the next subsection. These 4 versions allow us to study the separate effects of the two presented mechanisms, namely the novel regular simplex-based sampling strategy and the ellipsoidal trust-region mechanism. The different versions are defined in Table5 and are related to the values of μ1 and μ2 from Algorithm 3 and (29), respectively. To study the effects of the regular simplex-based sampling strategy apart from the ellipsoidal trust-region framework, version V1 can be compared with the ORBIT and NMSMAX algorithm. Note that in this case, the interpolation set management is purely input space-oriented. To study the effect of the ellipsoidal trust-region framework, version V2 can be compared with respect to V1. Recall that in this case, the shape of the trust-region is adapted to the output behavior of the objective function. To study the effect of the novel sampling strategy in combination with the ellipsoidal trust-region, versions V2 to V4 can be compared. Note that when parameter μ2 (29) is equal to 1, it corresponds with a spherical trust-region framework. It is not possible to study the behavior of the ellipsoidal trust-region framework in combination with the sampling strategy of ORBIT, given that this strategy does not assure that a better-than-linear model is available every iteration.Table 5
Tested versions of the algorithm.
m10.50.751m20.95V2V3V41V1


### 6.1. Numerical Implementation

All described procedures were implemented inMatlab R2017 computational environment and run on a Dell Laptop with an Intel i7 CPU with 4 cores. All code is available at the GitHub repository in [59]. Some of the subproblems were solved using the available implementations of routine in Matlab. These are documented next.As mentioned in Section4, once an affine subset was determined using Algorithm 3, the available set of interpolation points was reconsidered to enhance the modelling capabilities of the surrogate model. Additional points were added as long as they did not corrupt the well-poisedness of (14) and if the total amount of interpolation points did not exceed the limit qM. A similar procedure is performed by ORBIT and was adopted from that source code.The interpolation framework used here is that of universal kriging. We employed the DACE toolbox to obtain the modelling parameters as described in Section3. This framework requires the user to define the lower and upper limits of the hyperparameter β, β1, and β2, respectively. More specifically, we used Gaussian RBFs for the following reasons: they generally perform well and exhibit nice interpolation properties that result in a reasonably smooth optimization surface. Other kernels could be considered as well, however, their influence should be limited.Subproblem (3) was solved using the Matlab function fminsearchcon, an extension of the function fminsearch made available on MathWorks File exchange [60]. Since absolute CPU time was not considered, the solution procedure was executed until fully converged.


### 6.2. Discussion of Results

#### 6.2.1. Benchmark Problems

By assessing the results displayed in Table3, one can observe that, on the whole, the proposed algorithm outperforms both NMSMAX and ORBIT, disregarding some coincidental draws.Considering that V1 outperformed ORBIT in all but one occasion, we can conclude that the regular simplex-based sampling strategy offers a significant advantage. One might, however, argue that the increased implementational and computational complexity of Algorithm 3 does not justify the advantageous optimization setting in case that less challenging optimization problems are considered. In none of the benchmark problems did V1 outperform V2, V3, or V4, indicating that engaging an ellipsoidal trust-region framework is clearly a beneficial endeavor regardless, given the low implementational cost. However, the downside is that to obtain a curved model, a minimum ofn+2 points are required, making the advanced regular simplex sampling strategy indispensable.Whether V2, V3, or V4 outperformed one another depends on the objective function at hand. However, on the whole, we could state that V2 is the better candidate for lower-dimensional problems, whilst V3 is the better candidate for higher dimensional problems, as will be confirmed in the next section. Also, considering the convergence history of the different algorithms, it was noted that the ORBIT algorithm occasionally stagnated on a local plateau of the objective function. Since this behavior was also exhibited by V1 (e.g., GULF, TRID; although V1 recovered each time), it demonstrates that compensating for the output behavior of the objective function by the input space-oriented sampling strategy and by altering the shape of the trust-regions is quintessential to recover from such apparent plateaus.


#### 6.2.2. Extended Powell Singular and Extended Rosenbrock

To demonstrate and comment on the mechanics of the elaborated methods in greater detail, convergence results on the 8-dimensional extended Powell singular and extended Rosenbrock are presented in Figure4 and discussed below.Figure 4
Convergence history represented by the evolution of the best relative function value in the function of the total number of function evaluations. (a) Extended Powell Singular,n=8. (b) Extended Rosenbrock, n=8.
(a)(b)Consider the convergence history for the extended Powell singular function in Figure4(a). We will initially discuss the history of the algorithms V1, ORBIT, and NMSMAX. All exhibit approximately the same convergence behavior during the initial stage of the optimization (up to 100 function evaluations). After 100 function evaluations, the NMSMAX algorithm clearly starts to diverge from this trend. It is because of the fact that the NMSMAX algorithm is bounded to operate on its vertices. At the common point of approximately 100 function evaluations, it reaches a region where the simplex size needs to shrink sufficiently to attain a size that is small enough to accurately grasp the local curvature of the objective function. A similar phenomenon is observed for the ORBIT algorithm, which starts to diverge from the global trend approximately after 200 function evaluations. The phenomenon occurs later for the ORBIT algorithm since it is not bounded to operate on the vertices of its standard simplex, and the solution of (3) is only limited to the convex set, ℬxk,Δk. Therefore, it can maintain a larger trust-region throughout the iterations than NMSMAX. That algorithm V1 can maintain a larger trust-region longer than ORBIT, and NMSMAX is a result of the improved geometry of the interpolation points. Apart from that, algorithm V1 too fails to reach the benchmark condition within the admissible budget. Now, consider the behavior of algorithm V2. The sampling strategy is equivalent to that of V1, however, the modelling space is now normalized with respect to the estimated curvature of the objective function. It corresponds with an ellipsoidal trust-region geometry in the original optimization space. The effect is hardly observable during the initial stage but clearly demonstrates its usefulness after approximately 500 function evaluations, where V2 diverges from V1 and successfully reaches the benchmark condition. It is argued that because of the transformation, the modeling and optimization steps of the algorithm are now executed in a space, where the objective function behaves more as if it was a quadratic function. The effect of the parameter μ1 can be assessed by comparing V2, V3, and V4. The parameter determines how strictly the volume of an exact regular simplex has to be respected. It results in a larger number of interpolation points to be added in every iteration. Therefore, algorithm V3 and V4 diverge from the common trend at about the same moment V1 and ORBIT do. In this region of the optimization space, the size of the trust-region (simplex) starts to shrink for all algorithms. It means that less function evaluations can be recycled from one trust-region to another. Algorithms V3 and V4 set a stricter limit on the volume of the simplex and are, therefore, forced to add more points to complement the recycled set. As a result, they consume more function evaluations than V1 and V2 to cross the same region. Nonetheless, this increase in function evaluations benefits the overall modelling procedure, and V2 is surpassed by V3 eventually.These observations are confirmed by the convergence history for the extended Rosenbrock function, see Figure4(b). Here, the results are directly affected by the local behavior of the objective function and all algorithms start to diverge at the starting point. In this case, V1 does reach the benchmark condition, again illustrating the improved scattering of points with respect to state-of-the-art. We argue this is a result of the regular simplex based sampling strategy. The beneficial effect of the ellipsoidal trust-region framework is again clearly present as V2 to V4 converge faster than V1. In conclusion, one can observe that the same effect is present for V3 as was the case for the extended Powell singular function. Function evaluations are consumed at a higher rate in the initial stage. This phenomenon ultimately benefits the modelling procedure in the final optimization stage, resulting in a steeper convergence history and even the surpassing of V2 with respect to convergence speed.


## 6.1. Numerical Implementation

All described procedures were implemented inMatlab R2017 computational environment and run on a Dell Laptop with an Intel i7 CPU with 4 cores. All code is available at the GitHub repository in [59]. Some of the subproblems were solved using the available implementations of routine in Matlab. These are documented next.As mentioned in Section4, once an affine subset was determined using Algorithm 3, the available set of interpolation points was reconsidered to enhance the modelling capabilities of the surrogate model. Additional points were added as long as they did not corrupt the well-poisedness of (14) and if the total amount of interpolation points did not exceed the limit qM. A similar procedure is performed by ORBIT and was adopted from that source code.The interpolation framework used here is that of universal kriging. We employed the DACE toolbox to obtain the modelling parameters as described in Section3. This framework requires the user to define the lower and upper limits of the hyperparameter β, β1, and β2, respectively. More specifically, we used Gaussian RBFs for the following reasons: they generally perform well and exhibit nice interpolation properties that result in a reasonably smooth optimization surface. Other kernels could be considered as well, however, their influence should be limited.Subproblem (3) was solved using the Matlab function fminsearchcon, an extension of the function fminsearch made available on MathWorks File exchange [60]. Since absolute CPU time was not considered, the solution procedure was executed until fully converged.


## 6.2. Discussion of Results

### 6.2.1. Benchmark Problems

By assessing the results displayed in Table3, one can observe that, on the whole, the proposed algorithm outperforms both NMSMAX and ORBIT, disregarding some coincidental draws.Considering that V1 outperformed ORBIT in all but one occasion, we can conclude that the regular simplex-based sampling strategy offers a significant advantage. One might, however, argue that the increased implementational and computational complexity of Algorithm 3 does not justify the advantageous optimization setting in case that less challenging optimization problems are considered. In none of the benchmark problems did V1 outperform V2, V3, or V4, indicating that engaging an ellipsoidal trust-region framework is clearly a beneficial endeavor regardless, given the low implementational cost. However, the downside is that to obtain a curved model, a minimum ofn+2 points are required, making the advanced regular simplex sampling strategy indispensable.Whether V2, V3, or V4 outperformed one another depends on the objective function at hand. However, on the whole, we could state that V2 is the better candidate for lower-dimensional problems, whilst V3 is the better candidate for higher dimensional problems, as will be confirmed in the next section. Also, considering the convergence history of the different algorithms, it was noted that the ORBIT algorithm occasionally stagnated on a local plateau of the objective function. Since this behavior was also exhibited by V1 (e.g., GULF, TRID; although V1 recovered each time), it demonstrates that compensating for the output behavior of the objective function by the input space-oriented sampling strategy and by altering the shape of the trust-regions is quintessential to recover from such apparent plateaus.


### 6.2.2. Extended Powell Singular and Extended Rosenbrock

To demonstrate and comment on the mechanics of the elaborated methods in greater detail, convergence results on the 8-dimensional extended Powell singular and extended Rosenbrock are presented in Figure4 and discussed below.Figure 4
Convergence history represented by the evolution of the best relative function value in the function of the total number of function evaluations. (a) Extended Powell Singular,n=8. (b) Extended Rosenbrock, n=8.
(a)(b)Consider the convergence history for the extended Powell singular function in Figure4(a). We will initially discuss the history of the algorithms V1, ORBIT, and NMSMAX. All exhibit approximately the same convergence behavior during the initial stage of the optimization (up to 100 function evaluations). After 100 function evaluations, the NMSMAX algorithm clearly starts to diverge from this trend. It is because of the fact that the NMSMAX algorithm is bounded to operate on its vertices. At the common point of approximately 100 function evaluations, it reaches a region where the simplex size needs to shrink sufficiently to attain a size that is small enough to accurately grasp the local curvature of the objective function. A similar phenomenon is observed for the ORBIT algorithm, which starts to diverge from the global trend approximately after 200 function evaluations. The phenomenon occurs later for the ORBIT algorithm since it is not bounded to operate on the vertices of its standard simplex, and the solution of (3) is only limited to the convex set, ℬxk,Δk. Therefore, it can maintain a larger trust-region throughout the iterations than NMSMAX. That algorithm V1 can maintain a larger trust-region longer than ORBIT, and NMSMAX is a result of the improved geometry of the interpolation points. Apart from that, algorithm V1 too fails to reach the benchmark condition within the admissible budget. Now, consider the behavior of algorithm V2. The sampling strategy is equivalent to that of V1, however, the modelling space is now normalized with respect to the estimated curvature of the objective function. It corresponds with an ellipsoidal trust-region geometry in the original optimization space. The effect is hardly observable during the initial stage but clearly demonstrates its usefulness after approximately 500 function evaluations, where V2 diverges from V1 and successfully reaches the benchmark condition. It is argued that because of the transformation, the modeling and optimization steps of the algorithm are now executed in a space, where the objective function behaves more as if it was a quadratic function. The effect of the parameter μ1 can be assessed by comparing V2, V3, and V4. The parameter determines how strictly the volume of an exact regular simplex has to be respected. It results in a larger number of interpolation points to be added in every iteration. Therefore, algorithm V3 and V4 diverge from the common trend at about the same moment V1 and ORBIT do. In this region of the optimization space, the size of the trust-region (simplex) starts to shrink for all algorithms. It means that less function evaluations can be recycled from one trust-region to another. Algorithms V3 and V4 set a stricter limit on the volume of the simplex and are, therefore, forced to add more points to complement the recycled set. As a result, they consume more function evaluations than V1 and V2 to cross the same region. Nonetheless, this increase in function evaluations benefits the overall modelling procedure, and V2 is surpassed by V3 eventually.These observations are confirmed by the convergence history for the extended Rosenbrock function, see Figure4(b). Here, the results are directly affected by the local behavior of the objective function and all algorithms start to diverge at the starting point. In this case, V1 does reach the benchmark condition, again illustrating the improved scattering of points with respect to state-of-the-art. We argue this is a result of the regular simplex based sampling strategy. The beneficial effect of the ellipsoidal trust-region framework is again clearly present as V2 to V4 converge faster than V1. In conclusion, one can observe that the same effect is present for V3 as was the case for the extended Powell singular function. Function evaluations are consumed at a higher rate in the initial stage. This phenomenon ultimately benefits the modelling procedure in the final optimization stage, resulting in a steeper convergence history and even the surpassing of V2 with respect to convergence speed.


## 6.2.1. Benchmark Problems

By assessing the results displayed in Table3, one can observe that, on the whole, the proposed algorithm outperforms both NMSMAX and ORBIT, disregarding some coincidental draws.Considering that V1 outperformed ORBIT in all but one occasion, we can conclude that the regular simplex-based sampling strategy offers a significant advantage. One might, however, argue that the increased implementational and computational complexity of Algorithm 3 does not justify the advantageous optimization setting in case that less challenging optimization problems are considered. In none of the benchmark problems did V1 outperform V2, V3, or V4, indicating that engaging an ellipsoidal trust-region framework is clearly a beneficial endeavor regardless, given the low implementational cost. However, the downside is that to obtain a curved model, a minimum ofn+2 points are required, making the advanced regular simplex sampling strategy indispensable.Whether V2, V3, or V4 outperformed one another depends on the objective function at hand. However, on the whole, we could state that V2 is the better candidate for lower-dimensional problems, whilst V3 is the better candidate for higher dimensional problems, as will be confirmed in the next section. Also, considering the convergence history of the different algorithms, it was noted that the ORBIT algorithm occasionally stagnated on a local plateau of the objective function. Since this behavior was also exhibited by V1 (e.g., GULF, TRID; although V1 recovered each time), it demonstrates that compensating for the output behavior of the objective function by the input space-oriented sampling strategy and by altering the shape of the trust-regions is quintessential to recover from such apparent plateaus.


## 6.2.2. Extended Powell Singular and Extended Rosenbrock

To demonstrate and comment on the mechanics of the elaborated methods in greater detail, convergence results on the 8-dimensional extended Powell singular and extended Rosenbrock are presented in Figure4 and discussed below.Figure 4
Convergence history represented by the evolution of the best relative function value in the function of the total number of function evaluations. (a) Extended Powell Singular,n=8. (b) Extended Rosenbrock, n=8.
(a)(b)Consider the convergence history for the extended Powell singular function in Figure4(a). We will initially discuss the history of the algorithms V1, ORBIT, and NMSMAX. All exhibit approximately the same convergence behavior during the initial stage of the optimization (up to 100 function evaluations). After 100 function evaluations, the NMSMAX algorithm clearly starts to diverge from this trend. It is because of the fact that the NMSMAX algorithm is bounded to operate on its vertices. At the common point of approximately 100 function evaluations, it reaches a region where the simplex size needs to shrink sufficiently to attain a size that is small enough to accurately grasp the local curvature of the objective function. A similar phenomenon is observed for the ORBIT algorithm, which starts to diverge from the global trend approximately after 200 function evaluations. The phenomenon occurs later for the ORBIT algorithm since it is not bounded to operate on the vertices of its standard simplex, and the solution of (3) is only limited to the convex set, ℬxk,Δk. Therefore, it can maintain a larger trust-region throughout the iterations than NMSMAX. That algorithm V1 can maintain a larger trust-region longer than ORBIT, and NMSMAX is a result of the improved geometry of the interpolation points. Apart from that, algorithm V1 too fails to reach the benchmark condition within the admissible budget. Now, consider the behavior of algorithm V2. The sampling strategy is equivalent to that of V1, however, the modelling space is now normalized with respect to the estimated curvature of the objective function. It corresponds with an ellipsoidal trust-region geometry in the original optimization space. The effect is hardly observable during the initial stage but clearly demonstrates its usefulness after approximately 500 function evaluations, where V2 diverges from V1 and successfully reaches the benchmark condition. It is argued that because of the transformation, the modeling and optimization steps of the algorithm are now executed in a space, where the objective function behaves more as if it was a quadratic function. The effect of the parameter μ1 can be assessed by comparing V2, V3, and V4. The parameter determines how strictly the volume of an exact regular simplex has to be respected. It results in a larger number of interpolation points to be added in every iteration. Therefore, algorithm V3 and V4 diverge from the common trend at about the same moment V1 and ORBIT do. In this region of the optimization space, the size of the trust-region (simplex) starts to shrink for all algorithms. It means that less function evaluations can be recycled from one trust-region to another. Algorithms V3 and V4 set a stricter limit on the volume of the simplex and are, therefore, forced to add more points to complement the recycled set. As a result, they consume more function evaluations than V1 and V2 to cross the same region. Nonetheless, this increase in function evaluations benefits the overall modelling procedure, and V2 is surpassed by V3 eventually.These observations are confirmed by the convergence history for the extended Rosenbrock function, see Figure4(b). Here, the results are directly affected by the local behavior of the objective function and all algorithms start to diverge at the starting point. In this case, V1 does reach the benchmark condition, again illustrating the improved scattering of points with respect to state-of-the-art. We argue this is a result of the regular simplex based sampling strategy. The beneficial effect of the ellipsoidal trust-region framework is again clearly present as V2 to V4 converge faster than V1. In conclusion, one can observe that the same effect is present for V3 as was the case for the extended Powell singular function. Function evaluations are consumed at a higher rate in the initial stage. This phenomenon ultimately benefits the modelling procedure in the final optimization stage, resulting in a steeper convergence history and even the surpassing of V2 with respect to convergence speed.


## 7. Conclusion

This paper is about reducing the number of objective function evaluations in derivative-free numerical optimization with surrogate models. Here, the radial basis functions are considered to be surrogate models with the interpolation framework provided by universal kriging. In each iteration of an interpolation-based optimization algorithm, a trust region can be used, wherein an optimization is performed using the surrogate model. The interpolation framework consists of using an interpolation set, of which, a subset of interpolation points needs to be well-poised so that the quality of the surrogate model is maximal. Interpolation set management aims at maximizing this range of model validity.Currently, at leastq=n+1 interpolation points are considered to guarantee a valid surrogate model. We advanced on interpolation-based optimization by relying on a subset of at least q=n+2 interpolation points within the interpolation set management to improve the scattering of the interpolation points, but more importantly, to have surrogate models that exhibit curvature. For that purpose, we defined fully linear surrogate models in Definition 1 that can be incorporated in the trust-region Algorithm 2. We implemented RBF models using a flexible multivariate interpolation framework for which universal kriging was used (see Section 3). In Lemma 1, we provide for a given RBF model with linear tail, information on the boundedness of the tolerance of the RBF model, objective function, and information on their gradients, i.e., constants Кf and Кg in Definition 1. In this manner, we are able to predict whether an interpolation model can be classified as being fully linear. Since we consider q=n+2 interpolation points that represent a regular simplex (see Section 4.3), we propose in algorithms 3 and 4 how the interpolation set management needs to be adapted, i.e., how the input points that need to be evaluated in the objective function need to be chosen.By leveraging on the curvature of the surrogate model that is valid in the presented interpolation set management, we introduced an ellipsoidal trust-region framework that allows morphed trust-regions that adapt to the objective function geometry. This adaptation is implemented by means of curvature normalization with coordinate transformation (24) and filtered Hessian approximation (29). In this manner, we allow to incorporate the output behavior of the objective function in the shape of the trust-region.The proposed interpolation set management with proper input of interpolation points that follow a well-poised geometry, together with the output behavior of the objective function that is incorporated in the ellipsoidal trust-region, have shown their benefits with respect to the reduction of the number of objective function evaluations when performing elaborate numerical optimizations on various test functions withn≤8. To show the impact of internal parameters of the algorithms, i.e., μ1 in Algorithm 3, which is related to the input set, and μ2 in (29), which is related to the ellipsoidal trust-region shape, the results of 4 different versions of the presented algorithm are shown in Figure 4 and Table 3. The results clearly illustrate the impact of the used input interpolation points and the ellipsoidal trust region. Moreover, the convergence histories show faster convergence compared to algorithms incorporating at least n+1 interpolation points and compared to those that use a spherical trust-region [61].



---

*Source: 8362294-2022-03-02.xml*